| 1 | ! ================================================================================================================================= |
|---|
| 2 | ! MODULE : init_top |
|---|
| 3 | ! |
|---|
| 4 | ! CONTACT : orchidee-help _at_ ipsl.jussieu.fr |
|---|
| 5 | ! |
|---|
| 6 | ! LICENCE : IPSL (2006) |
|---|
| 7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
|---|
| 8 | ! |
|---|
| 9 | !>\BRIEF This module computes the parameters for TOPMODEL. |
|---|
| 10 | !! |
|---|
| 11 | !!\n DESCRIPTION : contains. |
|---|
| 12 | !! |
|---|
| 13 | !! RECENT CHANGE(S) : None |
|---|
| 14 | !! |
|---|
| 15 | !! REFERENCE(S) : None |
|---|
| 16 | !! |
|---|
| 17 | !! SVN : |
|---|
| 18 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/shushi.peng/ORCHIDEE/src_sechiba/init_top.f90 $ |
|---|
| 19 | !! $Date: 2015-03-10 10:16:04 +0100 (Tue, 10 Mar 2015) $ |
|---|
| 20 | !! $Revision: 2535 $ |
|---|
| 21 | !! \n |
|---|
| 22 | !_ ================================================================================================================================ |
|---|
| 23 | |
|---|
| 24 | !----------------------------------------------------------------- |
|---|
| 25 | !--------------- special set of characters for SCCS information |
|---|
| 26 | !----------------------------------------------------------------- |
|---|
| 27 | ! %Z% Lib:%F%, Version:%I%, Date:%D%, Last modified:%E% |
|---|
| 28 | !----------------------------------------------------------------- |
|---|
| 29 | ! ###################### |
|---|
| 30 | MODULE init_top |
|---|
| 31 | ! ###################### |
|---|
| 32 | ! |
|---|
| 33 | ! |
|---|
| 34 | ! |
|---|
| 35 | USE ioipsl |
|---|
| 36 | USE constantes |
|---|
| 37 | ! USE constantes_soil |
|---|
| 38 | USE pft_parameters |
|---|
| 39 | !pss:! USE constantes_veg |
|---|
| 40 | !! USE reqdprec |
|---|
| 41 | USE gammad_inc |
|---|
| 42 | |
|---|
| 43 | IMPLICIT NONE |
|---|
| 44 | INTEGER, PARAMETER :: IDOUBLE=KIND(1.D0) ! Double precision |
|---|
| 45 | CONTAINS |
|---|
| 46 | ! |
|---|
| 47 | ! ######spl |
|---|
| 48 | SUBROUTINE init_top_main(kjpindex, lalo, veget_max, PWWILT, PWSAT, PD_TOP, PM, PTI_MIN, PTI_MAX, & |
|---|
| 49 | PTI_MOY, PTI_STD, PTI_SKEW, PTAB_FSAT, PTAB_WTOP, PTAB_FWET, PTAB_WTOP_WET, ZPAS) |
|---|
| 50 | ! |
|---|
| 51 | ! |
|---|
| 52 | ! ##################################################################### |
|---|
| 53 | ! |
|---|
| 54 | !!**** *INIT_TOP* |
|---|
| 55 | !! |
|---|
| 56 | !! PURPOSE |
|---|
| 57 | !! ======= |
|---|
| 58 | ! |
|---|
| 59 | ! Calculates the new array of the Datin-Saulnier TOPMODEL framework fonction for xsat and compute each |
|---|
| 60 | ! satured fraction for each xsat value of the grids cells but also the active TOPMODEL-layer array, |
|---|
| 61 | ! the normalized mean deficit array. |
|---|
| 62 | ! For calculate new array, we use the incomplete gamma function. (see gamma_inc.f for more detail) |
|---|
| 63 | ! |
|---|
| 64 | ! |
|---|
| 65 | !------------------------------------------------------------------------------- |
|---|
| 66 | ! |
|---|
| 67 | !* 0. DECLARATIONS |
|---|
| 68 | ! =================== |
|---|
| 69 | ! |
|---|
| 70 | !* 0.1 declarations of arguments |
|---|
| 71 | ! |
|---|
| 72 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size |
|---|
| 73 | REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees) |
|---|
| 74 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI -> infty) |
|---|
| 75 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PWWILT |
|---|
| 76 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PWSAT |
|---|
| 77 | ! PWWILT = the wilting point volumetric |
|---|
| 78 | ! water content (m3 m-3) |
|---|
| 79 | ! PWSAT = saturation volumetric water content |
|---|
| 80 | ! of the soil (m3 m-3) |
|---|
| 81 | ! |
|---|
| 82 | |
|---|
| 83 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PTI_MIN |
|---|
| 84 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PTI_MAX |
|---|
| 85 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PTI_STD |
|---|
| 86 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PTI_SKEW |
|---|
| 87 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PM |
|---|
| 88 | REAL(r_std), INTENT (in) :: PD_TOP |
|---|
| 89 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: PTI_MOY |
|---|
| 90 | ! PTI_MOY = ti mean |
|---|
| 91 | ! PTI_MIN = ti min |
|---|
| 92 | ! PTI_MAX = ti max |
|---|
| 93 | ! PTI_STD = ti standard deviation |
|---|
| 94 | ! PTI_SKEW = ti skewness |
|---|
| 95 | ! PM = exponential decay factor of the local deficit |
|---|
| 96 | ! PD_TOP = Topmodel active soil depth |
|---|
| 97 | ! |
|---|
| 98 | REAL(r_std), DIMENSION(kjpindex,1000), INTENT(OUT) :: PTAB_FSAT, PTAB_WTOP,PTAB_FWET, PTAB_WTOP_WET |
|---|
| 99 | REAL(r_std), DIMENSION(kjpindex), INTENT (OUT) :: ZPAS |
|---|
| 100 | ! PTAB_FWET = Wetland fraction array |
|---|
| 101 | ! PTAB_FSAT = Satured fraction array |
|---|
| 102 | ! PTAB_WTOP = Active TOPMODEL-layer array |
|---|
| 103 | ! |
|---|
| 104 | !* 0.2 declarations of local variables |
|---|
| 105 | ! |
|---|
| 106 | REAL(KIND=IDOUBLE), DIMENSION(kjpindex) :: ZXSAT |
|---|
| 107 | ! ZXSAT = Initial satured index for all index |
|---|
| 108 | |
|---|
| 109 | |
|---|
| 110 | |
|---|
| 111 | ! |
|---|
| 112 | REAL(KIND=IDOUBLE) :: ZXI, ZPHI, ZNU, ZTI_MOY, ZTI_MAX, ZTI_MIN |
|---|
| 113 | ! ZTI_MOY = ti mean after regression |
|---|
| 114 | ! ZXI = ti pdf parameter |
|---|
| 115 | ! ZPHI = ti pdf parameter |
|---|
| 116 | ! ZNU = ti pdf parameter |
|---|
| 117 | ! |
|---|
| 118 | REAL(KIND=IDOUBLE) :: ZFTOT, ZYMAX, ZYMIN |
|---|
| 119 | REAL(KIND=IDOUBLE) :: ZGYMAX, ZGYMIN |
|---|
| 120 | ! ZFTOT = total fraction of a grid cell |
|---|
| 121 | ! ZYMAX = yi maximum variable |
|---|
| 122 | ! ZYMIN = yi minimum variable |
|---|
| 123 | ! ZGYMAX = incomplete gamma function for ymax (GAMSTAR result) |
|---|
| 124 | ! ZGYMIN = incomplete gamma function for ymin (GAMSTAR result) |
|---|
| 125 | ! |
|---|
| 126 | REAL(KIND=IDOUBLE) :: ZXSAT_IND, ZYSAT, ZY0, ZDMOY, ZXMOY, ZFMED, ZF0 |
|---|
| 127 | ! ZXSAT_IND = Satured index for all index |
|---|
| 128 | ! ZYSAT = changing variable of satured index |
|---|
| 129 | ! ZY0 = changing variable of dry index |
|---|
| 130 | ! ZDMOY = grid cell average deficit (= Dbar/M) |
|---|
| 131 | ! ZXMOY = ti mean value on wet (1-fsat-f0) fraction |
|---|
| 132 | ! ZF0 = dry fraction |
|---|
| 133 | ! ZFMED = wet (1-fsat-f0) fraction |
|---|
| 134 | ! |
|---|
| 135 | REAL(KIND=IDOUBLE) :: ZG, ZGYSAT, ZGY0 |
|---|
| 136 | ! ZG = GAM result |
|---|
| 137 | ! ZGYSAT = the incomplete gamma function for ysat (GAMSTAR result) |
|---|
| 138 | ! ZGY0 = the incomplete gamma function for y0 (GAMSTAR result) |
|---|
| 139 | ! |
|---|
| 140 | REAL(KIND=IDOUBLE) :: ZD0, ZGAM_PHI |
|---|
| 141 | ! ZD0 = Normalized TOPMODEL maximum deficit D0/M coefficient |
|---|
| 142 | ! ZPAS = pas for calculate the new xsat values array |
|---|
| 143 | ! |
|---|
| 144 | ! |
|---|
| 145 | INTEGER :: IFLG, IFLGST |
|---|
| 146 | ! IFLG = incomplete gamma function error flag (GAM result) |
|---|
| 147 | ! IFLGST = incomplete gamma function error flag (GAMSTAR result) |
|---|
| 148 | ! (see gamma_inc.f for more detail) |
|---|
| 149 | ! |
|---|
| 150 | INTEGER :: INI, I, IND, JSI_MIN, JSI_MAX, IPAS, qq, I_zti_moy |
|---|
| 151 | ! |
|---|
| 152 | REAL(r_std),DIMENSION(kjpindex) :: sum_veget |
|---|
| 153 | |
|---|
| 154 | !------------------------------------------------------------------------------- |
|---|
| 155 | ! |
|---|
| 156 | ! 1 TOPMODEL SURFACE RUNOFF SCHEME |
|---|
| 157 | ! ================================ |
|---|
| 158 | ! |
|---|
| 159 | ! 1.1 initialisation of local variable |
|---|
| 160 | ! ---------------------------------------- |
|---|
| 161 | ! |
|---|
| 162 | |
|---|
| 163 | sum_veget=SUM(veget_max(:,2:13),2) |
|---|
| 164 | |
|---|
| 165 | |
|---|
| 166 | ! Grid cells number |
|---|
| 167 | ! |
|---|
| 168 | |
|---|
| 169 | INI = kjpindex |
|---|
| 170 | IPAS = SIZE(PTAB_FSAT,2) |
|---|
| 171 | ! |
|---|
| 172 | ! GAM result (not use here !) |
|---|
| 173 | ! |
|---|
| 174 | ZG = 0.0 |
|---|
| 175 | ! |
|---|
| 176 | ! ti_sat values |
|---|
| 177 | ! |
|---|
| 178 | ZXSAT(:) = PTI_MIN(:)-ZPAS(:) |
|---|
| 179 | |
|---|
| 180 | ! |
|---|
| 181 | ! 1.2 Algorithme |
|---|
| 182 | ! ------------------ |
|---|
| 183 | ! |
|---|
| 184 | !grid cells loops |
|---|
| 185 | ! |
|---|
| 186 | DO I=1,INI |
|---|
| 187 | |
|---|
| 188 | !la boucle suivante permet de calculer la fraction a wfc (1er pas de boucle) et la fraction wetland (2nd pas de temps de boucle) |
|---|
| 189 | DO I_zti_moy=1,2 |
|---|
| 190 | |
|---|
| 191 | ! |
|---|
| 192 | ! IF(PTI_SKEW(I)/=-99.99)THEN |
|---|
| 193 | IF((PTI_SKEW(I)/=-99.99) .AND. (PTI_STD(I)/=-99.99) .AND. (PM(I) .GE. 0.01))THEN |
|---|
| 194 | ! write(*,*) I, PTI_SKEW(I) |
|---|
| 195 | ! |
|---|
| 196 | ! 1.2.0 first initialisation |
|---|
| 197 | ! -------------------------- |
|---|
| 198 | ! |
|---|
| 199 | ZXI = 0.0 |
|---|
| 200 | ZPHI = 0.0 |
|---|
| 201 | ZNU = 0.0 |
|---|
| 202 | ! |
|---|
| 203 | ! Wolock and McCabe (2000) linear regression equation between the mean |
|---|
| 204 | ! topographic index computed with a 1000 meter DEM and a 100 meter DEM. |
|---|
| 205 | ! |
|---|
| 206 | ! |
|---|
| 207 | ! ZTI_MOY=0.961*PTI_MOY(I)-1.957 |
|---|
| 208 | |
|---|
| 209 | |
|---|
| 210 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 211 | ZTI_MOY=PTI_MOY(I) |
|---|
| 212 | ZTI_MIN=PTI_MIN(I) |
|---|
| 213 | ZTI_MAX=PTI_MAX(I) |
|---|
| 214 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 215 | ZTI_MOY=PTI_MOY(I)+ SHIFT_fsat_fwet !shift de la distribution d indice topo |
|---|
| 216 | ZTI_MIN=PTI_MIN(I)+ SHIFT_fsat_fwet |
|---|
| 217 | ZTI_MAX=PTI_MAX(I)+ SHIFT_fsat_fwet |
|---|
| 218 | !! non necessaire de modifier le skewness et la std: SKWENESS(X+a)=SKWNESS(X) et STD(X+a)=STD(X) |
|---|
| 219 | ENDIF |
|---|
| 220 | |
|---|
| 221 | |
|---|
| 222 | ! Calculate topographic index pdf parameters |
|---|
| 223 | ! |
|---|
| 224 | |
|---|
| 225 | ZXI = PTI_SKEW(I)*PTI_STD(I)/2. |
|---|
| 226 | if (ZXI .NE. 0.0) then |
|---|
| 227 | ZPHI = (PTI_STD(I)/ZXI)**2 |
|---|
| 228 | ZNU = PTI_MOY(I)-ZPHI*ZXI |
|---|
| 229 | |
|---|
| 230 | !vire deserts |
|---|
| 231 | if ((ZPHI .LT. 100.0) .AND. ((PWSAT(I)-PWWILT(I)) .GT. min_sechiba) & |
|---|
| 232 | & .AND. (PWSAT(I) .GT. min_sechiba) .AND. (sum_veget(I) .GT. 0.01)) then |
|---|
| 233 | |
|---|
| 234 | ZD0 = (PWSAT(I)-PWWILT(I))*PD_TOP/PM(I) |
|---|
| 235 | |
|---|
| 236 | ! Initialise |
|---|
| 237 | ! |
|---|
| 238 | ZGYMAX = 0.0 |
|---|
| 239 | ZGYMIN = 0.0 |
|---|
| 240 | ZFTOT = 0.0 |
|---|
| 241 | ZYMIN = 0.0 |
|---|
| 242 | ZYMAX = 0.0 |
|---|
| 243 | ! |
|---|
| 244 | ! variable changing yi ---> (ti-nu)/xi |
|---|
| 245 | ! |
|---|
| 246 | ZYMIN = (ZTI_MIN-ZNU)/ZXI |
|---|
| 247 | ZYMAX = (ZTI_MAX-ZNU)/ZXI |
|---|
| 248 | ! |
|---|
| 249 | ! Supress numerical artifact |
|---|
| 250 | ! |
|---|
| 251 | ZYMIN = MAX(0.0,ZYMIN) |
|---|
| 252 | ! |
|---|
| 253 | ! Errors flags indicating a number of error condition in G and GSTAR |
|---|
| 254 | ! (see gamma_inc.f for more detail) |
|---|
| 255 | ! |
|---|
| 256 | IFLG =0 |
|---|
| 257 | IFLGST =0 |
|---|
| 258 | ! |
|---|
| 259 | ! Computation of F(0 --> ymin) |
|---|
| 260 | ! |
|---|
| 261 | CALL DGAM(ZPHI,ZYMIN,2.,ZG,ZGYMIN,IFLG,IFLGST) |
|---|
| 262 | ! |
|---|
| 263 | ! if the incomplete gamma function don't work, print why |
|---|
| 264 | ! |
|---|
| 265 | IF (IFLGST/=0) PRINT*,'MAILLE =',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'YMIN= ',ZYMIN |
|---|
| 266 | ! |
|---|
| 267 | ! Computation of F(0 --> ymax) |
|---|
| 268 | ! |
|---|
| 269 | CALL DGAM(ZPHI,ZYMAX,2.,ZG,ZGYMAX,IFLG,IFLGST) |
|---|
| 270 | ! |
|---|
| 271 | ! if the incomplete gamma function don't work, print why |
|---|
| 272 | ! |
|---|
| 273 | IF (IFLGST/=0) PRINT*,'MAILLE =',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'YMAX= ',ZYMAX |
|---|
| 274 | ! |
|---|
| 275 | ! FTOT = F(0 --> ymax) - F(0 --> ymin) |
|---|
| 276 | ! |
|---|
| 277 | ZFTOT=ZGYMAX-ZGYMIN |
|---|
| 278 | if (ZFTOT .NE. 0.0) then |
|---|
| 279 | ! |
|---|
| 280 | ! initialization of loop control variables |
|---|
| 281 | ! |
|---|
| 282 | JSI_MAX = 0 |
|---|
| 283 | ! |
|---|
| 284 | ! Define the new limits for the satured index loop |
|---|
| 285 | ! |
|---|
| 286 | JSI_MIN = 1 |
|---|
| 287 | JSI_MAX = IPAS |
|---|
| 288 | ZPAS(I) = (ZTI_MAX-ZTI_MIN)/(IPAS-1) |
|---|
| 289 | ! |
|---|
| 290 | ! 1.2.2 Calculate all topmodel arrays |
|---|
| 291 | ! ----------------------------------- |
|---|
| 292 | ! |
|---|
| 293 | ! Satured index loop |
|---|
| 294 | ! |
|---|
| 295 | DO IND=JSI_MIN,JSI_MAX |
|---|
| 296 | ! |
|---|
| 297 | ! initialize of loops variables |
|---|
| 298 | ! |
|---|
| 299 | ZXSAT_IND = 0.0 |
|---|
| 300 | ZYSAT = 0.0 |
|---|
| 301 | ZY0 = 0.0 |
|---|
| 302 | ZDMOY = 0.0 |
|---|
| 303 | ZXMOY = 0.0 |
|---|
| 304 | ZFMED = 0.0 |
|---|
| 305 | ! |
|---|
| 306 | ! Initialize of incomplete gamma function flags and variables |
|---|
| 307 | ! |
|---|
| 308 | IFLG = 0 |
|---|
| 309 | IFLGST = 0 |
|---|
| 310 | ZGYSAT = 0.0 |
|---|
| 311 | ZGY0 = 0.0 |
|---|
| 312 | ! |
|---|
| 313 | ! calculate xsat for all new index |
|---|
| 314 | ! |
|---|
| 315 | ZXSAT_IND=ZTI_MIN+(IND-1)*ZPAS(I) |
|---|
| 316 | ! |
|---|
| 317 | ! Changing variable to compute incomplete gamma function |
|---|
| 318 | ! |
|---|
| 319 | ZYSAT=(ZXSAT_IND-ZNU)/ZXI |
|---|
| 320 | ZY0 =((ZXSAT_IND-ZD0)-ZNU)/ZXI |
|---|
| 321 | ! |
|---|
| 322 | ! Calculate Y0 and ysat and assume ymin < y0 < ymax ! |
|---|
| 323 | |
|---|
| 324 | ZYSAT=MAX(ZYMIN,MIN(ZYSAT,ZYMAX)) |
|---|
| 325 | ZY0 =MAX(ZYMIN,MIN(ZY0,ZYMAX)) |
|---|
| 326 | ! |
|---|
| 327 | ! call incomplete gamma function for xsat |
|---|
| 328 | ! |
|---|
| 329 | CALL DGAM(ZPHI,ZYSAT,2.,ZG,ZGYSAT,IFLG,IFLGST) |
|---|
| 330 | ! |
|---|
| 331 | ! if the incomplete gamma function don't works, print why |
|---|
| 332 | ! |
|---|
| 333 | IF (IFLGST/=0) print*,'MAILLE= ',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'YSAT= ',ZYSAT |
|---|
| 334 | ! |
|---|
| 335 | ! call incomplete gamma function for xsat |
|---|
| 336 | ! |
|---|
| 337 | CALL DGAM(ZPHI,ZY0,2.,ZG,ZGY0,IFLG,IFLGST) |
|---|
| 338 | ! |
|---|
| 339 | ! if the incomplete gamma function don't works, print why |
|---|
| 340 | ! |
|---|
| 341 | IF (IFLGST/=0) print*,'MAILLE= ',I,'FLGST= ',IFLGST,'PHI= ',ZPHI,'Y0= ',ZY0 |
|---|
| 342 | ! |
|---|
| 343 | ! compute satured fraction as FSAT = F(0 --> ymax) - F(0 --> ysat) |
|---|
| 344 | ! |
|---|
| 345 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 346 | PTAB_FSAT(I,IND)=(ZGYMAX-ZGYSAT)/ZFTOT |
|---|
| 347 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 348 | PTAB_FWET(I,IND)=(ZGYMAX-ZGYSAT)/ZFTOT |
|---|
| 349 | ENDIF |
|---|
| 350 | ! |
|---|
| 351 | ! Calculate FMED = F(0 --> ysat) - F(0 --> y0) |
|---|
| 352 | ! |
|---|
| 353 | ZFMED=MAX(0.0,((ZGYSAT-ZGY0)/ZFTOT)) |
|---|
| 354 | ! |
|---|
| 355 | ! compute driest fraction as F0 = 1-Fsat-Fwet |
|---|
| 356 | ! |
|---|
| 357 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 358 | ZF0=MAX(0.0,(1.0-PTAB_FSAT(I,IND)-ZFMED)) |
|---|
| 359 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 360 | ZF0=MAX(0.0,(1.0-PTAB_FWET(I,IND)-ZFMED)) |
|---|
| 361 | ENDIF |
|---|
| 362 | |
|---|
| 363 | ! |
|---|
| 364 | IF (ZFMED/=0.0) THEN |
|---|
| 365 | ! |
|---|
| 366 | ! Compute the new x mean, xmoy', over the wet fraction Fwet |
|---|
| 367 | ! |
|---|
| 368 | CALL gamma(ZPHI, ZGAM_PHI) |
|---|
| 369 | |
|---|
| 370 | ZXMOY = ZNU+ZXI*(ZPHI+(EXP(-ZY0)*(ZY0**(ZPHI/2))*(ZY0**(ZPHI/2)) & |
|---|
| 371 | -EXP(-ZYSAT)*(ZYSAT**(ZPHI/2))*(ZYSAT**(ZPHI/2)))/(ZFMED*ZGAM_PHI)) |
|---|
| 372 | |
|---|
| 373 | ! supress numerical artifacs |
|---|
| 374 | ! |
|---|
| 375 | ZXMOY =MAX((ZXSAT_IND-ZD0),MIN(ZXSAT_IND,ZXMOY)) |
|---|
| 376 | ! |
|---|
| 377 | ! Calculate the mean normalysed deficit as Dbar/M = (1-fsat-f0)*(xsat-xmoy')+f0*D0/M |
|---|
| 378 | ! |
|---|
| 379 | ZDMOY = ZFMED*(ZXSAT_IND-ZXMOY)+ZF0*ZD0 |
|---|
| 380 | ! |
|---|
| 381 | ENDIF |
|---|
| 382 | ! |
|---|
| 383 | ! supress numerical artifacs |
|---|
| 384 | ! |
|---|
| 385 | ZDMOY = MAX(0.0,MIN(ZDMOY,ZD0)) |
|---|
| 386 | ! |
|---|
| 387 | ! Solves Dbar = (Wsat-WT)*d_top with Dbar/M (=ZDMOY) = (Wsat-WT)*d_top/M |
|---|
| 388 | ! |
|---|
| 389 | !modifie le calcul de PTAB_WTOP du au fait que le deficit est calcule par rapport a wfc et non par rapport a wsat |
|---|
| 390 | |
|---|
| 391 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 392 | PTAB_WTOP(I,IND) = PWSAT(I)-(PM(I)*ZDMOY/PD_TOP) |
|---|
| 393 | !pss: should not minus PWWILT, otherwise not consistent with hydro_subgrid.f90 |
|---|
| 394 | ! PTAB_WTOP(I,IND) = PTAB_WTOP(I,IND) - PWWILT(I) |
|---|
| 395 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 396 | !ne modifie que une des deux fonctions (c.a.d PTAB_FSAT/PTAB_FWET et pas PTAB_WTOP) |
|---|
| 397 | !PTAB_WTOP_WET(I,IND) = PWSAT(I)-(PM(I)*ZDMOY/PD_TOP(I)) |
|---|
| 398 | !PTAB_WTOP_WET(I,IND) = PTAB_WTOP_WET(I,IND) - PWWILT(I) |
|---|
| 399 | PTAB_WTOP_WET(I,IND) = PTAB_WTOP(I,IND) |
|---|
| 400 | ENDIF |
|---|
| 401 | |
|---|
| 402 | !suprime eventuels pbs |
|---|
| 403 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 404 | IF ( PTAB_WTOP(I,IND) <= 0.0 ) PTAB_FSAT(I,IND) = 0.0 |
|---|
| 405 | !pss: normally not need to decide the over-saturation as zero saturated fraction |
|---|
| 406 | ! IF ( PTAB_WTOP(I,IND) >= PWSAT(I) ) PTAB_FSAT(I,IND) = 0.0 |
|---|
| 407 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 408 | IF ( PTAB_WTOP_WET(I,IND) <= 0.0 ) PTAB_FWET(I,IND) = 0.0 |
|---|
| 409 | !pss: normally not need to decide the over-saturation as zero saturated fraction |
|---|
| 410 | ! IF ( PTAB_WTOP_WET(I,IND) >= PWSAT(I) ) PTAB_FWET(I,IND) = 0.0 |
|---|
| 411 | ENDIF |
|---|
| 412 | |
|---|
| 413 | |
|---|
| 414 | ! |
|---|
| 415 | ENDDO |
|---|
| 416 | |
|---|
| 417 | else !if (ZFTOT .NE. 0.0) then |
|---|
| 418 | |
|---|
| 419 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 420 | |
|---|
| 421 | PTAB_WTOP(I,:) = 0.0 |
|---|
| 422 | PTAB_FSAT(I,:) = 0.0 |
|---|
| 423 | |
|---|
| 424 | |
|---|
| 425 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 426 | |
|---|
| 427 | PTAB_WTOP_WET(I,:) = 0.0 |
|---|
| 428 | PTAB_FWET(I,:) = 0.0 |
|---|
| 429 | |
|---|
| 430 | |
|---|
| 431 | ENDIF |
|---|
| 432 | |
|---|
| 433 | endif !if (ZFTOT .NE. 0.0) then |
|---|
| 434 | |
|---|
| 435 | |
|---|
| 436 | else ! if ((ZPHI .LT. 100.0) .AND. ((PWSAT(I)-PWWILT(I)) .GT. min_sechiba) .AND. (PWSAT(I) .GT. min_sechiba) .AND. (sum_veget(I) .GT. 0.02)) |
|---|
| 437 | |
|---|
| 438 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 439 | |
|---|
| 440 | PTAB_WTOP(I,:) = 0.0 |
|---|
| 441 | PTAB_FSAT(I,:) = 0.0 |
|---|
| 442 | |
|---|
| 443 | |
|---|
| 444 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 445 | |
|---|
| 446 | PTAB_WTOP_WET(I,:) = 0.0 |
|---|
| 447 | PTAB_FWET(I,:) = 0.0 |
|---|
| 448 | |
|---|
| 449 | |
|---|
| 450 | ENDIF |
|---|
| 451 | |
|---|
| 452 | !fin de la boucle if sur la condition ZPHI > 100.0 |
|---|
| 453 | endif ! if ((ZPHI .LT. 100.0) .AND.... |
|---|
| 454 | |
|---|
| 455 | else |
|---|
| 456 | |
|---|
| 457 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 458 | |
|---|
| 459 | PTAB_WTOP(I,:) = 0.0 |
|---|
| 460 | PTAB_FSAT(I,:) = 0.0 |
|---|
| 461 | |
|---|
| 462 | |
|---|
| 463 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 464 | |
|---|
| 465 | PTAB_WTOP_WET(I,:) = 0.0 |
|---|
| 466 | PTAB_FWET(I,:) = 0.0 |
|---|
| 467 | |
|---|
| 468 | |
|---|
| 469 | ENDIF |
|---|
| 470 | |
|---|
| 471 | |
|---|
| 472 | |
|---|
| 473 | |
|---|
| 474 | endif !if (ZXI .NE. 0.0) |
|---|
| 475 | |
|---|
| 476 | ELSE |
|---|
| 477 | |
|---|
| 478 | |
|---|
| 479 | IF (I_zti_moy .EQ. 1)THEN |
|---|
| 480 | |
|---|
| 481 | PTAB_WTOP(I,:) = 0.0 |
|---|
| 482 | PTAB_FSAT(I,:) = 0.0 |
|---|
| 483 | |
|---|
| 484 | |
|---|
| 485 | ELSEIF (I_zti_moy .EQ. 2)THEN |
|---|
| 486 | |
|---|
| 487 | PTAB_WTOP_WET(I,:) = 0.0 |
|---|
| 488 | PTAB_FWET(I,:) = 0.0 |
|---|
| 489 | |
|---|
| 490 | |
|---|
| 491 | ENDIF |
|---|
| 492 | |
|---|
| 493 | |
|---|
| 494 | |
|---|
| 495 | ENDIF !IF((PTI_SKEW(I)/=-99.99) .AND. (PTI_STD(I)/=-99.99) .AND. (PM(I).GE. 0.01)) |
|---|
| 496 | ! |
|---|
| 497 | ENDDO |
|---|
| 498 | |
|---|
| 499 | |
|---|
| 500 | ENDDO |
|---|
| 501 | ! |
|---|
| 502 | !------------------------------------------------------------------------------- |
|---|
| 503 | ! |
|---|
| 504 | END SUBROUTINE init_top_main |
|---|
| 505 | ! |
|---|
| 506 | !******************************************* |
|---|
| 507 | !* FUNCTION GAMMA(X) * |
|---|
| 508 | !* --------------------------------------- * |
|---|
| 509 | !* Returns the value of Gamma(x) in double * |
|---|
| 510 | !* precision as EXP(LN(GAMMA(X))) for X>0. * |
|---|
| 511 | !******************************************* |
|---|
| 512 | !!! real*8 Function gamma(xx) |
|---|
| 513 | SUBROUTINE gamma(xx,gamma_re) |
|---|
| 514 | IMPLICIT NONE |
|---|
| 515 | REAL(r_std), INTENT(IN) :: xx |
|---|
| 516 | REAL(r_std), INTENT(OUT) :: gamma_re |
|---|
| 517 | REAL cof(6),stp,half,one,fpf,x,tmp,ser |
|---|
| 518 | INTEGER j |
|---|
| 519 | DATA cof,stp /76.18009173d0,-86.50532033d0,24.01409822d0, & |
|---|
| 520 | -1.231739516d0,0.120858003d-2,-0.536382d-5,2.50662827465d0/ |
|---|
| 521 | DATA half,one,fpf /0.5d0,1.0d0,5.5d0/ |
|---|
| 522 | x=xx-one |
|---|
| 523 | ! PRINT*, 'x',x |
|---|
| 524 | tmp=x+fpf |
|---|
| 525 | ! PRINT*, 'tmp', tmp |
|---|
| 526 | tmp=(x+half)*DLOG(tmp)-tmp |
|---|
| 527 | ! PRINT*, 'tmp', tmp |
|---|
| 528 | ser=one |
|---|
| 529 | do j=1,6 |
|---|
| 530 | x=x+one |
|---|
| 531 | ser=ser+cof(j)/x |
|---|
| 532 | end do |
|---|
| 533 | ! PRINT*, 'ser', ser |
|---|
| 534 | gamma_re = DEXP(tmp+DLOG(stp*ser)) |
|---|
| 535 | ! return |
|---|
| 536 | END SUBROUTINE gamma |
|---|
| 537 | |
|---|
| 538 | END MODULE init_top |
|---|
| 539 | ! |
|---|
