PageRenderTime 67ms CodeModel.GetById 24ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/phys/module_sf_noahmpdrv.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1545 lines | 1028 code | 164 blank | 353 comment | 19 complexity | 8474cce8cbf9614dca04b90f4579e1e9 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. MODULE module_sf_noahmpdrv
    2. 

  2. 

  3. !-------------------------------
    4. 

  4.   USE module_sf_noahmplsm
    5. 

  5.   USE module_sf_urban
    6. 

  6.   USE module_model_constants, ONLY : R_D, CP, XLF, XLV, RHOWATER, KARMAN
    7. 

  7.   USE module_sf_noahdrv, ONLY : SOIL_VEG_GEN_PARM
    8. 

  8.   USE module_sf_noah_seaice
    9. 

  9.   USE module_sf_noahlsm_glacial_only
    10. 

  10.   USE MODULE_RA_GFDLETA, ONLY: CAL_MON_DAY
    11. 

  11. #ifdef WRF_CHEM
    12. 

  12.   USE module_data_gocart_dust
    13. 

  13. #endif
    14. 

  14. !-------------------------------
    15. 

  15. 

  16. !
    17. 

  17. CONTAINS
    18. 

  18. !
    19. 

  19.   SUBROUTINE noahmplsm(DZ8W,QV3D,P8W3D,T3D,TSK,                      &
    20. 

  20.        HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,GLW,SMSTAV,SMSTOT, &
    21. 

  21.        SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,VEGFRA,     &
    22. 

  22.        ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,XICE_THRESHOLD,ISICE,EMISS,EMBCK,   &
    23. 

  23.        SNOWC,QSFC,RAINBL,                            &
    24. 

  24.        num_soil_layers,DT,DZS,ITIMESTEP,             &
    25. 

  25.        SMOIS,TSLB,SNOW,CANWAT,                       &
    26. 

  26.        CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,qz0,       & !H
    27. 

  27.        myj,RIB,frpcpn,                               &
    28. 

  28.        SH2O,SNOWH,                                   & !H
    29. 

  29.        U_PHY,V_PHY,                                  & !I
    30. 

  30.        COSZ_URB2D, XLAT_URB2D,                       & !I
    31. 

  31.        SNOALB,                                       & !I
    32. 

  32.        SNOTIME,ACSNOM,ACSNOW,                        & !O
    33. 

  33.        idveg    ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz ,iopt_inf , &
    34. 

  34.        iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot,iopt_stc ,                     &
    35. 

  35.        isnowxy  ,tvxy     ,tgxy     ,canicexy ,                               &
    36. 

  36.        canliqxy ,eahxy    ,tahxy    ,cmxy     ,chxy     ,                     &
    37. 

  37.        fwetxy   ,sneqvoxy ,alboldxy ,qsnowxy  ,wslakexy ,zwtxy    ,waxy     , &
    38. 

  38.        wtxy     ,tsnoxy   ,zsnsoxy  ,snicexy  ,snliqxy  ,lfmassxy ,rtmassxy , &
    39. 

  39.        stmassxy ,woodxy   ,stblcpxy ,fastcpxy ,xlaixy   ,xsaixy   ,           &
    40. 

  40.        tradxy   ,tsxy     ,neexy    ,gppxy    ,nppxy    ,fvegxy   ,qinxy    , &
    41. 

  41.        runsfxy  ,runsbxy  ,ecanxy   ,edirxy   ,etranxy  ,fsaxy    ,firaxy   , &
    42. 

  42.        aparxy   ,psnxy    ,savxy    ,sagxy    ,                               &
    43. 

  43.        fsnoxy   ,YR       ,JULIAN   ,                                         &
    44. 

  44.        potevp,                                       & !O
    45. 

  45. 

  46. !jref:start
    47. 

  47.        qcxy     ,pblhxy   ,isurban  ,iz0tlnd  ,dx     , & !I
    48. 

  48.        chstarxy ,t2mvxy   ,t2mbxy   ,rssunxy  ,rsshaxy, bgapxy ,wgapxy ,gapxy,   & !O
    49. 

  49.        tgvxy    ,tgbxy    ,q2mvxy   ,q2mbxy  ,shdmaxxy, chvxy  ,chbxy  ,         &
    50. 

  50. !jref:end
    51. 

  51. 

  52.        ids,ide, jds,jde, kds,kde,                    &
    53. 

  53.        ims,ime, jms,jme, kms,kme,                    &
    54. 

  54.        its,ite, jts,jte, kts,kte                     )
    55. 

  55. !----------------------------------------------------------------
    56. 

  56.     IMPLICIT NONE
    57. 

  57. !----------------------------------------------------------------
    58. 

  58. !----------------------------------------------------------------
    59. 

  59. ! --- atmospheric (WRF generic) variables
    60. 

  60. !-- DT          time step (seconds)
    61. 

  61. !-- DZ8W        thickness of layers (m)
    62. 

  62. !-- T3D         temperature (K)
    63. 

  63. !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
    64. 

  64. !-- P3D         3D pressure (Pa)
    65. 

  65. !-- FLHC        exchange coefficient for heat (m/s)
    66. 

  66. !-- FLQC        exchange coefficient for moisture (m/s)
    67. 

  67. !-- PSFC        surface pressure (Pa)
    68. 

  68. !-- XLAND       land mask (1 for land, 2 for water)
    69. 

  69. !-- QGH         saturated mixing ratio at 2 meter
    70. 

  70. !-- GSW         downward short wave flux at ground surface (W/m^2)
    71. 

  71. !-- GLW         downward long wave flux at ground surface (W/m^2)
    72. 

  72. !-- History variables
    73. 

  73. !-- CANWAT      canopy moisture content (mm)
    74. 

  74. !-- TSK         surface temperature (K)
    75. 

  75. !-- TSLB        soil temp (k)
    76. 

  76. !-- SMOIS       total soil moisture content (volumetric fraction)
    77. 

  77. !-- SH2O        unfrozen soil moisture content (volumetric fraction)
    78. 

  78. !                note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
    79. 

  79. !-- SNOWH       actual snow depth (m)
    80. 

  80. !-- SNOW        liquid water-equivalent snow depth (m)
    81. 

  81. !-- ALBEDO      time-varying surface albedo including snow effect (unitless fraction)
    82. 

  82. !-- ALBBCK      background surface albedo (unitless fraction)
    83. 

  83. !-- CHS          surface exchange coefficient for heat and moisture (m s-1);
    84. 

  84. !-- CHS2        2m surface exchange coefficient for heat  (m s-1);
    85. 

  85. !-- CQS2        2m surface exchange coefficient for moisture (m s-1);
    86. 

  86. ! --- soil variables
    87. 

  87. !-- num_soil_layers   the number of soil layers
    88. 

  88. !-- ZS          depths of centers of soil layers   (m)
    89. 

  89. !-- DZS         thicknesses of soil layers (m)
    90. 

  90. !-- SLDPTH      thickness of each soil layer (m, same as DZS)
    91. 

  91. !-- TMN         soil temperature at lower boundary (K)
    92. 

  92. !-- SMCMAX      porosity, i.e. saturated value of soil moisture (volumetric)
    93. 

  93. !-- NROOT       number of root layers, a function of veg type, determined
    94. 

  94. !               in subroutine redprm.
    95. 

  95. !-- SMSTAV      Soil moisture availability for evapotranspiration (
    96. 

  96. !                   fraction between SMCWLT and SMCMXA)
    97. 

  97. !-- SMSTOT      Total soil moisture content frozen+unfrozen) in the soil column (mm)
    98. 

  98. ! --- snow variables
    99. 

  99. !-- SNOWC       fraction snow coverage (0-1.0)
    100. 

  100. ! --- vegetation variables
    101. 

  101. !-- SNOALB      upper bound on maximum albedo over deep snow
    102. 

  102. !-- Z0BRD       Background fixed roughness length (M)
    103. 

  103. !-- Z0          Background vroughness length (M) as function
    104. 

  104. !-- ZNT         Time varying roughness length (M) as function
    105. 

  105. !-- ALBD(IVGTPK,ISN) background albedo reading from a table
    106. 

  106. ! --- LSM output
    107. 

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

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

  109. !-- LH          upward moisture flux at the surface (W m-2)
    110. 

  110. !-- GRDFLX(I,J) ground heat flux (W m-2)
    111. 

  111. !-- FDOWN       radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
    112. 

  112. !----------------------------------------------------------------------------
    113. 

  113. !-- EC          canopy water evaporation ((W m-2)
    114. 

  114. !-- EDIR        direct soil evaporation (W m-2)
    115. 

  115. !-- ESNOW       sublimation from (or deposition to if <0) snowpack (W m-2)
    116. 

  116. !-- DEW         dewfall (or frostfall for t<273.15) (M)
    117. 

  117. ! ----------------------------------------------------------------------
    118. 

  118. !-- ETP         potential evaporation (W m-2)
    119. 

  119. ! ----------------------------------------------------------------------
    120. 

  120. !-- FLX1        precip-snow sfc (W m-2)
    121. 

  121. !-- FLX2        freezing rain latent heat flux (W m-2)
    122. 

  122. !-- FLX3        phase-change heat flux from snowmelt (W m-2)
    123. 

  123. ! ----------------------------------------------------------------------
    124. 

  124. !-- ACSNOM      snow melt (mm) (water equivalent)
    125. 

  125. !-- ACSNOW      accumulated snow fall (mm) (water equivalent)
    126. 

  126. !-- POTEVP      accumulated potential evaporation (W/m^2)
    127. 

  127. !-- RIB         Bulk Richardson number from SFCLAY routine
    128. 

  128. ! ----------------------------------------------------------------------
    129. 

  129. !-- RUNOFF1     surface runoff (m s-1), not infiltrating the surface
    130. 

  130. !-- RUNOFF2     subsurface runoff (m s-1), drainage out bottom of last
    131. 

  131. !                  soil layer (baseflow)
    132. 

  132. ! ----------------------------------------------------------------------
    133. 

  133. !-- RC          canopy resistance (s m-1)
    134. 

  134. !-- PC          plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
    135. 

  135. !-- EMISS       surface emissivity (between 0 and 1)
    136. 

  136. !-- EMBCK       Background surface emissivity (between 0 and 1)
    137. 

  137. !-- SHDMAX      Maximum vegetation fraction
    138. 

  138. 

  139. !-- ROVCP       R/CP
    140. 

  140. !               (R_d/R_v) (dimensionless)
    141. 

  141. !-- ids         start index for i in domain
    142. 

  142. !-- ide         end index for i in domain
    143. 

  143. !-- jds         start index for j in domain
    144. 

  144. !-- jde         end index for j in domain
    145. 

  145. !-- kds         start index for k in domain
    146. 

  146. !-- kde         end index for k in domain
    147. 

  147. !-- ims         start index for i in memory
    148. 

  148. !-- ime         end index for i in memory
    149. 

  149. !-- jms         start index for j in memory
    150. 

  150. !-- jme         end index for j in memory
    151. 

  151. !-- kms         start index for k in memory
    152. 

  152. !-- kme         end index for k in memory
    153. 

  153. !-- its         start index for i in tile
    154. 

  154. !-- ite         end index for i in tile
    155. 

  155. !-- jts         start index for j in tile
    156. 

  156. !-- jte         end index for j in tile
    157. 

  157. !-- kts         start index for k in tile
    158. 

  158. !-- kte         end index for k in tile
    159. 

  159. !
    160. 

  160. !-- SR          fraction of frozen precip (0.0 to 1.0)
    161. 

  161. !----------------------------------------------------------------
    162. 

  162. 

  163. ! IN only
    164. 

  164. 

  165.     INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &
    166. 

  166.          &                           ims,ime, jms,jme, kms,kme,  &
    167. 

  167.          &                           its,ite, jts,jte, kts,kte
    168. 

  168.     
    169. 

  169.     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
    170. 

  170.     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
    171. 

  171.     REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
    172. 

  172.     REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
    173. 

  173. 

  174. 

  175.     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    176. 

  176.          &   INTENT(IN   )    ::                            TMN, &
    177. 

  177.          &                                                XLAND, &
    178. 

  178.          &                                                 XICE, &
    179. 

  179.          &                                               VEGFRA, &
    180. 

  180.          &                                               SNOALB, &
    181. 

  181.          &                                                  GSW, &
    182. 

  182.          &                                               SWDOWN, &
    183. 

  183.          &                                                  GLW, &
    184. 

  184.          &                                                   Z0, &
    185. 

  185.          &                                               ALBBCK, &
    186. 

  186.          &                                               RAINBL, &
    187. 

  187.          &                                                EMBCK, &
    188. 

  188.          &                                                   SR
    189. 

  189. 

  190. 

  191.     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme )            , &
    192. 

  192.          INTENT(IN   )    ::                           QV3D, &
    193. 

  193.          p8w3D, &
    194. 

  194.          DZ8W, &
    195. 

  195.          T3D
    196. 

  196. !jref:start - changed to inout        
    197. 

  197.     REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    198. 

  198.          INTENT(INOUT   )               ::               QGH,  &
    199. 

  199.          CHS,   &
    200. 

  200.          CPM
    201. 

  201. !jref:end         
    202. 

  202. 

  203.     INTEGER, DIMENSION( ims:ime, jms:jme )                     , &
    204. 

  204.          INTENT(IN   )    ::                         IVGTYP, &
    205. 

  205.          ISLTYP
    206. 

  206. 

  207.     INTEGER, INTENT(IN)       ::     num_soil_layers,ITIMESTEP
    208. 

  208. !jref:start - xice_threshold
    209. 

  209.     REAL,     INTENT(IN   )   ::     DT,ROVCP,XICE_THRESHOLD
    210. 

  210.     INTEGER,  INTENT(IN   )   ::     ISICE
    211. 

  211. !jref:end    
    212. 

  212. 

  213.     REAL,     DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
    214. 

  214. 

  215. ! IN and OUT
    216. 

  216. 

  217.     REAL,     DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
    218. 

  218.          &    INTENT(INOUT)   ::                          SMOIS, &
    219. 

  219.          &                                                 SH2O, &
    220. 

  220.          &                                                 TSLB
    221. 

  221. 

  222.     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    223. 

  223.          &   INTENT(INOUT)    ::                            TSK, &
    224. 

  224.          &                                                  HFX, &
    225. 

  225.          &                                                  QFX, &
    226. 

  226.          &                                                   LH, &
    227. 

  227.          &                                               GRDFLX, &
    228. 

  228.          &                                                 QSFC, &
    229. 

  229.          &                                                 CQS2, &
    230. 

  230.          &                                                 CHS2, &
    231. 

  231.          &                                                 SNOW, &
    232. 

  232.          &                                                SNOWC, &
    233. 

  233.          &                                                SNOWH, &
    234. 

  234.          &                                               CANWAT, &
    235. 

  235.          &                                               SMSTAV, &
    236. 

  236.          &                                               SMSTOT, &
    237. 

  237.          &                                            SFCRUNOFF, &
    238. 

  238.          &                                             UDRUNOFF, &
    239. 

  239.          &                                               ACSNOM, &
    240. 

  240.          &                                               ACSNOW, &
    241. 

  241.          &                                                EMISS, &
    242. 

  242.          &                                               POTEVP, &
    243. 

  243.          &                                                  RIB, &
    244. 

  244.          &                                               ALBEDO, &
    245. 

  245.          &                                                  ZNT
    246. 

  246. 

  247.     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    248. 

  248.          INTENT(OUT)    ::                        CHKLOWQ
    249. 

  249.     REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) ::        QZ0
    250. 

  250. 

  251. !niuin: 
    252. 

  252. ! in
    253. 

  253.     INTEGER,  INTENT(IN) :: idveg     !dynamic vegetation (1 -> off ; 2 -> on) with opt_crs = 1      
    254. 

  254.     INTEGER,  INTENT(IN) :: iopt_crs  !canopy stomatal resistance (1-> Ball-Berry; 2->Jarvis)
    255. 

  255.     INTEGER,  INTENT(IN) :: iopt_btr  !soil moisture factor for stomatal resistance (1-> Noah; 2-> CLM; 3-> SSiB)
    256. 

  256.     INTEGER,  INTENT(IN) :: iopt_run  !runoff and groundwater (1->SIMGM; 2->SIMTOP; 3->Schaake96; 4->BATS)
    257. 

  257.     INTEGER,  INTENT(IN) :: iopt_sfc  !surface layer drag coeff (CH & CM) (1->M-O; 2->Chen97)
    258. 

  258.     INTEGER,  INTENT(IN) :: iopt_frz  !supercooled liquid water (1-> NY06; 2->Koren99)
    259. 

  259.     INTEGER,  INTENT(IN) :: iopt_inf  !frozen soil permeability (1-> NY06; 2->Koren99)
    260. 

  260.     INTEGER,  INTENT(IN) :: iopt_rad  !radiation transfer (1->gap=F(3D,cosz); 2->gap=0; 3->gap=1-Fveg)
    261. 

  261.     INTEGER,  INTENT(IN) :: iopt_alb  !snow surface albedo (1->BATS; 2->CLASS)
    262. 

  262.     INTEGER,  INTENT(IN) :: iopt_snf  !rainfall & snowfall (1-Jordan91; 2->BATS; 3->Noah)
    263. 

  263.     INTEGER,  INTENT(IN) :: iopt_tbot !lower boundary of soil temperature (1->zero-flux; 2->Noah)
    264. 

  264.     INTEGER,  INTENT(IN) :: iopt_stc  !snow/soil temperature time scheme
    265. 

  265. 

  266. ! in & out
    267. 

  267.     INTEGER, INTENT(IN) :: YR
    268. 

  268.     REAL,    INTENT(IN) :: JULIAN
    269. 

  269.     INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy     !actual no. of snow layers
    270. 

  270.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy        !vegetation canopy temperature
    271. 

  271.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy        !ground surface temperature
    272. 

  272.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy    !canopy-intercepted ice (mm)
    273. 

  273.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy    !canopy-intercepted liquid water (mm)
    274. 

  274.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy       !canopy air vapor pressure (pa)
    275. 

  275.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy       !canopy air temperature (k)
    276. 

  276.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy        !momentum drag coefficient
    277. 

  277.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy        !sensible heat exchange coefficient
    278. 

  278.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy      !wetted or snowed fraction of the canopy (-)
    279. 

  279.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy    !snow mass at last time step(mm h2o)
    280. 

  280.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy    !snow albedo at last time step (-)
    281. 

  281.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy     !snowfall on the ground [mm/s]
    282. 

  282.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy    !lake water storage [mm]
    283. 

  283. 

  284.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy       !water table depth [m]
    285. 

  285.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy        !water in the "aquifer" [mm]
    286. 

  286.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy        !groundwater storage [mm]
    287. 

  287.     REAL, DIMENSION(ims:ime,-2:num_soil_layers,jms:jme), INTENT(INOUT) :: zsnsoxy  !snow layer depth [m]
    288. 

  288.     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: tsnoxy   !snow temperature [K]
    289. 

  289.     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snicexy  !snow layer ice [mm]
    290. 

  290.     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snliqxy  !snow layer liquid water [mm]
    291. 

  291.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy    !leaf mass [g/m2]
    292. 

  292.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy    !mass of fine roots [g/m2]
    293. 

  293.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy    !stem mass [g/m2]
    294. 

  294.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy      !mass of wood (incl. woody roots) [g/m2]
    295. 

  295.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy    !stable carbon in deep soil [g/m2]
    296. 

  296.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy    !short-lived carbon, shallow soil [g/m2]
    297. 

  297.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xlaixy      !leaf area index
    298. 

  298.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy      !stem area index
    299. 

  299. !jref:start    
    300. 

  300.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: SNOTIME      !snow age time
    301. 

  301. !jref:end    
    302. 

  302. 

  303. !out
    304. 

  304.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: tradxy      !surface radiative temperature (k)
    305. 

  305.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: tsxy        !surface temperature (k)
    306. 

  306.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: neexy       !net ecosys exchange (g/m2/s CO2)
    307. 

  307.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: gppxy       !gross primary assimilation [g/m2/s C]
    308. 

  308.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: nppxy       !net primary productivity [g/m2/s C]
    309. 

  309.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: fvegxy      !greenness vegetation fraction [-]
    310. 

  310.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: qinxy       !groundwater recharge [mm/s]
    311. 

  311.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: runsfxy     !surface runoff [mm/s]
    312. 

  312.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: runsbxy     !subsurface runoff [mm/s]
    313. 

  313.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: ecanxy      !evaporation of intercepted water (mm/s)
    314. 

  314.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: edirxy      !soil surface evaporation rate (mm/s]
    315. 

  315.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: etranxy     !transpiration rate (mm/s)
    316. 

  316.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: fsaxy       !total absorbed solar radiation (w/m2)
    317. 

  317.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: firaxy      !total net longwave rad (w/m2) [+ to atm]
    318. 

  318.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: aparxy      !photosyn active energy by canopy (w/m2)
    319. 

  319.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: psnxy       !total photosynthesis (umol co2/m2/s) [+]
    320. 

  320.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: savxy       !solar rad absorbed by veg. (w/m2)
    321. 

  321.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: sagxy       !solar rad absorbed by ground (w/m2)
    322. 

  322.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: fsnoxy      !snow cover fraction (-)
    323. 

  323. 

  324. !jref:start 
    325. 

  325.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: chstarxy    !effective ch
    326. 

  326.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: t2mvxy      !2m temperature of vegetation part
    327. 

  327.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: t2mbxy      !2m temperature of bare ground part
    328. 

  328.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: q2mvxy      !2m mixing ratio of vegetation part
    329. 

  329.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: q2mbxy      !2m mixing ratio of bare ground part
    330. 

  330. 

  331.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN)   :: qcxy        !cloud water mixing ratio
    332. 

  332.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN)   :: pblhxy      !Planetary boundary layer from sfclay
    333. 

  333.     INTEGER                          , INTENT(IN)   :: isurban
    334. 

  334.     INTEGER                          , INTENT(IN)   :: iz0tlnd
    335. 

  335.     REAL                             , INTENT(IN)   :: dx
    336. 

  336.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: rssunxy        !sunlit leaf stomatal resistance (s/m)
    337. 

  337.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: rsshaxy        !shaded leaf stomatal resistance (s/m)
    338. 

  338.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: bgapxy        !between gap fraction
    339. 

  339.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: wgapxy        !within gap fraction
    340. 

  340.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: gapxy        !within gap fraction
    341. 

  341.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tgvxy
    342. 

  342.     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tgbxy        
    343. 

  343.     REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   shdmaxxy
    344. 

  344.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: chvxy        !sensible heat exchange coefficient vegetated
    345. 

  345.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: chbxy        !sensible heat exchange coefficient bare-ground
    346. 

  346. !jref:end 
    347. 

  347. 

  348. 

  349. !niuout
    350. 

  350. 

  351. ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
    352. 

  352. 

  353.     INTEGER :: YEARLEN
    354. 

  354. 

  355.     REAL  ::  ETP, SSOIL,EC, ESNOW,             &
    356. 

  356.          FLX1,FLX2,FLX3,DEW,FDOWN,RC,PC,FFROZP
    357. 

  357. 

  358. !niuin
    359. 

  359. !locals (prognostic):
    360. 

  360. 

  361.     INTEGER                             :: isnow        !actual no. of snow layers
    362. 

  362.     REAL, DIMENSION(-2:num_soil_layers) :: stc          !snow/soil tmperatures
    363. 

  363.     REAL, DIMENSION( 1:num_soil_layers) :: smc          !vol. soil moisture (m3/m3)
    364. 

  364.     REAL, DIMENSION( 1:num_soil_layers) :: smh2o        !vol. soil liquid water (m3/m3)
    365. 

  365.     REAL                                :: tv           !vegetation canopy temperature
    366. 

  366.     REAL                                :: tg           !ground surface temperature
    367. 

  367.     REAL                                :: canice       !canopy-intercepted ice (mm)
    368. 

  368.     REAL                                :: canliq       !canopy-intercepted liquid water (mm)
    369. 

  369.     REAL                                :: snowd        !snow depth (m)
    370. 

  370.     REAL                                :: swe          !snow water equivalent (mm)
    371. 

  371.     REAL                                :: eah          !canopy air vapor pressure (pa)
    372. 

  372.     REAL                                :: tah          !canopy air temperature (k)
    373. 

  373.     REAL                                :: cm           !momentum drag coefficient
    374. 

  374.     REAL                                :: ch           !sensible heat exchange coefficient
    375. 

  375.     REAL                                :: fwet         !wetted or snowed fraction of the canopy (-)
    376. 

  376.     REAL                                :: sneqvo       !snow mass at last time step(mm h2o)
    377. 

  377.     REAL                                :: albold       !snow albedo at last time step (-)
    378. 

  378.     REAL                                :: qsnow        !snowfall on the ground [mm/s]
    379. 

  379.     REAL                                :: wslake       !lake water storage [mm]
    380. 

  380. 

  381.     REAL                                :: zwt          !water table depth [m]
    382. 

  382.     REAL                                :: wa           !water in the "aquifer" [mm]
    383. 

  383.     REAL                                :: wt           !groundwater storage [mm]
    384. 

  384.     REAL, DIMENSION(-2:num_soil_layers) :: zsnso        !snow layer depth [m]
    385. 

  385.     REAL, DIMENSION(-2:              0) :: tsno         !snow temperature [K]
    386. 

  386.     REAL, DIMENSION(-2:              0) :: snice        !snow layer ice [mm]
    387. 

  387.     REAL, DIMENSION(-2:              0) :: snliq        !snow layer liquid water [mm]
    388. 

  388.     REAL                                :: lfmass       !leaf mass [g/m2]
    389. 

  389.     REAL                                :: rtmass       !mass of fine roots [g/m2]
    390. 

  390.     REAL                                :: stmass       !stem mass [g/m2]
    391. 

  391.     REAL                                :: wood         !mass of wood (incl. woody roots) [g/m2]
    392. 

  392.     REAL                                :: stblcp       !stable carbon in deep soil [g/m2]
    393. 

  393.     REAL                                :: fastcp       !short-lived carbon, shallow soil [g/m2]
    394. 

  394.     REAL                                :: plai         !leaf area index
    395. 

  395.     REAL                                :: psai         !stem area index
    396. 

  396. 

  397. !jref:start
    398. 

  398.     REAL                                :: chstar2
    399. 

  399.     REAL                                :: cqstar2
    400. 

  400.     REAL                                :: chstar       !effective ch
    401. 

  401.     REAL                                :: tstar
    402. 

  402.     REAL                                :: t2mv         !2m temperature of vegetation part
    403. 

  403.     REAL                                :: t2mb         !2m temperature of bare ground part
    404. 

  404.     REAL                                :: q2mv         !2m mixing ratio of vegetation part
    405. 

  405.     REAL                                :: q2mb         !2m mixing ratio of bare ground part
    406. 

  406.     REAL                                :: qc           !
    407. 

  407.     REAL                                :: t2m
    408. 

  408.     REAL                                :: pblh 
    409. 

  409.     REAL                                :: qsfc1d
    410. 

  410.     REAL, DIMENSION(ims:ime,jms:jme)    :: tstarxy      !effective skin temperature
    411. 

  411.     REAL, DIMENSION(ims:ime,jms:jme)    :: chstar2xy    !effective 2m exchange coefficients
    412. 

  412.     REAL                                :: rssun
    413. 

  413.     REAL                                :: rssha
    414. 

  414.     REAL                                :: bgap
    415. 

  415.     REAL                                :: wgap
    416. 

  416.     REAL                                :: gap
    417. 

  417.     REAL                                :: tgv
    418. 

  418.     REAL                                :: tgb
    419. 

  419.     REAL                                :: snowhk
    420. 

  420.     REAL                                :: snotime1
    421. 

  421.     REAL                                :: qv1d         !mixing ratio
    422. 

  422.     REAL                                :: dz8w1d
    423. 

  423.     REAL                                :: shdmax
    424. 

  424.     REAL                                :: chv          !sensible heat exchange coefficient vegetated
    425. 

  425.     REAL                                :: chb          !sensible heat exchange coefficient bare-ground
    426. 

  426. !jref:end
    427. 

  427. 

  428. !out (outputs)
    429. 

  429. 

  430.     REAL                                :: trad         !surface radiative temperature (k)
    431. 

  431.     REAL                                :: ts           !surface temperature (k)
    432. 

  432.     REAL                                :: nee          !net ecosys exchange (g/m2/s CO2)
    433. 

  433.     REAL                                :: gpp          !gross primary assimilation [g/m2/s C]
    434. 

  434.     REAL                                :: npp          !net primary productivity [g/m2/s C]
    435. 

  435.     REAL                                :: fveg         !greenness vegetation fraction [-]
    436. 

  436.     REAL                                :: qin          !groundwater recharge [mm/s]
    437. 

  437.     REAL                                :: runsf        !surface runoff [mm/s]
    438. 

  438.     REAL                                :: runsb        !subsurface runoff [mm/s]
    439. 

  439.     REAL                                :: ecan         !evaporation of intercepted water (mm/s)
    440. 

  440.     REAL                                :: esoil        !soil surface evaporation rate (mm/s]
    441. 

  441.     REAL                                :: etran        !transpiration rate (mm/s)
    442. 

  442.     REAL                                :: fsa          !total absorbed solar radiation (w/m2)
    443. 

  443.     REAL                                :: fira         !total net longwave rad (w/m2) [+ to atm]
    444. 

  444.     REAL                                :: fsh          !total sensible heat (w/m2) [+ to atm]
    445. 

  445.     REAL                                :: flh          !total latent heat (w/m2) [+ to atm]
    446. 

  446.     REAL                                :: apar         !photosyn active energy by canopy (w/m2)
    447. 

  447.     REAL                                :: psn          !total photosynthesis (umol co2/m2/s) [+]
    448. 

  448.     REAL                                :: sav          !solar rad absorbed by veg. (w/m2)
    449. 

  449.     REAL                                :: sag          !solar rad absorbed by ground (w/m2)
    450. 

  450.     REAL                                :: fsno         !snow cover fraction (-)
    451. 

  451.     REAL                                :: salb         !surface albedo (-)
    452. 

  452.     REAL                                :: errwat
    453. 

  453.     REAL                                :: qmelt
    454. 

  454.     REAL                                :: ponding
    455. 

  455.     REAL                                :: ponding1
    456. 

  456.     REAL                                :: ponding2
    457. 

  457. 

  458. !local
    459. 

  459. 

  460.     real                                :: fsr          !total reflected solar radiation (w/m2)
    461. 

  461.     real                                :: fcev         !canopy evaporation heat (w/m2) [+ to atm]
    462. 

  462.     real                                :: fgev         !ground evaporation heat (w/m2) [+ to atm]
    463. 

  463.     real                                :: fctr         !transpiration heat flux (w/m2) [+ to atm]
    464. 

  464.     real, dimension(-2:              0) :: ficeold      !snow layer liquid water [mm]
    465. 

  465. 

  466.     INTEGER :: ILOC      !grid index
    467. 

  467.     INTEGER :: JLOC      !grid index
    468. 

  468.     INTEGER :: ISC       !soil color index
    469. 

  469.     INTEGER :: IST       !surface type 1-soil; 2-lake
    470. 

  470. 

  471. !niuout
    472. 

  472. 

  473.     LOGICAL,    INTENT(IN   )    ::     myj,frpcpn
    474. 

  474. 

  475. ! DECLARATIONS - LOGICAL
    476. 

  476. ! ----------------------------------------------------------------------
    477. 

  477.     LOGICAL, PARAMETER :: LOCAL=.false.
    478. 

  478.     LOGICAL :: FRZGRA, SNOWNG
    479. 

  479. 

  480.     LOGICAL :: IPRINT
    481. 

  481. 

  482. ! ----------------------------------------------------------------------
    483. 

  483. ! DECLARATIONS - INTEGER
    484. 

  484. ! ----------------------------------------------------------------------
    485. 

  485.     INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
    486. 

  486.     INTEGER :: NROOT
    487. 

  487.     INTEGER :: KZ ,K
    488. 

  488.     INTEGER :: NS
    489. 

  489. ! ----------------------------------------------------------------------
    490. 

  490. ! DECLARATIONS - REAL
    491. 

  491. ! ----------------------------------------------------------------------
    492. 

  492. 

  493.     REAL  :: DQSDT2, LWDN, PRCP, PSFC, UU, VV, CO2AIR, O2AIR,         &
    494. 

  494.          &   Q2SAT,Q2SATI,SFCPRS,SFCTMP,SHDFAC,SNOALB1,               &
    495. 

  495.          &   SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ETA, ETA_KINEMATIC,         &
    496. 

  496.          &   EMBRD, FOLN, LAT,                                        &
    497. 

  497.          &   Z0K,RUNOFF1,RUNOFF2,SOLNET,E2SAT,SFCTSNO
    498. 

  498. 

  499.     REAL :: RIBB
    500. 

  500.     REAL ::  FDTW
    501. 

  501. 

  502.     REAL  :: EMISSI
    503. 

  503. 

  504.     REAL  :: SNCOVR,SNEQV,CHK,TH2
    505. 

  505. 

  506.     REAL  :: SMCMAX,SNOMLT,SOILM,SOILW,Q1,T1
    507. 

  507. 

  508.     REAL  :: Z0BRD
    509. 

  509. !
    510. 

  510.     REAL  :: COSZ
    511. 

  511. !
    512. 

  512. !niu      REAL, DIMENSION(1:num_soil_layers)::  SLDPTH, STC,SMC,SWC
    513. 

  513.     REAL, DIMENSION(1:num_soil_layers)::  SLDPTH,SWC
    514. 

  514. 

  515. !jref:start    
    516. 

  516.     REAL, DIMENSION(1:num_soil_layers):: STCNEW
    517. 

  517. !jref:end     
    518. 

  518. !
    519. 

  519.     REAL, DIMENSION(1:num_soil_layers) ::     ZSOIL, RTDIS
    520. 

  520.     REAL, PARAMETER  :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65,   &
    521. 

  521.          T0=273.16E0, ELWV=2.50E6,  A23M4=A2*(A3-A4)
    522. 

  522. 

  523.     ! Used for calculating the 2-m Potential Temperature:
    524. 

  524.     REAL, PARAMETER                                       :: CAPA=R_D/CP
    525. 

  525.     REAL                                                  :: APELM
    526. 

  526.     REAL                                                  :: APES
    527. 

  527.     REAL                                                  :: SFCTH2
    528. 

  528. 

  529. ! ----------------------------------------------------------------------
    530. 

  530. ! ----------------------------------------------------------------------
    531. 

  531. 

  532. ! MEK JUL2007
    533. 

  533.     FDTW=DT/(XLV*RHOWATER)
    534. 

  534. ! debug printout
    535. 

  535.     IPRINT=.false.
    536. 

  536. 

  537. !      SLOPETYP=2
    538. 

  538.     SLOPETYP=1
    539. 

  539. !      SHDMIN=0.00
    540. 

  540. 

  541.     YEARLEN = 365
    542. 

  542.     if (mod(YR,4) == 0) then
    543. 

  543.        YEARLEN = 366
    544. 

  544.        if (mod(YR,100) == 0) then
    545. 

  545.           YEARLEN = 365
    546. 

  546.           if (mod(YR,400) == 0) then
    547. 

  547.              YEARLEN = 366
    548. 

  548.           endif
    549. 

  549.        endif
    550. 

  550.     endif
    551. 

  551. 

  552.     NSOIL=num_soil_layers
    553. 

  553. 

  554.     DO NS=1,NSOIL
    555. 

  555.        SLDPTH(NS)=DZS(NS)
    556. 

  556.     ENDDO
    557. 

  557. 

  558.     call noahmp_options(idveg     ,iopt_crs  ,iopt_btr  ,iopt_run  ,iopt_sfc  ,iopt_frz , &
    559. 

  559.          iopt_inf  ,iopt_rad  ,iopt_alb  ,iopt_snf  ,iopt_tbot, iopt_stc )
    560. 

  560. 

  561.     ISC = 4                  ! soil color: assuming a middle color category ?????????
    562. 

  562. 

  563.     ZSOIL(1) = -SLDPTH(1)    ! move out of  x-y do loops
    564. 

  564.     DO KZ = 2, NSOIL
    565. 

  565.        ZSOIL(KZ)         = -SLDPTH(KZ) + ZSOIL(KZ-1)
    566. 

  566.     END DO
    567. 

  567.     FOLN                 = 1.0
    568. 

  568. 

  569. !niuout
    570. 

  570. 

  571.     DO J=jts,jte
    572. 

  572. 

  573.        IF(ITIMESTEP.EQ.1)THEN
    574. 

  574.           DO I=its,ite
    575. 

  575. 

  576. !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
    577. 

  577.              IF((XLAND(I,J)-1.5).GE.0.)THEN
    578. 

  578. ! check sea-ice point
    579. 

  579.                 IF(XICE(I,J).EQ.1..and.IPRINT)PRINT*,' sea-ice at water point, I=',I,  &
    580. 

  580.                      'J=',J
    581. 

  581. !***   Open Water Case
    582. 

  582.                 SMSTAV(I,J)=1.0
    583. 

  583.                 SMSTOT(I,J)=1.0
    584. 

  584.                 DO NS=1,NSOIL
    585. 

  585.                    SMOIS(I,NS,J)=1.0
    586. 

  586.                    TSLB(I,NS,J)=273.16                                          !STEMP
    587. 

  587.                 ENDDO
    588. 

  588.              ELSE
    589. 

  589.                 IF(XICE(I,J).EQ.1.)THEN
    590. 

  590. !***        SEA-ICE CASE
    591. 

  591.                    SMSTAV(I,J)=1.0
    592. 

  592.                    SMSTOT(I,J)=1.0
    593. 

  593.                    DO NS=1,NSOIL
    594. 

  594.                       SMOIS(I,NS,J)=1.0
    595. 

  595.                    ENDDO
    596. 

  596.                 ENDIF
    597. 

  597.              ENDIF
    598. 

  598. !
    599. 

  599.           ENDDO
    600. 

  600.        ENDIF                                                               ! end of initialization over ocean
    601. 

  601. 

  602. 

  603. !-----------------------------------------------------------------------
    604. 

  604.        DO I=its,ite
    605. 

  605. ! surface pressure
    606. 

  606.           PSFC=P8w3D(i,1,j)
    607. 

  607. ! pressure in middle of lowest layer
    608. 

  608.           SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
    609. 

  609. ! convert from mixing ratio to specific humidity
    610. 

  610.           Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
    611. 

  611. !
    612. 

  612. !         Q2SAT=QGH(I,j)
    613. 

  613.           Q2SAT=QGH(I,J)/(1.0+QGH(I,J))        ! Q2SAT is sp humidity
    614. 

  614. ! add check on myj=.true.
    615. 

  615. !        IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
    616. 

  616.           IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
    617. 

  617.              CHKLOWQ(I,J)=0.
    618. 

  618.           ELSE
    619. 

  619.              CHKLOWQ(I,J)=1.
    620. 

  620.           ENDIF
    621. 

  621. 

  622.           SFCTMP=T3D(i,1,j)
    623. 

  623.           ZLVL=0.5*DZ8W(i,1,j)
    624. 

  624. 

  625. !        TH2=SFCTMP+(0.0097545*ZLVL)
    626. 

  626. ! calculate SFCTH2 via Exner function vs lapse-rate (above)
    627. 

  627.           APES=(1.E5/PSFC)**CAPA
    628. 

  628.           APELM=(1.E5/SFCPRS)**CAPA
    629. 

  629.           SFCTH2=SFCTMP*APELM
    630. 

  630.           TH2=SFCTH2/APES
    631. 

  631. !
    632. 

  632.           EMISSI = EMISS(I,J)
    633. 

  633. !          LWDN=GLW(I,J)*EMISSI
    634. 

  634.           LWDN=GLW(I,J)
    635. 

  635. ! SOLDN is total incoming solar
    636. 

  636.           SOLDN=SWDOWN(I,J)
    637. 

  637. ! GSW is net downward solar
    638. 

  638. !        SOLNET=GSW(I,J)
    639. 

  639. ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
    640. 

  640.           SOLNET=SOLDN*(1.-ALBEDO(I,J))
    641. 

  641.           PRCP=RAINBL(i,j)/DT
    642. 

  642.           VEGTYP=IVGTYP(I,J)
    643. 

  643.           SOILTYP=ISLTYP(I,J)
    644. 

  644.           SHDFAC=VEGFRA(I,J)/100.
    645. 

  645.           T1=TSK(I,J)
    646. 

  646.           CHK=CHS(I,J)
    647. 

  647.           SNOALB1=SNOALB(I,J)     !NEW
    648. 

  648. 

  649. ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
    650. 

  650. ! SR from e.g. Ferrier microphysics
    651. 

  651. ! otherwise define from 1st atmos level temperature
    652. 

  652.           IF(FRPCPN) THEN
    653. 

  653.              FFROZP=SR(I,J)
    654. 

  654.           ELSE
    655. 

  655.              IF (SFCTMP <=  273.15) THEN
    656. 

  656.                 FFROZP = 1.0
    657. 

  657.              ELSE
    658. 

  658.                 FFROZP = 0.0
    659. 

  659.              ENDIF
    660. 

  660.           ENDIF
    661. 

  661. !***
    662. 

  662.           IF((XLAND(I,J)-1.5).GE.0.)THEN                                  ! begining of land/sea if block
    663. 

  663. ! Open water points
    664. 

  664.           ELSE
    665. 

  665. ! Land or sea-ice case
    666. 

  666. 

  667.              IF (XICE(I,J) .GT. 0.5) THEN
    668. 

  668.                 ICE=1
    669. 

  669.              ELSE
    670. 

  670.                 ICE=0
    671. 

  671.              ENDIF
    672. 

  672.              DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
    673. 

  673. 

  674.              IF(SNOW(I,J).GT.0.0)THEN
    675. 

  675. ! snow on surface (use ice saturation properties)
    676. 

  676.                 SFCTSNO=SFCTMP
    677. 

  677.                 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
    678. 

  678.                 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
    679. 

  679.                 Q2SATI=Q2SATI/(1.0+Q2SATI)    ! spec. hum.
    680. 

  680.                 IF(T1 .GT. 273.15)THEN
    681. 

  681. ! warm ground temps, weight the saturation between ice and water according to SNOWC
    682. 

  682.                    Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
    683. 

  683.                    DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
    684. 

  684.                 ELSE
    685. 

  685. ! cold ground temps, use ice saturation only
    686. 

  686.                    Q2SAT=Q2SATI
    687. 

  687.                    DQSDT2=Q2SATI*6174./(SFCTSNO**2)
    688. 

  688.                 ENDIF
    689. 

  689. ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
    690. 

  690.                 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
    691. 

  691.              ENDIF
    692. 

  692. 

  693.              IF(ICE.EQ.0)THEN
    694. 

  694.                 TBOT=TMN(I,J)
    695. 

  695.              ELSE
    696. 

  696.                 TBOT=271.16
    697. 

  697.              ENDIF
    698. 

  698.              IF(VEGTYP.EQ.25) SHDFAC=0.0000
    699. 

  699.              IF(VEGTYP.EQ.26) SHDFAC=0.0000
    700. 

  700.              IF(VEGTYP.EQ.27) SHDFAC=0.0000
    701. 

  701.              IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
    702. 

  702.                 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
    703. 

  703.                 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
    704. 

  704.                 SOILTYP=7
    705. 

  705.              ENDIF
    706. 

  706. 

  707. !-------------------------------------------
    708. 

  708.              ALBBRD=ALBBCK(I,J)
    709. 

  709.              Z0BRD=Z0(I,J)
    710. 

  710.              EMBRD=EMBCK(I,J)
    711. 

  711. !jref:start - check if this is correct!! Maybe snowd
    712. 

  712.              RIBB=RIB(I,J)
    713. 

  713.              SNOTIME1 = SNOTIME(I,J)
    714. 

  714. !jref:end             
    715. 

  715. 

  716. !FEI: temporaray arrays above need to be changed later by using SI
    717. 

  717. 

  718. !niu          DO 70 NS=1,NSOIL
    719. 

  719. !niu            SMC(NS)=SMOIS(I,NS,J)
    720. 

  720. !niu            STC(NS)=TSLB(I,NS,J)                                          !STEMP
    721. 

  721. !niu            SWC(NS)=SH2O(I,NS,J)
    722. 

  722. !niu   70     CONTINUE
    723. 

  723. !
    724. 

  724. 

  725.              IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
    726. 

  726.                   IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
    727. 

  727.                 VEGTYP = ISURBAN
    728. 

  728.              ENDIF
    729. 

  729. 

  730.              IST    = 1 
    731. 

  731.              IF(VEGTYP == 16) IST = 2 ! lake points
    732. 

  732. 

  733.              CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,SLDPTH,ZSOIL,NSOIL,ISURBAN)
    734. 

  734. 

  735.              UU                   = U_PHY(I,1,J)
    736. 

  736.              VV                   = V_PHY(I,1,J)
    737. 

  737.              CO2AIR               = 395.E-06 * SFCPRS       !partial pressure co2 (pa)
    738. 

  738.              O2AIR                = 0.209    * SFCPRS       !partial pressure  o2 (pa)
    739. 

  739.              COSZ                 = COSZ_URB2D(I,J)
    740. 

  740.              LAT                  = XLAT_URB2D(I,J)
    741. 

  741. 

  742.              isnow                = isnowxy (i,j)
    743. 

  743.              stc  (isnow+1:    0) = tsnoxy  (i,isnow+1:    0,j)
    744. 

  744.              stc  (      1:nsoil) = tslb    (i,      1:nsoil,j)
    745. 

  745.              smc (       1:nsoil) = smois   (i,      1:nsoil,j)
    746. 

  746.              smh2o(      1:nsoil) = sh2o    (i,      1:nsoil,j)
    747. 

  747.              tv                   = tvxy    (i,j)
    748. 

  748.              tg                   = tgxy    (i,j)
    749. 

  749.              canliq               = canliqxy(i,j)
    750. 

  750.              canice               = canicexy(i,j)
    751. 

  751.              snowd                = snowh   (i,j)
    752. 

  752.              swe                  = snow    (i,j)
    753. 

  753.              eah                  = eahxy   (i,j)
    754. 

  754.              tah                  = tahxy   (i,j)
    755. 

  755.              cm                   = cmxy    (i,j)
    756. 

  756.              ch                   = chxy    (i,j)
    757. 

  757. !jref:start             
    758. 

  758.              chstar               = chs     (i,j)
    759. 

  759.              chstar2              = chs2    (i,j)
    760. 

  760.              cqstar2              = cqs2    (i,j)
    761. 

  761.              tstar                = T1
    762. 

  762.              qc                   = qcxy    (i,j)
    763. 

  763.              pblh                 = pblhxy  (i,j)
    764. 

  764.              qsfc1d               = qsfc    (i,j)
    765. 

  765.              t2mv                 = t2mvxy  (i,j)
    766. 

  766.              t2mb                 = t2mbxy  (i,j)
    767. 

  767.              q2mv                 = q2mvxy  (i,j)
    768. 

  768.              q2mb                 = q2mbxy  (i,j)
    769. 

  769.              qv1d                 = qv3d    (i,1,j) ! seaice/glacial needs mixing ratio (q2k = specific hum).             
    770. 

  770.              dz8w1d               = dz8w    (i,1,j)
    771. 

  771.              shdmax               = shdmaxxy (i,j)/100. !fraction
    772. 

  772. !jref:end             
    773. 

  773. 

  774.              fwet                 = fwetxy  (i,j)
    775. 

  775.              sneqvo               = sneqvoxy(i,j)
    776. 

  776.              albold               = alboldxy(i,j)
    777. 

  777.              qsnow                = qsnowxy (i,j)
    778. 

  778.              wslake               = wslakexy(i,j)
    779. 

  779. 

  780.              zwt                  = zwtxy   (i,j)
    781. 

  781.              wa                   = waxy    (i,j)
    782. 

  782.              wt                   = wtxy    (i,j)
    783. 

  783.              zsnso(isnow+1:nsoil) = zsnsoxy (i,isnow+1:nsoil,j)
    784. 

  784.              snice(isnow+1:    0) = snicexy (i,isnow+1:    0,j)
    785. 

  785.              snliq(isnow+1:    0) = snliqxy (i,isnow+1:    0,j)
    786. 

  786. 

  787.              lfmass               = lfmassxy(i,j)
    788. 

  788.              rtmass               = rtmassxy(i,j)
    789. 

  789.              stmass               = stmassxy(i,j)
    790. 

  790.              wood                 = woodxy  (i,j)
    791. 

  791.              stblcp               = stblcpxy(i,j)
    792. 

  792.              fastcp               = fastcpxy(i,j)
    793. 

  793.              plai                 = xlaixy  (i,j)
    794. 

  794.              psai                 = xsaixy  (i,j)
    795. 

  795. 

  796.              ficeold(isnow+1:0)   = snicexy(i,isnow+1:0,j) &
    797. 

  797.                   /(snicexy(i,isnow+1:0,j)+snliqxy(i,isnow+1:0,j))
    798. 

  798. 

  799. ! glacial, seaice split - jref
    800. 

  800. 

  801.     IF ( XICE(I,J) >= XICE_THRESHOLD ) THEN
    802. 

  802. 

  803.        SH2O  (i,1:nsoil,j) = 1.0
    804. 

  804.        XLAIXY(i,j)         = 0.01
    805. 

  805. 

  806.        cycle ! Skip any processing at sea-ice points
    807. 

  807. 

  808.     ELSE IF ( VEGTYP == ISICE ) THEN
    809. 

  809. 

  810.        SNCOVR = SNOWC(I,J)
    811. 

  811.        swe = swe*0.001 !jref mm -> m
    812. 

  812.        if ( (swe.ne.0..AND.snowd.eq.0.).or.(snowd.le.swe) )THEN
    813. 

  813.              snowd= 5.*swe
    814. 

  814.        endif
    815. 

  815. 

  816.        CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH,   &    !C
    817. 

  817.             &    LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,              &    !F
    818. 

  818.             &    TH2,Q2SAT,DQSDT2,                                &    !I
    819. 

  819.             &    ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, &    !S
    820. 

  820.             &    tstar,STC(1:NSOIL),snowd,swe,salb,chstar,        &    !H
    821. 

  821.             &    ETA,fsh, ETA_KINEMATIC,FDOWN,                    &    !O
    822. 

  822.             &    ESNOW,DEW,                                       &    !O
    823. 

  823.             &    ETP,SSOIL,                                       &    !O
    824. 

  824.             &    FLX1,FLX2,FLX3,                                  &    !O
    825. 

  825.             &    SNOMLT,SNCOVR,                                   &    !O
    826. 

  826.             &    runsf,                                         &    !O
    827. 

  827.             &    Q1,                                              &    !D
    828. 

  828.             &    SNOTIME1,                                        &
    829. 

  829.             &    RIBB)
    830. 

  830. 

  831.        tgb    = sfctmp ! Bare ground temperature will be the surface temperature over glacial points. 
    832. 

  832.        tgv    = 0.0    ! Temperature under vegetation undefined over glacial points.
    833. 

  833.        swe    = swe*1000.
    834. 

  834.        plai   = 0.01  ! Should make this zero?
    835. 

  835.        smc    = 1.00
    836. 

  836.        smh2o  = 1.00  ! Something else?
    837. 

  837.        runsb  = 0.00
    838. 

  838.        fgev   = ETA
    839. 

  839.        fcev   = 0.
    840. 

  840.        fctr   = 0.
    841. 

  841.        soilm  = 1.0   ! Something else?
    842. 

  842. !       SMAV    = 1.00  ! Something else?
    843. 

  843.        SNOWC(I,J) = 1.0
    844. 

  844.        
    845. 

  845.        QFX(I,J) = eta_kinematic
    846. 

  846.        POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
    847. 

  847. 

  848.        CHS2(I,J) = CQS2(I,J)
    849. 

  849.        IF ( Q1 .GT. QSFC(I,J) ) THEN
    850. 

  850.           CQS2(I,J) = CHS(I,J)
    851. 

  851.        ENDIF
    852. 

  852. 

  853.     ELSE
    854. 

  854. !jref:end
    855. 

  855. 

  856.              nee = -1.E36
    857. 

  857.              npp = -1.E36
    858. 

  858. 

  859. #if 0
    860. 

  860.              if ( I == 15 .and. J == 5 ) then
    861. 

  861.                 ! Intent (IN) or Intent (INOUT), but not Intent (OUT)
    862. 

  862.                 write(*,'("Before call to NOAHMP_SFLX, at point ", I8, I8)') i, j
    863. 

  863.                 write(*,'(10x, "ICE        = ",  I10   )') ICE
    864. 

  864.                 write(*,'(10x, "IST        = ",  I10   )') IST
    865. 

  865.                 write(*,'(10x, "VEGTYP     = ",  I10   )') VEGTYP
    866. 

  866.                 write(*,'(10x, "ISC        = ",  I10   )') ISC
    867. 

  867.                 write(*,'(10x, "NSOIL      = ",  I10   )') NSOIL
    868. 

  868.                 write(*,'(10x, "ZSOIL      = ", 7F20.10)') ZSOIL
    869. 

  869.                 write(*,'(10x, "DT         = ",  F20.10)') DT
    870. 

  870.                 write(*,'(10x, "QV1D       = ",  F20.10)') QV1D
    871. 

  871.                 write(*,'(10x, "SFCTMP     = ",  F20.10)') SFCTMP
    872. 

  872.                 write(*,'(10x, "UU         = ",  F20.10)') UU
    873. 

  873.                 write(*,'(10x, "VV         = ",  F20.10)') VV
    874. 

  874.                 write(*,'(10x, "SOLDN      = ",  F20.10)') SOLDN
    875. 

  875.                 write(*,'(10x, "LWDN       = ",  F20.10)') LWDN
    876. 

  876.                 write(*,'(10x, "PRCP       = ",  F20.10)') PRCP
    877. 

  877.                 write(*,'(10x, "ZLVL       = ",  F20.10)') ZLVL
    878. 

  878.                 write(*,'(10x, "CO2AIR     = ",  F20.10)') CO2AIR
    879. 

  879.                 write(*,'(10x, "O2AIR      = ",  F20.10)') O2AIR
    880. 

  880.                 write(*,'(10x, "COSZ       = ",  F20.10)') COSZ
    881. 

  881.                 write(*,'(10x, "TBOT       = ",  F20.10)') TBOT
    882. 

  882.                 write(*,'(10x, "FOLN       = ",  F20.10)') FOLN
    883. 

  883.                 write(*,'(10x, "SFCPRS     = ",  F20.10)') SFCPRS
    884. 

  884.                 write(*,'(10x, "SHDFAC     = ",  F20.10)') SHDFAC
    885. 

  885.                 write(*,'(10x, "LAT        = ",  F20.10)') LAT
    886. 

  886.                 write(*,'(10x, "DZ8W1D     = ",  F20.10)') DZ8W1D
    887. 

  887.                 write(*,'(10x, "EAH        = ",  F20.10)') EAH
    888. 

  888.                 write(*,'(10x, "TAH        = ",  F20.10)') TAH
    889. 

  889.                 write(*,'(10x, "FWET       = ",  F20.10)') FWET
    890. 

  890.                 write(*,'(10x, "FICEOLD    = ", 7F20.10)') FICEOLD
    891. 

  891.                 write(*,'(10x, "QSNOW      = ",  F20.10)') QSNOW
    892. 

  892.                 write(*,'(10x, "SNEQVO     = ",  F20.10)') SNEQVO
    893. 

  893.                 write(*,'(10x, "ISNOW      = ",  F20.10)') ISNOW
    894. 

  894.                 write(*,'(10x, "ZSNSO      = ", 7F20.10)') ZSNSO
    895. 

  895.                 write(*,'(10x, "CANLIQ     = ",  F20.10)') CANLIQ
    896. 

  896.                 write(*,'(10x, "CANICE     = ",  F20.10)') CANICE
    897. 

  897.                 write(*,'(10x, "SNOWD      = ",  F20.10)') SNOWD
    898. 

  898.                 write(*,'(10x, "SWE        = ",  F20.10)') SWE
    899. 

  899.                 write(*,'(10x, "SNICE      = ", 7F20.10)') SNICE
    900. 

  900.                 write(*,'(10x, "SNLIQ      = ", 7F20.10)') SNLIQ
    901. 

  901.                 write(*,'(10x, "TV         = ",  F20.10)') TV
    902. 

  902.                 write(*,'(10x, "TG         = ",  F20.10)') TG
    903. 

  903.                 write(*,'(10x, "STC        = ", 7F20.10)') STC
    904. 

  904.                 write(*,'(10x, "SMH2O      = ", 7F20.10)') SMH2O
    905. 

  905.                 write(*,'(10x, "SMC        = ", 7F20.10)') SMC
    906. 

  906.                 write(*,'(10x, "ZWT        = ",  F20.10)') ZWT
    907. 

  907.                 write(*,'(10x, "WA         = ",  F20.10)') WA
    908. 

  908.                 write(*,'(10x, "WT         = ",  F20.10)') WT
    909. 

  909.                 write(*,'(10x, "WSLAKE     = ",  F20.10)') WSLAKE
    910. 

  910.                 write(*,'(10x, "LFMASS     = ",  F20.10)') LFMASS
    911. 

  911.                 write(*,'(10x, "RTMASS     = ",  F20.10)') RTMASS
    912. 

  912.                 write(*,'(10x, "STMASS     = ",  F20.10)') STMASS
    913. 

  913.                 write(*,'(10x, "WOOD       = ",  F20.10)') WOOD
    914. 

  914.                 write(*,'(10x, "STBLCP     = ",  F20.10)') STBLCP
    915. 

  915.                 write(*,'(10x, "FASTCP     = ",  F20.10)') FASTCP
    916. 

  916.                 write(*,'(10x, "PLAI       = ",  F20.10)') PLAI
    917. 

  917.                 write(*,'(10x, "PSAI       = ",  F20.10)') PSAI
    918. 

  918.                 write(*,'(10x, "ALBOLD     = ",  F20.10)') ALBOLD
    919. 

  919.                 write(*,'(10x, "CM         = ",  F20.10)') CM
    920. 

  920.                 write(*,'(10x, "CH         = ",  F20.10)') CH
    921. 

  921.                 write(*,'(10x, "DX         = ",  F20.10)') DX
    922. 

  922.                 write(*,'(10x, "ISURBAN    = ",  I10   )') ISURBAN
    923. 

  923.                 write(*,'(10x, "IZ0TLND    = ",  I10   )') IZ0TLND
    924. 

  924.                 write(*,'(10x, "QC         = ",  F20.10)') QC
    925. 

  925.                 write(*,'(10x, "PBLH       = ",  F20.10)') PBLH
    926. 

  926.                 write(*,'(10x, "QSFC1D     = ",  F20.10)') QSFC1D
    927. 

  927.                 write(*,'(10x, "PSFC       = ",  F20.10)') PSFC
    928. 

  928.              endif
    929. 

  929. #endif
    930. 

  930. 

  931.              CALL NOAHMP_SFLX (&
    932. 

  932.                   I       , J       , LAT     , YEARLEN , JULIAN  , COSZ    , & ! IN : Time/Space-related
    933. 

  933.                   DT      , DX      , DZ8W1D  , NSOIL   , ZSOIL   , 3       , & ! IN : Model configuration 
    934. 

  934.                   SHDFAC  , SHDMAX  , VEGTYP  , ISURBAN , ICE     , IST     , & ! IN : Vegetation/Soil characteristics
    935. 

  935.                   ISC     ,                                                   & ! IN : Vegetation/Soil characteristics
    936. 

  936.                   IZ0TLND ,                                                   & ! IN : User options
    937. 

  937.                   SFCTMP  , SFCPRS  , PSFC    , UU      , VV      , QV1D    , & ! IN : Forcing
    938. 

  938.                   QC      , SOLDN   , LWDN    , PRCP    , TBOT    , CO2AIR  , & ! IN : Forcing
    939. 

  939.                   O2AIR   , FOLN    , FICEOLD , PBLH   ,                      & ! IN : Forcing
    940. 

  940.                   ZLVL    , ALBOLD  , SNEQVO  ,                               & ! IN/OUT : 
    941. 

  941.                   STC     , SMH2O   , SMC     , TAH     , EAH     , FWET    , & ! IN/OUT : 
    942. 

  942.                   CANLIQ  , CANICE  , TV      , TG      , QSFC1D  , QSNOW   , & ! IN/OUT : 
    943. 

  943.                   ISNOW   , ZSNSO   , SNOWD   , SWE     , SNICE   , SNLIQ   , & ! IN/OUT : 
    944. 

  944.                   ZWT     , WA      , WT      , WSLAKE  , LFMASS  , RTMASS  , & ! IN/OUT : 
    945. 

  945.                   STMASS  , WOOD    , STBLCP  , FASTCP  , PLAI    , PSAI    , & ! IN/OUT : 
    946. 

  946.                   CM      , CH      , CHSTAR  ,                               & ! IN/OUT : 
    947. 

  947.                   FSA     , FSR     , FIRA    , FSH     , SSOIL   , FCEV    , & ! OUT : 
    948. 

  948.                   FGEV    , FCTR    , ECAN    , ETRAN   , ESOIL   , TRAD    , & ! OUT : 
    949. 

  949.                   TS      , TGB     , TGV     , T2MV    , T2MB    , TSTAR   , & ! OUT : 
    950. 

  950.                   Q1      , Q2MV    , Q2MB    , RUNSF   , RUNSB   , APAR    , & ! OUT : 
    951. 

  951.                   PSN     , SAV     , SAG     , FSNO    , NEE     , GPP     , & ! OUT : 
    952. 

  952.                   NPP     , FVEG    , SALB    , QMELT   , PONDING , PONDING1, & ! OUT : 
    953. 

  953.                   PONDING2, RSSUN   , RSSHA   , BGAP    , WGAP    , GAP     , & ! OUT : 
    954. 

  954.                   ERRWAT  , CHV     , CHB     , EMISSI)                         ! OUT :
    955. 

  955. 

  956. #if 0
    957. 

  957.             if ( I == 15 .and. J == 5 ) then
    958. 

  958.                 ! Intent (OUT) or Intent (INOUT), but not Intent (IN)
    959. 

  959.                 write(*,'("After  call to NOAHMP_SFLX, at point ", I8, I8)') i, j
    960. 

  960.                 write(*,'(10x, "ZLVL       = ", 7F20.10)') ZLVL
    961. 

  961.                 write(*,'(10x, "EAH        = ",  F20.10)') EAH
    962. 

  962.                 write(*,'(10x, "TAH        = ",  F20.10)') TAH
    963. 

  963.                 write(*,'(10x, "FWET       = ",  F20.10)') FWET
    964. 

  964.                 write(*,'(10x, "QSNOW      = ",  F20.10)') QSNOW
    965. 

  965.                 write(*,'(10x, "SNEQVO     = ",  F20.10)') SNEQVO
    966. 

  966.                 write(*,'(10x, "ISNOW      = ",  F20.10)') ISNOW
    967. 

  967.                 write(*,'(10x, "ZSNSO      = ", 7F20.10)') ZSNSO
    968. 

  968.                 write(*,'(10x, "CANLIQ     = ",  F20.10)') CANLIQ
    969. 

  969.                 write(*,'(10x, "CANICE     = ",  F20.10)') CANICE
    970. 

  970.                 write(*,'(10x, "SNOWD      = ",  F20.10)') SNOWD
    971. 

  971.                 write(*,'(10x, "SWE        = ",  F20.10)') SWE
    972. 

  972.                 write(*,'(10x, "SNICE      = ", 3F20.10)') SNICE
    973. 

  973.                 write(*,'(10x, "SNLIQ      = ", 3F20.10)') SNLIQ
    974. 

  974.                 write(*,'(10x, "TV         = ",  F20.10)') TV
    975. 

  975.                 write(*,'(10x, "TG         = ",  F20.10)') TG
    976. 

  976.                 write(*,'(10x, "STC        = ", 7F20.10)') STC
    977. 

  977.                 write(*,'(10x, "SMH2O      = ", 7F20.10)') SMH2O
    978. 

  978.                 write(*,'(10x, "SMC        = ", 7F20.10)') SMC
    979. 

  979.                 write(*,'(10x, "ZWT        = ",  F20.10)') ZWT
    980. 

  980.                 write(*,'(10x, "WA         = ",  F20.10)') WA
    981. 

  981.                 write(*,'(10x, "WT         = ",  F20.10)') WT
    982. 

  982.                 write(*,'(10x, "WSLAKE     = ",  F20.10)') WSLAKE
    983. 

  983.                 write(*,'(10x, "LFMASS     = ",  F20.10)') LFMASS
    984. 

  984.                 write(*,'(10x, "RTMASS     = ",  F20.10)') RTMASS
    985. 

  985.                 write(*,'(10x, "STMASS     = ",  F20.10)') STMASS
    986. 

  986.                 write(*,'(10x, "WOOD       = ",  F20.10)') WOOD
    987. 

  987.                 write(*,'(10x, "STBLCP     = ",  F20.10)') STBLCP
    988. 

  988.                 write(*,'(10x, "FASTCP     = ",  F20.10)') FASTCP
    989. 

  989.                 write(*,'(10x, "PLAI       = ",  F20.10)') PLAI
    990. 

  990.                 write(*,'(10x, "PSAI       = ",  F20.10)') PSAI
    991. 

  991.                 write(*,'(10x, "ALBOLD     = ",  F20.10)') ALBOLD
    992. 

  992.                 write(*,'(10x, "CM         = ",  F20.10)') CM
    993. 

  993.                 write(*,'(10x, "CH         = ",  F20.10)') CH
    994. 

  994.                 write(*,'(10x, "FSA        = ",  F20.10)') FSA
    995. 

  995.                 write(*,'(10x, "FSR        = ",  F20.10)') FSR
    996. 

  996.                 write(*,'(10x, "FIRA       = ",  F20.10)') FIRA
    997. 

  997.                 write(*,'(10x, "FSH        = ",  F20.10)') FSH
    998. 

  998.                 write(*,'(10x, "SSOIL      = ",  F20.10)') SSOIL
    999. 

  999.                 write(*,'(10x, "FCEV       = ",  F20.10)') FCEV
    1000. 

  1000.                 write(*,'(10x, "FGEV       = ",  F20.10)') FGEV
    1001. 

  1001.                 write(*,'(10x, "FCTR       = ",  F20.10)') FCTR
    1002. 

  1002.                 write(*,'(10x, "TRAD       = ",  F20.10)') TRAD
    1003. 

  1003.                 write(*,'(10x, "ECAN       = ",  F20.10)') ECAN
    1004. 

  1004.                 write(*,'(10x, "ETRAN      = ",  F20.10)') ETRAN
    1005. 

  1005.                 write(*,'(10x, "ESOIL      = ",  F20.10)') ESOIL
    1006. 

  1006.                 write(*,'(10x, "RUNSF      = ",  F20.10)') RUNSF
    1007. 

  1007.                 write(*,'(10x, "RUNSB      = ",  F20.10)') RUNSB
    1008. 

  1008.                 write(*,'(10x, "APAR       = ",  F20.10)') APAR
    1009. 

  1009.                 write(*,'(10x, "PSN        = ",  F20.10)') PSN
    1010. 

  1010.                 write(*,'(10x, "SAV        = ",  F20.10)') SAV
    1011. 

  1011.                 write(*,'(10x, "SAG        = ",  F20.10)') SAG
    1012. 

  1012.                 write(*,'(10x, "FSNO       = ",  F20.10)') FSNO
    1013. 

  1013.                 write(*,'(10x, "NEE        = ",  F20.10)') NEE
    1014. 

  1014.                 write(*,'(10x, "GPP        = ",  F20.10)') GPP
    1015. 

  1015.                 write(*,'(10x, "NPP        = ",  F20.10)') NPP
    1016. 

  1016.                 write(*,'(10x, "TS         = ",  F20.10)') TS
    1017. 

  1017.                 write(*,'(10x, "FVEG       = ",  F20.10)') FVEG
    1018. 

  1018.                 write(*,'(10x, "SALB       = ",  F20.10)') SALB
    1019. 

  1019.                 write(*,'(10x, "ERRWAT     = ",  F20.10)') ERRWAT
    1020. 

  1020.                 write(*,'(10x, "QMELT      = ",  F20.10)') QMELT
    1021. 

  1021.                 write(*,'(10x, "PONDING    = ",  F20.10)') PONDING
    1022. 

  1022.                 write(*,'(10x, "PONDING1   = ",  F20.10)') PONDING1
    1023. 

  1023.                 write(*,'(10x, "PONDING2   = ",  F20.10)') PONDING2
    1024. 

  1024.                 write(*,'(10x, "QSFC1D     = ",  F20.10)') QSFC1D
    1025. 

  1025.                 write(*,'(10x, "CHSTAR     = ",  F20.10)') CHSTAR
    1026. 

  1026.                 write(*,'(10x, "TSTAR      = ",  F20.10)') TSTAR
    1027. 

  1027.                 write(*,'(10x, "T2MV       = ",  F20.10)') T2MV
    1028. 

  1028.                 write(*,'(10x, "T2MB       = ",  F20.10)') T2MB
    1029. 

  1029.                 write(*,'(10x, "RSSUN      = ",  F20.10)') RSSUN
    1030. 

  1030.                 write(*,'(10x, "RSSHA      = ",  F20.10)') RSSHA
    1031. 

  1031.                 write(*,'(10x, "BGAP       = ",  F20.10)') BGAP
    1032. 

  1032.                 write(*,'(10x, "WGAP       = ",  F20.10)') WGAP
    1033. 

  1033.                 write(*,'(10x, "GAP        = ",  F20.10)') GAP
    1034. 

  1034.                 write(*,'(10x, "TGV        = ",  F20.10)') TGV
    1035. 

  1035.                 write(*,'(10x, "TGB        = ",  F20.10)') TGB
    1036. 

  1036.                 write(*,'(10x, "Q1         = ",  F20.10)') Q1
    1037. 

  1037.              endif
    1038. 

  1038. #endif
    1039. 

  1039. 

  1040.              !Q1           = eah * 0.622 / (SFCPRS - 0.378*eah)
    1041. 

  1041. 

  1042.              chs2     (i,j)               = chstar2
    1043. 

  1043.              cqs2     (i,j)               = cqstar2
    1044. 

  1044.              QFX      (I,J)               = ecan + esoil + etran
    1045. 

  1045.              SNOWC    (I,J)               = fsno
    1046. 

  1046. 

  1047.    ENDIF ! glacial, seaice split ends 
    1048. 

  1048. !jref:end 
    1049. 

  1049. 

  1050.              isnowxy  (i,j)               = isnow
    1051. 

  1051.              canliqxy (i,j)               = canliq
    1052. 

  1052.              canicexy (i,j)               = canice
    1053. 

  1053.              snowh    (i,j)               = snowd
    1054. 

  1054.              snow     (i,j)               = swe
    1055. 

  1055.              zsnsoxy  (i,isnow+1:nsoil,j) = zsnso (isnow+1:nsoil)
    1056. 

  1056.              tslb     (i,      1:nsoil,j) = stc   (      1:nsoil)
    1057. 

  1057.              tsnoxy   (i,isnow+1:    0,j) = stc   (isnow+1:    0)
    1058. 

  1058.              smois    (i,      1:nsoil,j) = smc   (      1:nsoil)
    1059. 

  1059.              sh2o     (i,      1:nsoil,j) = smh2o (      1:nsoil)
    1060. 

  1060.              snicexy  (i,isnow+1:    0,j) = snice (isnow+1:    0)
    1061. 

  1061.              snliqxy  (i,isnow+1:    0,j) = snliq (isnow+1:    0)
    1062. 

  1062.              tvxy     (i,j)               = tv
    1063. 

  1063.              tgxy     (i,j)               = tg
    1064. 

  1064.              zwtxy    (i,j)               = zwt
    1065. 

  1065.              waxy     (i,j)               = wa
    1066. 

  1066.              wtxy     (i,j)               = wt
    1067. 

  1067.              lfmassxy (i,j)               = lfmass
    1068. 

  1068.              rtmassxy (i,j)               = rtmass
    1069. 

  1069.              stmassxy (i,j)               = stmass
    1070. 

  1070.              woodxy   (i,j)               = wood
    1071. 

  1071.              stblcpxy (i,j)               = stblcp
    1072. 

  1072.              fastcpxy (i,j)               = fastcp
    1073. 

  1073.              xlaixy   (i,j)               = plai
    1074. 

  1074.              xsaixy   (i,j)               = psai
    1075. 

  1075.              emiss    (i,j)               = emissi
    1076. 

  1076.              eahxy    (i,j)               = eah
    1077. 

  1077.              tahxy    (i,j)               = tah
    1078. 

  1078.              fwetxy   (i,j)               = fwet
    1079. 

  1079.              sneqvoxy (i,j)               = sneqvo
    1080. 

  1080.              alboldxy (i,j)               = albold
    1081. 

  1081.              qsnowxy  (i,j)               = qsnow
    1082. 

  1082.              wslakexy (i,j)               = wslake
    1083. 

  1083.              cmxy     (i,j)               = cm
    1084. 

  1084. !jref:start             
    1085. 

  1085.              chxy     (i,j)               = chstar
    1086. 

  1086.              rssunxy  (i,j)               = rssun
    1087. 

  1087.              rsshaxy  (i,j)               = rssha
    1088. 

  1088.              bgapxy   (i,j)               = bgap
    1089. 

  1089.              wgapxy   (i,j)               = wgap
    1090. 

  1090.              gapxy    (i,j)               = gap
    1091. 

  1091.              tgvxy    (i,j)               = tgv
    1092. 

  1092.              tgbxy    (i,j)               = tgb
    1093. 

  1093.              chvxy    (i,j)               = chv
    1094. 

  1094.              chbxy    (i,j)               = chb
    1095. 

  1095. !jref:end             
    1096. 

  1096. 

  1097.              !for output
    1098. 

  1098. 

  1099.              runsfxy  (i,j)               = runsf
    1100. 

  1100.              runsbxy  (i,j)               = runsb
    1101. 

  1101.              ecanxy   (i,j)               = ecan
    1102. 

  1102.              edirxy   (i,j)               = esoil
    1103. 

  1103.              etranxy  (i,j)               = etran
    1104. 

  1104.              aparxy   (i,j)               = apar
    1105. 

  1105.              psnxy    (i,j)               = psn
    1106. 

  1106.              savxy    (i,j)               = sav
    1107. 

  1107.              sagxy    (i,j)               = sag
    1108. 

  1108.              fsnoxy   (i,j)               = fsno
    1109. 

  1109.              fsaxy    (i,j)               = fsa
    1110. 

  1110.              firaxy   (i,j)               = fira
    1111. 

  1111.              hfx      (i,j)               = fsh
    1112. 

  1112.              lh       (i,j)               = fcev + fgev + fctr
    1113. 

  1113.              grdflx   (i,j)               = ssoil
    1114. 

  1114.              tradxy   (i,j)               = trad
    1115. 

  1115.              tsxy     (i,j)               = ts
    1116. 

  1116.              neexy    (i,j)               = nee
    1117. 

  1117.              gppxy    (i,j)               = gpp
    1118. 

  1118.              nppxy    (i,j)               = npp
    1119. 

  1119.              fvegxy   (i,j)               = fveg
    1120. 

  1120. !jref:4/21/2011
    1121. 

  1121.              t2mvxy   (i,j)               = t2mv
    1122. 

  1122.              t2mbxy   (i,j)               = t2mb
    1123. 

  1123.              q2mvxy   (i,j)               = q2mv
    1124. 

  1124.              q2mbxy   (i,j)               = q2mb
    1125. 

  1125.              chstarxy (i,j)               = chstar
    1126. 

  1126.              chs      (i,j)               = chstar
    1127. 

  1127.              tstarxy  (i,j)               = tstar
    1128. 

  1128. !jref:4/21/2011 
    1129. 

  1129. 

  1130.              CANWAT(I,J)  = canliqxy (i,j) + canicexy (i,j)
    1131. 

  1131.              IF ( SALB > -999 ) THEN
    1132. 

  1132.                 ALBEDO(I,J)  = salb
    1133. 

  1133.              ENDIF
    1134. 

  1134.              TSK(I,J)     = tradxy   (i,j)
    1135. 

  1135. !KWM          TSK(I,J)     = tstarxy  (i,j)
    1136. 

  1136. !niu          POTEVP(I,J)  = ???
    1137. 

  1137. !jref         CHS2(I,J)    = chxy     (i,j)
    1138. 

  1138. 

  1139.              !IF (Q1.GT.QSFC(I,J)) THEN
    1140. 

  1140.              !   CQS2(I,J) = CHS(I,J)
    1141. 

  1141.              !END IF
    1142. 

  1142.              QSFC(I,J)    = Q1/(1.0-Q1)
    1143. 

  1143. !jref: specific humidity to mixing ratio
    1144. 

  1144.              q2mvxy(i,j)  = q2mvxy(i,j)/(1.0-q2mvxy(i,j))
    1145. 

  1145. !             IF (VEGTYP == ISURBAN) write(*,*) "IN SFCDRV: q2mb=",q2mb,"q2mbxy(i,j)=",q2mbxy(i,j)
    1146. 

  1146.              q2mbxy(i,j)  = q2mbxy(i,j)/(1.0-q2mbxy(i,j))
    1147. 

  1147. 

  1148. !***  DIAGNOSTICS                   
    1149. 

  1149. !jref:start - THESE SHOULD BE LOOKED AT!!!
    1150. 

  1150.              SNOTIME(I,J) = SNOTIME1
    1151. 

  1151.              SMSTAV(I,J)=SOILW
    1152. 

  1152.              SMSTOT(I,J)=SOILM*1000.
    1153. 

  1153. !         Convert the water unit into mm
    1154. 

  1154.              SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+runsfxy(i,j)*DT*1000.0
    1155. 

  1155.              UDRUNOFF(I,J)=UDRUNOFF(I,J)+runsbxy(i,j)*DT*1000.0
    1156. 

  1156. !jref        SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
    1157. 

  1157. !jref        UDRUNOFF(I,J)=UDRUNOFF(I,J)+(RUNOFF2+RUNOFF3)*DT*1000.0
    1158. 

  1158. !jref:end             
    1159. 

  1159. ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
    1160. 

  1160.              IF(FFROZP.GT.0.5)THEN
    1161. 

  1161.                 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
    1162. 

  1162.              ENDIF
    1163. 

  1163.              IF(SNOW(I,J).GT.0.)THEN
    1164. 

  1164. !KWM                ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
    1165. 

  1165.              ENDIF
    1166. 

  1166. 

  1167.           ENDIF                                                           ! endif of land-sea test
    1168. 

  1168. !jref:start  make sure exchange coeff and TSK include water points         
    1169. 

  1169. !          IF((XLAND(I,J)-1.5).GE.0.)THEN                                  ! begining of land/sea if block
    1170. 

  1170. !             chstar2xy(i,j) = CHS2(i,j)
    1171. 

  1171. !             chstarxy(i,j) = CHS(i,j)    
    1172. 

  1172. !             tstarxy(i,j)  = T1 !TSK(i,j) test with T1
    1173. 

  1173. !          ENDIF
    1174. 

  1174. !jref:end
    1175. 

  1175. 

  1176. 

  1177.        ENDDO
    1178. 

  1178. 

  1179.     ENDDO                                                                ! of J loop
    1180. 

  1180. !------------------------------------------------------
    1181. 

  1181.   END SUBROUTINE noahmplsm
    1182. 

  1182. !------------------------------------------------------
    1183. 

  1183. 

  1184.   SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP ,           &
    1185. 

  1185.        TSLB , SMOIS , SH2O , DZS , FNDSOILW , FNDSNOWH ,                      &
    1186. 

  1186.        TSK, isnowxy , tvxy     ,tgxy     ,canicexy ,                          &
    1187. 

  1187.        canliqxy ,eahxy    ,tahxy    ,cmxy     ,chxy     ,                     &
    1188. 

  1188.        fwetxy   ,sneqvoxy ,alboldxy ,qsnowxy  ,wslakexy ,zwtxy    ,waxy     , &
    1189. 

  1189.        wtxy     ,tsnoxy   ,zsnsoxy  ,snicexy  ,snliqxy  ,lfmassxy ,rtmassxy , &
    1190. 

  1190.        stmassxy ,woodxy   ,stblcpxy ,fastcpxy ,xsaixy   , &
    1191. 

  1191. !jref:start
    1192. 

  1192.        t2mvxy   ,t2mbxy   ,chstarxy ,            &
    1193. 

  1193. !jref:end       
    1194. 

  1194.        num_soil_layers, restart,                 &
    1195. 

  1195.        allowed_to_read ,                         &
    1196. 

  1196.        ids,ide, jds,jde, kds,kde,                &
    1197. 

  1197.        ims,ime, jms,jme, kms,kme,                &
    1198. 

  1198.        its,ite, jts,jte, kts,kte                 )
    1199. 

  1199. 

  1200. ! Initializing Canopy air temperature to 287 K seems dangerous to me [KWM].
    1201. 

  1201. 

  1202.     INTEGER, INTENT(IN   )    ::     ids,ide, jds,jde, kds,kde,  &
    1203. 

  1203.          &                           ims,ime, jms,jme, kms,kme,  &
    1204. 

  1204.          &                           its,ite, jts,jte, kts,kte
    1205. 

  1205. 

  1206.     INTEGER, INTENT(IN)       ::     num_soil_layers
    1207. 

  1207. 

  1208.     LOGICAL, INTENT(IN)       ::     restart,                    &
    1209. 

  1209.          &                           allowed_to_read
    1210. 

  1210. 

  1211.     REAL,    DIMENSION( num_soil_layers), INTENT(IN)    ::     DZS  ! Thickness of the soil layers [m]
    1212. 

  1212. 

  1213.     REAL,    DIMENSION( ims:ime, num_soil_layers, jms:jme ) ,    &
    1214. 

  1214.          &   INTENT(INOUT)    ::     SMOIS,                      &
    1215. 

  1215.          &                           SH2O,                       &
    1216. 

  1216.          &                           TSLB
    1217. 

  1217. 

  1218.     REAL,    DIMENSION( ims:ime, jms:jme ) ,                     &
    1219. 

  1219.          &   INTENT(INOUT)    ::     SNOW,                       &
    1220. 

  1220.          &                           SNOWH,                      &
    1221. 

  1221.          &                           CANWAT
    1222. 

  1222. 

  1223.     INTEGER, DIMENSION( ims:ime, jms:jme ),                      &
    1224. 

  1224.          &   INTENT(IN)       ::     ISLTYP
    1225. 

  1225. 

  1226.     LOGICAL, INTENT(IN)       ::     FNDSOILW,                   &
    1227. 

  1227.          &                           FNDSNOWH
    1228. 

  1228. 

  1229.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: TSK         !skin temperature (k)
    1230. 

  1230.     INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy     !actual no. of snow layers
    1231. 

  1231.     REAL, DIMENSION(ims:ime,-2:num_soil_layers,jms:jme), INTENT(INOUT) :: zsnsoxy  !snow layer depth [m]
    1232. 

  1232.     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: tsnoxy   !snow temperature [K]
    1233. 

  1233.     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snicexy  !snow layer ice [mm]
    1234. 

  1234.     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snliqxy  !snow layer liquid water [mm]
    1235. 

  1235.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy        !vegetation canopy temperature
    1236. 

  1236.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy        !ground surface temperature
    1237. 

  1237.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy    !canopy-intercepted ice (mm)
    1238. 

  1238.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy    !canopy-intercepted liquid water (mm)
    1239. 

  1239.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy       !canopy air vapor pressure (pa)
    1240. 

  1240.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy       !canopy air temperature (k)
    1241. 

  1241.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy        !momentum drag coefficient
    1242. 

  1242.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy        !sensible heat exchange coefficient
    1243. 

  1243.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy      !wetted or snowed fraction of the canopy (-)
    1244. 

  1244.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy    !snow mass at last time step(mm h2o)
    1245. 

  1245.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy    !snow albedo at last time step (-)
    1246. 

  1246.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy     !snowfall on the ground [mm/s]
    1247. 

  1247.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy    !lake water storage [mm]
    1248. 

  1248.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy       !water table depth [m]
    1249. 

  1249.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy        !water in the "aquifer" [mm]
    1250. 

  1250.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy        !groundwater storage [mm]
    1251. 

  1251.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy    !leaf mass [g/m2]
    1252. 

  1252.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy    !mass of fine roots [g/m2]
    1253. 

  1253.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy    !stem mass [g/m2]
    1254. 

  1254.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy      !mass of wood (incl. woody roots) [g/m2]
    1255. 

  1255.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy    !stable carbon in deep soil [g/m2]
    1256. 

  1256.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy    !short-lived carbon, shallow soil [g/m2]
    1257. 

  1257.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy      !stem area index
    1258. 

  1258. 

  1259. !jref:start
    1260. 

  1260.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mvxy        !2m temperature vegetation part (k)
    1261. 

  1261.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mbxy        !2m temperature bare ground part (k)
    1262. 

  1262.     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chstarxy    !effective exchange coefficient
    1263. 

  1263. !jref:end
    1264. 

  1264. 

  1265. 

  1266.     REAL, DIMENSION(1:num_soil_layers)  :: ZSOIL      ! Depth of the soil layer bottom (m) from 
    1267. 

  1267.     !                                                   the surface (negative)
    1268. 

  1268. 

  1269.     REAL                      :: BX, SMCMAX, PSISAT, FREE
    1270. 

  1270. 

  1271.     REAL, PARAMETER           :: BLIM  = 5.5
    1272. 

  1272.     REAL, PARAMETER           :: HLICE = 3.335E5
    1273. 

  1273.     REAL, PARAMETER           :: GRAV = 9.81
    1274. 

  1274.     REAL, PARAMETER           :: T0 = 273.15
    1275. 

  1275. 

  1276.     INTEGER                   :: errflag
    1277. 

  1277. 

  1278.     character(len=80) :: err_message
    1279. 

  1279.     character(len=4)  :: MMINSL
    1280. 

  1280.     character(len=*), intent(in) :: MMINLU
    1281. 

  1281.     MMINSL='STAS'
    1282. 

  1282. 

  1283.     call read_mp_veg_parameters(trim(MMINLU))
    1284. 

  1284. 

  1285.     !
    1286. 

  1286.     ! initialize three Noah LSM related tables
    1287. 

  1287.     !
    1288. 

  1288.     IF ( allowed_to_read ) THEN
    1289. 

  1289.        CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
    1290. 

  1290.        CALL  SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
    1291. 

  1291.     ENDIF
    1292. 

  1292. 

  1293.     IF( .NOT. restart ) THEN
    1294. 

  1294. 

  1295.        itf=min0(ite,ide-1)
    1296. 

  1296.        jtf=min0(jte,jde-1)
    1297. 

  1297. 

  1298.        errflag = 0
    1299. 

  1299.        DO j = jts,jtf
    1300. 

  1300.           DO i = its,itf
    1301. 

  1301.              IF ( ISLTYP( i,j ) .LT. 1 ) THEN
    1302. 

  1302.                 errflag = 1
    1303. 

  1303.                 WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
    1304. 

  1304.                 CALL wrf_message(err_message)
    1305. 

  1305.              ENDIF
    1306. 

  1306.           ENDDO
    1307. 

  1307.        ENDDO
    1308. 

  1308.        IF ( errflag .EQ. 1 ) THEN
    1309. 

  1309.           CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
    1310. 

  1310.                "of ISLTYP. Is this field in the input?" )
    1311. 

  1311.        ENDIF
    1312. 

  1312. #ifdef WRF_CHEM
    1313. 

  1313.        !
    1314. 

  1314.        ! need this parameter for dust parameterization in wrf/chem
    1315. 

  1315.        !
    1316. 

  1316.        do I=1,NSLTYPE
    1317. 

  1317.           porosity(i)=maxsmc(i)
    1318. 

  1318.        enddo
    1319. 

  1319. #endif
    1320. 

  1320. 

  1321. ! initialize soil liquid water content SH2O
    1322. 

  1322. 

  1323. !  IF(.NOT.FNDSOILW) THEN
    1324. 

  1324. 

  1325. ! If no SWC, do the following
    1326. 

  1326. !         PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
    1327. 

  1327.        DO J = jts , jtf
    1328. 

  1328.           DO I = its , itf
    1329. 

  1329.              BX = BB(ISLTYP(I,J))
    1330. 

  1330.              SMCMAX = MAXSMC(ISLTYP(I,J))
    1331. 

  1331.              PSISAT = SATPSI(ISLTYP(I,J))
    1332. 

  1332.              IF ( ( bx > 0.0 ) .AND. ( smcmax > 0.0 ) .AND. ( psisat > 0.0 ) ) THEN
    1333. 

  1333.                 DO NS=1, num_soil_layers
    1334. 

  1334. 

  1335.                    IF ( TSLB(I,NS,J) < 273.149 ) THEN
    1336. 

  1336.                       ! SH2O <= SMOIS for T < 273.149K (-0.001C)
    1337. 

  1337. 

  1338.                       ! First guess of SH2O following explicit solution for 
    1339. 

  1339.                       ! Flerchinger Eqn from Koren et al, JGR, 1999, Eqn 17 
    1340. 

  1340.                       ! (KCOUNT=0 in function FRH2O).
    1341. 

  1341.                       BX = BB(ISLTYP(I,J))
    1342. 

  1342.                       SMCMAX = MAXSMC(ISLTYP(I,J))
    1343. 

  1343.                       PSISAT = SATPSI(ISLTYP(I,J))
    1344. 

  1344.                       IF ( BX >  BLIM ) BX = BLIM
    1345. 

  1345.                       FK=(( (HLICE/(GRAV*(-PSISAT))) *                              &
    1346. 

  1346.                            ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
    1347. 

  1347.                       FK = MAX(FK, 0.02)
    1348. 

  1348.                       SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
    1349. 

  1349. 

  1350.                       ! Use iterative solution for liquid soil water content
    1351. 

  1351.                       ! using function FRH2O, with the initial guess for SH2O 
    1352. 

  1352.                       ! from the above explicit first guess.
    1353. 

  1353.                       CALL FRH2O ( FREE , TSLB(I,NS,J) , SMOIS(I,NS,J) , SH2O(I,NS,J) )
    1354. 

  1354.                       SH2O(I,NS,J) = FREE
    1355. 

  1355.                    ELSE
    1356. 

  1356.                       ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
    1357. 

  1357.                       SH2O(I,NS,J)=SMOIS(I,NS,J)
    1358. 

  1358.                    ENDIF
    1359. 

  1359.                 END DO
    1360. 

  1360.              ELSE
    1361. 

  1361.                 DO NS=1, num_soil_layers
    1362. 

  1362.                    SH2O(I,NS,J)=SMOIS(I,NS,J)
    1363. 

  1363.                 END DO
    1364. 

  1364.              ENDIF
    1365. 

  1365.           ENDDO
    1366. 

  1366.        ENDDO
    1367. 

  1367. !  ENDIF
    1368. 

  1368. 

  1369. 

  1370.        !
    1371. 

  1371.        ! initialize physical snow height SNOWH
    1372. 

  1372.        !
    1373. 

  1373.        IF(.NOT.FNDSNOWH)THEN
    1374. 

  1374.           ! If no SNOWH do the following
    1375. 

  1375.           CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
    1376. 

  1376.           DO J = jts,jtf
    1377. 

  1377.              DO I = its,itf
    1378. 

  1378.                 SNOWH(I,J)=SNOW(I,J)*0.005               ! SNOW in mm and SNOWH in m
    1379. 

  1379.              ENDDO
    1380. 

  1380.           ENDDO
    1381. 

  1381.        ENDIF
    1382. 

  1382. 

  1383.        DO J = jts,jtf
    1384. 

  1384.           DO I = its,itf
    1385. 

  1385.              tvxy       (I,J) = TSK(I,J)
    1386. 

  1386.              tgxy       (I,J) = TSK(I,J)
    1387. 

  1387.              CANWAT     (I,J) = 0.0
    1388. 

  1388.              canliqxy   (I,J) = CANWAT(I,J)
    1389. 

  1389.              canicexy   (I,J) = 0.
    1390. 

  1390.              eahxy      (I,J) = 2000. 
    1391. 

  1391.              tahxy      (I,J) = 287.
    1392. 

  1392. !jref:start
    1393. 

  1393.              t2mvxy     (I,J) = TSK(I,J)
    1394. 

  1394.              t2mbxy     (I,J) = TSK(I,J)
    1395. 

  1395.              chstarxy   (I,J) = 0.0     
    1396. 

  1396. !jref:end
    1397. 

  1397. 

  1398.              cmxy       (I,J) = 0.0
    1399. 

  1399.              chxy       (I,J) = 0.0
    1400. 

  1400.              fwetxy     (I,J) = 0.0
    1401. 

  1401.              sneqvoxy   (I,J) = 0.0
    1402. 

  1402.              alboldxy   (I,J) = 0.65
    1403. 

  1403.              qsnowxy    (I,J) = 0.0
    1404. 

  1404.              wslakexy   (I,J) = 0.0
    1405. 

  1405. 

  1406.              waxy       (I,J) = 4900.                                       !???
    1407. 

  1407.              wtxy       (I,J) = waxy(i,j)                                   !???
    1408. 

  1408.              zwtxy      (I,J) = (25. + 2.0) - waxy(i,j)/1000/0.2            !???
    1409. 

  1409. 

  1410.              lfmassxy   (I,J) = 50.         !
    1411. 

  1411.              stmassxy   (I,J) = 50.0        !
    1412. 

  1412.              rtmassxy   (I,J) = 500.0       !
    1413. 

  1413.              woodxy     (I,J) = 500.0       !
    1414. 

  1414.              stblcpxy   (I,J) = 1000.0      !
    1415. 

  1415.              fastcpxy   (I,J) = 1000.0      !
    1416. 

  1416.              xsaixy     (I,J) = 0.1         !
    1417. 

  1417. 

  1418.           enddo
    1419. 

  1419.        enddo
    1420. 

  1420. 

  1421.        ! Given the soil layer thicknesses (in DZS), initialize the soil layer
    1422. 

  1422.        ! depths from the surface.
    1423. 

  1423.        ZSOIL(1)         = -DZS(1)          ! negative
    1424. 

  1424.        DO NS=2, num_soil_layers
    1425. 

  1425.           ZSOIL(NS)       = ZSOIL(NS-1) - DZS(NS)
    1426. 

  1426.        END DO
    1427. 

  1427. 

  1428.        ! Initialize snow/soil layer arrays ZSNSOXY, TSNOXY, SNICEXY, SNLIQXY, 
    1429. 

  1429.        ! and ISNOWXY
    1430. 

  1430.        CALL snow_init ( ims , ime , jms , jme , its , itf , jts , jtf , 3 , &
    1431. 

  1431.             &           num_soil_layers , zsoil , snow , tgxy , snowh ,     &
    1432. 

  1432.             &           zsnsoxy , tsnoxy , snicexy , snliqxy , isnowxy )
    1433. 

  1433. 

  1434.     ENDIF
    1435. 

  1435.   END SUBROUTINE NOAHMP_INIT
    1436. 

  1436. 

  1437. !------------------------------------------------------------------------------------------
    1438. 

  1438. !------------------------------------------------------------------------------------------
    1439. 

  1439. 

  1440.   SUBROUTINE SNOW_INIT ( ims , ime , jms , jme , its , itf , jts , jtf ,                  &
    1441. 

  1441.        &                 NSNOW , NSOIL , ZSOIL , SWE , TGXY , SNODEP ,                    &
    1442. 

  1442.        &                 ZSNSOXY , TSNOXY , SNICEXY ,SNLIQXY , ISNOWXY )
    1443. 

  1443. !------------------------------------------------------------------------------------------
    1444. 

  1444. !   Initialize snow arrays for Noah-MP LSM, based in input SNOWDEP, NSNOW
    1445. 

  1445. !   ISNOWXY is an index array, indicating the index of the top snow layer.  Valid indices
    1446. 

  1446. !           for snow layers range from 0 (no snow) and -1 (shallow snow) to (-NSNOW)+1 (deep snow).
    1447. 

  1447. !   TSNOXY holds the temperature of the snow layer.  Snow layers are initialized with 
    1448. 

  1448. !          temperature = ground temperature [?].  Snow-free levels in the array have value 0.0
    1449. 

  1449. !   SNICEXY is the frozen content of a snow layer.  Initial estimate based on SNODEP and SWE
    1450. 

  1450. !   SNLIQXY is the liquid content of a snow layer.  Initialized to 0.0
    1451. 

  1451. !   ZNSNOXY is the layer depth from the surface.  
    1452. 

  1452. !------------------------------------------------------------------------------------------
    1453. 

  1453.     IMPLICIT NONE
    1454. 

  1454. !------------------------------------------------------------------------------------------
    1455. 

  1455.     INTEGER, INTENT(IN)                              :: ims, ime, jms, jme
    1456. 

  1456.     INTEGER, INTENT(IN)                              :: its, itf, jts, jtf
    1457. 

  1457.     INTEGER, INTENT(IN)                              :: NSNOW
    1458. 

  1458.     INTEGER, INTENT(IN)                              :: NSOIL
    1459. 

  1459.     REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SWE 
    1460. 

  1460.     REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SNODEP
    1461. 

  1461.     REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: TGXY
    1462. 

  1462.     REAL,    INTENT(IN), DIMENSION(1:NSOIL)          :: ZSOIL
    1463. 

  1463. 

  1464.     INTEGER, INTENT(OUT), DIMENSION(ims:ime, jms:jme)                :: ISNOWXY ! Top snow layer index
    1465. 

  1465.     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:NSOIL,jms:jme) :: ZSNSOXY ! Snow/soil layer depth from surface [m]
    1466. 

  1466.     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: TSNOXY  ! Snow layer temperature [K]
    1467. 

  1467.     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: SNICEXY ! Snow layer ice content [mm]
    1468. 

  1468.     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: SNLIQXY ! snow layer liquid content [mm]
    1469. 

  1469. 

  1470. ! Local variables:
    1471. 

  1471. !   DZSNO   holds the thicknesses of the various snow layers.
    1472. 

  1472. !   DZSNOSO holds the thicknesses of the various soil/snow layers.
    1473. 

  1473.     INTEGER                           :: I,J,IZ
    1474. 

  1474.     REAL,   DIMENSION(-NSNOW+1:    0) :: DZSNO
    1475. 

  1475.     REAL,   DIMENSION(-NSNOW+1:NSOIL) :: DZSNSO
    1476. 

  1476. 

  1477. !------------------------------------------------------------------------------------------
    1478. 

  1478. 

  1479.     DO J = jts , jtf
    1480. 

  1480.        DO I = its , itf
    1481. 

  1481.           IF ( SNODEP(I,J) < 0.025 ) THEN
    1482. 

  1482.              ISNOWXY(I,J) = 0
    1483. 

  1483.              DZSNO(-NSNOW+1:0) = 0.
    1484. 

  1484.           ELSE
    1485. 

  1485.              IF ( ( SNODEP(I,J) >= 0.025 ) .AND. ( SNODEP(I,J) <= 0.05 ) ) THEN
    1486. 

  1486.                 ISNOWXY(I,J)    = -1
    1487. 

  1487.                 DZSNO(0)  = SNODEP(I,J)
    1488. 

  1488.              ELSE IF ( ( SNODEP(I,J) > 0.05 ) .AND. ( SNODEP(I,J) <= 0.10 ) ) THEN
    1489. 

  1489.                 ISNOWXY(I,J)    = -2
    1490. 

  1490.                 DZSNO(-1) = SNODEP(I,J)/2.
    1491. 

  1491.                 DZSNO( 0) = SNODEP(I,J)/2.
    1492. 

  1492.              ELSE IF ( (SNODEP(I,J) > 0.10 ) .AND. ( SNODEP(I,J) <= 0.25 ) ) THEN
    1493. 

  1493.                 ISNOWXY(I,J)    = -2
    1494. 

  1494.                 DZSNO(-1) = 0.05
    1495. 

  1495.                 DZSNO( 0) = SNODEP(I,J) - DZSNO(-1)
    1496. 

  1496.              ELSE IF ( ( SNODEP(I,J) > 0.25 ) .AND. ( SNODEP(I,J) <= 0.35 ) ) THEN
    1497. 

  1497.                 ISNOWXY(I,J)    = -3
    1498. 

  1498.                 DZSNO(-2) = 0.05
    1499. 

  1499.                 DZSNO(-1) = 0.5*(SNODEP(I,J)-DZSNO(-2))
    1500. 

  1500.                 DZSNO( 0) = 0.5*(SNODEP(I,J)-DZSNO(-2))
    1501. 

  1501.              ELSE IF ( SNODEP(I,J) > 0.35 ) THEN
    1502. 

  1502.                 ISNOWXY(I,J)     = -3
    1503. 

  1503.                 DZSNO(-2) = 0.05
    1504. 

  1504.                 DZSNO(-1) = 0.10
    1505. 

  1505.                 DZSNO( 0) = SNODEP(I,J) - DZSNO(-1) - DZSNO(-2)
    1506. 

  1506.              ELSE
    1507. 

  1507.                 CALL wrf_error_fatal("Problem with the logic assigning snow layers.")
    1508. 

  1508.              END IF
    1509. 

  1509.           END IF
    1510. 

  1510. 

  1511.           TSNOXY (I,-NSNOW+1:0,J) = 0.
    1512. 

  1512.           SNICEXY(I,-NSNOW+1:0,J) = 0.
    1513. 

  1513.           SNLIQXY(I,-NSNOW+1:0,J) = 0.
    1514. 

  1514.           DO IZ = ISNOWXY(I,J)+1 , 0
    1515. 

  1515.              TSNOXY(I,IZ,J)  = TGXY(I,J)  ! [k]
    1516. 

  1516.              SNLIQXY(I,IZ,J) = 0.00
    1517. 

  1517.              SNICEXY(I,IZ,J) = 1.00 * DZSNO(IZ) * (SWE(I,J)/SNODEP(I,J))  ! [kg/m3]
    1518. 

  1518.           END DO
    1519. 

  1519. 

  1520.           ! Assign local variable DZSNSO, the soil/snow layer thicknesses, for snow layers
    1521. 

  1521.           DO IZ = ISNOWXY(I,J)+1 , 0
    1522. 

  1522.              DZSNSO(IZ) = -DZSNO(IZ)
    1523. 

  1523.           END DO
    1524. 

  1524. 

  1525.           ! Assign local variable DZSNSO, the soil/snow layer thicknesses, for soil layers
    1526. 

  1526.           DZSNSO(1) = ZSOIL(1)
    1527. 

  1527.           DO IZ = 2 , NSOIL
    1528. 

  1528.              DZSNSO(IZ) = (ZSOIL(IZ) - ZSOIL(IZ-1))
    1529. 

  1529.           END DO
    1530. 

  1530. 

  1531.           ! Assign ZSNSOXY, the layer depths, for soil and snow layers
    1532. 

  1532.           ZSNSOXY(I,ISNOWXY(I,J)+1,J) = DZSNSO(ISNOWXY(I,J)+1)
    1533. 

  1533.           DO IZ = ISNOWXY(I,J)+2 , NSOIL
    1534. 

  1534.              ZSNSOXY(I,IZ,J) = ZSNSOXY(I,IZ-1,J) + DZSNSO(IZ)
    1535. 

  1535.           ENDDO
    1536. 

  1536. 

  1537.        END DO
    1538. 

  1538.     END DO
    1539. 

  1539. 

  1540.   END SUBROUTINE SNOW_INIT
    1541. 

  1541. !
    1542. 

  1542. !------------------------------------------------------------------------------------------
    1543. 

  1543. !------------------------------------------------------------------------------------------
    1544. 

  1544. !
    1545. 

  1545. END MODULE module_sf_noahmpdrv
    1546. 

  1546.