PageRenderTime 65ms CodeModel.GetById 20ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_sf_noahdrv.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1834 lines | 1144 code | 180 blank | 510 comment | 22 complexity | afe644bb40ac0e37a1f4953dbdf3283e MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. MODULE module_sf_noahdrv
    2. 

  2. 

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

  4.   USE module_sf_noahlsm
    5. 

  5.   USE module_sf_urban
    6. 

  6.   USE module_sf_bep
    7. 

  7.   USE module_sf_bep_bem
    8. 

  8. #ifdef WRF_CHEM
    9. 

  9.   USE module_data_gocart_dust
    10. 

  10. #endif
    11. 

  11. !-------------------------------
    12. 

  12. 

  13. !
    14. 

  14. CONTAINS
    15. 

  15. !
    16. 

  16. !----------------------------------------------------------------
    17. 

  17. ! Urban related variable are added to arguments - urban
    18. 

  18. !----------------------------------------------------------------
    19. 

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

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

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

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

  23.                   SNOWC,QSFC,RAINBL,MMINLU,                     &
    24. 

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

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

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

  27.                   myj,frpcpn,                                   &
    28. 

  28.                   SH2O,SNOWH,                                   & !H
    29. 

  29.                   U_PHY,V_PHY,                                  & !I
    30. 

  30.                   SNOALB,SHDMIN,SHDMAX,                         & !I
    31. 

  31.                   SNOTIME,                                      & !?
    32. 

  32.                   ACSNOM,ACSNOW,                                & !O
    33. 

  33.                   SNOPCX,                                       & !O
    34. 

  34.                   POTEVP,                                       & !O
    35. 

  35.                   SMCREL,                                       & !O
    36. 

  36.                   XICE_THRESHOLD,                               &
    37. 

  37.                   RDLAI2D,USEMONALB,                            &
    38. 

  38.                   RIB,                                          & !?
    39. 

  39.                   NOAHRES,                                      &
    40. 

  40.                   ids,ide, jds,jde, kds,kde,                    &
    41. 

  41.                   ims,ime, jms,jme, kms,kme,                    &
    42. 

  42.                   its,ite, jts,jte, kts,kte,                    &
    43. 

  43.                   sf_urban_physics,                             &
    44. 

  44.                   CMR_SFCDIF,CHR_SFCDIF,CMC_SFCDIF,CHC_SFCDIF,  &
    45. 

  45. !Optional Urban
    46. 

  46.                   TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & !H urban
    47. 

  47.                   UC_URB2D,                                     & !H urban
    48. 

  48.                   XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D,  & !H urban
    49. 

  49.                   TRL_URB3D,TBL_URB3D,TGL_URB3D,                & !H urban
    50. 

  50.                   SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,TS_URB2D,  & !H urban
    51. 

  51.                   PSIM_URB2D,PSIH_URB2D,U10_URB2D,V10_URB2D,    & !O urban
    52. 

  52.                   GZ1OZ0_URB2D,  AKMS_URB2D,                    & !O urban
    53. 

  53.                   TH2_URB2D,Q2_URB2D, UST_URB2D,                & !O urban
    54. 

  54.                   DECLIN_URB,COSZ_URB2D,OMG_URB2D,              & !I urban
    55. 

  55.                   XLAT_URB2D,                                   & !I urban
    56. 

  56.                   num_roof_layers, num_wall_layers,             & !I urban
    57. 

  57.                   num_road_layers, DZR, DZB, DZG,               & !I urban
    58. 

  58.                   FRC_URB2D,UTYPE_URB2D,                        & !O
    59. 

  59.                   num_urban_layers,                             & !I multi-layer urban
    60. 

  60.                   trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,      & !H multi-layer urban
    61. 

  61.                   tlev_urb3d,qlev_urb3d,                        & !H multi-layer urban
    62. 

  62.                   tw1lev_urb3d,tw2lev_urb3d,                    & !H multi-layer urban
    63. 

  63.                   tglev_urb3d,tflev_urb3d,                      & !H multi-layer urban
    64. 

  64.                   sf_ac_urb3d,lf_ac_urb3d,cm_ac_urb3d,          & !H multi-layer urban
    65. 

  65.                   sfvent_urb3d,lfvent_urb3d,                    & !H multi-layer urban
    66. 

  66.                   sfwin1_urb3d,sfwin2_urb3d,                    & !H multi-layer urban
    67. 

  67.                   sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,    & !H multi-layer urban
    68. 

  68.                   th_phy,rho,p_phy,ust,                         & !I multi-layer urban
    69. 

  69.                   gmt,julday,xlong,xlat,                        & !I multi-layer urban
    70. 

  70.                   a_u_bep,a_v_bep,a_t_bep,a_q_bep,              & !O multi-layer urban
    71. 

  71.                   a_e_bep,b_u_bep,b_v_bep,                      & !O multi-layer urban
    72. 

  72.                   b_t_bep,b_q_bep,b_e_bep,dlg_bep,              & !O multi-layer urban
    73. 

  73.                   dl_u_bep,sf_bep,vl_bep                         )   !O multi-layer urban         
    74. 

  74. 

  75. !----------------------------------------------------------------
    76. 

  76.     IMPLICIT NONE
    77. 

  77. !----------------------------------------------------------------
    78. 

  78. !----------------------------------------------------------------
    79. 

  79. ! --- atmospheric (WRF generic) variables
    80. 

  80. !-- DT          time step (seconds)
    81. 

  81. !-- DZ8W        thickness of layers (m)
    82. 

  82. !-- T3D         temperature (K)
    83. 

  83. !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
    84. 

  84. !-- P3D         3D pressure (Pa)
    85. 

  85. !-- FLHC        exchange coefficient for heat (m/s)
    86. 

  86. !-- FLQC        exchange coefficient for moisture (m/s)
    87. 

  87. !-- PSFC        surface pressure (Pa)
    88. 

  88. !-- XLAND       land mask (1 for land, 2 for water)
    89. 

  89. !-- QGH         saturated mixing ratio at 2 meter
    90. 

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

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

  92. !-- History variables
    93. 

  93. !-- CANWAT      canopy moisture content (mm)
    94. 

  94. !-- TSK         surface temperature (K)
    95. 

  95. !-- TSLB        soil temp (k)
    96. 

  96. !-- SMOIS       total soil moisture content (volumetric fraction)
    97. 

  97. !-- SH2O        unfrozen soil moisture content (volumetric fraction)
    98. 

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

  99. !-- SNOWH       actual snow depth (m)
    100. 

  100. !-- SNOW        liquid water-equivalent snow depth (m)
    101. 

  101. !-- ALBEDO      time-varying surface albedo including snow effect (unitless fraction)
    102. 

  102. !-- ALBBCK      background surface albedo (unitless fraction)
    103. 

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

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

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

  106. ! --- soil variables
    107. 

  107. !-- num_soil_layers   the number of soil layers
    108. 

  108. !-- ZS          depths of centers of soil layers   (m)
    109. 

  109. !-- DZS         thicknesses of soil layers (m)
    110. 

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

  111. !-- TMN         soil temperature at lower boundary (K)
    112. 

  112. !-- SMCWLT      wilting point (volumetric)
    113. 

  113. !-- SMCDRY      dry soil moisture threshold where direct evap from
    114. 

  114. !               top soil layer ends (volumetric)
    115. 

  115. !-- SMCREF      soil moisture threshold below which transpiration begins to
    116. 

  116. !                   stress (volumetric)
    117. 

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

  118. !-- NROOT       number of root layers, a function of veg type, determined
    119. 

  119. !               in subroutine redprm.
    120. 

  120. !-- SMSTAV      Soil moisture availability for evapotranspiration (
    121. 

  121. !                   fraction between SMCWLT and SMCMXA)
    122. 

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

  123. ! --- snow variables
    124. 

  124. !-- SNOWC       fraction snow coverage (0-1.0)
    125. 

  125. ! --- vegetation variables
    126. 

  126. !-- SNOALB      upper bound on maximum albedo over deep snow
    127. 

  127. !-- SHDMIN      minimum areal fractional coverage of annual green vegetation
    128. 

  128. !-- SHDMAX      maximum areal fractional coverage of annual green vegetation
    129. 

  129. !-- XLAI        leaf area index (dimensionless)
    130. 

  130. !-- Z0BRD       Background fixed roughness length (M)
    131. 

  131. !-- Z0          Background vroughness length (M) as function
    132. 

  132. !-- ZNT         Time varying roughness length (M) as function
    133. 

  133. !-- ALBD(IVGTPK,ISN) background albedo reading from a table
    134. 

  134. ! --- LSM output
    135. 

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

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

  137. !-- LH          upward moisture flux at the surface (W m-2)
    138. 

  138. !-- GRDFLX(I,J) ground heat flux (W m-2)
    139. 

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

  140. !----------------------------------------------------------------------------
    141. 

  141. !-- EC          canopy water evaporation ((W m-2)
    142. 

  142. !-- EDIR        direct soil evaporation (W m-2)
    143. 

  143. !-- ET          plant transpiration from a particular root layer (W m-2)
    144. 

  144. !-- ETT         total plant transpiration (W m-2)
    145. 

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

  146. !-- DRIP        through-fall of precip and/or dew in excess of canopy
    147. 

  147. !                 water-holding capacity (m)
    148. 

  148. !-- DEW         dewfall (or frostfall for t<273.15) (M)
    149. 

  149. !-- SMAV        Soil Moisture Availability for each layer, as a fraction
    150. 

  150. !                 between SMCWLT and SMCMAX (dimensionless fraction)
    151. 

  151. ! ----------------------------------------------------------------------
    152. 

  152. !-- BETA        ratio of actual/potential evap (dimensionless)
    153. 

  153. !-- ETP         potential evaporation (W m-2)
    154. 

  154. ! ----------------------------------------------------------------------
    155. 

  155. !-- FLX1        precip-snow sfc (W m-2)
    156. 

  156. !-- FLX2        freezing rain latent heat flux (W m-2)
    157. 

  157. !-- FLX3        phase-change heat flux from snowmelt (W m-2)
    158. 

  158. ! ----------------------------------------------------------------------
    159. 

  159. !-- ACSNOM      snow melt (mm) (water equivalent)
    160. 

  160. !-- ACSNOW      accumulated snow fall (mm) (water equivalent)
    161. 

  161. !-- SNOPCX      snow phase change heat flux (W/m^2)
    162. 

  162. !-- POTEVP      accumulated potential evaporation (W/m^2)
    163. 

  163. !-- RIB         Documentation needed!!!
    164. 

  164. ! ----------------------------------------------------------------------
    165. 

  165. !-- RUNOFF1     surface runoff (m s-1), not infiltrating the surface
    166. 

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

  167. !                  soil layer (baseflow)
    168. 

  168. !  important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
    169. 

  169. !-- RUNOFF3     numerical trunctation in excess of porosity (smcmax)
    170. 

  170. !                  for a given soil layer at the end of a time step (m s-1).
    171. 

  171. !SFCRUNOFF     Surface Runoff (mm)
    172. 

  172. !UDRUNOFF      Total Underground Runoff (mm), which is the sum of RUNOFF2 and RUNOFF3
    173. 

  173. ! ----------------------------------------------------------------------
    174. 

  174. !-- RC          canopy resistance (s m-1)
    175. 

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

  176. !-- RSMIN       minimum canopy resistance (s m-1)
    177. 

  177. !-- RCS         incoming solar rc factor (dimensionless)
    178. 

  178. !-- RCT         air temperature rc factor (dimensionless)
    179. 

  179. !-- RCQ         atmos vapor pressure deficit rc factor (dimensionless)
    180. 

  180. !-- RCSOIL      soil moisture rc factor (dimensionless)
    181. 

  181. 

  182. !-- EMISS       surface emissivity (between 0 and 1)
    183. 

  183. !-- EMBCK       Background surface emissivity (between 0 and 1)
    184. 

  184. 

  185. !-- ROVCP       R/CP
    186. 

  186. !               (R_d/R_v) (dimensionless)
    187. 

  187. !-- ids         start index for i in domain
    188. 

  188. !-- ide         end index for i in domain
    189. 

  189. !-- jds         start index for j in domain
    190. 

  190. !-- jde         end index for j in domain
    191. 

  191. !-- kds         start index for k in domain
    192. 

  192. !-- kde         end index for k in domain
    193. 

  193. !-- ims         start index for i in memory
    194. 

  194. !-- ime         end index for i in memory
    195. 

  195. !-- jms         start index for j in memory
    196. 

  196. !-- jme         end index for j in memory
    197. 

  197. !-- kms         start index for k in memory
    198. 

  198. !-- kme         end index for k in memory
    199. 

  199. !-- its         start index for i in tile
    200. 

  200. !-- ite         end index for i in tile
    201. 

  201. !-- jts         start index for j in tile
    202. 

  202. !-- jte         end index for j in tile
    203. 

  203. !-- kts         start index for k in tile
    204. 

  204. !-- kte         end index for k in tile
    205. 

  205. !
    206. 

  206. !-- SR          fraction of frozen precip (0.0 to 1.0)
    207. 

  207. !----------------------------------------------------------------
    208. 

  208. 

  209. ! IN only
    210. 

  210. 

  211.    INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &
    212. 

  212.                                     ims,ime, jms,jme, kms,kme,  &
    213. 

  213.                                     its,ite, jts,jte, kts,kte
    214. 

  214. 

  215.    INTEGER,  INTENT(IN   )   ::  sf_urban_physics               !urban
    216. 

  216.    INTEGER,  INTENT(IN   )   ::  isurban
    217. 

  217.    INTEGER,  INTENT(IN   )   ::  isice
    218. 

  218. 

  219.    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    220. 

  220.             INTENT(IN   )    ::                            TMN, &
    221. 

  221.                                                          XLAND, &
    222. 

  222.                                                           XICE, &
    223. 

  223.                                                         VEGFRA, &
    224. 

  224.                                                         SHDMIN, &
    225. 

  225.                                                         SHDMAX, &
    226. 

  226.                                                         SNOALB, &
    227. 

  227.                                                            GSW, &
    228. 

  228.                                                         SWDOWN, & !added 10 jan 2007
    229. 

  229.                                                            GLW, &
    230. 

  230.                                                         RAINBL, &
    231. 

  231.                                                         EMBCK,  &
    232. 

  232.                                                         SR
    233. 

  233. 

  234.    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    235. 

  235.             INTENT(INOUT)    ::                         ALBBCK, &
    236. 

  236.                                                             Z0
    237. 

  237.    CHARACTER(LEN=*), INTENT(IN   )    ::                 MMINLU
    238. 

  238. 

  239.    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme )            , &
    240. 

  240.             INTENT(IN   )    ::                           QV3D, &
    241. 

  241.                                                          p8w3D, &
    242. 

  242.                                                           DZ8W, &
    243. 

  243.                                                           T3D
    244. 

  244.    REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    245. 

  245.              INTENT(IN   )               ::               QGH,  &
    246. 

  246.                                                           CPM
    247. 

  247. 

  248.    INTEGER, DIMENSION( ims:ime, jms:jme )                     , &
    249. 

  249.             INTENT(IN   )    ::                         IVGTYP, &
    250. 

  250.                                                         ISLTYP
    251. 

  251. 

  252.    INTEGER, INTENT(IN)       ::     num_soil_layers,ITIMESTEP
    253. 

  253. 

  254.    REAL,     INTENT(IN   )   ::     DT,ROVCP
    255. 

  255. 

  256.    REAL,     DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
    257. 

  257. 

  258. ! IN and OUT
    259. 

  259. 

  260.    REAL,     DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
    261. 

  261.              INTENT(INOUT)   ::                          SMOIS, & ! total soil moisture
    262. 

  262.                                                          SH2O,  & ! new soil liquid
    263. 

  263.                                                          TSLB     ! TSLB     STEMP
    264. 

  264. 

  265.    REAL,     DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
    266. 

  266.              INTENT(OUT)     ::                         SMCREL
    267. 

  267. 

  268.    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    269. 

  269.             INTENT(INOUT)    ::                            TSK, & !was TGB (temperature)
    270. 

  270.                                                            HFX, &
    271. 

  271.                                                            QFX, &
    272. 

  272.                                                             LH, &
    273. 

  273.                                                         GRDFLX, &
    274. 

  274.                                                           QSFC,&
    275. 

  275.                                                           CQS2,&
    276. 

  276.                                                           CHS,   &
    277. 

  277.                                                           CHS2,&
    278. 

  278.                                                           SNOW, &
    279. 

  279.                                                          SNOWC, &
    280. 

  280.                                                          SNOWH, & !new
    281. 

  281.                                                         CANWAT, &
    282. 

  282.                                                         SMSTAV, &
    283. 

  283.                                                         SMSTOT, &
    284. 

  284.                                                      SFCRUNOFF, &
    285. 

  285.                                                       UDRUNOFF, &
    286. 

  286.                                                         ACSNOM, &
    287. 

  287.                                                         ACSNOW, &
    288. 

  288.                                                        SNOTIME, &
    289. 

  289.                                                         SNOPCX, &
    290. 

  290.                                                         EMISS,  &
    291. 

  291.                                                           RIB,  &
    292. 

  292.                                                         POTEVP, &
    293. 

  293.                                                         ALBEDO, &
    294. 

  294.                                                            ZNT
    295. 

  295.    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    296. 

  296.             INTENT(OUT)      ::                         NOAHRES
    297. 

  297. 

  298.    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    299. 

  299.                INTENT(OUT)    ::                        CHKLOWQ
    300. 

  300.    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LAI
    301. 

  301.    REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) ::        QZ0
    302. 

  302. 

  303.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMR_SFCDIF
    304. 

  304.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHR_SFCDIF
    305. 

  305.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMC_SFCDIF
    306. 

  306.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHC_SFCDIF
    307. 

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

  308. 

  309.       REAL, DIMENSION(1:num_soil_layers) ::  ET
    310. 

  310. 

  311.       REAL, DIMENSION(1:num_soil_layers) ::  SMAV
    312. 

  312. 

  313.       REAL  ::  BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT,        &
    314. 

  314.                 FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI,  &
    315. 

  315. !                RCS,RCT,RCQ,RCSOIL
    316. 

  316.                 RCS,RCT,RCQ,RCSOIL,FFROZP
    317. 

  317. 

  318.     LOGICAL,    INTENT(IN   )    ::     myj,frpcpn
    319. 

  319. 

  320. ! DECLARATIONS - LOGICAL
    321. 

  321. ! ----------------------------------------------------------------------
    322. 

  322.       LOGICAL, PARAMETER :: LOCAL=.false.
    323. 

  323.       LOGICAL :: FRZGRA, SNOWNG
    324. 

  324. 

  325.       LOGICAL :: IPRINT
    326. 

  326. 

  327. ! ----------------------------------------------------------------------
    328. 

  328. ! DECLARATIONS - INTEGER
    329. 

  329. ! ----------------------------------------------------------------------
    330. 

  330.       INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
    331. 

  331.       INTEGER :: NROOT
    332. 

  332.       INTEGER :: KZ ,K
    333. 

  333.       INTEGER :: NS
    334. 

  334. ! ----------------------------------------------------------------------
    335. 

  335. ! DECLARATIONS - REAL
    336. 

  336. ! ----------------------------------------------------------------------
    337. 

  337. 

  338.       REAL  :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN,                    &
    339. 

  339.                Q2SAT,Q2SATI,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1,         &
    340. 

  340.                SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
    341. 

  341.                EMBRD,                                                    &
    342. 

  342.                Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT,SOLNET,E2SAT,SFCTSNO,   &
    343. 

  343. ! mek, WRF testing, expanded diagnostics
    344. 

  344.                SOLUP,LWUP,RNET,RES,Q1SFC,TAIRV,SATFLG
    345. 

  345. ! MEK MAY 2007
    346. 

  346.       REAL ::  FDTLIW
    347. 

  347. ! MEK JUL2007 for pot. evap.
    348. 

  348.       REAL :: RIBB
    349. 

  349.       REAL ::  FDTW
    350. 

  350. 

  351.       REAL  :: EMISSI
    352. 

  352. 

  353.       REAL  :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
    354. 

  354. 

  355.       REAL  :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
    356. 

  356.       REAL  :: SNOTIME1    ! LSTSNW1 INITIAL NUMBER OF TIMESTEPS SINCE LAST SNOWFALL
    357. 

  357. 

  358.       REAL  :: DUMMY,Z0BRD
    359. 

  359. !
    360. 

  360.       REAL  :: COSZ, SOLARDIRECT
    361. 

  361. !
    362. 

  362.       REAL, DIMENSION(1:num_soil_layers)::  SLDPTH, STC,SMC,SWC
    363. 

  363. !
    364. 

  364.       REAL, DIMENSION(1:num_soil_layers) ::     ZSOIL, RTDIS
    365. 

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

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

  367. ! MEK MAY 2007
    368. 

  368.       REAL, PARAMETER  :: ROW=1.E3,ELIW=XLF,ROWLIW=ROW*ELIW
    369. 

  369. 

  370. ! ----------------------------------------------------------------------
    371. 

  371. ! DECLARATIONS START - urban
    372. 

  372. ! ----------------------------------------------------------------------
    373. 

  373. 

  374. ! input variables surface_driver --> lsm
    375. 

  375.      INTEGER, INTENT(IN) :: num_roof_layers
    376. 

  376.      INTEGER, INTENT(IN) :: num_wall_layers
    377. 

  377.      INTEGER, INTENT(IN) :: num_road_layers
    378. 

  378.      REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
    379. 

  379.      REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
    380. 

  380.      REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
    381. 

  381.      REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
    382. 

  382.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
    383. 

  383.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
    384. 

  384.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
    385. 

  385.      REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
    386. 

  386.      REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
    387. 

  387.      REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: TH_PHY
    388. 

  388.      REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: P_PHY
    389. 

  389.      REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: RHO
    390. 

  390.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
    391. 

  391. 

  392.      LOGICAL, intent(in) :: rdlai2d
    393. 

  393.      LOGICAL, intent(in) :: USEMONALB
    394. 

  394. 

  395. ! input variables lsm --> urban
    396. 

  396.      INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
    397. 

  397.      REAL :: TA_URB       ! potential temp at 1st atmospheric level [K]
    398. 

  398.      REAL :: QA_URB       ! mixing ratio at 1st atmospheric level  [kg/kg]
    399. 

  399.      REAL :: UA_URB       ! wind speed at 1st atmospheric level    [m/s]
    400. 

  400.      REAL :: U1_URB       ! u at 1st atmospheric level             [m/s]
    401. 

  401.      REAL :: V1_URB       ! v at 1st atmospheric level             [m/s]
    402. 

  402.      REAL :: SSG_URB      ! downward total short wave radiation    [W/m/m]
    403. 

  403.      REAL :: LLG_URB      ! downward long wave radiation           [W/m/m]
    404. 

  404.      REAL :: RAIN_URB     ! precipitation                          [mm/h]
    405. 

  405.      REAL :: RHOO_URB     ! air density                            [kg/m^3]
    406. 

  406.      REAL :: ZA_URB       ! first atmospheric level                [m]
    407. 

  407.      REAL :: DELT_URB     ! time step                              [s]
    408. 

  408.      REAL :: SSGD_URB     ! downward direct short wave radiation   [W/m/m]
    409. 

  409.      REAL :: SSGQ_URB     ! downward diffuse short wave radiation  [W/m/m]
    410. 

  410.      REAL :: XLAT_URB     ! latitude                               [deg]
    411. 

  411.      REAL :: COSZ_URB     ! cosz
    412. 

  412.      REAL :: OMG_URB      ! hour angle
    413. 

  413.      REAL :: ZNT_URB      ! roughness length                       [m]
    414. 

  414.      REAL :: TR_URB
    415. 

  415.      REAL :: TB_URB
    416. 

  416.      REAL :: TG_URB
    417. 

  417.      REAL :: TC_URB
    418. 

  418.      REAL :: QC_URB
    419. 

  419.      REAL :: UC_URB
    420. 

  420.      REAL :: XXXR_URB
    421. 

  421.      REAL :: XXXB_URB
    422. 

  422.      REAL :: XXXG_URB
    423. 

  423.      REAL :: XXXC_URB
    424. 

  424.      REAL, DIMENSION(1:num_roof_layers) :: TRL_URB  ! roof layer temp [K]
    425. 

  425.      REAL, DIMENSION(1:num_wall_layers) :: TBL_URB  ! wall layer temp [K]
    426. 

  426.      REAL, DIMENSION(1:num_road_layers) :: TGL_URB  ! road layer temp [K]
    427. 

  427.      LOGICAL  :: LSOLAR_URB
    428. 

  428. ! state variable surface_driver <--> lsm <--> urban
    429. 

  429.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
    430. 

  430.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
    431. 

  431.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
    432. 

  432.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
    433. 

  433.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
    434. 

  434.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
    435. 

  435.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
    436. 

  436.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
    437. 

  437.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
    438. 

  438.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
    439. 

  439.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
    440. 

  440.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
    441. 

  441.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
    442. 

  442.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
    443. 

  443. !
    444. 

  444.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
    445. 

  445. 

  446.      REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
    447. 

  447.      REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
    448. 

  448.      REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
    449. 

  449. 

  450. ! output variable lsm --> surface_driver
    451. 

  451.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
    452. 

  452.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
    453. 

  453.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
    454. 

  454.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
    455. 

  455.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
    456. 

  456.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
    457. 

  457.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
    458. 

  458. !
    459. 

  459.      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
    460. 

  460. !
    461. 

  461.      REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
    462. 

  462.      REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: FRC_URB2D
    463. 

  463.      INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
    464. 

  464. 

  465. 

  466. ! output variables urban --> lsm
    467. 

  467.      REAL :: TS_URB     ! surface radiative temperature    [K]
    468. 

  468.      REAL :: QS_URB     ! surface humidity                 [-]
    469. 

  469.      REAL :: SH_URB     ! sensible heat flux               [W/m/m]
    470. 

  470.      REAL :: LH_URB     ! latent heat flux                 [W/m/m]
    471. 

  471.      REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic  [kg/m/m/s]
    472. 

  472.      REAL :: SW_URB     ! upward short wave radiation flux [W/m/m]
    473. 

  473.      REAL :: ALB_URB    ! time-varying albedo            [fraction]
    474. 

  474.      REAL :: LW_URB     ! upward long wave radiation flux  [W/m/m]
    475. 

  475.      REAL :: G_URB      ! heat flux into the ground        [W/m/m]
    476. 

  476.      REAL :: RN_URB     ! net radiation                    [W/m/m]
    477. 

  477.      REAL :: PSIM_URB   ! shear f for momentum             [-]
    478. 

  478.      REAL :: PSIH_URB   ! shear f for heat                 [-]
    479. 

  479.      REAL :: GZ1OZ0_URB   ! shear f for heat                 [-]
    480. 

  480.      REAL :: U10_URB    ! wind u component at 10 m         [m/s]
    481. 

  481.      REAL :: V10_URB    ! wind v component at 10 m         [m/s]
    482. 

  482.      REAL :: TH2_URB    ! potential temperature at 2 m     [K]
    483. 

  483.      REAL :: Q2_URB     ! humidity at 2 m                  [-]
    484. 

  484.      REAL :: CHS_URB
    485. 

  485.      REAL :: CHS2_URB
    486. 

  486.      REAL :: UST_URB
    487. 

  487. ! Variables for multi-layer UCM (Martilli et al. 2002)
    488. 

  488.    REAL, OPTIONAL, INTENT(IN  )   ::                                   GMT 
    489. 

  489.    INTEGER, OPTIONAL, INTENT(IN  ) ::                               JULDAY
    490. 

  490.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )        ::XLAT, XLONG
    491. 

  491.    INTEGER, INTENT(IN  ) ::                               NUM_URBAN_LAYERS
    492. 

  492.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: trb_urb4d
    493. 

  493.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1_urb4d
    494. 

  494.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2_urb4d
    495. 

  495.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tgb_urb4d
    496. 

  496.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tlev_urb3d
    497. 

  497.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: qlev_urb3d
    498. 

  498.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
    499. 

  499.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
    500. 

  500.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tglev_urb3d
    501. 

  501.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tflev_urb3d
    502. 

  502.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
    503. 

  503.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
    504. 

  504.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
    505. 

  505.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
    506. 

  506.    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
    507. 

  507.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
    508. 

  508.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
    509. 

  509.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
    510. 

  510.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
    511. 

  511.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
    512. 

  512.    REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d
    513. 

  513.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep   !Implicit momemtum component X-direction
    514. 

  514.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep   !Implicit momemtum component Y-direction
    515. 

  515.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep   !Implicit component pot. temperature
    516. 

  516.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep   !Implicit momemtum component X-direction
    517. 

  517.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep   !Implicit component TKE
    518. 

  518.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep   !Explicit momentum component X-direction
    519. 

  519.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep   !Explicit momentum component Y-direction
    520. 

  520.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep   !Explicit component pot. temperature
    521. 

  521.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep   !Implicit momemtum component Y-direction
    522. 

  522.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep   !Explicit component TKE
    523. 

  523.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep    !Fraction air volume in grid cell
    524. 

  524.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep   !Height above ground
    525. 

  525.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep  !Fraction air at the face of grid cell
    526. 

  526.    REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep  !Length scale
    527. 

  527. 

  528. ! Local variables for multi-layer UCM (Martilli et al. 2002)
    529. 

  529.    REAL,    DIMENSION( ims:ime, jms:jme ) :: HFX_RURAL,LH_RURAL,GRDFLX_RURAL,RN_RURAL
    530. 

  530.    REAL,    DIMENSION( ims:ime, jms:jme ) :: QFX_RURAL,QSFC_RURAL,UMOM_RURAL,VMOM_RURAL
    531. 

  531.    REAL,    DIMENSION( ims:ime, jms:jme ) :: ALB_RURAL,EMISS_RURAL,UST_RURAL,TSK_RURAL
    532. 

  532. !   REAL,    DIMENSION( ims:ime, jms:jme ) :: GRDFLX_URB
    533. 

  533. !   REAL,    DIMENSION( ims:ime, jms:jme ) :: QFX_URB,QSFC_URB,UMOM_URB,VMOM_URB
    534. 

  534.    REAL,    DIMENSION( ims:ime, jms:jme ) :: HFX_URB,UMOM_URB,VMOM_URB
    535. 

  535.    REAL,    DIMENSION( ims:ime, jms:jme ) :: QFX_URB
    536. 

  536. !   REAL,    DIMENSION( ims:ime, jms:jme ) :: ALBEDO_URB,EMISS_URB,UMOM,VMOM,UST
    537. 

  537.    REAL, DIMENSION(ims:ime,jms:jme) ::EMISS_URB
    538. 

  538.    REAL, DIMENSION(ims:ime,jms:jme) :: RL_UP_URB
    539. 

  539.    REAL, DIMENSION(ims:ime,jms:jme) ::RS_ABS_URB
    540. 

  540.    REAL, DIMENSION(ims:ime,jms:jme) ::GRDFLX_URB
    541. 

  541.    REAL :: SIGMA_SB,RL_UP_RURAL,RL_UP_TOT,RS_ABS_TOT,UMOM,VMOM
    542. 

  542.    REAL :: r1,r2,r3
    543. 

  543.    REAL :: CMR_URB, CHR_URB, CMC_URB, CHC_URB
    544. 

  544. ! ----------------------------------------------------------------------
    545. 

  545. ! DECLARATIONS END - urban
    546. 

  546. ! ----------------------------------------------------------------------
    547. 

  547. 

  548.   REAL, PARAMETER  :: CAPA=R_D/CP
    549. 

  549.   REAL :: APELM,APES,SFCTH2,PSFC
    550. 

  550.   real, intent(in) :: xice_threshold
    551. 

  551.   character(len=80) :: message_text
    552. 

  552. 

  553. ! MEK MAY 2007
    554. 

  554.       FDTLIW=DT/ROWLIW
    555. 

  555. ! MEK JUL2007
    556. 

  556.       FDTW=DT/(XLV*RHOWATER)
    557. 

  557. ! debug printout
    558. 

  558.          IPRINT=.false.
    559. 

  559. 

  560. !      SLOPETYP=2
    561. 

  561.       SLOPETYP=1
    562. 

  562. !      SHDMIN=0.00
    563. 

  563. 

  564. 

  565.       NSOIL=num_soil_layers
    566. 

  566. 

  567.      DO NS=1,NSOIL
    568. 

  568.      SLDPTH(NS)=DZS(NS)
    569. 

  569.      ENDDO
    570. 

  570. 

  571.      JLOOP : DO J=jts,jte
    572. 

  572. 

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

  574.         DO 50 I=its,ite
    575. 

  575. !*** initialize soil conditions for IHOP 31 May case
    576. 

  576. !         IF((XLAND(I,J)-1.5) < 0.)THEN
    577. 

  577. !            if (I==108.and.j==85) then
    578. 

  578. !                  DO NS=1,NSOIL
    579. 

  579. !                      SMOIS(I,NS,J)=0.10
    580. 

  580. !                      SH2O(I,NS,J)=0.10
    581. 

  581. !                  enddo
    582. 

  582. !             endif
    583. 

  583. !         ENDIF
    584. 

  584. 

  585. !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
    586. 

  586.           IF((XLAND(I,J)-1.5).GE.0.)THEN
    587. 

  587. ! check sea-ice point
    588. 

  588. #if 0
    589. 

  589.             IF( XICE(I,J).GE. XICE_THRESHOLD .and. IPRINT ) PRINT*, ' sea-ice at water point, I=',I,'J=',J
    590. 

  590. #endif
    591. 

  591. !***   Open Water Case
    592. 

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

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

  594.             DO NS=1,NSOIL
    595. 

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

  596.               TSLB(I,NS,J)=273.16                                          !STEMP
    597. 

  597.               SMCREL(I,NS,J)=1.0
    598. 

  598.             ENDDO
    599. 

  599.           ELSE
    600. 

  600.             IF ( XICE(I,J) .GE. XICE_THRESHOLD ) THEN
    601. 

  601. !***        SEA-ICE CASE
    602. 

  602.               SMSTAV(I,J)=1.0
    603. 

  603.               SMSTOT(I,J)=1.0
    604. 

  604.               DO NS=1,NSOIL
    605. 

  605.                 SMOIS(I,NS,J)=1.0
    606. 

  606.                 SMCREL(I,NS,J)=1.0
    607. 

  607.               ENDDO
    608. 

  608.             ENDIF
    609. 

  609.           ENDIF
    610. 

  610. !
    611. 

  611.    50   CONTINUE
    612. 

  612.       ENDIF                                                               ! end of initialization over ocean
    613. 

  613. 

  614. !-----------------------------------------------------------------------
    615. 

  615.       ILOOP : DO I=its,ite
    616. 

  616. ! surface pressure
    617. 

  617.         PSFC=P8w3D(i,1,j)
    618. 

  618. ! pressure in middle of lowest layer
    619. 

  619.         SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
    620. 

  620. ! convert from mixing ratio to specific humidity
    621. 

  621.          Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
    622. 

  622. !
    623. 

  623. !         Q2SAT=QGH(I,j)
    624. 

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

  625. ! add check on myj=.true.
    626. 

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

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

  628.           SATFLG=0.
    629. 

  629.           CHKLOWQ(I,J)=0.
    630. 

  630.         ELSE
    631. 

  631.           SATFLG=1.0
    632. 

  632.           CHKLOWQ(I,J)=1.
    633. 

  633.         ENDIF
    634. 

  634. 

  635.         SFCTMP=T3D(i,1,j)
    636. 

  636.         ZLVL=0.5*DZ8W(i,1,j)
    637. 

  637. 

  638. !        TH2=SFCTMP+(0.0097545*ZLVL)
    639. 

  639. ! calculate SFCTH2 via Exner function vs lapse-rate (above)
    640. 

  640.          APES=(1.E5/PSFC)**CAPA
    641. 

  641.          APELM=(1.E5/SFCPRS)**CAPA
    642. 

  642.          SFCTH2=SFCTMP*APELM
    643. 

  643.          TH2=SFCTH2/APES
    644. 

  644. !
    645. 

  645.          EMISSI = EMISS(I,J)
    646. 

  646.          LWDN=GLW(I,J)*EMISSI
    647. 

  647. ! SOLDN is total incoming solar
    648. 

  648.         SOLDN=SWDOWN(I,J)
    649. 

  649. ! GSW is net downward solar
    650. 

  650. !        SOLNET=GSW(I,J)
    651. 

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

  652.         SOLNET=SOLDN*(1.-ALBEDO(I,J))
    653. 

  653.         PRCP=RAINBL(i,j)/DT
    654. 

  654.         VEGTYP=IVGTYP(I,J)
    655. 

  655.         SOILTYP=ISLTYP(I,J)
    656. 

  656.         SHDFAC=VEGFRA(I,J)/100.
    657. 

  657.         T1=TSK(I,J)
    658. 

  658.         CHK=CHS(I,J)
    659. 

  659.         SHMIN=SHDMIN(I,J)/100. !NEW
    660. 

  660.         SHMAX=SHDMAX(I,J)/100. !NEW
    661. 

  661. ! convert snow water equivalent from mm to meter
    662. 

  662.         SNEQV=SNOW(I,J)*0.001
    663. 

  663. ! snow depth in meters
    664. 

  664.         SNOWHK=SNOWH(I,J)
    665. 

  665.         SNCOVR=SNOWC(I,J)
    666. 

  666. 

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

  668. ! SR from e.g. Ferrier microphysics
    669. 

  669. ! otherwise define from 1st atmos level temperature
    670. 

  670.        IF(FRPCPN) THEN
    671. 

  671.           FFROZP=SR(I,J)
    672. 

  672.         ELSE
    673. 

  673.           IF (SFCTMP <=  273.15) THEN
    674. 

  674.             FFROZP = 1.0
    675. 

  675. 	  ELSE
    676. 

  676. 	    FFROZP = 0.0
    677. 

  677. 	  ENDIF
    678. 

  678.         ENDIF
    679. 

  679. !***
    680. 

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

  681. ! Open water points
    682. 

  682.           TSK_RURAL(I,J)=TSK(I,J)
    683. 

  683.           HFX_RURAL(I,J)=HFX(I,J)
    684. 

  684.           QFX_RURAL(I,J)=QFX(I,J)
    685. 

  685.           LH_RURAL(I,J)=LH(I,J)
    686. 

  686.           EMISS_RURAL(I,J)=EMISS(I,J)
    687. 

  687.           GRDFLX_RURAL(I,J)=GRDFLX(I,J)
    688. 

  688.         ELSE
    689. 

  689. ! Land or sea-ice case
    690. 

  690. 

  691.           IF (XICE(I,J) >= XICE_THRESHOLD) THEN
    692. 

  692.              ! Sea-ice point
    693. 

  693.              ICE = 1
    694. 

  694.           ELSE IF ( VEGTYP == ISICE ) THEN
    695. 

  695.              ! Land-ice point
    696. 

  696.              ICE = -1
    697. 

  697.           ELSE
    698. 

  698.              ! Neither sea ice or land ice.
    699. 

  699.              ICE=0
    700. 

  700.           ENDIF
    701. 

  701.           DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
    702. 

  702. 

  703.           IF(SNOW(I,J).GT.0.0)THEN
    704. 

  704. ! snow on surface (use ice saturation properties)
    705. 

  705.             SFCTSNO=SFCTMP
    706. 

  706.             E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
    707. 

  707.             Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
    708. 

  708.             Q2SATI=Q2SATI/(1.0+Q2SATI)    ! spec. hum.
    709. 

  709.             IF (T1 .GT. 273.14) THEN
    710. 

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

  711.               Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
    712. 

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

  713.             ELSE
    714. 

  714. ! cold ground temps, use ice saturation only
    715. 

  715.               Q2SAT=Q2SATI
    716. 

  716.               DQSDT2=Q2SATI*6174./(SFCTSNO**2)
    717. 

  717.             ENDIF
    718. 

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

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

  720.           ENDIF
    721. 

  721. 

  722.           IF(ICE.EQ.1)THEN
    723. 

  723.              ! Sea-ice point has deep-level temperature of -2 C
    724. 

  724.              TBOT=271.16
    725. 

  725.           ELSE
    726. 

  726.              ! Land-ice or land points have the usual deep-soil temperature.
    727. 

  727.              TBOT=TMN(I,J)
    728. 

  728.           ENDIF
    729. 

  729.           IF(VEGTYP.EQ.25) SHDFAC=0.0000
    730. 

  730.           IF(VEGTYP.EQ.26) SHDFAC=0.0000
    731. 

  731.           IF(VEGTYP.EQ.27) SHDFAC=0.0000
    732. 

  732.           IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
    733. 

  733. #if 0
    734. 

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

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

  736. #endif
    737. 

  737.             SOILTYP=7
    738. 

  738.           ENDIF
    739. 

  739.           SNOALB1 = SNOALB(I,J)
    740. 

  740.           CMC=CANWAT(I,J)
    741. 

  741. 

  742. !-------------------------------------------
    743. 

  743. !*** convert snow depth from mm to meter
    744. 

  744. !
    745. 

  745. !          IF(RDMAXALB) THEN
    746. 

  746. !           SNOALB=ALBMAX(I,J)*0.01
    747. 

  747. !         ELSE
    748. 

  748. !           SNOALB=MAXALB(IVGTPK)*0.01
    749. 

  749. !         ENDIF
    750. 

  750. 

  751. !        SNOALB1=0.80
    752. 

  752. !        SHMIN=0.00
    753. 

  753.         ALBBRD=ALBBCK(I,J)
    754. 

  754.         Z0BRD=Z0(I,J)
    755. 

  755.         EMBRD=EMBCK(I,J)
    756. 

  756.         SNOTIME1 = SNOTIME(I,J)
    757. 

  757.         RIBB=RIB(I,J)
    758. 

  758. !FEI: temporaray arrays above need to be changed later by using SI
    759. 

  759. 

  760.           DO NS=1,NSOIL
    761. 

  761.             SMC(NS)=SMOIS(I,NS,J)
    762. 

  762.             STC(NS)=TSLB(I,NS,J)                                          !STEMP
    763. 

  763.             SWC(NS)=SH2O(I,NS,J)
    764. 

  764.           ENDDO
    765. 

  765. !
    766. 

  766.           if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
    767. 

  767.             SNOWHK= 5.*SNEQV
    768. 

  768.           endif
    769. 

  769. !
    770. 

  770. 

  771. !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
    772. 

  772. ! the "NATURAL" category in the VEGPARM.TBL
    773. 

  773. 	
    774. 

  774. 	   IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
    775. 

  775.                 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
    776. 

  776.                   IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
    777. 

  777. 		 VEGTYP = NATURAL
    778. 

  778.                  SHDFAC = SHDTBL(NATURAL)
    779. 

  779.                  ALBEDOK =0.2         !  0.2
    780. 

  780.                  ALBBRD  =0.2         !0.2
    781. 

  781.                  EMISSI = 0.98                                 !for VEGTYP=5
    782. 

  782. 		 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
    783. 

  783.                    if(sf_urban_physics.eq.1)then
    784. 

  784.            T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
    785. 

  785.                    elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
    786. 

  786.                 r1= (tsk(i,j)**4.)
    787. 

  787.                 r2= frc_urb2d(i,j)*(ts_urb2d(i,j)**4.)
    788. 

  788.                 r3= (1.-frc_urb2d(i,j))
    789. 

  789.                 t1= ((r1-r2)/r3)**.25
    790. 

  790.                    endif
    791. 

  791. 	         ELSE
    792. 

  792. 		 T1 = TSK(I,J)
    793. 

  793.                  ENDIF
    794. 

  794.                 ENDIF
    795. 

  795.            ELSE
    796. 

  796.                  IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
    797. 

  797.                   IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
    798. 

  798.                   VEGTYP = ISURBAN
    799. 

  799.            	 ENDIF
    800. 

  800.            ENDIF
    801. 

  801. 

  802. #if 0
    803. 

  803.           IF(IPRINT) THEN
    804. 

  804. !
    805. 

  805.        print*, 'BEFORE SFLX, in Noahlsm_driver'
    806. 

  806.        print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL,   &
    807. 

  807.        'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
    808. 

  808.         LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN,      &
    809. 

  809.         'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K,   &
    810. 

  810.          'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
    811. 

  811.          'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
    812. 

  812.          'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
    813. 

  813.           TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
    814. 

  814.           STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
    815. 

  815.           'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT,      &
    816. 

  816.           'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC,   &
    817. 

  817.           'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
    818. 

  818.           'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
    819. 

  819.           'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
    820. 

  820.           'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
    821. 

  821.           'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS,  &
    822. 

  822.           'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW,     &
    823. 

  823.           'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
    824. 

  824.           'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
    825. 

  825.            endif
    826. 

  826. #endif
    827. 

  827. 

  828.           IF (rdlai2d) THEN
    829. 

  829.              xlai = lai(i,j)
    830. 

  830.           endif
    831. 

  831. 

  832.     IF ( XICE(I,J) >= XICE_THRESHOLD ) THEN
    833. 

  833. 

  834.        ! Sea-ice case
    835. 

  835. 

  836.        DO NS = 1, NSOIL
    837. 

  837.           SH2O(I,NS,J) = 1.0
    838. 

  838.        ENDDO
    839. 

  839.        LAI(I,J) = 0.01
    840. 

  840. 

  841.        CYCLE ILOOP
    842. 

  842. 

  843.     ELSE
    844. 

  844. 

  845.        ! Land and glacial land.
    846. 

  846. 

  847.        CALL SFLX (I,J,ICE,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH,  &    !C
    848. 

  848.                  LOCAL,                                           &    !L
    849. 

  849.                  LUTYPE, SLTYPE,                                  &    !CL
    850. 

  850.                  LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY,         &    !F
    851. 

  851.                  DUMMY,DUMMY, DUMMY,                              &    !F PRCPRAIN not used
    852. 

  852.                  TH2,Q2SAT,DQSDT2,                                &    !I
    853. 

  853.                  VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX,      &    !I
    854. 

  854.                  ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, &    !S
    855. 

  855.                  CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,&    !H
    856. 

  856.                  ETA,SHEAT, ETA_KINEMATIC,FDOWN,                  &    !O
    857. 

  857.                  EC,EDIR,ET,ETT,ESNOW,DRIP,DEW,                   &    !O
    858. 

  858.                  BETA,ETP,SSOIL,                                  &    !O
    859. 

  859.                  FLX1,FLX2,FLX3,                                  &    !O
    860. 

  860.                  SNOMLT,SNCOVR,                                   &    !O
    861. 

  861.                  RUNOFF1,RUNOFF2,RUNOFF3,                         &    !O
    862. 

  862.                  RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL,             &    !O
    863. 

  863.                  SOILW,SOILM,Q1,SMAV,                             &    !D
    864. 

  864.                  RDLAI2D,USEMONALB,                               &
    865. 

  865.                  SNOTIME1,                                        &
    866. 

  866.                  RIBB,                                            &
    867. 

  867.                  SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)
    868. 

  868.     ENDIF
    869. 

  869. 

  870.        lai(i,j) = xlai
    871. 

  871. 

  872. #if 0
    873. 

  873.           IF(IPRINT) THEN
    874. 

  874. 

  875.        print*, 'AFTER SFLX, in Noahlsm_driver'
    876. 

  876.        print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL,   &
    877. 

  877.        'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
    878. 

  878.         LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN,      &
    879. 

  879.         'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K,   &
    880. 

  880.          'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
    881. 

  881.           'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
    882. 

  882.          'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
    883. 

  883.           TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
    884. 

  884.           STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
    885. 

  885.           'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT,      &
    886. 

  886.           'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC,   &
    887. 

  887.           'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
    888. 

  888.           'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
    889. 

  889.           'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
    890. 

  890.           'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
    891. 

  891.           'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS,  &
    892. 

  892.           'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW,     &
    893. 

  893.           'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
    894. 

  894.           'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
    895. 

  895.            endif
    896. 

  896. #endif
    897. 

  897. 

  898. !***  UPDATE STATE VARIABLES
    899. 

  899.           CANWAT(I,J)=CMC
    900. 

  900.           SNOW(I,J)=SNEQV*1000.
    901. 

  901. !          SNOWH(I,J)=SNOWHK*1000.
    902. 

  902.           SNOWH(I,J)=SNOWHK                   ! SNOWHK in meters
    903. 

  903.           ALBEDO(I,J)=ALBEDOK
    904. 

  904.           ALB_RURAL(I,J)=ALBEDOK
    905. 

  905.           ALBBCK(I,J)=ALBBRD
    906. 

  906.           Z0(I,J)=Z0BRD
    907. 

  907.           EMISS(I,J) = EMISSI
    908. 

  908.           EMISS_RURAL(I,J) = EMISSI
    909. 

  909. ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
    910. 

  910.           ZNT(I,J)=Z0K
    911. 

  911.           TSK(I,J)=T1
    912. 

  912.           TSK_RURAL(I,J)=T1
    913. 

  913.           HFX(I,J)=SHEAT
    914. 

  914.           HFX_RURAL(I,J)=SHEAT
    915. 

  915. ! MEk Jul07 add potential evap accum
    916. 

  916.         POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
    917. 

  917.           QFX(I,J)=ETA_KINEMATIC
    918. 

  918.           QFX_RURAL(I,J)=ETA_KINEMATIC
    919. 

  919.           LH(I,J)=ETA
    920. 

  920.           LH_RURAL(I,J)=ETA
    921. 

  921.           GRDFLX(I,J)=SSOIL
    922. 

  922.           GRDFLX_RURAL(I,J)=SSOIL
    923. 

  923.           SNOWC(I,J)=SNCOVR
    924. 

  924.           CHS2(I,J)=CQS2(I,J)
    925. 

  925.           SNOTIME(I,J) = SNOTIME1
    926. 

  926. !      prevent diagnostic ground q (q1) from being greater than qsat(tsk)
    927. 

  927. !      as happens over snow cover where the cqs2 value also becomes irrelevant
    928. 

  928. !      by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
    929. 

  929.           IF (Q1 .GT. QSFC(I,J)) THEN
    930. 

  930.             CQS2(I,J) = CHS(I,J)
    931. 

  931.           ENDIF
    932. 

  932. !          QSFC(I,J)=Q1
    933. 

  933. ! Convert QSFC back to mixing ratio
    934. 

  934.            QSFC(I,J)= Q1/(1.0-Q1)
    935. 

  935. !
    936. 

  936.            QSFC_RURAL(I,J)= Q1/(1.0-Q1)
    937. 

  937. ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
    938. 

  938. 

  939.           DO 80 NS=1,NSOIL
    940. 

  940.            SMOIS(I,NS,J)=SMC(NS)
    941. 

  941.            TSLB(I,NS,J)=STC(NS)                                        !  STEMP
    942. 

  942.            SH2O(I,NS,J)=SWC(NS)
    943. 

  943.    80     CONTINUE
    944. 

  944. !       ENDIF
    945. 

  945. 

  946.      !
    947. 

  947.      ! Residual of surface energy balance equation terms
    948. 

  948.      !
    949. 

  949.      noahres(i,j) =     ( solnet + lwdn ) - sheat + ssoil - eta - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
    950. 

  950. 

  951. 

  952.         IF (SF_URBAN_PHYSICS == 1 ) THEN                                              ! Beginning of UCM CALL if block
    953. 

  953. !--------------------------------------
    954. 

  954. ! URBAN CANOPY MODEL START - urban
    955. 

  955. !--------------------------------------
    956. 

  956. ! Input variables lsm --> urban
    957. 

  957. 

  958. 

  959.           IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
    960. 

  960.               IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33 ) THEN
    961. 

  961. 

  962. ! Call urban
    963. 

  963. 

  964. !
    965. 

  965.             UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
    966. 

  966. 

  967.             TA_URB    = SFCTMP           ! [K]
    968. 

  968.             QA_URB    = Q2K              ! [kg/kg]
    969. 

  969.             UA_URB    = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
    970. 

  970.             U1_URB    = U_PHY(I,1,J)
    971. 

  971.             V1_URB    = V_PHY(I,1,J)
    972. 

  972.             IF(UA_URB < 1.) UA_URB=1.    ! [m/s]
    973. 

  973.             SSG_URB   = SOLDN            ! [W/m/m]
    974. 

  974.             SSGD_URB  = 0.8*SOLDN        ! [W/m/m]
    975. 

  975.             SSGQ_URB  = SSG_URB-SSGD_URB ! [W/m/m]
    976. 

  976.             LLG_URB   = GLW(I,J)         ! [W/m/m]
    977. 

  977.             RAIN_URB  = RAINBL(I,J)      ! [mm]
    978. 

  978.             RHOO_URB  = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
    979. 

  979.             ZA_URB    = ZLVL             ! [m]
    980. 

  980.             DELT_URB  = DT               ! [sec]
    981. 

  981.             XLAT_URB  = XLAT_URB2D(I,J)  ! [deg]
    982. 

  982.             COSZ_URB  = COSZ_URB2D(I,J)  !
    983. 

  983.             OMG_URB   = OMG_URB2D(I,J)   !
    984. 

  984.             ZNT_URB   = ZNT(I,J)
    985. 

  985. 

  986.             LSOLAR_URB = .FALSE.
    987. 

  987. 

  988.             TR_URB = TR_URB2D(I,J)
    989. 

  989.             TB_URB = TB_URB2D(I,J)
    990. 

  990.             TG_URB = TG_URB2D(I,J)
    991. 

  991.             TC_URB = TC_URB2D(I,J)
    992. 

  992.             QC_URB = QC_URB2D(I,J)
    993. 

  993.             UC_URB = UC_URB2D(I,J)
    994. 

  994. 

  995.             DO K = 1,num_roof_layers
    996. 

  996.               TRL_URB(K) = TRL_URB3D(I,K,J)
    997. 

  997.             END DO
    998. 

  998.             DO K = 1,num_wall_layers
    999. 

  999.               TBL_URB(K) = TBL_URB3D(I,K,J)
    1000. 

  1000.             END DO
    1001. 

  1001.             DO K = 1,num_road_layers
    1002. 

  1002.               TGL_URB(K) = TGL_URB3D(I,K,J)
    1003. 

  1003.             END DO
    1004. 

  1004. 

  1005.             XXXR_URB = XXXR_URB2D(I,J)
    1006. 

  1006.             XXXB_URB = XXXB_URB2D(I,J)
    1007. 

  1007.             XXXG_URB = XXXG_URB2D(I,J)
    1008. 

  1008.             XXXC_URB = XXXC_URB2D(I,J)
    1009. 

  1009. !
    1010. 

  1010. !
    1011. 

  1011. !      Limits to avoid dividing by small number
    1012. 

  1012.             if (CHS(I,J) < 1.0E-02) then
    1013. 

  1013.                CHS(I,J)  = 1.0E-02
    1014. 

  1014.             endif
    1015. 

  1015.             if (CHS2(I,J) < 1.0E-02) then
    1016. 

  1016.                CHS2(I,J)  = 1.0E-02
    1017. 

  1017.             endif
    1018. 

  1018.             if (CQS2(I,J) < 1.0E-02) then
    1019. 

  1019.                CQS2(I,J)  = 1.0E-02
    1020. 

  1020.             endif
    1021. 

  1021. !
    1022. 

  1022.             CHS_URB  = CHS(I,J)
    1023. 

  1023.             CHS2_URB = CHS2(I,J)
    1024. 

  1024.             IF (PRESENT(CMR_SFCDIF)) THEN
    1025. 

  1025.                CMR_URB = CMR_SFCDIF(I,J)
    1026. 

  1026.                CHR_URB = CHR_SFCDIF(I,J)
    1027. 

  1027.                CMC_URB = CMC_SFCDIF(I,J)
    1028. 

  1028.                CHC_URB = CHC_SFCDIF(I,J)
    1029. 

  1029.             ENDIF
    1030. 

  1030. !
    1031. 

  1031. ! Call urban
    1032. 

  1032. 

  1033.             CALL urban(LSOLAR_URB,                                      & ! I
    1034. 

  1034.                        num_roof_layers,num_wall_layers,num_road_layers, & ! C
    1035. 

  1035.                        DZR,DZB,DZG,                                     & ! C
    1036. 

  1036.                        UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
    1037. 

  1037.                        SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB,     & ! I
    1038. 

  1038.                        ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB,              & ! I
    1039. 

  1039.                        XLAT_URB,DELT_URB,ZNT_URB,                       & ! I
    1040. 

  1040.                        CHS_URB, CHS2_URB,                               & ! I
    1041. 

  1041.                        TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB,   & ! H
    1042. 

  1042.                        TRL_URB,TBL_URB,TGL_URB,                         & ! H
    1043. 

  1043.                        XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB,          & ! H
    1044. 

  1044.                        TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB,    & ! O
    1045. 

  1045.                        SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
    1046. 

  1046.                        GZ1OZ0_URB,                                      & !O
    1047. 

  1047.                        CMR_URB, CHR_URB, CMC_URB, CHC_URB,              &
    1048. 

  1048.                        U10_URB, V10_URB, TH2_URB, Q2_URB,               & ! O
    1049. 

  1049.                        UST_URB)                                           !O
    1050. 

  1050. 

  1051. #if 0
    1052. 

  1052.           IF(IPRINT) THEN
    1053. 

  1053. 

  1054.        print*, 'AFTER CALL URBAN'
    1055. 

  1055.        print*,'num_roof_layers',num_roof_layers, 'num_wall_layers',  &
    1056. 

  1056.         num_wall_layers,                                             &
    1057. 

  1057.        'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
    1058. 

  1058.         TA_URB,                                                      &
    1059. 

  1059.         'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB',    &
    1060. 

  1060.          V1_URB,                                                     &
    1061. 

  1061.          'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB,  &
    1062. 

  1062.         'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB,   &
    1063. 

  1063.         'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
    1064. 

  1064.         'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB,   &
    1065. 

  1065.          'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
    1066. 

  1066.          TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB,   &
    1067. 

  1067.           'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB,   &
    1068. 

  1068.          'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
    1069. 

  1069.          'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB',   &
    1070. 

  1070.          LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
    1071. 

  1071.          'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB',   &
    1072. 

  1072.           RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB,          &
    1073. 

  1073.          'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB,      &
    1074. 

  1074.           'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
    1075. 

  1075.            endif
    1076. 

  1076. #endif
    1077. 

  1077. 

  1078.             TS_URB2D(I,J) = TS_URB
    1079. 

  1079. 

  1080.             ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK   ![-]
    1081. 

  1081.             HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT         ![W/m/m]
    1082. 

  1082.             QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
    1083. 

  1083.                      + (1-FRC_URB2D(I,J))*ETA_KINEMATIC                ![kg/m/m/s]
    1084. 

  1084.             LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA            ![W/m/m]
    1085. 

  1085.             GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL       ![W/m/m]
    1086. 

  1086.             TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1            ![K]
    1087. 

  1087.             Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1            ![-]
    1088. 

  1088. ! Convert QSFC back to mixing ratio
    1089. 

  1089.             QSFC(I,J)= Q1/(1.0-Q1)
    1090. 

  1090.             UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J)      ![m/s]
    1091. 

  1091. 

  1092. #if 0
    1093. 

  1093.     IF(IPRINT)THEN
    1094. 

  1094. 

  1095.     print*, ' FRC_URB2D', FRC_URB2D,                        &
    1096. 

  1096.     'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
    1097. 

  1097.     'ALBEDO(I,J)',  ALBEDO(I,J),                  &
    1098. 

  1098.     'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J),  &
    1099. 

  1099.     'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC',  &
    1100. 

  1100.      ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J),                  &
    1101. 

  1101.     'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J),        &
    1102. 

  1102.     'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
    1103. 

  1103.     'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J),          &
    1104. 

  1104.     'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
    1105. 

  1105.      endif
    1106. 

  1106. #endif
    1107. 

  1107. 

  1108. 

  1109. 

  1110. ! Renew Urban State Varialbes
    1111. 

  1111. 

  1112.             TR_URB2D(I,J) = TR_URB
    1113. 

  1113.             TB_URB2D(I,J) = TB_URB
    1114. 

  1114.             TG_URB2D(I,J) = TG_URB
    1115. 

  1115.             TC_URB2D(I,J) = TC_URB
    1116. 

  1116.             QC_URB2D(I,J) = QC_URB
    1117. 

  1117.             UC_URB2D(I,J) = UC_URB
    1118. 

  1118. 

  1119.             DO K = 1,num_roof_layers
    1120. 

  1120.               TRL_URB3D(I,K,J) = TRL_URB(K)
    1121. 

  1121.             END DO
    1122. 

  1122.             DO K = 1,num_wall_layers
    1123. 

  1123.               TBL_URB3D(I,K,J) = TBL_URB(K)
    1124. 

  1124.             END DO
    1125. 

  1125.             DO K = 1,num_road_layers
    1126. 

  1126.               TGL_URB3D(I,K,J) = TGL_URB(K)
    1127. 

  1127.             END DO
    1128. 

  1128.             XXXR_URB2D(I,J) = XXXR_URB
    1129. 

  1129.             XXXB_URB2D(I,J) = XXXB_URB
    1130. 

  1130.             XXXG_URB2D(I,J) = XXXG_URB
    1131. 

  1131.             XXXC_URB2D(I,J) = XXXC_URB
    1132. 

  1132. 

  1133.             SH_URB2D(I,J)    = SH_URB
    1134. 

  1134.             LH_URB2D(I,J)    = LH_URB
    1135. 

  1135.             G_URB2D(I,J)     = G_URB
    1136. 

  1136.             RN_URB2D(I,J)    = RN_URB
    1137. 

  1137.             PSIM_URB2D(I,J)  = PSIM_URB
    1138. 

  1138.             PSIH_URB2D(I,J)  = PSIH_URB
    1139. 

  1139.             GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
    1140. 

  1140.             U10_URB2D(I,J)   = U10_URB
    1141. 

  1141.             V10_URB2D(I,J)   = V10_URB
    1142. 

  1142.             TH2_URB2D(I,J)   = TH2_URB
    1143. 

  1143.             Q2_URB2D(I,J)    = Q2_URB
    1144. 

  1144.             UST_URB2D(I,J)   = UST_URB
    1145. 

  1145.             AKMS_URB2D(I,J)  = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
    1146. 

  1146.             IF (PRESENT(CMR_SFCDIF)) THEN
    1147. 

  1147.                CMR_SFCDIF(I,J) = CMR_URB
    1148. 

  1148.                CHR_SFCDIF(I,J) = CHR_URB
    1149. 

  1149.                CMC_SFCDIF(I,J) = CMC_URB
    1150. 

  1150.                CHC_SFCDIF(I,J) = CHC_URB
    1151. 

  1151.             ENDIF
    1152. 

  1152.           END IF
    1153. 

  1153. 

  1154.          ENDIF                                   ! end of UCM CALL if block
    1155. 

  1155. !--------------------------------------
    1156. 

  1156. ! Urban Part End - urban
    1157. 

  1157. !--------------------------------------
    1158. 

  1158. 

  1159. !***  DIAGNOSTICS
    1160. 

  1160.           SMSTAV(I,J)=SOILW
    1161. 

  1161.           SMSTOT(I,J)=SOILM*1000.
    1162. 

  1162.           DO NS=1,NSOIL
    1163. 

  1163.           SMCREL(I,NS,J)=SMAV(NS)
    1164. 

  1164.           ENDDO
    1165. 

  1165. !         Convert the water unit into mm
    1166. 

  1166.           SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
    1167. 

  1167.           UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
    1168. 

  1168. ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
    1169. 

  1169.           IF(FFROZP.GT.0.5)THEN
    1170. 

  1170.             ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
    1171. 

  1171.           ENDIF
    1172. 

  1172.           IF(SNOW(I,J).GT.0.)THEN
    1173. 

  1173.             ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
    1174. 

  1174. ! accumulated snow-melt energy
    1175. 

  1175.             SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
    1176. 

  1176.           ENDIF
    1177. 

  1177. 

  1178.         ENDIF                                                           ! endif of land-sea test
    1179. 

  1179. 

  1180.       ENDDO ILOOP                                                       ! of I loop
    1181. 

  1181.    ENDDO JLOOP                                                          ! of J loop
    1182. 

  1182. 

  1183.       IF (SF_URBAN_PHYSICS == 2) THEN
    1184. 

  1184. 

  1185. 

  1186.          do j=jts,jte
    1187. 

  1187.          do i=its,ite
    1188. 

  1188.             EMISS_URB(i,j)=0.
    1189. 

  1189.             RL_UP_URB(i,j)=0.
    1190. 

  1190.             RS_ABS_URB(i,j)=0.
    1191. 

  1191.             GRDFLX_URB(i,j)=0.
    1192. 

  1192.             b_q_bep(i,kts:kte,j)=0.
    1193. 

  1193.          end do
    1194. 

  1194.          end do
    1195. 

  1195.        CALL BEP(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy,   &
    1196. 

  1196.                 th_phy,rho,p_phy,swdown,glw,                           &
    1197. 

  1197.                 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
    1198. 

  1198.                 num_urban_layers,                                      &
    1199. 

  1199.                 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,               &
    1200. 

  1200.                 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,             &
    1201. 

  1201.                 a_u_bep,a_v_bep,a_t_bep,                               &
    1202. 

  1202.                 a_e_bep,b_u_bep,b_v_bep,                               &
    1203. 

  1203.                 b_t_bep,b_e_bep,dlg_bep,                               &
    1204. 

  1204.                 dl_u_bep,sf_bep,vl_bep,                                &
    1205. 

  1205.                 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb,             &
    1206. 

  1206.                 ids,ide, jds,jde, kds,kde,                             &
    1207. 

  1207.                 ims,ime, jms,jme, kms,kme,                             &
    1208. 

  1208.                 its,ite, jts,jte, kts,kte )
    1209. 

  1209. 

  1210.        ENDIF
    1211. 

  1211. 

  1212.        
    1213. 

  1213.        IF (SF_URBAN_PHYSICS == 3) THEN
    1214. 

  1214. 

  1215. 

  1216.          do j=jts,jte
    1217. 

  1217.          do i=its,ite
    1218. 

  1218.             EMISS_URB(i,j)=0.
    1219. 

  1219.             RL_UP_URB(i,j)=0.
    1220. 

  1220.             RS_ABS_URB(i,j)=0.
    1221. 

  1221.             GRDFLX_URB(i,j)=0.
    1222. 

  1222.             b_q_bep(i,kts:kte,j)=0.
    1223. 

  1223.          end do
    1224. 

  1224.          end do
    1225. 

  1225.           
    1226. 

  1226.        CALL BEP_BEM(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
    1227. 

  1227.                 th_phy,rho,p_phy,swdown,glw,                           &
    1228. 

  1228.                 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
    1229. 

  1229.                 num_urban_layers,                                      &
    1230. 

  1230.                 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,               &
    1231. 

  1231.                 tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d,       &
    1232. 

  1232.                 tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d,       &
    1233. 

  1233.                 cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d,                 &
    1234. 

  1234.                 sfwin1_urb3d,sfwin2_urb3d,                             &
    1235. 

  1235.                 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,             &
    1236. 

  1236.                 a_u_bep,a_v_bep,a_t_bep,                               &
    1237. 

  1237.                 a_e_bep,b_u_bep,b_v_bep,                               &
    1238. 

  1238.                 b_t_bep,b_e_bep,b_q_bep,dlg_bep,                       &
    1239. 

  1239.                 dl_u_bep,sf_bep,vl_bep,                                &
    1240. 

  1240.                 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb,qv3d,        &
    1241. 

  1241.                 ids,ide, jds,jde, kds,kde,                             &
    1242. 

  1242.                 ims,ime, jms,jme, kms,kme,                             &
    1243. 

  1243.                 its,ite, jts,jte, kts,kte )
    1244. 

  1244. 

  1245.        ENDIF
    1246. 

  1246. 

  1247.     if((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then         !Bep begin
    1248. 

  1248. ! fix the value of the Stefan-Boltzmann constant
    1249. 

  1249.          sigma_sb=5.67e-08
    1250. 

  1250.          do j=jts,jte
    1251. 

  1251.          do i=its,ite
    1252. 

  1252.             UMOM_URB(I,J)=0.
    1253. 

  1253.             VMOM_URB(I,J)=0.
    1254. 

  1254.             HFX_URB(I,J)=0.
    1255. 

  1255.             QFX_URB(I,J)=0.
    1256. 

  1256.          do k=kts,kte
    1257. 

  1257.             a_u_bep(i,k,j)=a_u_bep(i,k,j)*frc_urb2d(i,j)
    1258. 

  1258.             a_v_bep(i,k,j)=a_v_bep(i,k,j)*frc_urb2d(i,j)
    1259. 

  1259.             a_t_bep(i,k,j)=a_t_bep(i,k,j)*frc_urb2d(i,j)
    1260. 

  1260.             a_q_bep(i,k,j)=0.        
    1261. 

  1261.             a_e_bep(i,k,j)=0.
    1262. 

  1262.             b_u_bep(i,k,j)=b_u_bep(i,k,j)*frc_urb2d(i,j)
    1263. 

  1263.             b_v_bep(i,k,j)=b_v_bep(i,k,j)*frc_urb2d(i,j)
    1264. 

  1264.             b_t_bep(i,k,j)=b_t_bep(i,k,j)*frc_urb2d(i,j)
    1265. 

  1265.             b_q_bep(i,k,j)=b_q_bep(i,k,j)*frc_urb2d(i,j)
    1266. 

  1266.             b_e_bep(i,k,j)=b_e_bep(i,k,j)*frc_urb2d(i,j)
    1267. 

  1267.             HFX_URB(I,J)=HFX_URB(I,J)+B_T_BEP(I,K,J)*RHO(I,K,J)*CP*                &
    1268. 

  1268.                           DZ8W(I,K,J)*VL_BEP(I,K,J)
    1269. 

  1269.             QFX_URB(I,J)=QFX_URB(I,J)+B_Q_BEP(I,K,J)*               &
    1270. 

  1270.                           DZ8W(I,K,J)*VL_BEP(I,K,J)
    1271. 

  1271.             UMOM_URB(I,J)=UMOM_URB(I,J)+ (A_U_BEP(I,K,J)*U_PHY(I,K,J)+             &
    1272. 

  1272.                           B_U_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
    1273. 

  1273.             VMOM_URB(I,J)=VMOM_URB(I,J)+ (A_V_BEP(I,K,J)*V_PHY(I,K,J)+             &
    1274. 

  1274.                           B_V_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
    1275. 

  1275.             vl_bep(i,k,j)=(1.-frc_urb2d(i,j))+vl_bep(i,k,j)*frc_urb2d(i,j)
    1276. 

  1276.             sf_bep(i,k,j)=(1.-frc_urb2d(i,j))+sf_bep(i,k,j)*frc_urb2d(i,j)
    1277. 

  1277.          end do
    1278. 

  1278.             a_u_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/   &
    1279. 

  1279.                            ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_u_bep(i,1,j)
    1280. 

  1280.             a_v_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/   &
    1281. 

  1281.                            ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_v_bep(i,1,j)
    1282. 

  1282.             b_t_bep(i,1,j)=(1.-frc_urb2d(i,j))*hfx_rural(i,j)/dz8w(i,1,j)/rho(i,1,j)/CP+ &
    1283. 

  1283.                             b_t_bep(i,1,j)
    1284. 

  1284.             b_q_bep(i,1,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)/dz8w(i,1,j)/rho(i,1,j)+b_q_bep(i,1,j)
    1285. 

  1285.             umom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*u_phy(i,1,j)/               &
    1286. 

  1286.                          ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+umom_urb(i,j)
    1287. 

  1287.             vmom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*v_phy(i,1,j)/               &
    1288. 

  1288.                          ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+vmom_urb(i,j)
    1289. 

  1289.             sf_bep(i,1,j)=1.
    1290. 

  1290.            
    1291. 

  1291. ! compute upward longwave radiation from the rural part and total
    1292. 

  1292. !           rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
    1293. 

  1293. !           rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)   
    1294. 

  1294. !           emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
    1295. 

  1295. ! using the emissivity and the total longwave upward radiation estimate the averaged skin temperature
    1296. 

  1296.            IF (FRC_URB2D(I,J).GT.0.) THEN
    1297. 

  1297.               rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
    1298. 

  1298.               rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)   
    1299. 

  1299.               emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
    1300. 

  1300.               ts_urb2d(i,j)=(max(0.,(-rl_up_urb(i,j)-(1.-emiss_urb(i,j))*glw(i,j))/emiss_urb(i,j)/sigma_sb))**0.25
    1301. 

  1301.               tsk(i,j)=(max(0., (-1.*rl_up_tot-(1.-emiss(i,j))*glw(i,j) )/emiss(i,j)/sigma_sb))**.25
    1302. 

  1302.            rs_abs_tot=(1.-frc_urb2d(i,j))*swdown(i,j)*(1.-albedo(i,j))+frc_urb2d(i,j)*rs_abs_urb(i,j)  
    1303. 

  1303.           if(swdown(i,j).gt.0.)then
    1304. 

  1304.            albedo(i,j)=1.-rs_abs_tot/swdown(i,j)
    1305. 

  1305.           else
    1306. 

  1306.            albedo(i,j)=alb_rural(i,j)
    1307. 

  1307.           endif
    1308. 

  1308. ! rename *_urb to sh_urb2d,lh_urb2d,g_urb2d,rn_urb2d
    1309. 

  1309.          grdflx(i,j)= (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+frc_urb2d(i,j)*grdflx_urb(i,j) 
    1310. 

  1310.          qfx(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)+qfx_urb(i,j)
    1311. 

  1311. !        lh(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)*xlv
    1312. 

  1312.          lh(i,j)=qfx(i,j)*xlv
    1313. 

  1313.          HFX(I,J) = HFX_URB(I,J)+(1-FRC_URB2D(I,J))*HFX_RURAL(I,J)         ![W/m/m]
    1314. 

  1314.             SH_URB2D(I,J)    = HFX_URB(I,J)/FRC_URB2D(I,J)
    1315. 

  1315.             LH_URB2D(I,J)    = qfx_urb(i,j)*xlv
    1316. 

  1316.             G_URB2D(I,J)     = grdflx_urb(i,j)
    1317. 

  1317.             RN_URB2D(I,J)    = rs_abs_urb(i,j)+emiss_urb(i,j)*glw(i,j)-rl_up_urb(i,j)
    1318. 

  1318.             ust(i,j)=(umom**2.+vmom**2.)**.25
    1319. 

  1319. !          if(tsk(i,j).gt.350)write(*,*)'tsk too big!',i,j,tsk(i,j)
    1320. 

  1320. !          if(tsk(i,j).lt.260)write(*,*)'tsk too small!',i,j,tsk(i,j),rl_up_tot,rl_up_urb(i,j),rl_up_rural
    1321. 

  1321. !                    print*,'ivgtyp,i,j,sigma_sb',ivgtyp(i,j),i,j,sigma_sb
    1322. 

  1322. !                    print*,'hfx,lh,qfx,grdflx,ts_urb2d',hfx(i,j),lh(i,j),qfx(i,j),grdflx(i,j),ts_urb2d(i,j)
    1323. 

  1323. !                    print*,'tsk,albedo,emiss',tsk(i,j),albedo(i,j),emiss(i,j)
    1324. 

  1324. !                if(i.eq.56.and.j.eq.29)then
    1325. 

  1325. !                    print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
    1326. 

  1326. !                    print*,'emiss_rural,emiss_urb',emiss_rural(i,j),emiss_urb(i,j)
    1327. 

  1327. !                    print*,'rl_up_rural,rl_up_urb(i,j)',rl_up_rural,rl_up_urb(i,j)
    1328. 

  1328. !                    print*,'tsk_rural,ts_urb2d(i,j),tsk',tsk_rural(i,j),ts_urb2d(i,j),tsk(i,j)
    1329. 

  1329. !                    print*,'reconstruction fei',((emiss(i,j)*tsk(i,j)**4.-frc_urb2d(i,j)*emiss_urb(i,j)*ts_urb2d(i,j)**4.)/(emiss_rural(i,j)*(1.-frc_urb2d(i,j))))**.25
    1330. 

  1330. !                    print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
    1331. 

  1331. !                    print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
    1332. 

  1332. !                    print*,'qfx',qfx(i,j)
    1333. 

  1333. !                    print*,'ts_urb2d',ts_urb2d(i,j)
    1334. 

  1334. !                    print*,'ust',ust(i,j)
    1335. 

  1335. !                    print*,'swdown,glw',swdown(i,j),glw(i,j)
    1336. 

  1336. !                endif
    1337. 

  1337.             else
    1338. 

  1338.               SH_URB2D(I,J)    = 0.
    1339. 

  1339.               LH_URB2D(I,J)    = 0.
    1340. 

  1340.               G_URB2D(I,J)     = 0.
    1341. 

  1341.               RN_URB2D(I,J)    = 0.
    1342. 

  1342.             endif
    1343. 

  1343. !          IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
    1344. 

  1344. !                  IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
    1345. 

  1345. !                    print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
    1346. 

  1346. !                    print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
    1347. 

  1347. !                    print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
    1348. 

  1348. !                    print*,'qfx',qfx(i,j)
    1349. 

  1349. !                    print*,'ts_urb2d',ts_urb2d(i,j)
    1350. 

  1350. !                    print*,'ust',ust(i,j)
    1351. 

  1351. !          endif
    1352. 

  1352.         enddo
    1353. 

  1353.         enddo
    1354. 

  1354. 

  1355. 

  1356.        endif                                                           !Bep end
    1357. 

  1357. 

  1358. !------------------------------------------------------
    1359. 

  1359.    END SUBROUTINE lsm
    1360. 

  1360. !------------------------------------------------------
    1361. 

  1361. 

  1362.   SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV,    &
    1363. 

  1363.                      SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW,        &
    1364. 

  1364.                      ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
    1365. 

  1365.                      MMINLU,                                   &
    1366. 

  1366.                      SNOALB, FNDSOILW, FNDSNOWH, RDMAXALB,     &
    1367. 

  1367.                      num_soil_layers, restart,                 &
    1368. 

  1368.                      allowed_to_read ,                         &
    1369. 

  1369.                      ids,ide, jds,jde, kds,kde,                &
    1370. 

  1370.                      ims,ime, jms,jme, kms,kme,                &
    1371. 

  1371.                      its,ite, jts,jte, kts,kte                 )
    1372. 

  1372. 

  1373.    INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &
    1374. 

  1374.                                     ims,ime, jms,jme, kms,kme,  &
    1375. 

  1375.                                     its,ite, jts,jte, kts,kte
    1376. 

  1376. 

  1377.    INTEGER, INTENT(IN)       ::     num_soil_layers
    1378. 

  1378. 

  1379.    LOGICAL , INTENT(IN) :: restart , allowed_to_read
    1380. 

  1380. 

  1381.    REAL,    DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
    1382. 

  1382. 

  1383.    REAL,    DIMENSION( ims:ime, num_soil_layers, jms:jme )    , &
    1384. 

  1384.             INTENT(INOUT)    ::                          SMOIS, &  !Total soil moisture
    1385. 

  1385.                                                          SH2O,  &  !liquid soil moisture
    1386. 

  1386.                                                          TSLB      !STEMP
    1387. 

  1387. 

  1388.    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
    1389. 

  1389.             INTENT(INOUT)    ::                           SNOW, &
    1390. 

  1390.                                                          SNOWH, &
    1391. 

  1391.                                                          SNOWC, &
    1392. 

  1392.                                                         SNOALB, &
    1393. 

  1393.                                                         CANWAT, &
    1394. 

  1394.                                                         SMSTAV, &
    1395. 

  1395.                                                         SMSTOT, &
    1396. 

  1396.                                                      SFCRUNOFF, &
    1397. 

  1397.                                                       UDRUNOFF, &
    1398. 

  1398.                                                         ACSNOW, &
    1399. 

  1399.                                                         VEGFRA, &
    1400. 

  1400.                                                         ACSNOM
    1401. 

  1401. 

  1402.    INTEGER, DIMENSION( ims:ime, jms:jme )                     , &
    1403. 

  1403.             INTENT(IN)       ::                         IVGTYP, &
    1404. 

  1404.                                                         ISLTYP
    1405. 

  1405.    CHARACTER(LEN=*),  INTENT(IN)      ::                MMINLU
    1406. 

  1406. 

  1407.    LOGICAL, INTENT(IN)       ::                      FNDSOILW , &
    1408. 

  1408.                                                      FNDSNOWH
    1409. 

  1409.    LOGICAL, INTENT(IN)       ::                      RDMAXALB
    1410. 

  1410. 

  1411. 

  1412.    INTEGER                   :: L
    1413. 

  1413.    REAL                      :: BX, SMCMAX, PSISAT, FREE
    1414. 

  1414.    REAL, PARAMETER           :: BLIM = 5.5, HLICE = 3.335E5,    &
    1415. 

  1415.                                 GRAV = 9.81, T0 = 273.15
    1416. 

  1416.    INTEGER                   :: errflag
    1417. 

  1417. 

  1418.    character*256 :: MMINSL
    1419. 

  1419.         MMINSL='STAS'
    1420. 

  1420. !
    1421. 

  1421. 

  1422. ! initialize three Noah LSM related tables
    1423. 

  1423.    IF ( allowed_to_read ) THEN
    1424. 

  1424.      CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
    1425. 

  1425.      CALL  SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
    1426. 

  1426.    ENDIF
    1427. 

  1427. 

  1428. #ifdef WRF_CHEM
    1429. 

  1429. !
    1430. 

  1430. ! need this parameter for dust parameterization in wrf/chem
    1431. 

  1431. !
    1432. 

  1432.    do I=1,NSLTYPE
    1433. 

  1433.       porosity(i)=maxsmc(i)
    1434. 

  1434.       drypoint(i)=drysmc(i)
    1435. 

  1435.    enddo
    1436. 

  1436. #endif
    1437. 

  1437. 

  1438.    IF(.not.restart)THEN
    1439. 

  1439. 

  1440.    itf=min0(ite,ide-1)
    1441. 

  1441.    jtf=min0(jte,jde-1)
    1442. 

  1442. 

  1443.    errflag = 0
    1444. 

  1444.    DO j = jts,jtf
    1445. 

  1445.      DO i = its,itf
    1446. 

  1446.        IF ( ISLTYP( i,j ) .LT. 1 ) THEN
    1447. 

  1447.          errflag = 1
    1448. 

  1448.          WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
    1449. 

  1449.          CALL wrf_message(err_message)
    1450. 

  1450.        ENDIF
    1451. 

  1451.        IF(.not.RDMAXALB) THEN
    1452. 

  1452.           SNOALB(i,j)=MAXALB(IVGTYP(i,j))*0.01
    1453. 

  1453.        ENDIF
    1454. 

  1454.      ENDDO
    1455. 

  1455.    ENDDO
    1456. 

  1456.    IF ( errflag .EQ. 1 ) THEN
    1457. 

  1457.       CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
    1458. 

  1458.                             "of ISLTYP. Is this field in the input?" )
    1459. 

  1459.    ENDIF
    1460. 

  1460. 

  1461. ! initialize soil liquid water content SH2O
    1462. 

  1462. 

  1463. !  IF(.NOT.FNDSOILW) THEN
    1464. 

  1464. 

  1465. ! If no SWC, do the following
    1466. 

  1466. !         PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
    1467. 

  1467.         DO J = jts,jtf
    1468. 

  1468.         DO I = its,itf
    1469. 

  1469.           BX = BB(ISLTYP(I,J))
    1470. 

  1470.           SMCMAX = MAXSMC(ISLTYP(I,J))
    1471. 

  1471.           PSISAT = SATPSI(ISLTYP(I,J))
    1472. 

  1472.          if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
    1473. 

  1473.           DO NS=1, num_soil_layers
    1474. 

  1474. ! ----------------------------------------------------------------------
    1475. 

  1475. !SH2O  <= SMOIS for T < 273.149K (-0.001C)
    1476. 

  1476.              IF (TSLB(I,NS,J) < 273.149) THEN
    1477. 

  1477. ! ----------------------------------------------------------------------
    1478. 

  1478. ! first guess following explicit solution for Flerchinger Eqn from Koren
    1479. 

  1479. ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
    1480. 

  1480. ! ISLTPK is soil type
    1481. 

  1481.               BX = BB(ISLTYP(I,J))
    1482. 

  1482.               SMCMAX = MAXSMC(ISLTYP(I,J))
    1483. 

  1483.               PSISAT = SATPSI(ISLTYP(I,J))
    1484. 

  1484.               IF ( BX >  BLIM ) BX = BLIM
    1485. 

  1485.               FK=(( (HLICE/(GRAV*(-PSISAT))) *                              &
    1486. 

  1486.                  ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
    1487. 

  1487.               IF (FK < 0.02) FK = 0.02
    1488. 

  1488.               SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
    1489. 

  1489. ! ----------------------------------------------------------------------
    1490. 

  1490. ! now use iterative solution for liquid soil water content using
    1491. 

  1491. ! FUNCTION FRH2O with the initial guess for SH2O from above explicit
    1492. 

  1492. ! first guess.
    1493. 

  1493.               CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J),    &
    1494. 

  1494.                  SMCMAX,BX,PSISAT)
    1495. 

  1495.               SH2O(I,NS,J) = FREE
    1496. 

  1496.              ELSE             ! of IF (TSLB(I,NS,J)
    1497. 

  1497. ! ----------------------------------------------------------------------
    1498. 

  1498. ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
    1499. 

  1499.               SH2O(I,NS,J)=SMOIS(I,NS,J)
    1500. 

  1500. ! ----------------------------------------------------------------------
    1501. 

  1501.              ENDIF            ! of IF (TSLB(I,NS,J)
    1502. 

  1502.           END DO              ! of DO NS=1, num_soil_layers
    1503. 

  1503.          else                 ! of if ((bx > 0.0)
    1504. 

  1504.           DO NS=1, num_soil_layers
    1505. 

  1505.            SH2O(I,NS,J)=SMOIS(I,NS,J)
    1506. 

  1506.           END DO
    1507. 

  1507.          endif                ! of if ((bx > 0.0)
    1508. 

  1508.         ENDDO                 ! DO I = its,itf
    1509. 

  1509.         ENDDO                 ! DO J = jts,jtf
    1510. 

  1510. !  ENDIF                       ! of IF(.NOT.FNDSOILW)THEN
    1511. 

  1511. 

  1512. ! initialize physical snow height SNOWH
    1513. 

  1513. 

  1514.         IF(.NOT.FNDSNOWH)THEN
    1515. 

  1515. ! If no SNOWH do the following
    1516. 

  1516.           CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
    1517. 

  1517.           DO J = jts,jtf
    1518. 

  1518.           DO I = its,itf
    1519. 

  1519.             SNOWH(I,J)=SNOW(I,J)*0.005               ! SNOW in mm and SNOWH in m
    1520. 

  1520.           ENDDO
    1521. 

  1521.           ENDDO
    1522. 

  1522.         ENDIF
    1523. 

  1523. 

  1524. ! initialize canopy water to ZERO
    1525. 

  1525. 

  1526. !          GO TO 110
    1527. 

  1527. !         print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
    1528. 

  1528.           DO J = jts,jtf
    1529. 

  1529.           DO I = its,itf
    1530. 

  1530.             CANWAT(I,J)=0.0
    1531. 

  1531.           ENDDO
    1532. 

  1532.           ENDDO
    1533. 

  1533.  110      CONTINUE
    1534. 

  1534. 

  1535.    ENDIF
    1536. 

  1536. !------------------------------------------------------------------------------
    1537. 

  1537.   END SUBROUTINE lsminit
    1538. 

  1538. !------------------------------------------------------------------------------
    1539. 

  1539. 

  1540. 

  1541. 

  1542. !-----------------------------------------------------------------
    1543. 

  1543.         SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
    1544. 

  1544. !-----------------------------------------------------------------
    1545. 

  1545. 

  1546.         USE module_wrf_error
    1547. 

  1547.         IMPLICIT NONE
    1548. 

  1548. 

  1549.         CHARACTER(LEN=*), INTENT(IN) :: MMINLU, MMINSL
    1550. 

  1550.         integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
    1551. 

  1551.         integer :: ierr
    1552. 

  1552.         INTEGER , PARAMETER :: OPEN_OK = 0
    1553. 

  1553. 

  1554.         character*128 :: mess , message
    1555. 

  1555.         logical, external :: wrf_dm_on_monitor
    1556. 

  1556. 

  1557. 

  1558. !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
    1559. 

  1559. !             ALBBCK: SFC albedo (in percentage)
    1560. 

  1560. !                 Z0: Roughness length (m)
    1561. 

  1561. !             SHDFAC: Green vegetation fraction (in percentage)
    1562. 

  1562. !  Note: The ALBEDO, Z0, and SHDFAC values read from the following table
    1563. 

  1563. !          ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
    1564. 

  1564. !          the monthly green vegetation data
    1565. 

  1565. !             CMXTBL: MAX CNPY Capacity (m)
    1566. 

  1566. !             NROTBL: Rooting depth (layer)
    1567. 

  1567. !              RSMIN: Mimimum stomatal resistance (s m-1)
    1568. 

  1568. !              RSMAX: Max. stomatal resistance (s m-1)
    1569. 

  1569. !                RGL: Parameters used in radiation stress function
    1570. 

  1570. !                 HS: Parameter used in vapor pressure deficit functio
    1571. 

  1571. !               TOPT: Optimum transpiration air temperature. (K)
    1572. 

  1572. !             CMCMAX: Maximum canopy water capacity
    1573. 

  1573. !             CFACTR: Parameter used in the canopy inteception calculati
    1574. 

  1574. !               SNUP: Threshold snow depth (in water equivalent m) that
    1575. 

  1575. !                     implies 100% snow cover
    1576. 

  1576. !                LAI: Leaf area index (dimensionless)
    1577. 

  1577. !             MAXALB: Upper bound on maximum albedo over deep snow
    1578. 

  1578. !
    1579. 

  1579. !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
    1580. 

  1580. !
    1581. 

  1581. 

  1582.        IF ( wrf_dm_on_monitor() ) THEN
    1583. 

  1583. 

  1584.         OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
    1585. 

  1585.         IF(ierr .NE. OPEN_OK ) THEN
    1586. 

  1586.           WRITE(message,FMT='(A)') &
    1587. 

  1587.           'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
    1588. 

  1588.           CALL wrf_error_fatal ( message )
    1589. 

  1589.         END IF
    1590. 

  1590. 

  1591. 

  1592.         LUMATCH=0
    1593. 

  1593. 

  1594.         FIND_LUTYPE : DO WHILE (LUMATCH == 0)
    1595. 

  1595.            READ (19,*,END=2002)
    1596. 

  1596.            READ (19,*,END=2002)LUTYPE
    1597. 

  1597.            READ (19,*)LUCATS,IINDEX
    1598. 

  1598. 

  1599.            IF(LUTYPE.EQ.MMINLU)THEN
    1600. 

  1600.               WRITE( mess , * ) 'LANDUSE TYPE = ' // TRIM ( LUTYPE ) // ' FOUND', LUCATS,' CATEGORIES'
    1601. 

  1601.               CALL wrf_message( mess )
    1602. 

  1602.               LUMATCH=1
    1603. 

  1603.            ELSE
    1604. 

  1604.               call wrf_message ( "Skipping over LUTYPE = " // TRIM ( LUTYPE ) )
    1605. 

  1605.               DO LC = 1, LUCATS+12
    1606. 

  1606.                  read(19,*)
    1607. 

  1607.               ENDDO
    1608. 

  1608.            ENDIF
    1609. 

  1609.         ENDDO FIND_LUTYPE
    1610. 

  1610. ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
    1611. 

  1611.         IF ( SIZE(SHDTBL)       < LUCATS .OR. &
    1612. 

  1612.              SIZE(NROTBL)       < LUCATS .OR. &
    1613. 

  1613.              SIZE(RSTBL)        < LUCATS .OR. &
    1614. 

  1614.              SIZE(RGLTBL)       < LUCATS .OR. &
    1615. 

  1615.              SIZE(HSTBL)        < LUCATS .OR. &
    1616. 

  1616.              SIZE(SNUPTBL)      < LUCATS .OR. &
    1617. 

  1617.              SIZE(MAXALB)       < LUCATS .OR. &
    1618. 

  1618.              SIZE(LAIMINTBL)    < LUCATS .OR. &
    1619. 

  1619.              SIZE(LAIMAXTBL)    < LUCATS .OR. &
    1620. 

  1620.              SIZE(Z0MINTBL)     < LUCATS .OR. &
    1621. 

  1621.              SIZE(Z0MAXTBL)     < LUCATS .OR. &
    1622. 

  1622.              SIZE(ALBEDOMINTBL) < LUCATS .OR. &
    1623. 

  1623.              SIZE(ALBEDOMAXTBL) < LUCATS .OR. &
    1624. 

  1624.              SIZE(EMISSMINTBL ) < LUCATS .OR. &
    1625. 

  1625.              SIZE(EMISSMAXTBL ) < LUCATS ) THEN
    1626. 

  1626.            CALL wrf_error_fatal('Table sizes too small for value of LUCATS in module_sf_noahdrv.F')
    1627. 

  1627.         ENDIF
    1628. 

  1628. 

  1629.         IF(LUTYPE.EQ.MMINLU)THEN
    1630. 

  1630.           DO LC=1,LUCATS
    1631. 

  1631.               READ (19,*)IINDEX,SHDTBL(LC),                        &
    1632. 

  1632.                         NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
    1633. 

  1633.                         SNUPTBL(LC),MAXALB(LC), LAIMINTBL(LC),     &
    1634. 

  1634.                         LAIMAXTBL(LC),EMISSMINTBL(LC),             &
    1635. 

  1635.                         EMISSMAXTBL(LC), ALBEDOMINTBL(LC),         &
    1636. 

  1636.                         ALBEDOMAXTBL(LC), Z0MINTBL(LC), Z0MAXTBL(LC)
    1637. 

  1637.           ENDDO
    1638. 

  1638. !
    1639. 

  1639.           READ (19,*)
    1640. 

  1640.           READ (19,*)TOPT_DATA
    1641. 

  1641.           READ (19,*)
    1642. 

  1642.           READ (19,*)CMCMAX_DATA
    1643. 

  1643.           READ (19,*)
    1644. 

  1644.           READ (19,*)CFACTR_DATA
    1645. 

  1645.           READ (19,*)
    1646. 

  1646.           READ (19,*)RSMAX_DATA
    1647. 

  1647.           READ (19,*)
    1648. 

  1648.           READ (19,*)BARE
    1649. 

  1649.           READ (19,*)
    1650. 

  1650.           READ (19,*)NATURAL
    1651. 

  1651.         ENDIF
    1652. 

  1652. !
    1653. 

  1653.  2002   CONTINUE
    1654. 

  1654. 

  1655.         CLOSE (19)
    1656. 

  1656.         IF (LUMATCH == 0) then
    1657. 

  1657.            CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
    1658. 

  1658.         ENDIF
    1659. 

  1659.       ENDIF
    1660. 

  1660. 

  1661.       CALL wrf_dm_bcast_string  ( LUTYPE  , 4 )
    1662. 

  1662.       CALL wrf_dm_bcast_integer ( LUCATS  , 1 )
    1663. 

  1663.       CALL wrf_dm_bcast_integer ( IINDEX  , 1 )
    1664. 

  1664.       CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
    1665. 

  1665.       CALL wrf_dm_bcast_real    ( SHDTBL  , NLUS )
    1666. 

  1666.       CALL wrf_dm_bcast_real    ( NROTBL  , NLUS )
    1667. 

  1667.       CALL wrf_dm_bcast_real    ( RSTBL   , NLUS )
    1668. 

  1668.       CALL wrf_dm_bcast_real    ( RGLTBL  , NLUS )
    1669. 

  1669.       CALL wrf_dm_bcast_real    ( HSTBL   , NLUS )
    1670. 

  1670.       CALL wrf_dm_bcast_real    ( SNUPTBL , NLUS )
    1671. 

  1671.       CALL wrf_dm_bcast_real    ( LAIMINTBL    , NLUS )
    1672. 

  1672.       CALL wrf_dm_bcast_real    ( LAIMAXTBL    , NLUS )
    1673. 

  1673.       CALL wrf_dm_bcast_real    ( Z0MINTBL     , NLUS )
    1674. 

  1674.       CALL wrf_dm_bcast_real    ( Z0MAXTBL     , NLUS )
    1675. 

  1675.       CALL wrf_dm_bcast_real    ( EMISSMINTBL  , NLUS )
    1676. 

  1676.       CALL wrf_dm_bcast_real    ( EMISSMAXTBL  , NLUS )
    1677. 

  1677.       CALL wrf_dm_bcast_real    ( ALBEDOMINTBL , NLUS )
    1678. 

  1678.       CALL wrf_dm_bcast_real    ( ALBEDOMAXTBL , NLUS )
    1679. 

  1679.       CALL wrf_dm_bcast_real    ( MAXALB  , NLUS )
    1680. 

  1680.       CALL wrf_dm_bcast_real    ( TOPT_DATA    , 1 )
    1681. 

  1681.       CALL wrf_dm_bcast_real    ( CMCMAX_DATA  , 1 )
    1682. 

  1682.       CALL wrf_dm_bcast_real    ( CFACTR_DATA  , 1 )
    1683. 

  1683.       CALL wrf_dm_bcast_real    ( RSMAX_DATA  , 1 )
    1684. 

  1684.       CALL wrf_dm_bcast_integer ( BARE    , 1 )
    1685. 

  1685.       CALL wrf_dm_bcast_integer ( NATURAL , 1 )
    1686. 

  1686. 

  1687. !
    1688. 

  1688. !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
    1689. 

  1689. !
    1690. 

  1690.       IF ( wrf_dm_on_monitor() ) THEN
    1691. 

  1691.         OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
    1692. 

  1692.         IF(ierr .NE. OPEN_OK ) THEN
    1693. 

  1693.           WRITE(message,FMT='(A)') &
    1694. 

  1694.           'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
    1695. 

  1695.           CALL wrf_error_fatal ( message )
    1696. 

  1696.         END IF
    1697. 

  1697. 

  1698.         WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ', TRIM ( MMINSL )
    1699. 

  1699.         CALL wrf_message( mess )
    1700. 

  1700. 

  1701.         LUMATCH=0
    1702. 

  1702. 

  1703.         READ (19,*)
    1704. 

  1704.         READ (19,2000,END=2003)SLTYPE
    1705. 

  1705.  2000   FORMAT (A4)
    1706. 

  1706.         READ (19,*)SLCATS,IINDEX
    1707. 

  1707.         IF(SLTYPE.EQ.MMINSL)THEN
    1708. 

  1708.             WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ', TRIM ( SLTYPE ) , ' FOUND', &
    1709. 

  1709.                   SLCATS,' CATEGORIES'
    1710. 

  1710.             CALL wrf_message ( mess )
    1711. 

  1711.           LUMATCH=1
    1712. 

  1712.         ENDIF
    1713. 

  1713. ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
    1714. 

  1714.         IF ( SIZE(BB    ) < SLCATS .OR. &
    1715. 

  1715.              SIZE(DRYSMC) < SLCATS .OR. &
    1716. 

  1716.              SIZE(F11   ) < SLCATS .OR. &
    1717. 

  1717.              SIZE(MAXSMC) < SLCATS .OR. &
    1718. 

  1718.              SIZE(REFSMC) < SLCATS .OR. &
    1719. 

  1719.              SIZE(SATPSI) < SLCATS .OR. &
    1720. 

  1720.              SIZE(SATDK ) < SLCATS .OR. &
    1721. 

  1721.              SIZE(SATDW ) < SLCATS .OR. &
    1722. 

  1722.              SIZE(WLTSMC) < SLCATS .OR. &
    1723. 

  1723.              SIZE(QTZ   ) < SLCATS  ) THEN
    1724. 

  1724.            CALL wrf_error_fatal('Table sizes too small for value of SLCATS in module_sf_noahdrv.F')
    1725. 

  1725.         ENDIF
    1726. 

  1726.         IF(SLTYPE.EQ.MMINSL)THEN
    1727. 

  1727.           DO LC=1,SLCATS
    1728. 

  1728.               READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
    1729. 

  1729.                         REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC),   &
    1730. 

  1730.                         WLTSMC(LC), QTZ(LC)
    1731. 

  1731.           ENDDO
    1732. 

  1732.         ENDIF
    1733. 

  1733. 

  1734.  2003   CONTINUE
    1735. 

  1735. 

  1736.         CLOSE (19)
    1737. 

  1737.       ENDIF
    1738. 

  1738. 

  1739.       CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
    1740. 

  1740.       CALL wrf_dm_bcast_string  ( SLTYPE  , 4 )
    1741. 

  1741.       CALL wrf_dm_bcast_string  ( MMINSL  , 4 )  ! since this is reset above, see oct2 ^
    1742. 

  1742.       CALL wrf_dm_bcast_integer ( SLCATS  , 1 )
    1743. 

  1743.       CALL wrf_dm_bcast_integer ( IINDEX  , 1 )
    1744. 

  1744.       CALL wrf_dm_bcast_real    ( BB      , NSLTYPE )
    1745. 

  1745.       CALL wrf_dm_bcast_real    ( DRYSMC  , NSLTYPE )
    1746. 

  1746.       CALL wrf_dm_bcast_real    ( F11     , NSLTYPE )
    1747. 

  1747.       CALL wrf_dm_bcast_real    ( MAXSMC  , NSLTYPE )
    1748. 

  1748.       CALL wrf_dm_bcast_real    ( REFSMC  , NSLTYPE )
    1749. 

  1749.       CALL wrf_dm_bcast_real    ( SATPSI  , NSLTYPE )
    1750. 

  1750.       CALL wrf_dm_bcast_real    ( SATDK   , NSLTYPE )
    1751. 

  1751.       CALL wrf_dm_bcast_real    ( SATDW   , NSLTYPE )
    1752. 

  1752.       CALL wrf_dm_bcast_real    ( WLTSMC  , NSLTYPE )
    1753. 

  1753.       CALL wrf_dm_bcast_real    ( QTZ     , NSLTYPE )
    1754. 

  1754. 

  1755.       IF(LUMATCH.EQ.0)THEN
    1756. 

  1756.           CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
    1757. 

  1757.           CALL wrf_message( 'MATCH SOILPARM TABLE'                 )
    1758. 

  1758.           CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
    1759. 

  1759.       ENDIF
    1760. 

  1760. 

  1761. !
    1762. 

  1762. !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
    1763. 

  1763. !
    1764. 

  1764.       IF ( wrf_dm_on_monitor() ) THEN
    1765. 

  1765.         OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
    1766. 

  1766.         IF(ierr .NE. OPEN_OK ) THEN
    1767. 

  1767.           WRITE(message,FMT='(A)') &
    1768. 

  1768.           'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
    1769. 

  1769.           CALL wrf_error_fatal ( message )
    1770. 

  1770.         END IF
    1771. 

  1771. 

  1772.         READ (19,*)
    1773. 

  1773.         READ (19,*)
    1774. 

  1774.         READ (19,*) NUM_SLOPE
    1775. 

  1775. 

  1776.           SLPCATS=NUM_SLOPE
    1777. 

  1777. ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
    1778. 

  1778.           IF ( SIZE(slope_data) < NUM_SLOPE ) THEN
    1779. 

  1779.             CALL wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
    1780. 

  1780.           ENDIF
    1781. 

  1781. 

  1782.           DO LC=1,SLPCATS
    1783. 

  1783.               READ (19,*)SLOPE_DATA(LC)
    1784. 

  1784.           ENDDO
    1785. 

  1785. 

  1786.           READ (19,*)
    1787. 

  1787.           READ (19,*)SBETA_DATA
    1788. 

  1788.           READ (19,*)
    1789. 

  1789.           READ (19,*)FXEXP_DATA
    1790. 

  1790.           READ (19,*)
    1791. 

  1791.           READ (19,*)CSOIL_DATA
    1792. 

  1792.           READ (19,*)
    1793. 

  1793.           READ (19,*)SALP_DATA
    1794. 

  1794.           READ (19,*)
    1795. 

  1795.           READ (19,*)REFDK_DATA
    1796. 

  1796.           READ (19,*)
    1797. 

  1797.           READ (19,*)REFKDT_DATA
    1798. 

  1798.           READ (19,*)
    1799. 

  1799.           READ (19,*)FRZK_DATA
    1800. 

  1800.           READ (19,*)
    1801. 

  1801.           READ (19,*)ZBOT_DATA
    1802. 

  1802.           READ (19,*)
    1803. 

  1803.           READ (19,*)CZIL_DATA
    1804. 

  1804.           READ (19,*)
    1805. 

  1805.           READ (19,*)SMLOW_DATA
    1806. 

  1806.           READ (19,*)
    1807. 

  1807.           READ (19,*)SMHIGH_DATA
    1808. 

  1808.           READ (19,*)
    1809. 

  1809.           READ (19,*)LVCOEF_DATA
    1810. 

  1810.         CLOSE (19)
    1811. 

  1811.       ENDIF
    1812. 

  1812. 

  1813.       CALL wrf_dm_bcast_integer ( NUM_SLOPE    ,  1 )
    1814. 

  1814.       CALL wrf_dm_bcast_integer ( SLPCATS      ,  1 )
    1815. 

  1815.       CALL wrf_dm_bcast_real    ( SLOPE_DATA   ,  NSLOPE )
    1816. 

  1816.       CALL wrf_dm_bcast_real    ( SBETA_DATA   ,  1 )
    1817. 

  1817.       CALL wrf_dm_bcast_real    ( FXEXP_DATA   ,  1 )
    1818. 

  1818.       CALL wrf_dm_bcast_real    ( CSOIL_DATA   ,  1 )
    1819. 

  1819.       CALL wrf_dm_bcast_real    ( SALP_DATA    ,  1 )
    1820. 

  1820.       CALL wrf_dm_bcast_real    ( REFDK_DATA   ,  1 )
    1821. 

  1821.       CALL wrf_dm_bcast_real    ( REFKDT_DATA  ,  1 )
    1822. 

  1822.       CALL wrf_dm_bcast_real    ( FRZK_DATA    ,  1 )
    1823. 

  1823.       CALL wrf_dm_bcast_real    ( ZBOT_DATA    ,  1 )
    1824. 

  1824.       CALL wrf_dm_bcast_real    ( CZIL_DATA    ,  1 )
    1825. 

  1825.       CALL wrf_dm_bcast_real    ( SMLOW_DATA   ,  1 )
    1826. 

  1826.       CALL wrf_dm_bcast_real    ( SMHIGH_DATA  ,  1 )
    1827. 

  1827.       CALL wrf_dm_bcast_real    ( LVCOEF_DATA  ,  1 )
    1828. 

  1828. 

  1829. 

  1830. !-----------------------------------------------------------------
    1831. 

  1831.       END SUBROUTINE SOIL_VEG_GEN_PARM
    1832. 

  1832. !-----------------------------------------------------------------
    1833. 

  1833. 

  1834. END MODULE module_sf_noahdrv
    1835. 

  1835.