|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######spl 00002 SUBROUTINE INIT_TOP (HISBA, HTOPREG, KLUOUT, PPATCH, PRUNOFFD, & 00003 PDZG, PWWILT, PWSAT, PTI_MIN, & 00004 PTI_MAX, PTI_MEAN, PTI_STD, PTI_SKEW, & 00005 PSOILWGHT, PTAB_FSAT, PTAB_WTOP, PM ) 00006 ! 00007 ! ##################################################################### 00008 ! 00009 !!**** *INIT_TOP* 00010 !! 00011 !! PURPOSE 00012 !! ======= 00013 ! 00014 ! Calculates the new array of the Datin-Saulnier TOPMODEL framework fonction for xsat and compute each 00015 ! satured fraction for each xsat value of the grids cells but also the active TOPMODEL-layer array, the 00016 ! driest fraction array and the normalized mean deficit array. 00017 ! For calculate new array, we use the incomplete gamma function. (see gamma_inc.f for more detail) 00018 ! 00019 ! Note that over land point where topographic index do not exist, a VIC 00020 ! distribution is used with Bcoef at least equal to 0.1. This value can be 00021 ! change in namelist 00022 ! 00023 !------------------------------------------------------------------------------- 00024 ! 00025 USE MODD_SURF_PAR,ONLY : XUNDEF 00026 USE MODD_ISBA_n, ONLY : NPATCH, NGROUND_LAYER, NSIZE_NATURE_P 00027 ! 00028 USE MODD_SGH_PAR, ONLY : X2, X4, XREGP, XREGA 00029 ! 00030 USE MODI_DGAM 00031 USE MODI_GAMMAS 00032 ! 00033 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00034 USE PARKIND1 ,ONLY : JPRB 00035 ! 00036 USE MODI_ABOR1_SFX 00037 ! 00038 !* 0. DECLARATIONS 00039 ! =================== 00040 ! 00041 IMPLICIT NONE 00042 ! 00043 !* 0.1 declarations of arguments 00044 ! 00045 CHARACTER(LEN=*), INTENT(IN) :: HISBA ! type of ISBA version: 00046 ! ! '2-L' (default) 00047 ! ! '3-L' 00048 ! ! 'DIF' 00049 ! 00050 CHARACTER(LEN=*), INTENT(IN) :: HTOPREG ! Wolock and McCabe (2000) linear regression for Topmodel 00051 ! 'DEF' = Reg 00052 ! 'NON' 00053 ! 00054 INTEGER, INTENT(IN) :: KLUOUT 00055 ! 00056 REAL, DIMENSION(:,:), INTENT(IN) :: PPATCH 00057 ! 00058 REAL, DIMENSION(:,:,:),INTENT(IN) :: PDZG ! soil layer thickness 00059 ! 00060 REAL, DIMENSION(:,:,:),INTENT(IN) :: PSOILWGHT ! ISBA-DIF: weights for vertical 00061 ! ! integration of soil water and properties 00062 ! 00063 REAL, DIMENSION(:,:), INTENT(IN) :: PRUNOFFD ! depth over which sub-grid runoff is computed 00064 ! 00065 REAL, DIMENSION(:,:), INTENT(IN) :: PWWILT, PWSAT 00066 ! PWWILT = the wilting point volumetric 00067 ! water content (m3 m-3) 00068 ! PWSAT = saturation volumetric water content 00069 ! of the soil (m3 m-3) 00070 ! 00071 REAL, DIMENSION(:), INTENT(IN) :: PTI_MIN, PTI_MAX, PTI_STD, 00072 PTI_SKEW, PTI_MEAN 00073 ! PTI_MEAN = ti mean 00074 ! PTI_MIN = ti min 00075 ! PTI_MAX = ti max 00076 ! PTI_STD = ti standard deviation 00077 ! PTI_SKEW = ti skewness 00078 ! 00079 REAL, DIMENSION(:,:), INTENT(INOUT) :: PTAB_FSAT, PTAB_WTOP 00080 ! PTAB_FSAT = Satured fraction array 00081 ! PTAB_WTOP = Active TOPMODEL-layer array 00082 ! 00083 REAL, DIMENSION(:), INTENT(INOUT) :: PM 00084 ! PM = exponential decay factor of the local deficit 00085 ! 00086 !* 0.2 declarations of local variables 00087 ! 00088 REAL, DIMENSION(SIZE(PM)) :: ZD_TOP, ZWSAT_AVG, ZWWILT_AVG 00089 ! ZD_TOP = Topmodel active layer 00090 ! 00091 REAL :: ZXI, ZPHI, ZNU, ZTI_MEAN, ZTI_MIN, ZTI_MAX, ZTI_STD, ZTI_SKEW 00092 ! ZTI_MEAN= ti mean after regression 00093 ! ZXI = ti pdf parameter 00094 ! ZPHI = ti pdf parameter 00095 ! ZNU = ti pdf parameter 00096 ! 00097 REAL :: ZFTOT, ZYMAX, ZYMIN 00098 REAL :: ZGYMAX, ZGYMIN 00099 ! ZFTOT = total fraction of a grid cell 00100 ! ZYMAX = yi maximum variable 00101 ! ZYMIN = yi minimum variable 00102 ! ZGYMAX = incomplete gamma function for ymax (GAMSTAR result) 00103 ! ZGYMIN = incomplete gamma function for ymin (GAMSTAR result) 00104 ! 00105 REAL :: ZXSAT_IND, ZYSAT, ZY0, ZDMOY, ZXMOY, ZFMED, ZF0 00106 ! ZXSAT_IND = Satured index for all index 00107 ! ZYSAT = changing variable of satured index 00108 ! ZY0 = changing variable of dry index 00109 ! ZDMOY = grid cell average deficit (= Dbar/M) 00110 ! ZXMOY = ti mean value on wet (1-fsat-f0) fraction 00111 ! ZF0 = dry fraction 00112 ! ZFMED = wet (1-fsat-f0) fraction 00113 ! 00114 REAL :: ZG, ZGYSAT, ZGY0 00115 ! ZG = GAM result 00116 ! ZGYSAT = the incomplete gamma function for ysat (GAMSTAR result) 00117 ! ZGY0 = the incomplete gamma function for y0 (GAMSTAR result) 00118 ! 00119 REAL :: ZD0, ZPAS, ZCOR 00120 ! ZD0 = Normalized TOPMODEL maximum deficit D0/M coefficient 00121 ! ZPAS = pas for calculate the new xsat values array 00122 ! 00123 INTEGER :: IFLG, IFLGST 00124 ! IFLG = incomplete gamma function error flag (GAM result) 00125 ! IFLGST = incomplete gamma function error flag (GAMSTAR result) 00126 ! (see gamma_inc.f for more detail) 00127 ! 00128 REAL :: ZNO, ZAR, ZTOT 00129 ! 00130 REAL :: ZFUP, ZFDOWN, ZWUP, ZWDOWN, ZSLOPEW 00131 ! 00132 INTEGER, DIMENSION (1):: ID 00133 ! 00134 INTEGER :: INI, I, IND, JSI_MIN, JSI_MAX, IPAS, 00135 JL, INL, JPATCH 00136 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00137 ! 00138 !------------------------------------------------------------------------------- 00139 ! 00140 ! 1 TOPMODEL SURFACE RUNOFF SCHEME 00141 ! ================================ 00142 ! 00143 ! 1.1 initialisation of local variable 00144 ! ---------------------------------------- 00145 ! 00146 ! Grid cells number 00147 ! 00148 IF (LHOOK) CALL DR_HOOK('INIT_TOP',0,ZHOOK_HANDLE) 00149 INI = SIZE(PM(:)) 00150 INL = SIZE(PWSAT,2) 00151 IPAS = SIZE(PTAB_FSAT,2) 00152 ! 00153 ! GAM result (not use here !) 00154 ! 00155 ZG = 0.0 00156 ! 00157 ZD_TOP (:) = 0.0 00158 ZWSAT_AVG (:) = 0.0 00159 ZWWILT_AVG(:) = 0.0 00160 ! 00161 ! soil properties for runoff (m) 00162 ! 00163 IF (HISBA == 'DIF') THEN 00164 ! 00165 DO JPATCH=1,NPATCH 00166 IF (NSIZE_NATURE_P(JPATCH) == 0 ) CYCLE 00167 DO JL=1,INL 00168 DO I=1,INI 00169 ZD_TOP (I) = ZD_TOP (I) + PPATCH(I,JPATCH)*PSOILWGHT(I,JL,JPATCH) 00170 ZWSAT_AVG (I) = ZWSAT_AVG (I) + PPATCH(I,JPATCH)*PSOILWGHT(I,JL,JPATCH)*PWSAT (I,JL) 00171 ZWWILT_AVG(I) = ZWWILT_AVG(I) + PPATCH(I,JPATCH)*PSOILWGHT(I,JL,JPATCH)*PWWILT(I,JL) 00172 ENDDO 00173 ENDDO 00174 ENDDO 00175 ! 00176 WHERE(ZD_TOP(:)>0.0) 00177 ZWSAT_AVG (:)=ZWSAT_AVG (:)/ZD_TOP(:) 00178 ZWWILT_AVG(:)=ZWWILT_AVG(:)/ZD_TOP(:) 00179 ENDWHERE 00180 ! 00181 ELSE 00182 ! 00183 DO JPATCH=1,NPATCH 00184 IF (NSIZE_NATURE_P(JPATCH) == 0 ) CYCLE 00185 DO I=1,INI 00186 ZD_TOP(I)=ZD_TOP(I)+PRUNOFFD(I,JPATCH)*PPATCH(I,JPATCH) 00187 ENDDO 00188 ENDDO 00189 ! 00190 ZWSAT_AVG (:) = PWSAT (:,1) 00191 ZWWILT_AVG(:) = PWWILT(:,1) 00192 ! 00193 ENDIF 00194 ! 00195 ! 00196 ! 1.2 Algorithme 00197 ! ------------------ 00198 ! 00199 !grid cells loops 00200 ! 00201 ZAR = 0.0 00202 ZTOT= 0.0 00203 ZNO = 0.0 00204 ! 00205 DO I=1,INI 00206 ! 00207 IF(PTI_MEAN(I)==XUNDEF)THEN 00208 ! 00209 ! *Case where the Topographics index are not defined. 00210 ! -------------------------------------------------------- 00211 ZNO=ZNO+1.0 00212 PTAB_FSAT(I,:)=0.0 00213 PTAB_WTOP(I,:)=XUNDEF 00214 ! 00215 PM(I) =XUNDEF 00216 ! 00217 ELSE 00218 ! 00219 ZTOT = ZTOT + 1.0 00220 ! 00221 ! 1.2.0 first initialisation 00222 ! -------------------------- 00223 ! 00224 ZXI = 0.0 00225 ZPHI = 0.0 00226 ZNU = 0.0 00227 ! 00228 ! Wolock and McCabe (2000) linear regression equation between the mean 00229 ! topographic index computed with a 1000 meter DEM and a 100 meter DEM. 00230 ! 00231 IF(HTOPREG=="DEF")THEN 00232 ZTI_MEAN=XREGP*PTI_MEAN(I)-XREGA 00233 IF (HISBA=='DIF'.OR.(PTI_MAX(I)-PTI_MIN(I))<0.2) THEN 00234 ZTI_MIN =XREGP*PTI_MIN (I)-XREGA 00235 ZTI_MAX =XREGP*PTI_MAX (I)-XREGA 00236 ZTI_STD =SQRT((XREGP**2)*(PTI_STD(I)**2)) 00237 ELSE 00238 ZTI_MIN =PTI_MIN (I) 00239 ZTI_MAX =PTI_MAX (I) 00240 ZTI_STD =PTI_STD (I) 00241 ENDIF 00242 ELSE 00243 ZTI_MEAN=PTI_MEAN(I) 00244 ZTI_MIN =PTI_MIN (I) 00245 ZTI_MAX =PTI_MAX (I) 00246 ZTI_STD =PTI_STD (I) 00247 ENDIF 00248 ! 00249 ! Calculate topographic index pdf parameters 00250 ! 00251 ! Numerical problem especialy over Greenland 00252 IF(PTI_SKEW(I)<=0.2)THEN 00253 ! 00254 ZTI_SKEW=0.2 00255 ! 00256 WRITE(KLUOUT,*)'TI_SKEW is too low or negatif (=',PTI_SKEW(I),'),' 00257 WRITE(KLUOUT,*)'then PHI is too big for the grid-cell',I,'So,GAMMA(PHI) -> +inf.' 00258 WRITE(KLUOUT,*)'The applied solution is to put TI_SKEW = 0.2' 00259 IF(ZTI_STD<1.0)THEN 00260 WRITE(KLUOUT,*)'In addition TI_STD is too low (=',ZTI_STD,'),' 00261 WRITE(KLUOUT,*)'The applied solution is to put TI_STD = 1.0' 00262 ZTI_STD=1.0 00263 ENDIF 00264 ! 00265 ZAR = ZAR +1.0 00266 ! 00267 ZXI = ZTI_SKEW*ZTI_STD/X2 00268 ZPHI = (ZTI_STD/ZXI)**X2 00269 ! 00270 ELSE 00271 ! 00272 ZXI = PTI_SKEW(I)*ZTI_STD/X2 00273 ZPHI = (ZTI_STD/ZXI)**X2 00274 ! 00275 ENDIF 00276 ! 00277 ZNU = ZTI_MEAN-ZPHI*ZXI 00278 ! 00279 ! Exponential decay factor of the local deficit 00280 ! 00281 PM(I) =(ZWSAT_AVG(I)-ZWWILT_AVG(I))*ZD_TOP(I)/X4 00282 ! 00283 ! 1.2.1 Calculate grid cell pdf total density FTOT = F(ymin --> ymax) 00284 ! ------------------------------------------------------------------- 00285 ! 00286 ! Normalized TOPMODEL maximum deficit D0/M coefficient 00287 ! 00288 ZD0 = (ZWSAT_AVG(I)-ZWWILT_AVG(I))*ZD_TOP(I)/PM(I) 00289 ! 00290 ! Initialise 00291 ! 00292 ZGYMAX = 0.0 00293 ZGYMIN = 0.0 00294 ZFTOT = 0.0 00295 ZYMIN = 0.0 00296 ZYMAX = 0.0 00297 ! 00298 ! variable changing yi ---> (ti-nu)/xi 00299 ! 00300 ZYMIN = (ZTI_MIN-ZNU)/ZXI 00301 ZYMAX = (ZTI_MAX-ZNU)/ZXI 00302 ! 00303 ! Supress numerical artifact 00304 ! 00305 ZCOR = ABS(MIN(0.0,ZYMIN)) 00306 ! 00307 ZYMIN = MAX(0.0,ZYMIN+ZCOR) 00308 ZYMAX = ZYMAX+ZCOR 00309 ! 00310 ! Errors flags indicating a number of error condition in G and GSTAR 00311 ! (see gamma_inc.f for more detail) 00312 ! 00313 IFLG =0 00314 IFLGST =0 00315 ! 00316 ! Computation of F(0 --> ymin) 00317 ! 00318 CALL DGAM(ZPHI,ZYMIN,10.,ZG,ZGYMIN,IFLG,IFLGST) 00319 ! 00320 ! if the incomplete gamma function don't work, print why 00321 ! 00322 IF (IFLGST/=0)THEN 00323 WRITE(KLUOUT,*)'GRID-CELL =',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'YMIN= ',ZYMIN 00324 CALL ABOR1_SFX('INIT_TOP: (1) PROBLEM WITH DGAM FUNCTION') 00325 ENDIF 00326 ! 00327 ! Computation of F(0 --> ymax) 00328 ! 00329 CALL DGAM(ZPHI,ZYMAX,10.,ZG,ZGYMAX,IFLG,IFLGST) 00330 ! 00331 ! if the incomplete gamma function don't work, print why 00332 ! 00333 IF (IFLGST/=0)THEN 00334 WRITE(KLUOUT,*)'GRID-CELL =',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'YMAX= ',ZYMAX 00335 CALL ABOR1_SFX('INIT_TOP: (2) PROBLEM WITH DGAM FUNCTION') 00336 ENDIF 00337 ! 00338 ! FTOT = F(0 --> ymax) - F(0 --> ymin) 00339 ! 00340 ZFTOT=ZGYMAX-ZGYMIN 00341 ! 00342 ! initialization water content and fraction 00343 ! 00344 PTAB_WTOP(I,1) = ZWSAT_AVG(I) 00345 PTAB_FSAT(I,1) = 1.0 00346 ! 00347 PTAB_WTOP(I,IPAS) = ZWWILT_AVG(I) 00348 PTAB_FSAT(I,IPAS) = 0.0 00349 ! 00350 ! Define the new limits for the satured index loop 00351 ! 00352 JSI_MIN = 2 00353 JSI_MAX = IPAS-1 00354 ZPAS = (ZTI_MAX-ZTI_MIN)/REAL(IPAS-1) 00355 ! 00356 ! 1.2.2 Calculate all topmodel arrays 00357 ! ----------------------------------- 00358 ! 00359 ! Satured index loop 00360 ! 00361 DO IND=JSI_MIN,JSI_MAX 00362 ! 00363 ! initialize of loops variables 00364 ! 00365 ZXSAT_IND = 0.0 00366 ZYSAT = 0.0 00367 ZY0 = 0.0 00368 ZDMOY = 0.0 00369 ZXMOY = 0.0 00370 ZFMED = 0.0 00371 ! 00372 ! Initialize of incomplete gamma function flags and variables 00373 ! 00374 IFLG = 0 00375 IFLGST = 0 00376 ZGYSAT = 0.0 00377 ZGY0 = 0.0 00378 ! 00379 ! calculate xsat for all new index 00380 ! 00381 ZXSAT_IND=ZTI_MIN+REAL(IND-1)*ZPAS 00382 ! 00383 ! Changing variable to compute incomplete gamma function 00384 ! 00385 ZYSAT=(ZXSAT_IND-ZNU)/ZXI 00386 ZY0 =((ZXSAT_IND-ZD0)-ZNU)/ZXI 00387 ! 00388 ! Calculate Y0 and ysat and assume ymin < y0 < ymax ! 00389 00390 ZYSAT=MAX(ZYMIN,MIN(ZYSAT+ZCOR,ZYMAX)) 00391 ZY0 =MAX(ZYMIN,MIN(ZY0 +ZCOR,ZYMAX)) 00392 ! 00393 ! call incomplete gamma function for xsat 00394 ! 00395 CALL DGAM(ZPHI,ZYSAT,10.,ZG,ZGYSAT,IFLG,IFLGST) 00396 ! 00397 ! if the incomplete gamma function don't works, print why 00398 ! 00399 IF (IFLGST/=0)THEN 00400 WRITE(KLUOUT,*)'GRID-CELL= ',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'YSAT= ',ZYSAT 00401 CALL ABOR1_SFX('INIT_TOP: (3) PROBLEM WITH DGAM FUNCTION') 00402 ENDIF 00403 ! 00404 ! call incomplete gamma function for xsat 00405 ! 00406 CALL DGAM(ZPHI,ZY0,10.,ZG,ZGY0,IFLG,IFLGST) 00407 ! 00408 ! if the incomplete gamma function don't works, print why 00409 ! 00410 IF (IFLGST/=0)THEN 00411 WRITE(KLUOUT,*)'GRID-CELL= ',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'Y0= ',ZY0 00412 CALL ABOR1_SFX('INIT_TOP: (4) PROBLEM WITH DGAM FUNCTION') 00413 ENDIF 00414 ! 00415 ! compute satured fraction as FSAT = F(0 --> ymax) - F(0 --> ysat) 00416 ! 00417 PTAB_FSAT(I,IND)=MAX(0.0,(ZGYMAX-ZGYSAT)/ZFTOT) 00418 ! 00419 ! Compute driest fraction 00420 ! 00421 ZF0=MAX(0.0,(ZGY0-ZGYMIN)/ZFTOT) 00422 ! 00423 ! Calculate FMED 00424 ! 00425 ZFMED=(1.0-PTAB_FSAT(I,IND)-ZF0) 00426 ! 00427 IF (ZFMED/=0.0) THEN 00428 ! 00429 ! Compute the new x mean, xmoy', over the wet fraction Fwet 00430 ! 00431 ZXMOY = ZNU+ZXI*(ZPHI-ZCOR+(EXP(-ZY0)*(ZY0**(ZPHI/X2))*(ZY0**(ZPHI/X2)) & 00432 - EXP(-ZYSAT)*(ZYSAT**(ZPHI/X2))*(ZYSAT**(ZPHI/X2)))/(ZFMED*GAMMAS(ZPHI))) 00433 ! 00434 ! supress numerical artifacs 00435 ! 00436 ZXMOY =MAX((ZXSAT_IND-ZD0),MIN(ZXSAT_IND,ZXMOY)) 00437 ! 00438 ! Calculate the mean normalysed deficit as Dbar/M = (1-fsat-f0)*(xsat-xmoy')+f0*D0/M 00439 ! 00440 ZDMOY = ZFMED*(ZXSAT_IND-ZXMOY)+ZF0*ZD0 00441 ! 00442 ENDIF 00443 ! 00444 ! supress numerical artifacs 00445 ! 00446 ZDMOY = MAX(0.0,MIN(ZDMOY,ZD0)) 00447 ! 00448 ! Solves Dbar = (Wsat-WT)*d_top with Dbar/M (=ZDMOY) = (Wsat-WT)*d_top/M 00449 ! 00450 PTAB_WTOP(I,IND) = ZWSAT_AVG(I)-(PM(I)*ZDMOY/ZD_TOP(I)) 00451 ! 00452 ENDDO 00453 ! 00454 ENDIF 00455 ! 00456 ENDDO 00457 ! 00458 ! supress numerical artifacs for boundaries conditions 00459 ! 00460 DO I=1,INI 00461 ! 00462 ! Upper boundary 00463 ! 00464 IF(PTAB_WTOP(I,2)==ZWSAT_AVG(I))THEN 00465 ! 00466 ZFUP=PTAB_FSAT(I,1) 00467 ZWUP=PTAB_WTOP(I,1) 00468 ! 00469 ID(:)=MAXLOC(PTAB_WTOP(I,:),PTAB_WTOP(I,:)<ZWSAT_AVG(I)) 00470 ! 00471 ZFDOWN=PTAB_FSAT(I,ID(1)) 00472 ZWDOWN=PTAB_WTOP(I,ID(1)) 00473 ! 00474 ZSLOPEW=(ZWUP-ZWDOWN)/(ZFUP-ZFDOWN) 00475 ! 00476 DO IND=2,ID(1)-1 00477 PTAB_WTOP(I,IND)=ZWDOWN+(PTAB_FSAT(I,IND)-ZFDOWN)*ZSLOPEW 00478 ENDDO 00479 ! 00480 ENDIF 00481 ! 00482 ! Lower boundary 00483 ! 00484 WHERE(PTAB_FSAT(I,:)<=0.0 ) 00485 PTAB_WTOP(I,:)=ZWWILT_AVG(I) 00486 ENDWHERE 00487 WHERE(PTAB_WTOP(I,:)<=ZWWILT_AVG(I)) 00488 PTAB_FSAT(I,:)=0.0 00489 ENDWHERE 00490 ! 00491 ENDDO 00492 ! 00493 WRITE(KLUOUT,*)'-------------------TOPMODEL SUM-UP-------------------------' 00494 WRITE(KLUOUT,*)'Number of grid-cells ',INI,'Number of Topmodel points',INT(ZTOT) 00495 IF(INI/=0) THEN 00496 WRITE(KLUOUT,*)'Percentage of non TOPMODEL grid-cells',(100.*ZNO/FLOAT(INI)) 00497 ENDIF 00498 IF(ZTOT>0.0)THEN 00499 WRITE(KLUOUT,*)'Percentage of arranged (TI-SKE=0.2) grid-cells',(100.*ZAR/ZTOT) 00500 ENDIF 00501 WRITE(KLUOUT,*)'-----------------------------------------------------------' 00502 ! 00503 IF (LHOOK) CALL DR_HOOK('INIT_TOP',1,ZHOOK_HANDLE) 00504 ! 00505 !------------------------------------------------------------------------------- 00506 ! 00507 END SUBROUTINE INIT_TOP
1.8.0
