PageRenderTime 68ms CodeModel.GetById 32ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_sf_gfs.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1767 lines | 1103 code | 59 blank | 605 comment | 7 complexity | 14ab5460849b4731c1df6fb8bf02c89f MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

    

  1. !!WRF:MODEL_LAYER:PHYSICS
    2. 

  2. !
    3. 

  3. MODULE module_sf_gfs
    4. 

  4. 

  5. 

  6. CONTAINS
    7. 

  7. 

  8. !-------------------------------------------------------------------
    9. 

  9.    SUBROUTINE SF_GFS(U3D,V3D,T3D,QV3D,P3D,                              &
    10. 

  10. 		 CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM,             &
    11. 

  11.                      ZNT,UST,PSIM,PSIH,                                 &
    12. 

  12. 		 XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,                    &
    13. 

  13.                      QGH,QSFC,U10,V10,                                  &
    14. 

  14.                      GZ1OZ0,WSPD,BR,ISFFLX,                             &
    15. 

  15.                      EP1,EP2,KARMAN,itimestep,                          &
    16. 

  16.                      ids,ide, jds,jde, kds,kde,                         &
    17. 

  17.                      ims,ime, jms,jme, kms,kme,                         &
    18. 

  18.                      its,ite, jts,jte, kts,kte                     )
    19. 

  19. !-------------------------------------------------------------------
    20. 

  20.       USE MODULE_GFS_MACHINE, ONLY : kind_phys
    21. 

  21.       USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
    22. 

  22. !-------------------------------------------------------------------
    23. 

  23.       IMPLICIT NONE
    24. 

  24. !-------------------------------------------------------------------
    25. 

  25. !-- U3D    	3D u-velocity interpolated to theta points (m/s)
    26. 

  26. !-- V3D    	3D v-velocity interpolated to theta points (m/s)
    27. 

  27. !-- T3D     	temperature (K)
    28. 

  28. !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
    29. 

  29. !-- P3D         3D pressure (Pa)
    30. 

  30. !-- CP		heat capacity at constant pressure for dry air (J/kg/K)
    31. 

  31. !-- ROVCP       R/CP
    32. 

  32. !-- R           gas constant for dry air (J/kg/K)
    33. 

  33. !-- XLV         latent heat of vaporization for water (J/kg)
    34. 

  34. !-- PSFC	surface pressure (Pa)
    35. 

  35. !-- ZNT		roughness length (m)
    36. 

  36. !-- UST		u* in similarity theory (m/s)
    37. 

  37. !-- PSIM        similarity stability function for momentum
    38. 

  38. !-- PSIH        similarity stability function for heat
    39. 

  39. !-- XLAND	land mask (1 for land, 2 for water)
    40. 

  40. !-- HFX		upward heat flux at the surface (W/m^2)
    41. 

  41. !-- QFX   	upward moisture flux at the surface (kg/m^2/s)
    42. 

  42. !-- LH          net upward latent heat flux at surface (W/m^2)
    43. 

  43. !-- TSK		surface temperature (K)
    44. 

  44. !-- FLHC        exchange coefficient for heat (m/s)
    45. 

  45. !-- FLQC        exchange coefficient for moisture (m/s)
    46. 

  46. !-- QGH         lowest-level saturated mixing ratio
    47. 

  47. !-- U10         diagnostic 10m u wind
    48. 

  48. !-- V10         diagnostic 10m v wind
    49. 

  49. !-- GZ1OZ0      log(z/z0) where z0 is roughness length
    50. 

  50. !-- WSPD        wind speed at lowest model level (m/s)
    51. 

  51. !-- BR          bulk Richardson number in surface layer
    52. 

  52. !-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
    53. 

  53. !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
    54. 

  54. !-- KARMAN      Von Karman constant
    55. 

  55. !-- ids         start index for i in domain
    56. 

  56. !-- ide         end index for i in domain
    57. 

  57. !-- jds         start index for j in domain
    58. 

  58. !-- jde         end index for j in domain
    59. 

  59. !-- kds         start index for k in domain
    60. 

  60. !-- kde         end index for k in domain
    61. 

  61. !-- ims         start index for i in memory
    62. 

  62. !-- ime         end index for i in memory
    63. 

  63. !-- jms         start index for j in memory
    64. 

  64. !-- jme         end index for j in memory
    65. 

  65. !-- kms         start index for k in memory
    66. 

  66. !-- kme         end index for k in memory
    67. 

  67. !-- its         start index for i in tile
    68. 

  68. !-- ite         end index for i in tile
    69. 

  69. !-- jts         start index for j in tile
    70. 

  70. !-- jte         end index for j in tile
    71. 

  71. !-- kts         start index for k in tile
    72. 

  72. !-- kte         end index for k in tile
    73. 

  73. !-------------------------------------------------------------------
    74. 

  74. 

  75.       INTEGER, INTENT(IN) ::            ids,ide, jds,jde, kds,kde,      &
    76. 

  76.                                         ims,ime, jms,jme, kms,kme,      &
    77. 

  77.                                         its,ite, jts,jte, kts,kte,      &
    78. 

  78.                                         ISFFLX,itimestep
    79. 

  79. 

  80.       REAL,    INTENT(IN) ::                                            &
    81. 

  81.                                         CP,                             &
    82. 

  82.                                         EP1,                            &
    83. 

  83.                                         EP2,                            &
    84. 

  84.                                         KARMAN,                         &
    85. 

  85.                                         R,                              &
    86. 

  86.                                         ROVCP,                          &
    87. 

  87.                                         XLV
    88. 

  88. 

  89.       REAL,    DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) ::      &
    90. 

  90.                                         P3D,                            &
    91. 

  91.                                         QV3D,                           &
    92. 

  92.                                         T3D,                            &
    93. 

  93.                                         U3D,                            &
    94. 

  94.                                         V3D
    95. 

  95. 

  96.       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(IN) ::               &
    97. 

  97.                                         TSK,                            &
    98. 

  98.                                         PSFC,                           &
    99. 

  99.                                         XLAND
    100. 

  100. 

  101.       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(INOUT) ::            &
    102. 

  102.                                         UST,                            &
    103. 

  103.                                         ZNT
    104. 

  104. 

  105.       REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) ::                 &
    106. 

  106.                                         BR,                             &
    107. 

  107.                                         CHS,                            &
    108. 

  108.                                         CHS2,                           &
    109. 

  109.                                         CPM,                            &
    110. 

  110.                                         CQS2,                           &
    111. 

  111.                                         FLHC,                           &
    112. 

  112.                                         FLQC,                           &
    113. 

  113.                                         GZ1OZ0,                         &
    114. 

  114.                                         HFX,                            &
    115. 

  115.                                         LH,                             &
    116. 

  116.                                         PSIM,                           &
    117. 

  117.                                         PSIH,                           &
    118. 

  118.                                         QFX,                            &
    119. 

  119.                                         QGH,                            &
    120. 

  120.                                         QSFC,                           &
    121. 

  121.                                         U10,                            &
    122. 

  122.                                         V10,                            &
    123. 

  123.                                         WSPD
    124. 

  124. 

  125. 

  126. !--------------------------- LOCAL VARS ------------------------------
    127. 

  127. 

  128.       REAL ::                           ESAT
    129. 

  129. 

  130.       REAL     (kind=kind_phys) ::                                      &
    131. 

  131.                                         RHOX
    132. 

  132. 

  133.       REAL     (kind=kind_phys), DIMENSION(its:ite) ::                  &
    134. 

  134.                                         CH,                             &
    135. 

  135.                                         CM,                             &
    136. 

  136.                                         DDVEL,                          &
    137. 

  137.                                         DRAIN,                          &
    138. 

  138.                                         EP,                             &
    139. 

  139.                                         EVAP,                           &
    140. 

  140.                                         FH,                             &
    141. 

  141.                                         FH2,                            &
    142. 

  142.                                         FM,                             &
    143. 

  143.                                         HFLX,                           &
    144. 

  144.                                         PH,                             &
    145. 

  145.                                         PM,                             &
    146. 

  146.                                         PRSL1,                          &
    147. 

  147.                                         PRSLKI,                         &
    148. 

  148.                                         PS,                             &
    149. 

  149.                                         Q1,                             &
    150. 

  150.                                         Q2M,                            &
    151. 

  151.                                         QSS,                            &
    152. 

  152.                                         QSURF,                          &
    153. 

  153.                                         RB,                             &
    154. 

  154.                                         RCL,                            &
    155. 

  155.                                         RHO1,                           &
    156. 

  156.                                         SLIMSK,                         &
    157. 

  157.                                         STRESS,                         &
    158. 

  158.                                         T1,                             &
    159. 

  159.                                         T2M,                            &
    160. 

  160.                                         THGB,                           &
    161. 

  161.                                         THX,                            &
    162. 

  162.                                         TSKIN,                          &
    163. 

  163.                                         SHELEG,                         &
    164. 

  164.                                         U1,                             &
    165. 

  165.                                         U10M,                           &
    166. 

  166.                                         USTAR,                          &
    167. 

  167.                                         V1,                             &
    168. 

  168.                                         V10M,                           &
    169. 

  169.                                         WIND,                           &
    170. 

  170.                                         Z0RL,                           &
    171. 

  171.                                         Z1
    172. 

  172. 

  173. 

  174.       INTEGER ::                                                        &
    175. 

  175.                                         I,                              &
    176. 

  176.                                         IM,                             &
    177. 

  177.                                         J,                              &
    178. 

  178.                                         K,                              &
    179. 

  179.                                         KM
    180. 

  180. 

  181. 

  182.    if(itimestep.eq.0) then
    183. 

  183.      CALL GFUNCPHYS
    184. 

  184.    endif
    185. 

  185. 

  186.    IM=ITE-ITS+1
    187. 

  187.    KM=KTE-KTS+1
    188. 

  188. 

  189.    DO J=jts,jte
    190. 

  190. 

  191.       DO i=its,ite
    192. 

  192.         DDVEL(I)=0.
    193. 

  193.         RCL(i)=1.
    194. 

  194.         PRSL1(i)=P3D(i,kts,j)*.001
    195. 

  195.         PS(i)=PSFC(i,j)*.001
    196. 

  196.         Q1(I) = QV3D(i,kts,j)
    197. 

  197. !        QSURF(I)=QSFC(I,J)
    198. 

  198.         QSURF(I)=0.
    199. 

  199.         SHELEG(I)=0.
    200. 

  200.         SLIMSK(i)=ABS(XLAND(i,j)-2.)
    201. 

  201.         TSKIN(i)=TSK(i,j)
    202. 

  202.         T1(I) = T3D(i,kts,j)
    203. 

  203.         U1(I) = U3D(i,kts,j)
    204. 

  204.         USTAR(I) = UST(i,j)
    205. 

  205.         V1(I) = V3D(i,kts,j)
    206. 

  206.         Z0RL(I) = ZNT(i,j)*100.
    207. 

  207.       ENDDO
    208. 

  208. 

  209.       DO i=its,ite
    210. 

  210.          PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
    211. 

  211.          THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
    212. 

  212.          THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
    213. 

  213.          RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
    214. 

  214.          Q1(I)=Q1(I)/(1.+Q1(I))
    215. 

  215.       ENDDO
    216. 

  216. 

  217. 

  218.       CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1,                                 &
    219. 

  219.                   SHELEG,TSKIN,QSURF,                                   &
    220. 

  220. !WRF              SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX,      &
    221. 

  221. !WRF              SLRAD,SNOWMT,DELT,                                    &
    222. 

  222.                   Z0RL,                                                 &
    223. 

  223. !WRF              TG3,GFLUX,F10M,                                       &
    224. 

  224.                   U10M,V10M,T2M,Q2M,                                    &
    225. 

  225. !WRF              ZSOIL,                                                &
    226. 

  226.                   CM,CH,RB,                                             &
    227. 

  227. !WRF              RHSCNPY,RHSMC,AIM,BIM,CIM,                            &
    228. 

  228.                   RCL,PRSL1,PRSLKI,SLIMSK,                              &
    229. 

  229.                   DRAIN,EVAP,HFLX,STRESS,EP,                            &
    230. 

  230.                   FM,FH,USTAR,WIND,DDVEL,                               &
    231. 

  231.                   PM,PH,FH2,QSS,Z1                                      )
    232. 

  232. 

  233. 

  234.       DO i=its,ite
    235. 

  235.         U10(i,j)=U10M(i)
    236. 

  236.         V10(i,j)=V10M(i)
    237. 

  237.         BR(i,j)=RB(i)
    238. 

  238.         CHS(I,J)=CH(I)*WIND(I)
    239. 

  239.         CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
    240. 

  240.         CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
    241. 

  241.         esat = fpvs(t1(i))
    242. 

  242.         QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
    243. 

  243.         QSFC(I,J)=qss(i)
    244. 

  244.         PSIH(i,j)=PH(i)
    245. 

  245.         PSIM(i,j)=PM(i)
    246. 

  246.         UST(i,j)=ustar(i)
    247. 

  247.         WSPD(i,j)=WIND(i)
    248. 

  248.         ZNT(i,j)=Z0RL(i)*.01
    249. 

  249.       ENDDO
    250. 

  250. 

  251.       DO i=its,ite
    252. 

  252.         FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
    253. 

  253.         FLQC(i,j)=RHO1(I)*CHS(I,J)
    254. 

  254.         GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
    255. 

  255.         CQS2(i,j)=CHS2(I,J)
    256. 

  256.       ENDDO
    257. 

  257. 

  258.       IF (ISFFLX.EQ.0) THEN
    259. 

  259.         DO i=its,ite
    260. 

  260.           HFX(i,j)=0.
    261. 

  261.           LH(i,j)=0.
    262. 

  262.           QFX(i,j)=0.
    263. 

  263.         ENDDO
    264. 

  264.       ELSE
    265. 

  265.         DO i=its,ite
    266. 

  266.           IF(XLAND(I,J)-1.5.GT.0.)THEN
    267. 

  267.             HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
    268. 

  268.           ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
    269. 

  269.             HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
    270. 

  270.             HFX(I,J)=AMAX1(HFX(I,J),-250.)
    271. 

  271.           ENDIF
    272. 

  272.           QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
    273. 

  273.           QFX(I,J)=AMAX1(QFX(I,J),0.)
    274. 

  274.           LH(I,J)=XLV*QFX(I,J)
    275. 

  275.         ENDDO
    276. 

  276.       ENDIF
    277. 

  277. 

  278. 

  279.     ENDDO
    280. 

  280. 

  281. 

  282.    END SUBROUTINE SF_GFS
    283. 

  283. 

  284. 

  285. !-------------------------------------------------------------------
    286. 

  286. 

  287.       SUBROUTINE PROGTM(IM,KM,PS,U1,V1,T1,Q1,                           &
    288. 

  288.      &                  SHELEG,TSKIN,QSURF,                             &
    289. 

  289. !WRF     &                  SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,   &
    290. 

  290. !WRF     &                  DLWFLX,SLRAD,SNOWMT,DELT,                   &
    291. 

  291.      &                  Z0RL,                                           &
    292. 

  292. !WRF     &                  TG3,GFLUX,F10M,                             &
    293. 

  293.      &                  U10M,V10M,T2M,Q2M,                              &
    294. 

  294. !WRF     &                  ZSOIL,                                      &
    295. 

  295.      &                  CM, CH, RB,                                     &
    296. 

  296. !WRF     &                  RHSCNPY,RHSMC,AIM,BIM,CIM,                  &
    297. 

  297.      &                  RCL,PRSL1,PRSLKI,SLIMSK,                        &
    298. 

  298.      &                  DRAIN,EVAP,HFLX,STRESS,EP,                      &
    299. 

  299.      &                  FM,FH,USTAR,WIND,DDVEL,                         &
    300. 

  300.      &                  PM,PH,FH2,QSS,Z1                                )
    301. 

  301. !
    302. 

  302. 

  303.       USE MODULE_GFS_MACHINE, ONLY : kind_phys
    304. 

  304.       USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
    305. 

  305.       USE MODULE_GFS_PHYSCONS, grav => con_g, SBC => con_sbc, HVAP => con_HVAP &
    306. 

  306.      &,             CP => con_CP, HFUS => con_HFUS, JCAL => con_JCAL    &
    307. 

  307.      &,             EPS => con_eps, EPSM1 => con_epsm1, t0c => con_t0c  &
    308. 

  308.      &,             RVRDM1 => con_FVirt, RD => con_RD
    309. 

  309.       implicit none
    310. 

  310. !
    311. 

  311. !     include 'constant.h'
    312. 

  312. !
    313. 

  313.       integer              IM, km
    314. 

  314. !
    315. 

  315.       real(kind=kind_phys), parameter :: cpinv=1.0/cp, HVAPI=1.0/HVAP
    316. 

  316.       real(kind=kind_phys) DELT
    317. 

  317.       INTEGER              SOILTYP(IM),  VEGTYPE(IM)
    318. 

  318.       real(kind=kind_phys) PS(IM),       U1(IM),      V1(IM),           &
    319. 

  319.      &                     T1(IM),       Q1(IM),      SHELEG(IM),       &
    320. 

  320.      &                     TSKIN(IM),    QSURF(IM),   SMC(IM,KM),       &
    321. 

  321.      &                     STC(IM,KM),   DM(IM),      SIGMAF(IM),       &
    322. 

  322.      &                     CANOPY(IM),   DLWFLX(IM),  SLRAD(IM),        &
    323. 

  323.      &                     SNOWMT(IM),   Z0RL(IM),    TG3(IM),          &
    324. 

  324.      &                     GFLUX(IM),    F10M(IM),    U10M(IM),         &
    325. 

  325.      &                     V10M(IM),     T2M(IM),     Q2M(IM),          &
    326. 

  326.      &                     ZSOIL(IM,KM), CM(IM),      CH(IM),           &
    327. 

  327.      &                     RB(IM),       RHSCNPY(IM), RHSMC(IM,KM),     &
    328. 

  328.      &                     AIM(IM,KM),   BIM(IM,KM),  CIM(IM,KM),       &
    329. 

  329.      &                     RCL(IM),      PRSL1(IM),   PRSLKI(IM),       &
    330. 

  330.      &                     SLIMSK(IM),   DRAIN(IM),   EVAP(IM),         &
    331. 

  331.      &                     HFLX(IM),     RNET(IM),    EP(IM),           &
    332. 

  332.      &                     FM(IM),       FH(IM),      USTAR(IM),        &
    333. 

  333.      &                     WIND(IM),     DDVEL(IM),   STRESS(IM)
    334. 

  334. !
    335. 

  335. !     Locals
    336. 

  336. !
    337. 

  337.       integer              k,i
    338. 

  338. !
    339. 

  339.       real(kind=kind_phys) CANFAC(IM),                                  &
    340. 

  340.      &                     DDZ(IM),     DDZ2(IM),    DELTA(IM),         &
    341. 

  341.      &                     DEW(IM),     DF1(IM),     DFT0(IM),          &
    342. 

  342.      &                     DFT2(IM),    DFT1(IM),                       &
    343. 

  343.      &                     DMDZ(IM),    DMDZ2(IM),   DTDZ1(IM),         &
    344. 

  344.      &                     DTDZ2(IM),   DTV(IM),     EC(IM),            &
    345. 

  345.      &                     EDIR(IM),    ETPFAC(IM),                     &
    346. 

  346.      &                     FACTSNW(IM), FH2(IM),     FM10(IM),          &
    347. 

  347.      &                     FX(IM),      GX(IM),                         &
    348. 

  348.      &                     HCPCT(IM),   HL1(IM),     HL12(IM),          &
    349. 

  349.      &                     HLINF(IM),   PARTLND(IM), PH(IM),            &
    350. 

  350.      &                     PH2(IM),     PM(IM),      PM10(IM),          &
    351. 

  351.      &                     PSURF(IM),   Q0(IM),      QS1(IM),           &
    352. 

  352.      &                     QSS(IM),     RAT(IM),     RCAP(IM),          &
    353. 

  353.      &                     RCH(IM),     RHO(IM),     RS(IM),            &
    354. 

  354.      &                     RSMALL(IM),  SLWD(IM),    SMCZ(IM),          &
    355. 

  355.      &                     SNET(IM),    SNOEVP(IM),  SNOWD(IM),         &
    356. 

  356.      &                     T1O(IM),     T2MO(IM),    TERM1(IM),         &
    357. 

  357.      &                     TERM2(IM),   THETA1(IM),  THV1(IM),          &
    358. 

  358.      &                     TREF(IM),    TSURF(IM),   TV1(IM),           &
    359. 

  359.      &                     TVS(IM),     TSURFO(IM),  TWILT(IM),         &
    360. 

  360.      &                     XX(IM),      XRCL(IM),    YY(IM),            &
    361. 

  361.      &                     Z0(IM),      Z0MAX(IM),   Z1(IM),            &
    362. 

  362.      &                     ZTMAX(IM),   ZZ(IM),      PS1(IM)
    363. 

  363. !
    364. 

  364.       real(kind=kind_phys) a0,    a0p,      a1,    a1p,     aa,  aa0,   &
    365. 

  365.      &                     aa1,   adtv,     alpha, arnu,    b1,  b1p,   &
    366. 

  366.      &                     b2,    b2p,      bb,    bb0,     bb1, bb2,   &
    367. 

  367.      &                     bfact, ca,       cc,    cc1,    cc2, cfactr, &
    368. 

  368.      &                     ch2o,  charnock, cice,  convrad, cq,  csoil, &
    369. 

  369.      &                     ctfil1,ctfil2,   delt2, df2,     dfsnow,     &
    370. 

  370.      &                     elocp, eth,      ff,  FMS,                   &
    371. 

  371. !WRF     &                     fhs,   funcdf,   funckt,g,      hl0, hl0inf, &
    372. 

  372.      &                     fhs,   g,        hl0,    hl0inf,             &
    373. 

  373.      &                     hl110, hlt,      hltinf,OLINF,   rcq, rcs,   &
    374. 

  374.      &                     rct,   restar,   rhoh2o,rnu,     RSI,        &
    375. 

  375.      &                     rss,   scanop,   sig2k, sigma,   smcdry,     &
    376. 

  376.      &                     t12,   t14,      tflx,  tgice,   topt,       &
    377. 

  377.      &                     val,   vis,      zbot,  snomin,  tem
    378. 

  378. !
    379. 

  379. !
    380. 

  380. 

  381.       PARAMETER (CHARNOCK=.014,CA=.4)!C CA IS THE VON KARMAN CONSTANT
    382. 

  382.       PARAMETER (G=grav,sigma=sbc)
    383. 

  383. 

  384.       PARAMETER (ALPHA=5.,A0=-3.975,A1=12.32,B1=-7.755,B2=6.041)
    385. 

  385.       PARAMETER (A0P=-7.941,A1P=24.75,B1P=-8.705,B2P=7.899,VIS=1.4E-5)
    386. 

  386.       PARAMETER (AA1=-1.076,BB1=.7045,CC1=-.05808)
    387. 

  387.       PARAMETER (BB2=-.1954,CC2=.009999)
    388. 

  388.       PARAMETER (ELOCP=HVAP/CP,DFSNOW=.31,CH2O=4.2E6,CSOIL=1.26E6)
    389. 

  389.       PARAMETER (SCANOP=.5,CFACTR=.5,ZBOT=-3.,TGICE=271.2)
    390. 

  390.       PARAMETER (CICE=1880.*917.,topt=298.)
    391. 

  391.       PARAMETER (RHOH2O=1000.,CONVRAD=JCAL*1.E4/60.)
    392. 

  392.       PARAMETER (CTFIL1=.5,CTFIL2=1.-CTFIL1)
    393. 

  393.       PARAMETER (RNU=1.51E-5,ARNU=.135*RNU)
    394. 

  394.       parameter (snomin=1.0e-9)
    395. 

  395. !
    396. 

  396.       LOGICAL FLAG(IM), FLAGSNW(IM)
    397. 

  397. !WRF      real(kind=kind_phys) KT1(IM),       KT2(IM),      KTSOIL,         &
    398. 

  398.       real(kind=kind_phys) KT1(IM),       KT2(IM),                      &
    399. 

  399.      &                     ET(IM,KM),                                   &
    400. 

  400.      &                     STSOIL(IM,KM), AI(IM,KM),    BI(IM,KM),      &
    401. 

  401.      &                     CI(IM,KM),     RHSTC(IM,KM)
    402. 

  402.       real(kind=kind_phys) rsmax(13), rgl(13),  rsmin(13), hs(13),      &
    403. 

  403.      &                     smmax(9),  smdry(9), smref(9),  smwlt(9)
    404. 

  404. 

  405. !
    406. 

  406. !  the 13 vegetation types are:
    407. 

  407. !
    408. 

  408. !  1  ...  broadleave-evergreen trees (tropical forest)
    409. 

  409. !  2  ...  broadleave-deciduous trees
    410. 

  410. !  3  ...  broadleave and needle leave trees (mixed forest)
    411. 

  411. !  4  ...  needleleave-evergreen trees
    412. 

  412. !  5  ...  needleleave-deciduous trees (larch)
    413. 

  413. !  6  ...  broadleave trees with groundcover (savanna)
    414. 

  414. !  7  ...  groundcover only (perenial)
    415. 

  415. !  8  ...  broadleave shrubs with perenial groundcover
    416. 

  416. !  9  ...  broadleave shrubs with bare soil
    417. 

  417. ! 10  ...  dwarf trees and shrubs with ground cover (trunda)
    418. 

  418. ! 11  ...  bare soil
    419. 

  419. ! 12  ...  cultivations (use parameters from type 7)
    420. 

  420. ! 13  ...  glacial
    421. 

  421. !
    422. 

  422.       data rsmax/13*5000./
    423. 

  423.       data rsmin/150.,100.,125.,150.,100.,70.,40.,                      &
    424. 

  424.      &           300.,400.,150.,999.,40.,999./
    425. 

  425.       data rgl/5*30.,65.,4*100.,999.,100.,999./
    426. 

  426.       data hs/41.69,54.53,51.93,47.35,47.35,54.53,36.35,                &
    427. 

  427.      &        3*42.00,999.,36.35,999./
    428. 

  428.       data smmax/.421,.464,.468,.434,.406,.465,.404,.439,.421/
    429. 

  429.       data smdry/.07,.14,.22,.08,.18,.16,.12,.10,.07/
    430. 

  430.       data smref/.283,.387,.412,.312,.338,.382,.315,.329,.283/
    431. 

  431.       data smwlt/.029,.119,.139,.047,.010,.103,.069,.066,.029/
    432. 

  432. !
    433. 

  433. !!!   save rsmax, rsmin, rgl, hs, smmax, smdry, smref, smwlt
    434. 

  434. !
    435. 

  435. 

  436. !WRF      DELT2 = DELT * 2.
    437. 

  437. !
    438. 

  438. !     ESTIMATE SIGMA ** K AT 2 M
    439. 

  439. !
    440. 

  440.       SIG2K = 1. - 4. * G * 2. / (CP * 280.)
    441. 

  441. !
    442. 

  442. !  INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIE
    443. 

  443. !  PSURF IS IN PASCALS
    444. 

  444. !  WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1
    445. 

  445. !  RHO IS DENSITY, QS1 IS SAT. HUM. AT LEVEL1 AND QSS IS SAT. HUM. AT
    446. 

  446. !  SURFACE
    447. 

  447. !  CONVERT SLRAD TO THE CIVILIZED UNIT FROM LANGLEY MINUTE-1 K-4
    448. 

  448. !  SURFACE ROUGHNESS LENGTH IS CONVERTED TO M FROM CM
    449. 

  449. !
    450. 

  450. !!
    451. 

  451. !     qs1 = fpvs(t1)
    452. 

  452. !     qss = fpvs(tskin)
    453. 

  453.       DO I=1,IM
    454. 

  454.         XRCL(I)  = SQRT(RCL(I))
    455. 

  455.         PSURF(I) = 1000. * PS(I)
    456. 

  456.         PS1(I)   = 1000. * PRSL1(I)
    457. 

  457. !       SLWD(I)  = SLRAD(I) * CONVRAD
    458. 

  458. !WRF        SLWD(I)  = SLRAD(I)
    459. 

  459. !
    460. 

  460. !  DLWFLX has been given a negative sign for downward longwave
    461. 

  461. !  snet is the net shortwave flux
    462. 

  462. !
    463. 

  463. !WRF        SNET(I) = -SLWD(I) - DLWFLX(I)
    464. 

  464.         WIND(I) = XRCL(I) * SQRT(U1(I) * U1(I) + V1(I) * V1(I))         &
    465. 

  465.      &              + MAX(0.0_kind_phys, MIN(DDVEL(I), 30.0_kind_phys))
    466. 

  466.         WIND(I) = MAX(WIND(I),1._kind_phys)
    467. 

  467.         Q0(I) = MAX(Q1(I),1.E-8_kind_phys)
    468. 

  468.         TSURF(I) = TSKIN(I)
    469. 

  469.         THETA1(I) = T1(I) * PRSLKI(I)
    470. 

  470.         TV1(I) = T1(I) * (1. + RVRDM1 * Q0(I))
    471. 

  471.         THV1(I) = THETA1(I) * (1. + RVRDM1 * Q0(I))
    472. 

  472.         TVS(I) = TSURF(I) * (1. + RVRDM1 * Q0(I))
    473. 

  473.         RHO(I) = PS1(I) / (RD * TV1(I))
    474. 

  474. !jfe    QS1(I) = 1000. * FPVS(T1(I))
    475. 

  475.         qs1(i) = fpvs(t1(i))
    476. 

  476.         QS1(I) = EPS * QS1(I) / (PS1(I) + EPSM1 * QS1(I))
    477. 

  477.         QS1(I) = MAX(QS1(I), 1.E-8_kind_phys)
    478. 

  478.         Q0(I) = min(QS1(I),Q0(I))
    479. 

  479. !jfe    QSS(I) = 1000. * FPVS(TSURF(I))
    480. 

  480.         qss(i) = fpvs(tskin(i))
    481. 

  481.         QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
    482. 

  482. !       RS = PLANTR
    483. 

  483.         RS(I) = 0.
    484. 

  484. !WRF        if(VEGTYPE(I).gt.0.) RS(I) = rsmin(VEGTYPE(I))
    485. 

  485.         Z0(I) = .01 * Z0RL(i)
    486. 

  486. !WRF        CANOPY(I)= MAX(CANOPY(I),0._kind_phys)
    487. 

  487.         DM(I) = 1.
    488. 

  488. !WRF
    489. 

  489.       GOTO 1111
    490. 

  490. !WRF
    491. 

  491.         FACTSNW(I) = 10.
    492. 

  492.         IF(SLIMSK(I).EQ.2.) FACTSNW(I) = 3.
    493. 

  493. !
    494. 

  494. !  SNOW DEPTH IN WATER EQUIVALENT IS CONVERTED FROM MM TO M UNIT
    495. 

  495. !
    496. 

  496.         SNOWD(I) = SHELEG(I) / 1000.
    497. 

  497.         FLAGSNW(I) = .FALSE.
    498. 

  498. !
    499. 

  499. !  WHEN SNOW DEPTH IS LESS THAN 1 MM, A PATCHY SNOW IS ASSUMED AND
    500. 

  500. !  SOIL IS ALLOWED TO INTERACT WITH THE ATMOSPHERE.
    501. 

  501. !  WE SHOULD EVENTUALLY MOVE TO A LINEAR COMBINATION OF SOIL AND
    502. 

  502. !  SNOW UNDER THE CONDITION OF PATCHY SNOW.
    503. 

  503. !
    504. 

  504.         IF(SNOWD(I).GT..001.OR.SLIMSK(I).EQ.2.) RS(I) = 0.
    505. 

  505.         IF(SNOWD(I).GT..001) FLAGSNW(I) = .TRUE.
    506. 

  506. !##DG  IF(LAT.EQ.LATD) THEN
    507. 

  507. !##DG    PRINT *, ' WIND,TV1,TVS,Q1,QS1,SNOW,SLIMSK=',
    508. 

  508. !##DG&   WIND,TV1,TVS,Q1,QS1,SNOWD,SLIMSK
    509. 

  509. !##DG    PRINT *, ' SNET, SLWD =', SNET, SLWD(I)
    510. 

  510. !##DG  ENDIF
    511. 

  511.         IF(SLIMSK(I).EQ.0.) THEN
    512. 

  512.           ZSOIL(I,1) = 0.
    513. 

  513.         ELSEIF(SLIMSK(I).EQ.1.) THEN
    514. 

  514.           ZSOIL(I,1) = -.10
    515. 

  515.         ELSE
    516. 

  516.           ZSOIL(I,1) = -3. / KM
    517. 

  517.         ENDIF
    518. 

  518. !WRF
    519. 

  519. 1111  CONTINUE
    520. 

  520. !WRF
    521. 

  521.       ENDDO
    522. 

  522. 

  523. !!
    524. 

  524. !WRF
    525. 

  525.       GOTO 2222
    526. 

  526. !WRF
    527. 

  527.       DO K = 2, KM
    528. 

  528.         DO I=1,IM
    529. 

  529.           IF(SLIMSK(I).EQ.0.) THEN
    530. 

  530.             ZSOIL(I,K) = 0.
    531. 

  531.           ELSEIF(SLIMSK(I).EQ.1.) THEN
    532. 

  532.             ZSOIL(I,K) = ZSOIL(I,K-1)                                   &
    533. 

  533.      &                   + (-2. - ZSOIL(I,1)) / (KM - 1)
    534. 

  534.           ELSE
    535. 

  535.             ZSOIL(I,K) = - 3. * FLOAT(K) / FLOAT(KM)
    536. 

  536.           ENDIF
    537. 

  537.         ENDDO
    538. 

  538.       ENDDO
    539. 

  539. !WRF
    540. 

  540. 2222  CONTINUE
    541. 

  541. !WRF
    542. 

  542. !!
    543. 

  543.       DO I=1,IM
    544. 

  544.         Z1(I) = -RD * TV1(I) * LOG(PS1(I)/PSURF(I)) / G
    545. 

  545.         DRAIN(I) = 0.
    546. 

  546.       ENDDO
    547. 

  547. 

  548. !!
    549. 

  549.       DO K = 1, KM
    550. 

  550.         DO I=1,IM
    551. 

  551.           ET(I,K) = 0.
    552. 

  552.           RHSMC(I,K) = 0.
    553. 

  553.           AIM(I,K) = 0.
    554. 

  554.           BIM(I,K) = 1.
    555. 

  555.           CIM(I,K) = 0.
    556. 

  556.           STSOIL(I,K) = STC(I,K)
    557. 

  557.         ENDDO
    558. 

  558.       ENDDO
    559. 

  559. 

  560.       DO I=1,IM
    561. 

  561.         EDIR(I) = 0.
    562. 

  562.         EC(I) = 0.
    563. 

  563.         EVAP(I) = 0.
    564. 

  564.         EP(I) = 0.
    565. 

  565.         SNOWMT(I) = 0.
    566. 

  566.         GFLUX(I) = 0.
    567. 

  567.         RHSCNPY(I) = 0.
    568. 

  568.         FX(I) = 0.
    569. 

  569.         ETPFAC(I) = 0.
    570. 

  570.         CANFAC(I) = 0.
    571. 

  571.       ENDDO
    572. 

  572. !
    573. 

  573. !  COMPUTE STABILITY DEPENDENT EXCHANGE COEFFICIENTS
    574. 

  574. !
    575. 

  575. !  THIS PORTION OF THE CODE IS PRESENTLY SUPPRESSED
    576. 

  576. !
    577. 

  577.       DO I=1,IM
    578. 

  578.         IF(SLIMSK(I).EQ.0.) THEN
    579. 

  579.           USTAR(I) = SQRT(G * Z0(I) / CHARNOCK)
    580. 

  580.         ENDIF
    581. 

  581. !
    582. 

  582. !  COMPUTE STABILITY INDICES (RB AND HLINF)
    583. 

  583. !
    584. 

  584. 

  585.         Z0MAX(I) = MIN(Z0(I),0.1 * Z1(I))
    586. 

  586.         ZTMAX(I) = Z0MAX(I)
    587. 

  587.         IF(SLIMSK(I).EQ.0.) THEN
    588. 

  588.           RESTAR = USTAR(I) * Z0MAX(I) / VIS
    589. 

  589.           RESTAR = MAX(RESTAR,.000001_kind_phys)
    590. 

  590. !         RESTAR = ALOG(RESTAR)
    591. 

  591. !         RESTAR = MIN(RESTAR,5.)
    592. 

  592. !         RESTAR = MAX(RESTAR,-5.)
    593. 

  593. !         RAT(I) = AA1 + BB1 * RESTAR + CC1 * RESTAR ** 2
    594. 

  594. !         RAT(I) = RAT(I) / (1. + BB2 * RESTAR
    595. 

  595. !    &                       + CC2 * RESTAR ** 2)
    596. 

  596. !  Rat taken from Zeng, Zhao and Dickinson 1997
    597. 

  597.           RAT(I) = 2.67 * restar ** .25 - 2.57
    598. 

  598.           RAT(I) = min(RAT(I),7._kind_phys)
    599. 

  599.           ZTMAX(I) = Z0MAX(I) * EXP(-RAT(I))
    600. 

  600.         ENDIF
    601. 

  601.       ENDDO
    602. 

  602. 

  603. !##DG  IF(LAT.EQ.LATD) THEN
    604. 

  604. !##DG    PRINT *, ' z0max, ztmax, restar, RAT(I) =',
    605. 

  605. !##DG &   z0max, ztmax, restar, RAT(I)
    606. 

  606. !##DG  ENDIF
    607. 

  607.       DO I = 1, IM
    608. 

  608.         DTV(I) = THV1(I) - TVS(I)
    609. 

  609.         ADTV = ABS(DTV(I))
    610. 

  610.         ADTV = MAX(ADTV,.001_kind_phys)
    611. 

  611.         DTV(I) = SIGN(1._kind_phys,DTV(I)) * ADTV
    612. 

  612.         RB(I) = G * DTV(I) * Z1(I) / (.5 * (THV1(I) + TVS(I))           &
    613. 

  613.      &          * WIND(I) * WIND(I))
    614. 

  614.         RB(I) = MAX(RB(I),-5000._kind_phys)
    615. 

  615. !        FM(I) = LOG((Z0MAX(I)+Z1(I)) / Z0MAX(I))
    616. 

  616. !        FH(I) = LOG((ZTMAX(I)+Z1(I)) / ZTMAX(I))
    617. 

  617.         FM(I) = LOG((Z1(I)) / Z0MAX(I))
    618. 

  618.         FH(I) = LOG((Z1(I)) / ZTMAX(I))
    619. 

  619.         HLINF(I) = RB(I) * FM(I) * FM(I) / FH(I)
    620. 

  620.         FM10(I) = LOG((Z0MAX(I)+10.) / Z0MAX(I))
    621. 

  621.         FH2(I) = LOG((ZTMAX(I)+2.) / ZTMAX(I))
    622. 

  622.       ENDDO
    623. 

  623. !##DG  IF(LAT.EQ.LATD) THEN
    624. 

  624. !##DG    PRINT *, ' DTV, RB(I), FM(I), FH(I), HLINF =',
    625. 

  625. !##DG &   dtv, rb, FM(I), FH(I), hlinf
    626. 

  626. !##DG  ENDIF
    627. 

  627. !
    628. 

  628. !  STABLE CASE
    629. 

  629. !
    630. 

  630.       DO I = 1, IM
    631. 

  631.         IF(DTV(I).GE.0.) THEN
    632. 

  632.           HL1(I) = HLINF(I)
    633. 

  633.         ENDIF
    634. 

  634.         IF(DTV(I).GE.0..AND.HLINF(I).GT..25) THEN
    635. 

  635.           HL0INF = Z0MAX(I) * HLINF(I) / Z1(I)
    636. 

  636.           HLTINF = ZTMAX(I) * HLINF(I) / Z1(I)
    637. 

  637.           AA = SQRT(1. + 4. * ALPHA * HLINF(I))
    638. 

  638.           AA0 = SQRT(1. + 4. * ALPHA * HL0INF)
    639. 

  639.           BB = AA
    640. 

  640.           BB0 = SQRT(1. + 4. * ALPHA * HLTINF)
    641. 

  641.           PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
    642. 

  642.           PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
    643. 

  643.           FMS = FM(I) - PM(I)
    644. 

  644.           FHS = FH(I) - PH(I)
    645. 

  645.           HL1(I) = FMS * FMS * RB(I) / FHS
    646. 

  646.         ENDIF
    647. 

  647.       ENDDO
    648. 

  648. !
    649. 

  649. !  SECOND ITERATION
    650. 

  650. !
    651. 

  651.       DO I = 1, IM
    652. 

  652.         IF(DTV(I).GE.0.) THEN
    653. 

  653.           HL0 = Z0MAX(I) * HL1(I) / Z1(I)
    654. 

  654.           HLT = ZTMAX(I) * HL1(I) / Z1(I)
    655. 

  655.           AA = SQRT(1. + 4. * ALPHA * HL1(I))
    656. 

  656.           AA0 = SQRT(1. + 4. * ALPHA * HL0)
    657. 

  657.           BB = AA
    658. 

  658.           BB0 = SQRT(1. + 4. * ALPHA * HLT)
    659. 

  659.           PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
    660. 

  660.           PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
    661. 

  661.           HL110 = HL1(I) * 10. / Z1(I)
    662. 

  662.           AA = SQRT(1. + 4. * ALPHA * HL110)
    663. 

  663.           PM10(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
    664. 

  664.           HL12(I) = HL1(I) * 2. / Z1(I)
    665. 

  665. !         AA = SQRT(1. + 4. * ALPHA * HL12(I))
    666. 

  666.           BB = SQRT(1. + 4. * ALPHA * HL12(I))
    667. 

  667.           PH2(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
    668. 

  668.         ENDIF
    669. 

  669.       ENDDO
    670. 

  670. !!
    671. 

  671. !##DG  IF(LAT.EQ.LATD) THEN
    672. 

  672. !##DG    PRINT *, ' HL1(I), PM, PH =',
    673. 

  673. !##DG &   HL1(I),  pm, ph
    674. 

  674. !##DG  ENDIF
    675. 

  675. !
    676. 

  676. !  UNSTABLE CASE
    677. 

  677. !
    678. 

  678. !
    679. 

  679. !  CHECK FOR UNPHYSICAL OBUKHOV LENGTH
    680. 

  680. !
    681. 

  681.       DO I=1,IM
    682. 

  682.         IF(DTV(I).LT.0.) THEN
    683. 

  683.           OLINF = Z1(I) / HLINF(I)
    684. 

  684.           IF(ABS(OLINF).LE.50. * Z0MAX(I)) THEN
    685. 

  685.             HLINF(I) = -Z1(I) / (50. * Z0MAX(I))
    686. 

  686.           ENDIF
    687. 

  687.         ENDIF
    688. 

  688.       ENDDO
    689. 

  689. !
    690. 

  690. !  GET PM AND PH
    691. 

  691. !
    692. 

  692.       DO I = 1, IM
    693. 

  693.         IF(DTV(I).LT.0..AND.HLINF(I).GE.-.5) THEN
    694. 

  694.           HL1(I) = HLINF(I)
    695. 

  695.           PM(I) = (A0 + A1 * HL1(I)) * HL1(I)                           &
    696. 

  696.      &            / (1. + B1 * HL1(I) + B2 * HL1(I) * HL1(I))
    697. 

  697.           PH(I) = (A0P + A1P * HL1(I)) * HL1(I)                         &
    698. 

  698.      &            / (1. + B1P * HL1(I) + B2P * HL1(I) * HL1(I))
    699. 

  699.           HL110 = HL1(I) * 10. / Z1(I)
    700. 

  700.           PM10(I) = (A0 + A1 * HL110) * HL110                           &
    701. 

  701.      &            / (1. + B1 * HL110 + B2 * HL110 * HL110)
    702. 

  702.           HL12(I) = HL1(I) * 2. / Z1(I)
    703. 

  703.           PH2(I) = (A0P + A1P * HL12(I)) * HL12(I)                      &
    704. 

  704.      &            / (1. + B1P * HL12(I) + B2P * HL12(I) * HL12(I))
    705. 

  705.         ENDIF
    706. 

  706.         IF(DTV(I).LT.0.AND.HLINF(I).LT.-.5) THEN
    707. 

  707.           HL1(I) = -HLINF(I)
    708. 

  708.           PM(I) = LOG(HL1(I)) + 2. * HL1(I) ** (-.25) - .8776
    709. 

  709.           PH(I) = LOG(HL1(I)) + .5 * HL1(I) ** (-.5) + 1.386
    710. 

  710.           HL110 = HL1(I) * 10. / Z1(I)
    711. 

  711.           PM10(I) = LOG(HL110) + 2. * HL110 ** (-.25) - .8776
    712. 

  712.           HL12(I) = HL1(I) * 2. / Z1(I)
    713. 

  713.           PH2(I) = LOG(HL12(I)) + .5 * HL12(I) ** (-.5) + 1.386
    714. 

  714.         ENDIF
    715. 

  715.       ENDDO
    716. 

  716. !
    717. 

  717. !  FINISH THE EXCHANGE COEFFICIENT COMPUTATION TO PROVIDE FM AND FH
    718. 

  718. !
    719. 

  719.       DO I = 1, IM
    720. 

  720. 

  721.         FM(I) = FM(I) - PM(I)
    722. 

  722.         FH(I) = FH(I) - PH(I)
    723. 

  723.         FM10(I) = FM10(I) - PM10(I)
    724. 

  724.         FH2(I) = FH2(I) - PH2(I)
    725. 

  725.         CM(I) = CA * CA / (FM(I) * FM(I))
    726. 

  726.         CH(I) = CA * CA / (FM(I) * FH(I))
    727. 

  727.         CQ = CH(I)
    728. 

  728.         STRESS(I) = CM(I) * WIND(I) * WIND(I)
    729. 

  729.         USTAR(I)  = SQRT(STRESS(I))
    730. 

  730. !       USTAR(I) = SQRT(CM(I) * WIND(I) * WIND(I))
    731. 

  731.       ENDDO
    732. 

  732. !##DG  IF(LAT.EQ.LATD) THEN
    733. 

  733. !##DG    PRINT *, ' FM, FH, CM, CH(I), USTAR =',
    734. 

  734. !##DG &   FM, FH, CM, ch, USTAR
    735. 

  735. !##DG  ENDIF
    736. 

  736. !
    737. 

  737. !  UPDATE Z0 OVER OCEAN
    738. 

  738. !
    739. 

  739.       DO I = 1, IM
    740. 

  740.         IF(SLIMSK(I).EQ.0.) THEN
    741. 

  741.           Z0(I) = (CHARNOCK / G) * USTAR(I) ** 2
    742. 

  742. !  NEW IMPLEMENTATION OF Z0
    743. 

  743. !         CC = USTAR(I) * Z0 / RNU
    744. 

  744. !         PP = CC / (1. + CC)
    745. 

  745. !         FF = G * ARNU / (CHARNOCK * USTAR(I) ** 3)
    746. 

  746. !         Z0 = ARNU / (USTAR(I) * FF ** PP)
    747. 

  747.           Z0(I) = MIN(Z0(I),.1_kind_phys)
    748. 

  748.           Z0(I) = MAX(Z0(I),1.E-7_kind_phys)
    749. 

  749.           Z0RL(I) = 100. * Z0(I)
    750. 

  750.         ENDIF
    751. 

  751.       ENDDO
    752. 

  752. 

  753. 	  GOTO 5555
    754. 

  754. !
    755. 

  755. !  RCP = RHO CP CH V
    756. 

  756. !
    757. 

  757.       DO I = 1, IM
    758. 

  758.         RCH(I) = RHO(I) * CP * CH(I) * WIND(I)
    759. 

  759.       ENDDO
    760. 

  760. 

  761. 

  762. !
    763. 

  763. !  SENSIBLE AND LATENT HEAT FLUX OVER OPEN WATER
    764. 

  764. !
    765. 

  765.       DO I = 1, IM
    766. 

  766.         IF(SLIMSK(I).EQ.0.) THEN
    767. 

  767.           EVAP(I) = ELOCP * RCH(I) * (QSS(I) - Q1(I))
    768. 

  768.           DM(I) = 1.
    769. 

  769.           QSURF(I) = QSS(I)
    770. 

  770.         ENDIF
    771. 

  771.       ENDDO
    772. 

  772. 

  773.         !
    774. 

  774.         !  COMPUTE SOIL/SNOW/ICE HEAT FLUX IN PREPARATION FOR SURFACE ENERGY
    775. 

  775.         !  BALANCE CALCULATION
    776. 

  776.         !
    777. 

  777. 	    DO I = 1, IM
    778. 

  778. 	      GFLUX(I) = 0.
    779. 

  779. 	      IF(SLIMSK(I).EQ.1.) THEN
    780. 

  780. 	        SMCZ(I) = .5 * (SMC(I,1) + .20)
    781. 

  781. 	        DFT0(I) = KTSOIL(SMCZ(I),SOILTYP(I))
    782. 

  782. 	      ELSEIF(SLIMSK(I).EQ.2.) THEN
    783. 

  783.         !  DF FOR ICE IS TAKEN FROM MAYKUT AND UNTERSTEINER
    784. 

  784.         !  DF IS IN SI UNIT OF W K-1 M-1
    785. 

  785. 	        DFT0(I) = 2.2
    786. 

  786. 	      ENDIF
    787. 

  787. 	    ENDDO
    788. 

  788.         !!
    789. 

  789. 	    DO I=1,IM
    790. 

  790. 	      IF(SLIMSK(I).NE.0.) THEN
    791. 

  791.         !         IF(SNOWD(I).GT..001) THEN
    792. 

  792. 	        IF(FLAGSNW(I)) THEN
    793. 

  793.         !
    794. 

  794.         !  WHEN SNOW COVERED, GROUND HEAT FLUX COMES FROM SNOW
    795. 

  795.         !
    796. 

  796. 		TFLX = MIN(T1(I), TSURF(I))
    797. 

  797. 		GFLUX(I) = -DFSNOW * (TFLX - STSOIL(I,1))                   &
    798. 

  798. 	   &                 / (FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
    799. 

  799. 	        ELSE
    800. 

  800. 		GFLUX(I) = DFT0(I) * (STSOIL(I,1) - TSURF(I))               &
    801. 

  801. 	   &                 / (-.5 * ZSOIL(I,1))
    802. 

  802. 	        ENDIF
    803. 

  803. 	        GFLUX(I) = MAX(GFLUX(I),-200._kind_phys)
    804. 

  804. 	        GFLUX(I) = MIN(GFLUX(I),+200._kind_phys)
    805. 

  805. 	      ENDIF
    806. 

  806. 	    ENDDO
    807. 

  807. 	    DO I = 1, IM
    808. 

  808. 	      FLAG(I) = SLIMSK(I).NE.0.
    809. 

  809. 	      PARTLND(I) = 1.
    810. 

  810. 	      IF(SNOWD(I).GT.0..AND.SNOWD(I).LE..001) THEN
    811. 

  811. 	        PARTLND(I) = 1. - SNOWD(I) / .001
    812. 

  812. 	      ENDIF
    813. 

  813. 	    ENDDO
    814. 

  814. 	    DO I = 1, IM
    815. 

  815. 	      SNOEVP(I) = 0.
    816. 

  816. 	      if(SNOWD(I).gt..001) PARTLND(I) = 0.
    817. 

  817. 	    ENDDO
    818. 

  818.         !
    819. 

  819.         !  COMPUTE POTENTIAL EVAPORATION FOR LAND AND SEA ICE
    820. 

  820.         !
    821. 

  821. 	    DO I = 1, IM
    822. 

  822. 	      IF(FLAG(I)) THEN
    823. 

  823. 	        T12 = T1(I) * T1(I)
    824. 

  824. 	        T14 = T12 * T12
    825. 

  825.         !
    826. 

  826.         !  RCAP = FNET - SIGMA T**4 + GFLX - RHO CP CH V (T1-THETA1)
    827. 

  827.         !
    828. 

  828. 	        RCAP(I) = -SLWD(I) - SIGMA * T14 + GFLUX(I)                   &
    829. 

  829. 	   &              - RCH(I) * (T1(I) - THETA1(I))
    830. 

  830.         !
    831. 

  831.         !  RSMALL = 4 SIGMA T**3 / RCH(I) + 1
    832. 

  832.         !
    833. 

  833. 	        RSMALL(I) = 4. * SIGMA * T1(I) * T12 / RCH(I) + 1.
    834. 

  834.         !
    835. 

  835.         !  DELTA = L / CP * DQS/DT
    836. 

  836.         !
    837. 

  837. 	        DELTA(I) = ELOCP * EPS * HVAP * QS1(I) / (RD * T12)
    838. 

  838.         !
    839. 

  839.         !  POTENTIAL EVAPOTRANSPIRATION ( WATTS / M**2 ) AND
    840. 

  840.         !  POTENTIAL EVAPORATION
    841. 

  841.         !
    842. 

  842. 	        TERM1(I) = ELOCP * RSMALL(I) * RCH(I)*(QS1(I)-Q0(I))
    843. 

  843. 	        TERM2(I) = RCAP(I) * DELTA(I)
    844. 

  844. 	        EP(I) = (ELOCP * RSMALL(I) * RCH(I) * (QS1(I) - Q0(I))        &
    845. 

  845. 	   &              + RCAP(I) * DELTA(I))
    846. 

  846. 	        EP(I) = EP(I) / (RSMALL(I) + DELTA(I))
    847. 

  847. 	      ENDIF
    848. 

  848. 	    ENDDO
    849. 

  849.         !
    850. 

  850.         !  ACTUAL EVAPORATION OVER LAND IN THREE PARTS : EDIR, ET, AND EC
    851. 

  851.         !
    852. 

  852.         !  DIRECT EVAPORATION FROM SOIL, THE UNIT GOES FROM M S-1 TO KG M-2 S-1
    853. 

  853.         !
    854. 

  854. 	    DO I = 1, IM
    855. 

  855. 	      FLAG(I) = SLIMSK(I).EQ.1..AND.EP(I).GT.0.
    856. 

  856. 	    ENDDO
    857. 

  857. 	    DO I = 1, IM
    858. 

  858. 	      IF(FLAG(I)) THEN
    859. 

  859. 	        DF1(I) = FUNCDF(SMC(I,1),SOILTYP(I))
    860. 

  860. 	        KT1(I) = FUNCKT(SMC(I,1),SOILTYP(I))
    861. 

  861. 	      endif
    862. 

  862. 	      if(FLAG(I).and.STC(I,1).lt.t0c) then
    863. 

  863. 	        DF1(I) = 0.
    864. 

  864. 	        KT1(I) = 0.
    865. 

  865. 	      endif
    866. 

  866. 	      IF(FLAG(I)) THEN
    867. 

  867.         !         TREF = .75 * THSAT(SOILTYP(I))
    868. 

  868. 	        TREF(I) = smref(SOILTYP(I))
    869. 

  869.         !         TWILT = TWLT(SOILTYP(I))
    870. 

  870. 	        TWILT(I) = smwlt(SOILTYP(I))
    871. 

  871. 	        smcdry = smdry(SOILTYP(I))
    872. 

  872.         !         FX(I) = -2. * DF1(I) * (SMC(I,1) - .23) / ZSOIL(I,1)
    873. 

  873.         !    &            - KT1(I)
    874. 

  874. 	        FX(I) = -2. * DF1(I) * (SMC(I,1) - smcdry) / ZSOIL(I,1)       &
    875. 

  875. 	   &            - KT1(I)
    876. 

  876. 	        FX(I) = MIN(FX(I), EP(I)/HVAP)
    877. 

  877. 	        FX(I) = MAX(FX(I),0._kind_phys)
    878. 

  878.         !
    879. 

  879.         !  SIGMAF IS THE FRACTION OF AREA COVERED BY VEGETATION
    880. 

  880.         !
    881. 

  881. 	        EDIR(I) = FX(I) * (1. - SIGMAF(I)) * PARTLND(I)
    882. 

  882. 	      ENDIF
    883. 

  883. 	    ENDDO
    884. 

  884.         !
    885. 

  885.         !  calculate stomatal resistance
    886. 

  886.         !
    887. 

  887. 	    DO I = 1, IM
    888. 

  888. 	      if(FLAG(I)) then
    889. 

  889.         !
    890. 

  890.         !  resistance due to PAR. We use net solar flux as proxy at the present time
    891. 

  891.         !
    892. 

  892. 	        ff = .55 * 2. * SNET(I) / rgl(VEGTYPE(I))
    893. 

  893. 	        rcs = (ff + RS(I)/rsmax(VEGTYPE(I))) / (1. + ff)
    894. 

  894. 	        rcs = max(rcs,.0001_kind_phys)
    895. 

  895. 	        rct = 1.
    896. 

  896. 	        rcq = 1.
    897. 

  897.         !
    898. 

  898.         !  resistance due to thermal effect
    899. 

  899.         !
    900. 

  900.         !         rct = 1. - .0016 * (topt - theta1) ** 2
    901. 

  901.         !         rct = max(rct,.0001)
    902. 

  902.         !
    903. 

  903.         !  resistance due to humidity
    904. 

  904.         !
    905. 

  905.         !         rcq = 1. / (1. + hs(VEGTYPE(I)) * (QS1(I) - Q0(I)))
    906. 

  906.         !         rcq = max(rcq,.0001)
    907. 

  907.         !
    908. 

  908.         !  compute resistance without the effect of soil moisture
    909. 

  909.         !
    910. 

  910. 	        RS(I) = RS(I) / (rcs * rct * rcq)
    911. 

  911. 	      endif
    912. 

  912. 	    ENDDO
    913. 

  913.         !
    914. 

  914.         !  TRANSPIRATION FROM ALL LEVELS OF THE SOIL
    915. 

  915.         !
    916. 

  916. 	    DO I = 1, IM
    917. 

  917. 	      IF(FLAG(I)) THEN
    918. 

  918. 	        CANFAC(I) = (CANOPY(I) / SCANOP) ** CFACTR
    919. 

  919. 	      endif
    920. 

  920. 	      IF(FLAG(I)) THEN
    921. 

  921. 	        ETPFAC(I) = SIGMAF(I)                                         &
    922. 

  922. 	   &           * (1. - CANFAC(I)) / HVAP
    923. 

  923. 	        GX(I) = (SMC(I,1) - TWILT(I)) / (TREF(I) - TWILT(I))
    924. 

  924. 	        GX(I) = MAX(GX(I),0._kind_phys)
    925. 

  925. 	        GX(I) = MIN(GX(I),1._kind_phys)
    926. 

  926.         !
    927. 

  927.         !  resistance due to soil moisture deficit
    928. 

  928.         !
    929. 

  929. 	        rss = GX(I) * (ZSOIL(I,1) / ZSOIL(I,km))
    930. 

  930. 	        rss = max(rss,.0001_kind_phys)
    931. 

  931. 	        RSI = RS(I) / rss
    932. 

  932.         !
    933. 

  933.         !  transpiration a la Monteith
    934. 

  934.         !
    935. 

  935. 	        eth = (TERM1(I) + TERM2(I)) /                                 &
    936. 

  936. 	   &          (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
    937. 

  937. 	        ET(I,1) = ETPFAC(I) * eth                                     &
    938. 

  938. 	   &            * PARTLND(I)
    939. 

  939. 	      ENDIF
    940. 

  940. 	    ENDDO
    941. 

  941.         !!
    942. 

  942. 	    DO K = 2, KM
    943. 

  943. 	      DO I=1,IM
    944. 

  944. 	        IF(FLAG(I)) THEN
    945. 

  945. 		GX(I) = (SMC(I,K) - TWILT(I)) / (TREF(I) - TWILT(I))
    946. 

  946. 		GX(I) = MAX(GX(I),0._kind_phys)
    947. 

  947. 		GX(I) = MIN(GX(I),1._kind_phys)
    948. 

  948.         !
    949. 

  949.         !  resistance due to soil moisture deficit
    950. 

  950.         !
    951. 

  951. 	        rss = GX(I) * ((ZSOIL(I,k) - ZSOIL(I,k-1))/ZSOIL(I,km))
    952. 

  952. 	        rss = max(rss,1.e-6_kind_phys)
    953. 

  953. 	        RSI = RS(I) / rss
    954. 

  954.         !
    955. 

  955.         !  transpiration a la Monteith
    956. 

  956.         !
    957. 

  957. 	        eth = (TERM1(I) + TERM2(I)) /                                 &
    958. 

  958. 	   &          (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
    959. 

  959. 		ET(I,K) = eth                                               &
    960. 

  960. 	   &               * ETPFAC(I) * PARTLND(I)
    961. 

  961. 	        ENDIF
    962. 

  962. 	      ENDDO
    963. 

  963. 	    ENDDO
    964. 

  964.         !!
    965. 

  965.          400  CONTINUE
    966. 

  966.         !
    967. 

  967.         !  CANOPY RE-EVAPORATION
    968. 

  968.         !
    969. 

  969. 	    DO I=1,IM
    970. 

  970. 	      IF(FLAG(I)) THEN
    971. 

  971. 	        EC(I) = SIGMAF(I) * CANFAC(I) * EP(I) / HVAP
    972. 

  972. 	        EC(I) = EC(I) * PARTLND(I)
    973. 

  973. 	        EC(I) = min(EC(I),CANOPY(I)/delt)
    974. 

  974. 	      ENDIF
    975. 

  975. 	    ENDDO
    976. 

  976.         !
    977. 

  977.         !  SUM UP TOTAL EVAPORATION
    978. 

  978.         !
    979. 

  979. 	    DO I = 1, IM
    980. 

  980. 	      IF(FLAG(I)) THEN
    981. 

  981. 	       EVAP(I) = EDIR(I) + EC(I)
    982. 

  982. 	      ENDIF
    983. 

  983. 	    ENDDO
    984. 

  984.         !!
    985. 

  985. 	    DO K = 1, KM
    986. 

  986. 	      DO I=1,IM
    987. 

  987. 	        IF(FLAG(I)) THEN
    988. 

  988. 		EVAP(I) = EVAP(I) + ET(I,K)
    989. 

  989. 	        ENDIF
    990. 

  990. 	      ENDDO
    991. 

  991. 	    ENDDO
    992. 

  992.         !!
    993. 

  993.         !
    994. 

  994.         !  RETURN EVAP UNIT FROM KG M-2 S-1 TO WATTS M-2
    995. 

  995.         !
    996. 

  996. 	    DO I=1,IM
    997. 

  997. 	      IF(FLAG(I)) THEN
    998. 

  998. 	        EVAP(I) = MIN(EVAP(I)*HVAP,EP(I))
    999. 

  999. 	      ENDIF
    1000. 

  1000. 	    ENDDO
    1001. 

  1001.         !##DG  IF(LAT.EQ.LATD) THEN
    1002. 

  1002.         !##DG    PRINT *, 'FX(I), SIGMAF, EDIR(I), ETPFAC=', FX(I)*HVAP,SIGMAF,
    1003. 

  1003.         !##DG&          EDIR(I)*HVAP,ETPFAC*HVAP
    1004. 

  1004.         !##DG    PRINT *, ' ET =', (ET(K)*HVAP,K=1,KM)
    1005. 

  1005.         !##DG    PRINT *, ' CANFAC(I), EC(I), EVAP', CANFAC(I),EC(I)*HVAP,EVAP
    1006. 

  1006.         !##DG  ENDIF
    1007. 

  1007.         !
    1008. 

  1008.         !  EVAPORATION OVER BARE SEA ICE
    1009. 

  1009.         !
    1010. 

  1010. 	    DO I = 1, IM
    1011. 

  1011.         !       IF(SLIMSK(I).EQ.2.AND.SNOWD(I).LE..001) THEN
    1012. 

  1012. 	      IF(SLIMSK(I).EQ.2.) THEN
    1013. 

  1013. 	        EVAP(I) = PARTLND(I) * EP(I)
    1014. 

  1014. 	      ENDIF
    1015. 

  1015. 	    ENDDO
    1016. 

  1016.         !
    1017. 

  1017.         !  TREAT DOWNWARD MOISTURE FLUX SITUATION
    1018. 

  1018.         !  (EVAP WAS PRESET TO ZERO SO NO UPDATE NEEDED)
    1019. 

  1019.         !  DEW IS CONVERTED FROM KG M-2 TO M TO CONFORM TO PRECIP UNIT
    1020. 

  1020.         !
    1021. 

  1021. 	    DO I = 1, IM
    1022. 

  1022. 	      FLAG(I) = SLIMSK(I).NE.0..AND.EP(I).LE.0.
    1023. 

  1023. 	      DEW(I) = 0.
    1024. 

  1024. 	    ENDDO
    1025. 

  1025. 	    DO I = 1, IM
    1026. 

  1026. 	      IF(FLAG(I)) THEN
    1027. 

  1027. 	        DEW(I) = -EP(I) * DELT / (HVAP * RHOH2O)
    1028. 

  1028. 	        EVAP(I) = EP(I)
    1029. 

  1029. 	        DEW(I) = DEW(I) * PARTLND(I)
    1030. 

  1030. 	        EVAP(I) = EVAP(I) * PARTLND(I)
    1031. 

  1031. 	        DM(I) = 1.
    1032. 

  1032. 	      ENDIF
    1033. 

  1033. 	    ENDDO
    1034. 

  1034.         !
    1035. 

  1035.         !  SNOW COVERED LAND AND SEA ICE
    1036. 

  1036.         !
    1037. 

  1037. 	    DO I = 1, IM
    1038. 

  1038. 	      FLAG(I) = SLIMSK(I).NE.0..AND.SNOWD(I).GT.0.
    1039. 

  1039. 	    ENDDO
    1040. 

  1040.         !
    1041. 

  1041.         !  CHANGE OF SNOW DEPTH DUE TO EVAPORATION OR SUBLIMATION
    1042. 

  1042.         !
    1043. 

  1043.         !  CONVERT EVAP FROM KG M-2 S-1 TO M S-1 TO DETERMINE THE REDUCTION OF S
    1044. 

  1044.         !
    1045. 

  1045. 	    DO I = 1, IM
    1046. 

  1046. 	      IF(FLAG(I)) THEN
    1047. 

  1047. 	        BFACT = SNOWD(I) / (DELT * EP(I) / (HVAP * RHOH2O))
    1048. 

  1048. 	        BFACT = MIN(BFACT,1._kind_phys)
    1049. 

  1049.         !
    1050. 

  1050.         !  THE EVAPORATION OF SNOW
    1051. 

  1051.         !
    1052. 

  1052. 	        IF(EP(I).LE.0.) BFACT = 1.
    1053. 

  1053. 	        IF(SNOWD(I).LE..001) THEN
    1054. 

  1054.         !           EVAP = (SNOWD(I)/.001)*BFACT*EP(I) + EVAP
    1055. 

  1055.         !           SNOEVP(I) = bfact * EP(I) * (1. - PARTLND(I))
    1056. 

  1056.         !           EVAP = EVAP + SNOEVP(I)
    1057. 

  1057. 		SNOEVP(I) = bfact * EP(I)
    1058. 

  1058.         !           EVAP   = EVAP + SNOEVP(I) * (1. - PARTLND(I))
    1059. 

  1059. 		EVAP(I)=EVAP(I)+SNOEVP(I)*(1.-PARTLND(I))
    1060. 

  1060. 	        ELSE
    1061. 

  1061.         !           EVAP(I) = BFACT * EP(I)
    1062. 

  1062. 		SNOEVP(I) = bfact * EP(I)
    1063. 

  1063. 		EVAP(I) = SNOEVP(I)
    1064. 

  1064. 	        ENDIF
    1065. 

  1065. 	        TSURF(I) = T1(I) +                                            &
    1066. 

  1066. 	   &          (RCAP(I) - GFLUX(I) - DFSNOW * (T1(I) - STSOIL(I,1))    &
    1067. 

  1067. 	   &           /(FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))           &
    1068. 

  1068.         !    &           + THETA1 - T1                                          &
    1069. 

  1069.         !    &           - BFACT * EP(I)) / (RSMALL(I) * RCH(I)                 &
    1070. 

  1070. 	   &           - SNOEVP(I)) / (RSMALL(I) * RCH(I)                     &
    1071. 

  1071. 	   &           + DFSNOW / (FACTSNW(I)* MAX(SNOWD(I),.001_kind_phys)))
    1072. 

  1072.         !         SNOWD(I) = SNOWD(I) - BFACT * EP(I) * DELT / (RHOH2O * HVAP)
    1073. 

  1073. 	        SNOWD(I) = SNOWD(I) - SNOEVP(I) * delt / (rhoh2o * hvap)
    1074. 

  1074. 	        SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
    1075. 

  1075. 	      ENDIF
    1076. 

  1076. 	    ENDDO
    1077. 

  1077.         !
    1078. 

  1078.         !  SNOW MELT (M)
    1079. 

  1079.         !
    1080. 

  1080.          500  CONTINUE
    1081. 

  1081. 	    DO I = 1, IM
    1082. 

  1082. 	      FLAG(I) = SLIMSK(I).NE.0.                                       &
    1083. 

  1083. 	   &            .AND.SNOWD(I).GT..0
    1084. 

  1084. 	    ENDDO
    1085. 

  1085. 	    DO I = 1, IM
    1086. 

  1086. 	      IF(FLAG(I).AND.TSURF(I).GT.T0C) THEN
    1087. 

  1087. 	        SNOWMT(I) = RCH(I) * RSMALL(I) * DELT                         &
    1088. 

  1088. 	   &              * (TSURF(I) - T0C) / (RHOH2O * HFUS)
    1089. 

  1089. 	        SNOWMT(I) = min(SNOWMT(I),SNOWD(I))
    1090. 

  1090. 	        SNOWD(I) = SNOWD(I) - SNOWMT(I)
    1091. 

  1091. 	        SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
    1092. 

  1092. 	        TSURF(I) = MAX(T0C,TSURF(I)                                   &
    1093. 

  1093. 	   &             -HFUS*SNOWMT(I)*RHOH2O/(RCH(I)*RSMALL(I)*DELT))
    1094. 

  1094. 	      ENDIF
    1095. 

  1095. 	    ENDDO
    1096. 

  1096.         !
    1097. 

  1097.         !  We need to re-evaluate evaporation because of snow melt
    1098. 

  1098.         !    the skin temperature is now bounded to 0 deg C
    1099. 

  1099.         !
    1100. 

  1100.         !     qss = fpvs(tsurf)
    1101. 

  1101. 	    DO I = 1, IM
    1102. 

  1102.         !       IF (SNOWD(I) .GT. 0.0) THEN
    1103. 

  1103. 	      IF (SNOWD(I) .GT. snomin) THEN
    1104. 

  1104.         !jfe      QSS(I) = 1000. * FPVS(TSURF(I))
    1105. 

  1105. 	        qss(i) = fpvs(tsurf(i))
    1106. 

  1106. 	        QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
    1107. 

  1107. 	        EVAP(I) = elocp * RCH(I) * (QSS(I) - Q0(I))
    1108. 

  1108. 	      ENDIF
    1109. 

  1109. 	    ENDDO
    1110. 

  1110.         !
    1111. 

  1111.         !  PREPARE TENDENCY TERMS FOR THE SOIL MOISTURE FIELD WITHOUT PRECIPITAT
    1112. 

  1112.         !  THE UNIT OF MOISTURE FLUX NEEDS TO BECOME M S-1 FOR SOIL MOISTURE
    1113. 

  1113.         !   HENCE THE FACTOR OF RHOH2O
    1114. 

  1114.         !
    1115. 

  1115. 	    DO I = 1, IM
    1116. 

  1116. 	      FLAG(I) = SLIMSK(I).EQ.1.
    1117. 

  1117. 	      if(FLAG(I)) then
    1118. 

  1118. 	        DF1(I) = FUNCDF(SMCZ(I),SOILTYP(I))
    1119. 

  1119. 	        KT1(I) = FUNCKT(SMCZ(I),SOILTYP(I))
    1120. 

  1120. 	      endif
    1121. 

  1121. 	      if(FLAG(I).and.STC(I,1).lt.t0c) then
    1122. 

  1122. 	        DF1(I) = 0.
    1123. 

  1123. 	        KT1(I) = 0.
    1124. 

  1124. 	      endif
    1125. 

  1125. 	      IF(FLAG(I)) THEN
    1126. 

  1126. 	        RHSCNPY(I) = -EC(I) + SIGMAF(I) * RHOH2O * DEW(I) / DELT
    1127. 

  1127. 	        SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
    1128. 

  1128. 	        DMDZ(I) = (SMC(I,1) - SMC(I,2)) / (-.5 * ZSOIL(I,2))
    1129. 

  1129. 	        RHSMC(I,1) = (DF1(I) * DMDZ(I) + KT1(I)                       &
    1130. 

  1130. 	   &        + (EDIR(I) + ET(I,1))) / (ZSOIL(I,1) * RHOH2O)
    1131. 

  1131. 	        RHSMC(I,1) = RHSMC(I,1) - (1. - SIGMAF(I)) * DEW(I) /         &
    1132. 

  1132. 	   &                 ( ZSOIL(I,1) * delt)
    1133. 

  1133. 	        DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
    1134. 

  1134.         !
    1135. 

  1135.         !  AIM, BIM, AND CIM ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
    1136. 

  1136.         !  IMPLICIT UPDATE OF THE SOIL MOISTURE
    1137. 

  1137.         !
    1138. 

  1138. 	        AIM(I,1) = 0.
    1139. 

  1139. 	        BIM(I,1) = DF1(I) * DDZ(I) / (-ZSOIL(I,1) * RHOH2O)
    1140. 

  1140. 	        CIM(I,1) = -BIM(I,1)
    1141. 

  1141. 	      ENDIF
    1142. 

  1142. 	    ENDDO
    1143. 

  1143.         !!
    1144. 

  1144. 	    DO K = 2, KM
    1145. 

  1145. 	      IF(K.LT.KM) THEN
    1146. 

  1146. 	        DO I=1,IM
    1147. 

  1147. 		IF(FLAG(I)) THEN
    1148. 

  1148. 		  DF2 = FUNCDF(SMCZ(I),SOILTYP(I))
    1149. 

  1149. 		  KT2(I) = FUNCKT(SMCZ(I),SOILTYP(I))
    1150. 

  1150. 		ENDIF
    1151. 

  1151. 		IF(FLAG(I).and.STC(I,k).lt.t0c) THEN
    1152. 

  1152. 		  df2 = 0.
    1153. 

  1153. 		  KT2(I) = 0.
    1154. 

  1154. 		ENDIF
    1155. 

  1155. 		IF(FLAG(I)) THEN
    1156. 

  1156. 		  DMDZ2(I) = (SMC(I,K) - SMC(I,K+1))                        &
    1157. 

  1157. 	   &                   / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
    1158. 

  1158. 		  SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
    1159. 

  1159. 		  RHSMC(I,K) = (DF2 * DMDZ2(I) + KT2(I)                     &
    1160. 

  1160. 	   &             - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K))               &
    1161. 

  1161. 	   &                     / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
    1162. 

  1162. 		  DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
    1163. 

  1163. 		  CIM(I,K) = -DF2 * DDZ2(I)                                 &
    1164. 

  1164. 	   &                / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
    1165. 

  1165. 		ENDIF
    1166. 

  1166. 	        ENDDO
    1167. 

  1167. 	      ELSE
    1168. 

  1168. 	        DO I = 1, IM
    1169. 

  1169. 		IF(FLAG(I)) THEN
    1170. 

  1170. 		  KT2(I) = FUNCKT(SMC(I,K),SOILTYP(I))
    1171. 

  1171. 		ENDIF
    1172. 

  1172. 		if(FLAG(I).and.STC(I,k).lt.t0c) KT2(I) = 0.
    1173. 

  1173. 		IF(FLAG(I)) THEN
    1174. 

  1174. 		  RHSMC(I,K) = (KT2(I)                                      &
    1175. 

  1175. 	   &             - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K))               &
    1176. 

  1176. 	   &                     / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
    1177. 

  1177. 		  DRAIN(I) = KT2(I)
    1178. 

  1178. 		  CIM(I,K) = 0.
    1179. 

  1179. 		ENDIF
    1180. 

  1180. 	        ENDDO
    1181. 

  1181. 	      ENDIF
    1182. 

  1182. 	      DO I = 1, IM
    1183. 

  1183. 	        IF(FLAG(I)) THEN
    1184. 

  1184. 		AIM(I,K) = -DF1(I) * DDZ(I)                                 &
    1185. 

  1185. 	   &                / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
    1186. 

  1186. 		BIM(I,K) = -(AIM(I,K) + CIM(I,K))
    1187. 

  1187. 		DF1(I) = DF2
    1188. 

  1188. 		KT1(I) = KT2(I)
    1189. 

  1189. 		DMDZ(I) = DMDZ2(I)
    1190. 

  1190. 		DDZ(I) = DDZ2(I)
    1191. 

  1191. 	        ENDIF
    1192. 

  1192. 	      ENDDO
    1193. 

  1193. 	    ENDDO
    1194. 

  1194.         !!
    1195. 

  1195.          600  CONTINUE
    1196. 

  1196.         !
    1197. 

  1197.         !  UPDATE SOIL TEMPERATURE AND SEA ICE TEMPERATURE
    1198. 

  1198.         !
    1199. 

  1199. 	    DO I=1,IM
    1200. 

  1200. 	      FLAG(I) = SLIMSK(I).NE.0.
    1201. 

  1201. 	    ENDDO
    1202. 

  1202.         !
    1203. 

  1203.         !  SURFACE TEMPERATURE IS PART OF THE UPDATE WHEN SNOW IS ABSENT
    1204. 

  1204.         !
    1205. 

  1205. 	    DO I=1,IM
    1206. 

  1206.         !       IF(FLAG(I).AND.SNOWD(I).LE..001) THEN
    1207. 

  1207. 	      IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
    1208. 

  1208. 	        YY(I) = T1(I) +                                               &
    1209. 

  1209.         !    &          (RCAP(I)-GFLUX(I) + THETA1 - T1(I)                      &
    1210. 

  1210. 	   &          (RCAP(I)-GFLUX(I)                                       &
    1211. 

  1211. 	   &           - EVAP(I)) / (RSMALL(I) * RCH(I))
    1212. 

  1212. 	        ZZ(I) = 1. + DFT0(I) / (-.5 * ZSOIL(I,1) * RCH(I) * RSMALL(I))
    1213. 

  1213. 	        XX(I) = DFT0(I) * (STSOIL(I,1) - YY(I)) /                     &
    1214. 

  1214. 	   &            (.5 * ZSOIL(I,1) * ZZ(I))
    1215. 

  1215. 	      ENDIF
    1216. 

  1216.         !       IF(FLAG(I).AND.SNOWD(I).GT..001) THEN
    1217. 

  1217. 	      IF(FLAG(I).AND.FLAGSNW(I)) THEN
    1218. 

  1218. 	        YY(I) = STSOIL(I,1)
    1219. 

  1219.         !
    1220. 

  1220.         !  HEAT FLUX FROM SNOW IS EXPLICIT IN TIME
    1221. 

  1221.         !
    1222. 

  1222. 	        ZZ(I) = 1.
    1223. 

  1223. 	        XX(I) = DFSNOW * (STSOIL(I,1) - TSURF(I))                     &
    1224. 

  1224. 	   &            / (-FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
    1225. 

  1225. 	      ENDIF
    1226. 

  1226. 	    ENDDO
    1227. 

  1227.         !
    1228. 

  1228.         !  COMPUTE THE FORCING AND THE IMPLICIT MATRIX ELEMENTS FOR UPDATE
    1229. 

  1229.         !
    1230. 

  1230.         !  CH2O IS THE HEAT CAPACITY OF WATER AND CSOIL IS THE HEAT CAPACITY OF
    1231. 

  1231.         !
    1232. 

  1232. 	    DO I = 1, IM
    1233. 

  1233. 	      IF(FLAG(I)) THEN
    1234. 

  1234. 	        SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
    1235. 

  1235. 	        DTDZ1(I) = (STSOIL(I,1) - STSOIL(I,2)) / (-.5 * ZSOIL(I,2))
    1236. 

  1236. 	        IF(SLIMSK(I).EQ.1.) THEN
    1237. 

  1237. 		DFT1(I) = KTSOIL(SMCZ(I),SOILTYP(I))
    1238. 

  1238. 		HCPCT(I) = SMC(I,1) * CH2O + (1. - SMC(I,1)) * CSOIL
    1239. 

  1239. 	        ELSE
    1240. 

  1240. 		DFT1(I) = DFT0(I)
    1241. 

  1241. 		HCPCT(I) = CICE
    1242. 

  1242. 	        ENDIF
    1243. 

  1243. 	        DFT2(I) = DFT1(I)
    1244. 

  1244. 	        DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
    1245. 

  1245.         !
    1246. 

  1246.         !  AI, BI, AND CI ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
    1247. 

  1247.         !  IMPLICIT UPDATE OF THE SOIL TEMPERATURE
    1248. 

  1248.         !
    1249. 

  1249. 	        AI(I,1) = 0.
    1250. 

  1250. 	        BI(I,1) = DFT1(I) * DDZ(I) / (-ZSOIL(I,1) * HCPCT(I))
    1251. 

  1251. 	        CI(I,1) = -BI(I,1)
    1252. 

  1252. 	        BI(I,1) = BI(I,1)                                             &
    1253. 

  1253. 	   &            + DFT0(I) / (.5 * ZSOIL(I,1) **2 * HCPCT(I) * ZZ(I))
    1254. 

  1254.         !         SS = DFT0(I) * (STSOIL(I,1) - YY(I))                          &
    1255. 

  1255.         !    &         / (.5 * ZSOIL(I,1) * ZZ(I))
    1256. 

  1256.         !         RHSTC(1) = (DFT1(I) * DTDZ1(I) - SS)
    1257. 

  1257. 	        RHSTC(I,1) = (DFT1(I) * DTDZ1(I) - XX(I))                     &
    1258. 

  1258. 	   &                 / (ZSOIL(I,1) * HCPCT(I))
    1259. 

  1259. 	      ENDIF
    1260. 

  1260. 	    ENDDO
    1261. 

  1261.         !!
    1262. 

  1262. 	    DO K = 2, KM
    1263. 

  1263. 	      DO I=1,IM
    1264. 

  1264. 	        IF(SLIMSK(I).EQ.1.) THEN
    1265. 

  1265. 		HCPCT(I) = SMC(I,K) * CH2O + (1. - SMC(I,K)) * CSOIL
    1266. 

  1266. 	        ELSEIF(SLIMSK(I).EQ.2.) THEN
    1267. 

  1267. 		HCPCT(I) = CICE
    1268. 

  1268. 	        ENDIF
    1269. 

  1269. 	      ENDDO
    1270. 

  1270. 	      IF(K.LT.KM) THEN
    1271. 

  1271. 	        DO I = 1, IM
    1272. 

  1272. 		IF(FLAG(I)) THEN
    1273. 

  1273. 		  DTDZ2(I) = (STSOIL(I,K) - STSOIL(I,K+1))                  &
    1274. 

  1274. 	   &                   / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
    1275. 

  1275. 		  SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
    1276. 

  1276. 		  IF(SLIMSK(I).EQ.1.) THEN
    1277. 

  1277. 		    DFT2(I) = KTSOIL(SMCZ(I),SOILTYP(I))
    1278. 

  1278. 		  ENDIF
    1279. 

  1279. 		  DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
    1280. 

  1280. 		  CI(I,K) = -DFT2(I) * DDZ2(I)                              &
    1281. 

  1281. 	   &                / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
    1282. 

  1282. 		ENDIF
    1283. 

  1283. 	        ENDDO
    1284. 

  1284. 	      ELSE
    1285. 

  1285.         !
    1286. 

  1286.         !  AT THE BOTTOM, CLIMATOLOGY IS ASSUMED AT 2M DEPTH FOR LAND AND
    1287. 

  1287.         !  FREEZING TEMPERATURE IS ASSUMED FOR SEA ICE AT Z(KM)
    1288. 

  1288. 	        DO I = 1, IM
    1289. 

  1289. 		IF(SLIMSK(I).EQ.1.) THEN
    1290. 

  1290. 		  DTDZ2(I) = (STSOIL(I,K) - TG3(I))                         &
    1291. 

  1291. 	   &              / (.5 * (ZSOIL(I,K-1) + ZSOIL(I,K)) - ZBOT)
    1292. 

  1292. 		  DFT2(I) = KTSOIL(SMC(I,K),SOILTYP(I))
    1293. 

  1293. 		  CI(I,K) = 0.
    1294. 

  1294. 		ENDIF
    1295. 

  1295. 		IF(SLIMSK(I).EQ.2.) THEN
    1296. 

  1296. 		  DTDZ2(I) = (STSOIL(I,K) - TGICE)                          &
    1297. 

  1297. 	   &                   / (.5 * ZSOIL(I,K-1) - .5 * ZSOIL(I,K))
    1298. 

  1298. 		  DFT2(I) = DFT1(I)
    1299. 

  1299. 		  CI(I,K) = 0.
    1300. 

  1300. 		ENDIF
    1301. 

  1301. 	        ENDDO
    1302. 

  1302. 	      ENDIF
    1303. 

  1303. 	      DO I = 1, IM
    1304. 

  1304. 	        IF(FLAG(I)) THEN
    1305. 

  1305. 		RHSTC(I,K) = (DFT2(I) * DTDZ2(I) - DFT1(I) * DTDZ1(I))      &
    1306. 

  1306. 	   &                 / ((ZSOIL(I,K) - ZSOIL(I,K-1)) * HCPCT(I))
    1307. 

  1307. 		AI(I,K) = -DFT1(I) * DDZ(I)                                 &
    1308. 

  1308. 	   &                / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
    1309. 

  1309. 		BI(I,K) = -(AI(I,K) + CI(I,K))
    1310. 

  1310. 		DFT1(I) = DFT2(I)
    1311. 

  1311. 		DTDZ1(I) = DTDZ2(I)
    1312. 

  1312. 		DDZ(I) = DDZ2(I)
    1313. 

  1313. 	        ENDIF
    1314. 

  1314. 	      ENDDO
    1315. 

  1315. 	    ENDDO
    1316. 

  1316.         !!
    1317. 

  1317.          700  CONTINUE
    1318. 

  1318.         !
    1319. 

  1319.         !  SOLVE THE TRI-DIAGONAL MATRIX
    1320. 

  1320.         !
    1321. 

  1321. 	    DO K = 1, KM
    1322. 

  1322. 	      DO I=1,IM
    1323. 

  1323. 	        IF(FLAG(I))  THEN
    1324. 

  1324. 		RHSTC(I,K) = RHSTC(I,K) * DELT2
    1325. 

  1325. 		AI(I,K) = AI(I,K) * DELT2
    1326. 

  1326. 		BI(I,K) = …
    1327.

Large files files are truncated, but you can click here to view the full file