cvsp_make_actlay.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\cvsp_make_actlay.f
00002 !
00080                      SUBROUTINE CVSP_MAKE_ACTLAY
00081 !                    ***************************
00082 !
00083 !
00084 !***********************************************************************
00085 ! SISYPHE   V6P3                                   12/02/2013
00086 !***********************************************************************
00087 !
00088 !
00089 !
00090 !
00091 !
00092 !
00093 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00094 !| ELAY0              |--| ALTHICKNESS INIT
00095 !| AVAIL(J,K,I)       |--| FRACTIONS
00096 !| ELAY%R(J)          |--| ALTHICKNESS REAL FOR POINT(J)
00097 !| ESTRAT%R(J)        |--| ASTHICKNESS REAL FOR POINT(J)
00098 !| ES(J,K)            |--| LAYERTHICKNESS FOR POINT(J) / LAYER(K)
00099 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00100 
00101       USE BIEF_DEF
00102       USE BIEF
00103       USE DECLARATIONS_SISYPHE
00104 
00105       IMPLICIT NONE
00106       INTEGER LNG,LU
00107       COMMON/INFO/LNG,LU
00108 !
00109 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00110 !
00111       LOGICAL DB
00112       INTEGER  I,J,K,M,JG, II
00113       DOUBLE PRECISION TEMP, Z_HIGH, Z_LOW, CVSP_INTEGRATE_VOLUME, AT
00114       DOUBLE PRECISION ASUM, NEW_ALT
00115 !
00116 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00117 !
00118       DOUBLE PRECISION, EXTERNAL :: CVSP_ALT
00119       DOUBLE PRECISION SUMES, SUMAV, ALT
00120       CHARACTER(LEN=30) DEBUGFILE
00121 !
00122 !-----------------------------------------------------------------------
00123 !
00124       AT = DT*LT/PERCOU
00125 !
00126       DO J=1,NPOIN
00127         NLAYER%I(J) = NOMBLAY
00128         NEW_ALT =  CVSP_ALT(J,ALT_MODEL)
00129 !
00130         IF(NEW_ALT.GT.(ZF%R(J)-ZR%R(J))) THEN
00131           NEW_ALT = (ZF%R(J)-ZR%R(J))
00132         ENDIF
00133 !
00134         SUMES = 0.D0
00135         DO K=1,1               !NLAYER%I(J)
00136           IF(K.EQ.1) THEN
00137             ALT = NEW_ALT
00138           ELSE
00139             ALT = ELAY0
00140           ENDIF
00141           SUMES = SUMES + ALT
00142 !
00143 !-----------------------------------------------------------------------
00144 ! ALT DOESN'T USE THE FULL DEPTH
00145 !-----------------------------------------------------------------------
00146 !
00147           IF((ZF%R(J)-ZR%R(J)).GE.NEW_ALT) THEN
00148 !
00149 !-----------------------------------------------------------------------
00150 ! SUMES < FULL DEPTH
00151 !-----------------------------------------------------------------------
00152 !
00153             IF(SUMES.LE.(ZF%R(J)-ZR%R(J))) THEN
00154               IF(K.EQ.1) THEN
00155                 ES(J,K) = NEW_ALT
00156               ELSE
00157                 IF(((ZF%R(J)-ZR%R(J))-SUMES).GT.2.D0*ELAY0) THEN
00158                   ES(J,K) = ELAY0
00159                 ENDIF
00160                 IF(((ZF%R(J)-ZR%R(J))-SUMES).LT.2.D0*ELAY0) THEN
00161                    ES(J,K) = (ZF%R(J)-ZR%R(J)) - SUMES
00162                   NLAYER%I(J) = K
00163                 ELSE
00164                   ES(J,K) = 0.D0
00165                 ENDIF
00166               ENDIF         ! K=1
00167 
00168               IF (K == 1) THEN
00169                 Z_HIGH = PRO_D(J,PRO_MAX(J),1)
00170               ELSE
00171                 Z_HIGH = PRO_D(J,PRO_MAX(J),1) - SUMES + ALT
00172                 IF (Z_HIGH.GT.PRO_D(J,PRO_MAX(J),1)) THEN
00173                   Z_HIGH=PRO_D(J,PRO_MAX(J),1)
00174                 ENDIF
00175               ENDIF
00176 
00177               Z_LOW =  Z_HIGH - ES(J,K)
00178 !
00179 !-----------------------------------------------------------------------
00180 ! THIS IS THE CORE!
00181 !-----------------------------------------------------------------------
00182 !
00183               ASUM = 0.D0
00184               TEMP = CVSP_INTEGRATE_VOLUME(J,1,Z_HIGH,Z_LOW,T1%R)
00185 !
00186 !-----------------------------------------------------------------------
00187 ! ASSIGN FRACTIONS
00188 !-----------------------------------------------------------------------
00189 !
00190               DO I=1,NSICLA
00191                 AVAIL(J,K,I) = T1%R(I) / ES(J,K)
00192                 ASUM = AVAIL(J,K,I) + ASUM
00193 
00194                 IF ((AVAIL(J,K,I)>1+ZERO)) THEN
00195                   WRITE(LU,*) "MAKE_AL_", J, K, I, AT,
00196      &                 AVAIL(J,K,I), T1%R(I), ES(J,K),Z_HIGH,Z_LOW
00197      &                 , ES(J,K)
00198                   WRITE(LU,*) "ES,ALT,ELAY0,NEWALT",ES(J,K),ALT
00199      &                 ,ELAY0,NEW_ALT
00200                   CALL PLANTE(1)
00201                 ENDIF
00202               ENDDO
00203 !
00204 !-----------------------------------------------------------------------
00205 ! TRUNCATION ERRORS NORMALIZED FOR ALL FRACTIONS
00206 !-----------------------------------------------------------------------
00207 !
00208               IF(ABS(ASUM - 1.D0).NE.0.D0) THEN
00209                 DO I=1,NSICLA
00210                   AVAIL(J,K,I) = AVAIL(J,K,I) / ASUM
00211                 END DO
00212               ENDIF
00213 
00214             ELSE
00215               WRITE(LU,*) 'UHM_NOW_OBSOLETE',J, K, I, AT
00216               CALL PLANTE(1)
00217               STOP
00218 !
00219               ES(J,K) = MAX(-1.D0*(SUMES-(ZF%R(J)-ZR%R(J))), 0.D0)
00220 !
00221               IF(ES(J,K).EQ.0.D0) THEN
00222                 DO I=1,NSICLA
00223                   AVAIL(J,K,I) = 1.D0 / NSICLA
00224                 ENDDO
00225               ELSE
00226                 IF(DB(JG,0)) CALL CVSP_P('./ERR/','IVES3_',JG)
00227 !
00228 !-----------------------------------------------------------------------
00229 ! ASSIGN Z-COORDINATE FOR UPPER AND LOWER BOUNDARY OF LAYER
00230 !-----------------------------------------------------------------------
00231 !
00232                 Z_LOW = PRO_D(J,PRO_MAX(J),1)
00233                 DO M=1,K
00234                   Z_HIGH = Z_LOW
00235                   Z_LOW = Z_LOW - ES(J,M)
00236                 ENDDO
00237 !
00238                 IF(DB(JG,0)) WRITE(LU,*) 'UHM_L_',K,I
00239                 TEMP=CVSP_INTEGRATE_VOLUME(J,1,Z_HIGH,Z_LOW,T1%R)
00240 !
00241                 DO I=1,NSICLA
00242                   AVAIL(J,K,I) = T1%R(I) / ES(J,K)
00243                   ASUM = AVAIL(J,K,I) + ASUM
00244 !
00245                   IF ((AVAIL(J,K,I)>1+ZERO)) THEN
00246                     WRITE(LU,*) "MAKE_AL_", J, K, I, AT,
00247      &                   AVAIL(J,K,I), T1%R(I), ES(J,K)
00248      &                   ,Z_HIGH,Z_LOW
00249                   ENDIF
00250                 END DO            ! NSICLA
00251 
00252                 IF(ABS(ASUM - 1) > ZERO) THEN
00253                   WRITE(LU,*) "MAKE_AL_ABS2", J, K, I, AT,
00254      &                 AVAIL(J,K,I), TEMP, ES(J,K),Z_HIGH,Z_LOW,ASUM
00255 !
00256                   DO I=1,NSICLA
00257                     AVAIL(J,K,I) = AVAIL(J,K,I) / ASUM
00258                   ENDDO
00259                 ENDIF
00260               ENDIF     ! ES(J,K).EQ.0.D0
00261             ENDIF    ! SUMES.LE.(ZF%R(J)-ZR%R(J))
00262 !
00263           ELSE
00264             ES(J,K) = 0.D0
00265             ES(J,1) = (ZF%R(J)-ZR%R(J))
00266             NLAYER%I(J) = 1
00267             DO I=1,NSICLA
00268               AVAIL(J,K,I) = 1.D0 / NSICLA
00269             ENDDO
00270 !
00271           ENDIF ! (ZF%R(J)-ZR%R(J)).GE.NEW_ALT
00272 !
00273         IF (K == 1) ELAY%R(J) = ES(J,K)
00274         IF (K == 2) ESTRAT%R(J) = ES(J,K)
00275 !
00276 !-----------------------------------------------------------------------
00277 ! STOPS IF THE NUMBER OF LAYERS IS REDUCED DURING CALCULATION
00278 !-----------------------------------------------------------------------
00279 !
00280         IF(NLAYER%I(J) == K) EXIT
00281 !
00282       ENDDO                     ! K
00283 !
00284 !-----------------------------------------------------------------------
00285 ! NOT ENOUGH SPACE FOR NLAYER
00286 !-----------------------------------------------------------------------
00287 !
00288       DO K=1,1                  !NLAYER%I(J)
00289         IF(ES(J,K).EQ.0.D0) THEN
00290           WRITE(LU,*) 'NOT ENOUGH SPACE FOR NLAYER',J
00291           WRITE(LU,*) 'POSSIBLE ERROR! RIGID BED? ',NLAYER%I(J)-1
00292           NLAYER%I(J) = NLAYER%I(J)-1
00293           EXIT
00294         ENDIF
00295       END DO
00296 !
00297       ENDDO                     ! J
00298 !
00299 !-----------------------------------------------------------------------
00300 ! CHECKS
00301 !-----------------------------------------------------------------------
00302 !
00303       DO J=1,NPOIN
00304         SUMES = 0.D0
00305         JG = J
00306         IF (NCSIZE > 1) JG = MESH%KNOLG%I(J)
00307         DO K=1,9
00308           SUMES = SUMES+ES(J,K)
00309           SUMAV = 0.D0
00310           DO I=1,NSICLA
00311             SUMAV=AVAIL(J,K,I)+SUMAV
00312             IF((AVAIL(J,K,I).GT.1.D0+1.D-12).OR.
00313      &          (AVAIL(J,K,I).LT.0.D0-1.D-12)) THEN
00314               WRITE(LU,*) 'ERROR AVAIL',J,K,I,AVAIL(J,K,I)
00315               CALL CVSP_P('./ERR/','ERR_AVAIL_VSP_',  JG)
00316               CALL LAYERS_P('./ERR/ERR_AVAIL_LAY_', JG)
00317               CALL PLANTE(1)
00318             ENDIF
00319           ENDDO              !I
00320 !
00321           IF((ABS(SUMAV-1.D0).LT.1.D-12).AND.
00322      &        (ABS(SUMAV-1.D0).GT.0.D0)) THEN
00323             AVAIL(J,K,1) = AVAIL(J,K,1) - (1.D0 - SUMAV)
00324           ENDIF
00325 !
00326           IF(ABS(SUMAV-1.D0).GT.1.D-12) THEN
00327             IF(ABS(SUMAV-1.D0).GT.1.D-4) THEN
00328               WRITE(LU,*) 'ERROR SUMAV TOO BAD:',J,K,SUMAV
00329               WRITE(UNIT=DEBUGFILE, FMT='(A,I4,A)')
00330      &             './ERR/','SUMAV_',JG,'_VSP_'
00331               WRITE(UNIT=DEBUGFILE, FMT='(A,I4,A)')
00332      &             './ERR/','SUMAV_',JG,'_LAY_'
00333               DO II=1,NSICLA
00334                 AVAIL(J,K,II) = AVAIL(J,K,II) / SUMAV
00335               ENDDO
00336             ENDIF
00337           ENDIF
00338         ENDDO                 !K
00339 !
00340         IF( ABS(SUMES-(ZF%R(J)-ZR%R(J))).LT.ZERO) THEN
00341           WRITE(LU,*) 'ERROR SUM ES',J,SUMES,ZF%R(J)-ZR%R(J)
00342           CALL CVSP_P('./ERR/','ERR_SUM_ES_VSP_', JG)
00343           CALL LAYERS_P('./ERR/ERR_SUM_ES_LAY_', JG)
00344           CALL PLANTE(1)
00345         ENDIF
00346 !
00347       ENDDO                    !J
00348 !
00349 !-----------------------------------------------------------------------
00350 !
00351       RETURN
00352       END SUBROUTINE CVSP_MAKE_ACTLAY

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