|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 !----------------------------------------------------------------- 00002 ! ########################### 00003 SUBROUTINE TOPODYN_LAT(PRW,PDEF,PKAPPA,PKAPPAC,GTOPD) 00004 !,PERROR) 00005 ! ########################### 00006 ! 00007 ! 00008 ! PURPOSE 00009 ! ------- 00010 ! to distribute laterally soil water following topodyn concept 00011 ! 00012 ! 00013 ! METHOD 00014 ! ------ 00015 ! 00016 ! EXTERNAL 00017 ! -------- 00018 ! none 00019 ! 00020 ! 00021 ! AUTHOR 00022 ! ------ 00023 ! 00024 ! G.-M. Saulnier * LTHE * 00025 ! K. Chancibault * CNRM * 00026 ! 00027 ! MODIFICATIONS 00028 ! ------------- 00029 ! 00030 ! Original 12/2003 00031 ! writing in fortran 90 12/2004 00032 !------------------------------------------------------------------------------------------ 00033 ! 00034 !* 0.0 DECLARATIONS 00035 ! ------------ 00036 USE MODD_TOPODYN, ONLY : NNCAT, NMESHT, XDMAXT, XDXT, XMPARA, NNMC, XCONN, NLINE,& 00037 XSLOP, XDAREA 00038 USE MODD_COUPLING_TOPD, ONLY : XWSTOPT, XDTOPT 00039 USE MODD_TOPD_PAR, ONLY : XSTEPK 00040 ! 00041 USE MODD_SURF_PAR, ONLY : XUNDEF 00042 ! 00043 USE MODI_FLOWDOWN 00044 USE MODI_ABOR1_SFX 00045 USE MODI_WRITE_FILE_VECMAP 00046 USE MODI_WRITE_FILE_MAP 00047 ! 00048 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00049 USE PARKIND1 ,ONLY : JPRB 00050 ! 00051 IMPLICIT NONE 00052 ! 00053 !* 0.1 declarations of arguments 00054 ! 00055 REAL, DIMENSION(:,:), INTENT(IN) :: PRW 00056 REAL, DIMENSION(:,:), INTENT(OUT) :: PDEF 00057 REAL, DIMENSION(:,:), INTENT(OUT) :: PKAPPA 00058 REAL, DIMENSION(:), INTENT(OUT) :: PKAPPAC 00059 LOGICAL, DIMENSION(:), INTENT(OUT) :: GTOPD 00060 ! 00061 !* 0.2 declarations of local variables 00062 ! 00063 LOGICAL :: GFOUND ! logical variable 00064 REAL :: ZSOMME 00065 REAL :: ZM ! XMPARA in m 00066 REAL :: ZDX ! XDXT in m 00067 REAL :: ZKVAL, ZKVALMIN, ZKVALMAX 00068 REAL :: ZDAV ! Averaged deficit (m) 00069 REAL :: ZDAV2 ! Averaged deficit on ZA-ZAS-ZAD (m) 00070 REAL :: ZNDMAXAV,ZNKAV ! temporary averaged maximal deficit and averaged similarity index 00071 REAL :: ZDMAXAV,ZKAV ! averaged maximal deficit and averaged similarity index 00072 REAL :: ZFUNC 00073 REAL :: ZDIF,ZDIFMIN ! difference calcul 00074 REAL :: ZNAS, ZNAD ! temporary saturated and dry relative catchment area 00075 REAL :: ZAS,ZAD ! saturated and dry relative catchment area 00076 REAL :: ZA ! total catchment area 00077 REAL :: ZTMP 00078 ! 00079 REAL, DIMENSION(NMESHT) :: ZDMAX ! XDMAXT in m 00080 REAL, DIMENSION(NMESHT) :: ZRW ! PRW in m 00081 REAL, DIMENSION(NMESHT) :: ZDINI ! initial deficit 00082 REAL, DIMENSION(NMESHT) :: ZMASK 00083 REAL, DIMENSION(NMESHT) :: ZKAPPA_PACK, ZDMAX_PACK 00084 !REAL, DIMENSION(NMESHT) :: ZTMP 00085 ! 00086 INTEGER :: J1, J2, JJ ! 00087 INTEGER :: INKAPPA ! number of steps in similarity index distribution 00088 INTEGER :: INPCON ! number of connected pixels 00089 INTEGER :: INAS ! number of saturated pixels 00090 INTEGER :: INAD ! number of dry pixels 00091 INTEGER :: I_DIM 00092 ! 00093 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00094 !----------------------------------------------------------------------------------- 00095 IF (LHOOK) CALL DR_HOOK('TOPODYN_LAT',0,ZHOOK_HANDLE) 00096 ! 00097 PKAPPA(:,: )= 0.0 00098 PKAPPAC(:) = 0. 00099 GTOPD(:) = .TRUE. 00100 PDEF(:,:) = 0. 00101 ! 00102 DO JJ = 1,NNCAT 00103 !* 0. Initialisation 00104 ! -------------- 00105 ZMASK(:) = 0.0 00106 ZRW(:) = 0.0 00107 ZDMAX(:) = 0.0 00108 INPCON = 0 00109 ZDAV = 0.0 00110 ZDAV2 = 0.0 00111 GFOUND = .FALSE. 00112 ZDIFMIN = 99999. 00113 ! 00114 ZRW(:) = PRW(JJ,:) 00115 ZDMAX(:) = XDMAXT(JJ,:) 00116 ZDX = XDXT(JJ) 00117 ZM = XMPARA(JJ) 00118 ZDINI(:) = ZDMAX(:) 00119 ! 00120 !* 0.2 definition of the catchment area concerned by lateral distribution 00121 ! ------------------------------------------------------------------ 00122 ! 00123 DO J1=1,NNMC(JJ) 00124 ! 00125 IF ( ZRW(J1).GT.0.0 ) THEN 00126 ZMASK(J1)=1.0 00127 ELSE 00128 ZMASK(J1)=0.0 00129 ENDIF 00130 ! 00131 ENDDO 00132 ! 00133 CALL FLOWDOWN(NNMC(JJ),ZMASK,XCONN(JJ,:,:),NLINE(JJ,:)) 00134 ! 00135 WHERE (ZMASK == 0.0) ZMASK = XUNDEF 00136 ! 00137 ! 00138 !* 1. Calcul of hydrological similarity and topographic indexes 00139 ! --------------------------------------------------------- 00140 !* 1.1 Calcul of averaged deficit and initialisation of indexes 00141 ! -------------------------------------------------------- 00142 ! 00143 ZA = NNMC(JJ) * ZDX**2 00144 ZTMP=0. 00145 ! 00146 DO J1=1,NNMC(JJ) 00147 ! 00148 IF (ZMASK(J1)/=XUNDEF) THEN 00149 ! 00150 PKAPPA(JJ,J1) = ZRW(J1) 00151 INPCON = INPCON + 1 00152 ZDINI(J1) = ZDMAX(J1) - ZRW(J1) 00153 ! 00154 IF ( ZDINI(J1) <0.0 ) THEN 00155 ! WRITE(*,*) J1,'Dini negatif !' 00156 ZTMP = ZTMP - ZDINI(J1) !we stock here water above saturation to be 00157 ! distributed among the others pixels 00158 ZDINI(J1) = 0. 00159 ENDIF 00160 ! 00161 ZDAV = ZDAV + ZDINI(J1) 00162 ! 00163 ELSE 00164 ! 00165 PKAPPA(JJ,J1) = XUNDEF 00166 ! 00167 ENDIF 00168 ! 00169 ENDDO 00170 ! 00171 IF (ZTMP>0.) THEN 00172 !write(*,*) COUNT(ZDINI(:)<0.),' pixels avec ZDINI negatif. Volume total :', ZTMP 00173 WHERE ( ZDINI(:)>0. ) ZDINI(:) = ZDINI(:)-ZTMP/(COUNT(ZDINI(:)>0.)) 00174 ENDIF 00175 ! 00176 IF (INPCON >= NNMC(JJ)/1000) THEN 00177 ! 00178 ZDAV = ZDAV / INPCON 00179 ZDAV = ZDAV / ZM 00180 ! 00181 !* 1.2 Propagation of indexes 00182 ! ---------------------- 00183 ! 00184 CALL FLOWDOWN(NNMC(JJ),PKAPPA(JJ,:),XCONN(JJ,:,:),NLINE(JJ,:)) 00185 ! 00186 !* 1.3 Distribution of indexes 00187 ! ---------------------- 00188 ! 00189 J2=1 00190 ! 00191 DO WHILE ( .NOT.GFOUND .AND. J2.LE.NNMC(JJ) ) 00192 ! 00193 IF (ZMASK(J2)/=XUNDEF) THEN 00194 ! 00195 GFOUND = .TRUE. 00196 ZKVAL = PKAPPA(JJ,J2) * XDAREA(JJ,J2) / XSLOP(JJ,J2) 00197 ZKVAL = LOG(ZKVAL) 00198 ZKVALMAX = ZKVAL 00199 ZKVALMIN = ZKVAL 00200 PKAPPA(JJ,J2) = ZKVAL 00201 ! 00202 ELSE 00203 ! 00204 PKAPPA(JJ,J2) = XUNDEF 00205 ! 00206 ENDIF 00207 ! 00208 J2 = J2 + 1 00209 ! 00210 ENDDO 00211 ! 00212 DO J1 = J2,NNMC(JJ) 00213 ! 00214 IF (ZMASK(J1)/=XUNDEF) THEN 00215 ! 00216 ZKVAL = PKAPPA(JJ,J1) * XDAREA(JJ,J1) / XSLOP(JJ,J1) 00217 ZKVAL = LOG(ZKVAL) 00218 ! 00219 IF (ZKVAL.GT.ZKVALMAX) THEN 00220 ZKVALMAX = ZKVAL 00221 ELSEIF (ZKVAL.LT.ZKVALMIN) THEN 00222 ZKVALMIN = ZKVAL 00223 ENDIF 00224 ! 00225 PKAPPA(JJ,J1) = ZKVAL 00226 ! 00227 ELSE 00228 ! 00229 PKAPPA(JJ,J1) = XUNDEF 00230 ! 00231 ENDIF 00232 ! 00233 ENDDO 00234 ! 00235 !* 1.4 Calcul of saturation index 00236 ! -------------------------- 00237 ! 00238 I_DIM = COUNT( ZMASK(1:NNMC(JJ))/=XUNDEF ) 00239 ZKAPPA_PACK(:) = XUNDEF 00240 ZDMAX_PACK (:) = XUNDEF 00241 ZKAPPA_PACK(1:I_DIM) = PACK(PKAPPA(JJ,1:NNMC(JJ)),ZMASK(1:NNMC(JJ))/=XUNDEF) 00242 ZDMAX_PACK (1:I_DIM) = PACK(ZDMAX (1:NNMC(JJ)),ZMASK(1:NNMC(JJ))/=XUNDEF) 00243 ! 00244 INKAPPA = INT((ZKVALMAX - ZKVALMIN) / XSTEPK) 00245 ! 00246 DO J1=1,INKAPPA 00247 ! 00248 ZKVAL = ZKVALMIN + (XSTEPK * (J1-1)) 00249 INAS = 0 00250 INAD = 0 00251 ZNDMAXAV = 0.0 00252 ZNKAV = 0.0 00253 ! 00254 DO J2=1,I_DIM 00255 ! 00256 IF ( ZKAPPA_PACK(J2).GE.ZKVAL ) THEN 00257 ! saturated pixel 00258 INAS = INAS + 1 00259 ELSEIF (ZKAPPA_PACK(J2).LE.( ZKVAL-(ZDMAX_PACK(J2)/ZM)) ) THEN 00260 ! dry pixel 00261 INAD = INAD + 1 00262 ZNDMAXAV = ZNDMAXAV + ZDMAX_PACK(J2) 00263 ELSE 00264 ZNKAV = ZNKAV + ZKAPPA_PACK(J2) 00265 ENDIF 00266 ! 00267 ENDDO 00268 ! 00269 IF (INAD == 0) THEN 00270 ZNDMAXAV = 0.0 00271 ELSE 00272 ZNDMAXAV = ZNDMAXAV / REAL(INAD) 00273 ENDIF 00274 ! 00275 IF ( INPCON == INAS .OR. INPCON == INAD .OR. INPCON == (INAD+INAS)) THEN 00276 ZNKAV = 0.0 00277 ELSE 00278 ZNKAV = ZNKAV / REAL(INPCON - INAD - INAS) 00279 ENDIF 00280 ! 00281 IF (INPCON /= 0) THEN 00282 ZNAS = REAL(INAS) / REAL(INPCON) 00283 ZNAD = REAL(INAD) / REAL(INPCON) 00284 ENDIF 00285 ! 00286 ZFUNC = (1 - ZNAS - ZNAD) * ( ZKVAL - ZNKAV ) 00287 IF (ZM /= 0.) ZFUNC = ZFUNC + (ZNAD * (ZNDMAXAV / ZM)) 00288 ! 00289 ZDIF = ABS( ZFUNC - ZDAV ) 00290 ! 00291 IF ( ZDIF.LT.ZDIFMIN ) THEN 00292 ! 00293 ZDIFMIN = ZDIF 00294 PKAPPAC(JJ) = ZKVAL 00295 ZAS = ZNAS 00296 ZAD = ZNAD 00297 ZDMAXAV = ZNDMAXAV 00298 ZKAV = ZNKAV 00299 ! 00300 ENDIF 00301 ! 00302 ENDDO 00303 ! 00304 !* 2. Local deficits calculation 00305 ! -------------------------- 00306 ! 00307 !* 2.1 New averaged deficit on A-Ad-As 00308 ! ------------------------------- 00309 ! 00310 ZDAV = ZDAV * ZM 00311 ! 00312 IF ( ZAS<1. .AND. ZAD<1. ) THEN 00313 ! 00314 ZDAV2 = (ZDAV - ZDMAXAV * ZAD) / (1 - ZAS - ZAD) 00315 ! 00316 ELSE 00317 ! 00318 IF (ZAS>=1.) WRITE(*,*) 'ALL THE AREA IS SATURATED' 00319 IF (ZAD>=1.) WRITE(*,*) 'ALL THE AREA HAS A MAXIMAL DEFICIT' 00320 WRITE(*,*) 'ALL THE AREA',ZAS,ZAD 00321 CALL ABOR1_SFX("TOPODYN_LAT: ALL THE AREA IS SATURATED OR HAS A MAXIMAL DEFICIT") 00322 ! 00323 ENDIF 00324 ! 00325 !* 2.2 Local deficits 00326 ! -------------- 00327 ! 00328 ZSOMME=0.0 00329 !ZTMP(:)=XUNDEF 00330 ! 00331 DO J1=1,NNMC(JJ) 00332 ! 00333 IF ( ZMASK(J1)/=XUNDEF ) THEN 00334 ! 00335 IF ( (PKAPPA(JJ,J1).GT.(PKAPPAC(JJ) - ZDMAX(J1)/ZM)) .AND. (PKAPPA(JJ,J1).LT.PKAPPAC(JJ)) ) THEN 00336 ! 00337 PDEF(JJ,J1) = ZM * (ZKAV - PKAPPA(JJ,J1)) + ZDAV2 00338 !ZTMP(J1) = 0.5 00339 IF (PDEF(JJ,J1) < 0.0) PDEF(JJ,J1) = 0.0 00340 ! 00341 ELSEIF ( PKAPPA(JJ,J1).GE.PKAPPAC(JJ) ) THEN 00342 ! 00343 PDEF(JJ,J1) = 0.0 00344 !ZTMP(J1) = 1. 00345 ! 00346 ELSEIF ( PKAPPA(JJ,J1).LE.(PKAPPAC(JJ) - ZDMAX(J1)/ZM) ) THEN 00347 ! 00348 PDEF(JJ,J1) = ZDMAX(J1) 00349 !ZTMP(J1) = -1.0 00350 ! 00351 ENDIF 00352 ! 00353 ! nouveau contenu en eau total (m) 00354 ZSOMME = ZSOMME + ( XWSTOPT(JJ,J1)*XDTOPT(JJ,J1) - PDEF(JJ,J1) ) 00355 ! 00356 ELSE 00357 ! 00358 PDEF(JJ,J1) = ZDMAX(J1) 00359 ! 00360 ENDIF 00361 ! 00362 ENDDO 00363 ! 00364 ! variation du contenu en eau total 00365 DO J1=1,NNMC(JJ) 00366 ! 00367 IF (PDEF(JJ,J1)<0.0) THEN 00368 WRITE(*,*) 'LAMBDA=',PKAPPA(JJ,J1),'LAMBDAC=',PKAPPAC(JJ) 00369 ENDIF 00370 ! 00371 ENDDO 00372 ! 00373 GTOPD(JJ)=.TRUE. 00374 ! 00375 ELSE 00376 ! 00377 ! 'Pas de redistribution laterale' 00378 GTOPD(JJ)=.FALSE. 00379 ! 00380 PKAPPA(JJ,:) = XUNDEF 00381 ! PDEF(JJ,:) = ZDMAX(:) - ZRW(:) 00382 PDEF(JJ,:) = ZDINI(:) 00383 PKAPPAC(JJ) = XUNDEF 00384 ! 00385 ENDIF 00386 ! 00387 ENDDO 00388 ! 00389 IF (LHOOK) CALL DR_HOOK('TOPODYN_LAT',1,ZHOOK_HANDLE) 00390 ! 00391 END SUBROUTINE TOPODYN_LAT 00392
1.8.0
