cvsp_integrate_volume.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\cvsp_integrate_volume.f
00002 !
00067              DOUBLE PRECISION FUNCTION CVSP_INTEGRATE_VOLUME
00068 !            ***********************************************
00069 !
00070      &(J,I,Z_HIGH,Z_LOW,A)
00071 !
00072 !***********************************************************************
00073 ! SISYPHE   V6P3                                   14/03/2013
00074 !***********************************************************************
00075 !
00076 !
00077 !
00078 !
00079 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00080 !| J              |<--| INDEX OF A POINT IN MESH
00081 !| I              |<--| INDEX OF A FRACTION IN VERTICAL SORTING PROFILE
00082 !| Z_HIGH         |<--| HIGHER DEPTH COORDINATE
00083 !| Z_LOW          |<--| LOWER  DEPTH COORDINATE
00084 !| A1...A10       |<--| INTEGRATED VOLUMES PER FRACTION
00085 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00086 !
00087       USE DECLARATIONS_SISYPHE
00088 !
00089       IMPLICIT NONE
00090       INTEGER LNG,LU
00091       COMMON/INFO/LNG,LU
00092 !
00093 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00094 !
00095       INTEGER,          INTENT(IN)  :: J
00096       INTEGER,          INTENT(IN)  :: I
00097       DOUBLE PRECISION, INTENT(IN)  :: Z_HIGH
00098       DOUBLE PRECISION, INTENT(IN)  :: Z_LOW
00099       DOUBLE PRECISION, INTENT(OUT) :: A(10)
00100 !
00101 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00102 !
00103       DOUBLE PRECISION TEMP,TEMPOLD, DHIG, DLOW, AT, SUMUP, FLOW, FHIG
00104       DOUBLE PRECISION TEMP2, TEMP2MAX, SUMUP2,TEMP3, TEMP3MAX, SUMUP3
00105       DOUBLE PRECISION CORRECT, CHSUM
00106       INTEGER L_CNT, MYCASE, F_CNT, REVCNT, HELPER, LASTCASE, JG, K
00107       LOGICAL RET,CVSP_CHECK_F
00108       LOGICAL, EXTERNAL:: DB
00109 !
00110 !-----------------------------------------------------------------------
00111 !
00112       JG = J
00113       IF (NCSIZE.GT.1) JG = MESH%KNOLG%I(J)
00114 !
00115       AT = DT*LT/PERCOU
00116       SUMUP = 0.D0
00117       SUMUP2  = 0.D0
00118       SUMUP3  = 0.D0
00119       MYCASE  = 0
00120 !
00121 !-----------------------------------------------------------------------
00122 ! DOING ALL FRACTIONS
00123 !-----------------------------------------------------------------------
00124 !
00125       DO F_CNT = 1, NSICLA
00126         TEMP = 0.D0
00127         TEMP2 = 0.D0
00128         TEMP2MAX = 0.D0
00129         TEMP3 = 0.D0
00130         TEMP3MAX = 0.D0
00131         CHSUM = 0.D0
00132 !
00133 !-----------------------------------------------------------------------
00134 ! GOING THROUGH ALL SECTIONS
00135 !-----------------------------------------------------------------------
00136 !
00137         DO L_CNT = 0,(PRO_MAX(J)-2)
00138           TEMPOLD = TEMP
00139           REVCNT = PRO_MAX(J)-L_CNT
00140 !
00141 !-----------------------------------------------------------------------
00142 ! DEPTH COORDINATES OF THE SECTION TO CHECK
00143 !-----------------------------------------------------------------------
00144 !
00145           DHIG = PRO_D(J,REVCNT,F_CNT)
00146           DLOW = PRO_D(J,REVCNT-1,F_CNT)
00147 !
00148 !-----------------------------------------------------------------------
00149 ! SORTING PROFILE SECTION TOTALLY INSIDE (CASE ZDDZ)
00150 !-----------------------------------------------------------------------
00151 !
00152           IF ( (DHIG <= Z_HIGH ) .AND.
00153      &        (DLOW >= Z_LOW  ) ) THEN
00154 
00155             FLOW = PRO_F(J,REVCNT-1,F_CNT)
00156             FHIG = PRO_F(J,REVCNT,F_CNT)
00157 !
00158             MYCASE  = MYCASE  + 1
00159             LASTCASE = 1
00160 !
00161             TEMP2 = 0.5D0*(FHIG+FLOW)*(DHIG-DLOW)  + TEMP2
00162             TEMP2MAX = Z_HIGH - DLOW
00163 !
00164             DO HELPER = 1, NSICLA
00165               CHSUM = PRO_F(J,REVCNT, HELPER) + CHSUM
00166             ENDDO
00167 !
00168             CHSUM = 1.D0 - CHSUM
00169 !
00170 !-----------------------------------------------------------------------
00171 ! SORTING PROFILE SECTION PARTIALLY LOWER (CASE ZDZD)
00172 !-----------------------------------------------------------------------
00173 !
00174           ELSEIF ((DHIG <= Z_HIGH) .AND.
00175      &            (DHIG >  Z_LOW ) .AND.
00176      &            (DLOW <  Z_LOW ) )  THEN
00177 
00178             FHIG = PRO_F(J,REVCNT,F_CNT)
00179             FLOW = PRO_F(J,REVCNT-1,F_CNT)
00180             FLOW = - ((FHIG-FLOW)/(DHIG-DLOW))*(DHIG-Z_LOW) + FHIG
00181 !
00182 !-----------------------------------------------------------------------
00183 ! CUT THE SECTION
00184 !-----------------------------------------------------------------------
00185 !
00186             DLOW = Z_LOW
00187 !
00188             MYCASE  = MYCASE  + 1000
00189             LASTCASE = 2
00190 !
00191             TEMP3 = 0.5D0*(FHIG+FLOW)*(DHIG-DLOW)  + TEMP3
00192             TEMP3MAX = DHIG - DLOW
00193 !
00194 !-----------------------------------------------------------------------
00195 ! SORTING PROFILE SECTION PARTIALLY HIGHER (CASE DZDZ)
00196 !-----------------------------------------------------------------------
00197 !
00198           ELSEIF ((DHIG > Z_HIGH) .AND.
00199      &            (DLOW >= Z_LOW) .AND.
00200      &            (DLOW <  Z_HIGH) ) THEN
00201 !
00202             FLOW = PRO_F(J,REVCNT-1,F_CNT)
00203 !
00204             FHIG = PRO_F(J,REVCNT,F_CNT)
00205             FHIG = ((FHIG-FLOW)/(DHIG-DLOW))*(Z_HIGH-DLOW) + FLOW
00206 !
00207 !-----------------------------------------------------------------------
00208 !INSERT
00209 !-----------------------------------------------------------------------
00210 !
00211 !CUT THE SECTION
00212             DHIG = Z_HIGH
00213 !
00214             MYCASE  = MYCASE + 1000000
00215             LASTCASE = 3
00216 !
00217 !-----------------------------------------------------------------------
00218 ! LAYER TOTALLY INSIDE ONE SECTION (CASE DZZD)
00219 !-----------------------------------------------------------------------
00220 !
00221           ELSEIF ((DHIG >= Z_HIGH) .AND.
00222      &            (DLOW <= Z_LOW) ) THEN
00223 !
00224             FHIG =
00225      &           - (PRO_F(J,REVCNT,F_CNT)-PRO_F(J,REVCNT-1,F_CNT)) /
00226      &           (DHIG-DLOW) *
00227      &           (DHIG-Z_HIGH)
00228      &           + PRO_F(J,REVCNT,F_CNT)
00229 !
00230             FLOW =
00231      &           - (PRO_F(J,REVCNT,F_CNT)-PRO_F(J,REVCNT-1,F_CNT)) /
00232      &           (DHIG-DLOW) *
00233      &           (DHIG-Z_LOW)
00234      &           + PRO_F(J,REVCNT,F_CNT)
00235 !
00236 !-----------------------------------------------------------------------
00237 ! CUT THE SECTION
00238 !-----------------------------------------------------------------------
00239 !
00240             DHIG = Z_HIGH
00241             DLOW = Z_LOW
00242 !
00243             LASTCASE = 8
00244             MYCASE  = MYCASE + 100000000
00245 !
00246 !-----------------------------------------------------------------------
00247 ! SECTION WITH 0 STRENGTH
00248 !-----------------------------------------------------------------------
00249 !
00250           ELSEIF (DHIG == DLOW) THEN
00251             FLOW = 0.D0
00252             FHIG = 0.D0
00253 !
00254             MYCASE  = MYCASE + 1000000000
00255             LASTCASE = 4
00256 !
00257 !A TRUE BUG!
00258           ELSEIF (Z_LOW > Z_HIGH) THEN
00259             FLOW = 0.D0
00260             FHIG = 0.D0
00261 !
00262             WRITE(LU,*) 'Z_LOW >= Z_HIGH', DHIG, DLOW, Z_HIGH, Z_LOW
00263             CALL CVSP_P('./','ZLOHI',JG)
00264             MYCASE  = MYCASE + 1000000000
00265             LASTCASE = 5
00266 !
00267 !A TRUE BUG!
00268           ELSEIF (DHIG < DLOW) THEN
00269             FLOW = 0.D0
00270             FHIG = 0.D0
00271 !
00272             WRITE(LU,*) 'DHIG <= DLOW', J,DHIG, DLOW, Z_HIGH, Z_LOW
00273             CALL CVSP_P('./','DLOHI',JG)
00274             MYCASE  = MYCASE + 1000000000
00275             LASTCASE = 6
00276 !
00277 !A SECTION THAT IS NOT INVOLVED / NOT A BUG!!
00278           ELSE
00279 !
00280             FLOW = 0.D0
00281             FHIG = 0.D0
00282 !
00283             MYCASE  = MYCASE + 1000000000
00284             LASTCASE = 7
00285 !
00286           ENDIF
00287 !
00288 !-----------------------------------------------------------------------
00289 ! TRAPEZOID FORMULA
00290 !-----------------------------------------------------------------------
00291 !
00292           TEMP = 0.5D0*(FHIG+FLOW)*(DHIG-DLOW) + TEMP
00293 
00294 !DEBUG!DEBUG!DEBUG!DEBUG!DEBUG!DEBUG
00295           IF (0.5D0*(FHIG+FLOW)*(DHIG-DLOW) < 0.D0) THEN
00296             WRITE(LU,FMT='(A,1X,2(I11,1X),11(G20.10,1X),1X,I11)')
00297      &           'INTEGRATE_VOL_ER_TMP:<0:'
00298      &           ,JG, I, AT, FHIG,FLOW,DHIG,DLOW, DHIG-DLOW, REVCNT,
00299      &           PRO_F(J,REVCNT-1,F_CNT),PRO_F(J,REVCNT,F_CNT),
00300      &           PRO_D(J,REVCNT-1,F_CNT),PRO_D(J,REVCNT,F_CNT),
00301      &           LASTCASE
00302             CALL CVSP_P('./','IVKT',JG)
00303             CALL PLANTE(1)
00304           ENDIF
00305 !DEBUG!DEBUG!DEBUG!DEBUG!DEBUG!DEBUG
00306 !
00307         ENDDO                 !SECTION
00308 !
00309 !-----------------------------------------------------------------------
00310 ! ADDING UP FRACTIONS FOR DEBUGGING PURPOSES
00311 !-----------------------------------------------------------------------
00312 !
00313         SUMUP = TEMP + SUMUP
00314         SUMUP2 = TEMP2 + SUMUP2
00315         SUMUP3 = TEMP3 + SUMUP3
00316 !
00317         A(F_CNT) = TEMP
00318 !
00319       ENDDO                    !FRACTION
00320 !
00321       CORRECT = (SUMUP / (Z_HIGH-Z_LOW))
00322       CVSP_INTEGRATE_VOLUME = A(I) / CORRECT
00323 !
00324 !DEBUG!DEBUG!DEBUG!DEBUG!DEBUG!DEBUG
00325       IF (CVSP_INTEGRATE_VOLUME < 0.D0) THEN
00326         CALL CVSP_P('./','IVK0',JG)
00327         WRITE(*,FMT='(A,2(I11),14(G20.10))')'INTEGRATE_VOL_ER:<0:'
00328      &       ,JG, I, AT, CVSP_INTEGRATE_VOLUME, MYCASE ,Z_HIGH,Z_LOW,
00329      &       SUMUP,SUMUP2,SUMUP3,CHSUM,A(1),A(2),A(3),A(4),A(5)
00330         CALL PLANTE(1)
00331       ENDIF
00332 !DEBUG!DEBUG!DEBUG!DEBUG!DEBUG!DEBUG
00333 !
00334 !-----------------------------------------------------------------------
00335 ! CHECKSUM OVER ALL FRACTIONS AND LAYERS
00336 !-----------------------------------------------------------------------
00337 !
00338       SUMUP  = SUMUP  - ABS(Z_HIGH-Z_LOW)
00339       SUMUP2 = SUMUP2 - ABS(TEMP2MAX)
00340       SUMUP3 = SUMUP3 - ABS(TEMP3MAX)
00341       CHSUM =  CHSUM / 5.D0
00342 !
00343 !-----------------------------------------------------------------------
00344 ! BEWARE: 1.D-5 IS UP TO 10G OF A 1000KG BUT HIGHER ACCURACY
00345 ! LEADS AGAIN TO FLOATING POINT TRUNCATION ERRORS ...
00346 !-----------------------------------------------------------------------
00347 !
00348       IF (ABS(SUMUP).GT.1.D-5) THEN
00349         CALL CVSP_P('./','IV_E',JG)
00350         WRITE(LU,*) 'INTEGRATE VOLUME ACCURRACY!!!', SUMUP, JG
00351         DO K = 1, PRO_MAX(J)
00352 !
00353 !-----------------------------------------------------------------------
00354 ! REMOVES NUMERIC INSTABILITIES
00355 !-----------------------------------------------------------------------
00356 !
00357           RET =  CVSP_CHECK_F(J,K,' IV_FIX:   ')
00358         ENDDO
00359       ENDIF
00360 !
00361 !-----------------------------------------------------------------------
00362 !
00363       RETURN
00364       END FUNCTION CVSP_INTEGRATE_VOLUME

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0