PageRenderTime 65ms CodeModel.GetById 21ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_sf_pxlsm.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1703 lines | 900 code | 218 blank | 585 comment | 8 complexity | f2952415ca981168cd0e25858a38680d MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. MODULE module_sf_pxlsm
    2. 

  2.  
    3. 

  3.   USE module_model_constants
    4. 

  4.   USE module_sf_pxlsm_data
    5. 

  5. 

  6.   INTEGER, PARAMETER   :: NSOLD=20
    7. 

  7.   REAL, PARAMETER      :: RD     = 287.04,   CPD   = 1004.67,             &
    8. 

  8.                           CPH2O  = 4.218E+3, CPICE = 2.106E+3,            &
    9. 

  9.                           LSUBF  = 3.335E+5, SIGMA = 5.67E-8,             &  
    10. 

  10.                           ROVCP  = RD / CPD
    11. 

  11.                           
    12. 

  12.   REAL, PARAMETER      :: CRANKP = 0.5                    ! CRANK-NIC PARAMETER
    13. 

  13.   REAL, PARAMETER      :: RIC    = 0.25                   ! critical Richardson number
    14. 

  14.   REAL, PARAMETER      :: DENW   = 1000.0                 ! water density in KG/M3                  
    15. 

  15.   REAL, PARAMETER      :: TAUINV = 1.0 / 86400.0          ! 1/1DAY(SEC)
    16. 

  16.   REAL, PARAMETER      :: T2TFAC = 1.0 / 10.0             ! Bottom soil temp response factor
    17. 

  17.   REAL, PARAMETER      :: PI = 3.1415926
    18. 

  18.   REAL, PARAMETER      :: PR0 = 0.95
    19. 

  19.   REAL, PARAMETER      :: CZO    = 0.032
    20. 

  20.   REAL, PARAMETER      :: OZO    = 1.E-4
    21. 

  21. 

  22. CONTAINS
    23. 

  23. !
    24. 

  24. 

  25. !-------------------------------------------------------------------------
    26. 

  26. !-------------------------------------------------------------------------
    27. 

  27.    SUBROUTINE pxlsm(U3D, V3D, DZ8W, QV3D,T3D,TH3D, RHO,         &   
    28. 

  28.                     PSFC, GSW, GLW, RAINBL, EMISS,              &
    29. 

  29.                     ITIMESTEP, NSOIL, DT, ANAL_INTERVAL,        &            
    30. 

  30.                     XLAND, XICE, ALBBCK, ALBEDO, SNOALB,        &
    31. 

  31.                     SMOIS, TSLB, MAVAIL, TA2, QA2,              &
    32. 

  32.                     ZS,DZS, PSIH,                               &
    33. 

  33.                     LANDUSEF,SOILCTOP,SOILCBOT,VEGFRA,VEGF_PX,  &
    34. 

  34.                     ISLTYP,RA,RS,LAI,NLCAT,NSCAT,               & 
    35. 

  35.                     HFX,QFX,LH,TSK,SST,ZNT,CANWAT,              &
    36. 

  36.                     GRDFLX,SHDMIN,SHDMAX,                       &
    37. 

  37.                     SNOWC,PBLH,RMOL,UST,CAPG,DTBL,              &
    38. 

  38.                     T2_NDG_OLD, T2_NDG_NEW,                     &     
    39. 

  39.                     Q2_NDG_OLD, Q2_NDG_NEW,                     & 
    40. 

  40.                     SN_NDG_OLD, SN_NDG_NEW, SNOW, SNOWH,SNOWNCV,&
    41. 

  41.                     T2OBS, Q2OBS, PXLSM_SMOIS_INIT,             &           
    42. 

  42.                     PXLSM_SOIL_NUDGE,                           &
    43. 

  43.                     ids,ide, jds,jde, kds,kde,                  &
    44. 

  44.                     ims,ime, jms,jme, kms,kme,                  &
    45. 

  45.                     its,ite, jts,jte, kts,kte                    )
    46. 

  46. 

  47. !-------------------------------------------------------------------------
    48. 

  48. !   THIS MODULE CONTAINS THE PLEIM-XIU LAND-SURFACE MODEL (PX-LSM).
    49. 

  49. !   IT IS DESIGNED TO SIMULATE CHARACTERISTICS OF THE LAND SURFACE AND
    50. 

  50. !   VEGETATION AND EXCHANGE WITH THE PLANETARY BOUNDARY LAYER (PBL). THE
    51. 

  51. !   SOIL MOISTURE MODEL IS BASED ON THE ISBA SCHEME DEVELOPED BY NOILHAN
    52. 

  52. !   AND PLANTON (1989) AND JACQUEMIN AND NOILHAN (1990) AND INCLUDES
    53. 

  53. !   PROGNOSTIC EQUATIONS FOR SOIL MOISTURE AND SOIL TEMPERATURE IN TWO
    54. 

  54. !   LAYERS (1 CM AND 1 M) AS WELL AS CANOPY WATER CONTENT.  SURFACE
    55. 

  55. !   MOISTURE FLUXES ARE MODELED BY 3 PATHWAYS: SOIL EVAPORATION, CANOPY
    56. 

  56. !   EVAPORATION, AND VEGETATIVE EVAPOTRANSPIRRATION.
    57. 

  57. !   EVAPOTRANSPIRATION DIRECTLY FROM THE ROOT ZONE SOIL LAYER IS MODELED
    58. 

  58. !   VIA A CANOPY RESISTANCE ANALOG ALGORITHM WHERE STOMATAL CONDUCTANCE
    59. 

  59. !   IS CONTROLLED BY SOLAR RADIATION, AIR TEMPERATURE, AIR HUMIDITY, AND
    60. 

  60. !   ROOT ZONE SOIL MOISTURE. REQUIRED VEGETATION CHARACTERISTICS DERIVED
    61. 

  61. !   FROM THE USGS LANDUSE DATA INCLUDE: LEAF AREA INDEX, FRACTIONAL VEGETATION
    62. 

  62. !   COVERAGE, ROUGHNESS LENGTH, AND MINIMUM STOMATAL RESISTANCE. AN INDIRECT
    63. 

  63. !   NUDGING SCHEME ADJUSTS SOIL MOISTURE ACCORDING TO DIFFERENCES BETWEEN
    64. 

  64. !   MODELED TEMPERATURE AND HUMIDITY AND ANALYSED SURFACE FIELDS.
    65. 

  65. !
    66. 

  66. ! References:
    67. 

  67. ! Pleim and Xiu, 1995: Development and testing of a surface flux and planetary 
    68. 

  68. !                      boundary layer model for application in mesoscale models.
    69. 

  69. !                      J. Appl. Meteoro., Vol. 34, 16-32.
    70. 

  70. ! Xiu and Pleim, 2001: Development of a land surface model. Part I: Application
    71. 

  71. !                      in a mesoscale meteorological model. J. Appl. Meteoro.,
    72. 

  72. !                      Vol. 40, 192-209.
    73. 

  73. ! Pleim and Xiu, 2003: Development of a land surface model. Part II: Data
    74. 

  74. !                      assimilation. J. Appl. Meteoro., Vol. 42, 1811-1822.
    75. 

  75. ! 
    76. 

  76. ! Pleim and Gilliam, 2009: An Indirect Data Assimilation Scheme for Deep Soil Temperature in the 
    77. 

  77. !                          Pleim–Xiu Land Surface Model. J. Appl. Meteor. Climatol., 48, 1362–1376.
    78. 

  78. !
    79. 

  79. ! Gilliam and Pleim, 2010: Performance assessment of new land-surface and planetary boundary layer 
    80. 

  80. !                          physics in the WRF-ARW. Journal of Applied Meteorology and Climatology, 49, 760-774. 
    81. 

  81. ! REVISION HISTORY:
    82. 

  82. !     AX     4/2005  - developed the initial WRF version based on the MM5 PX LSM
    83. 

  83. !     RG     2/2008  - Completed testing of the intial working version of PX LSM, released in WRFV3.0 early 2008
    84. 

  84. !     DW     8/2011  - Landuse specific versions of PX (USGS/NLCD/MODIS) were unified into a single code with
    85. 

  85. !                      landuse characteristics defined in module_sf_pxsfclay.F.
    86. 

  86. !     RG     12/2011 - Basic code clean, removed commented out debug statements, lined up columns, etc.
    87. 

  87. !     RG     01/2012 - Removed FIRSTIME Logic that computes PX Landuse characteristics at first time step only. This resulted 
    88. 

  88. !                      in different solutions when OpenMP was used and would not work with moving domains.
    89. 

  89. !                      
    90. 

  90. !--------------------------------------------------------------------------------------------------------------
    91. 

  91. !--------------------------------------------------------------------------------------------------------------
    92. 

  92. !   ARGUMENT LIST:
    93. 

  93. !
    94. 

  94. !... Inputs:
    95. 

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

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

  97. !-- DZ8W            dz between full levels (m)
    98. 

  98. !-- QV3D	    3D mixing ratio
    99. 

  99. !-- T3D             Temperature (K)
    100. 

  100. !-- TH3D            Theta (K)
    101. 

  101. !-- RHO             3D dry air density (kg/m^3)
    102. 

  102. 

  103. !-- PSFC		surface pressure (Pa)
    104. 

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

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

  106. !-- RAINBL		Timestep rainfall 
    107. 

  107. !-- EMISS		surface emissivity (between 0 and 1)
    108. 

  108. 

  109. !-- ITIMESTEP		time step number 
    110. 

  110. !-- NSOIL	        number of soil layers
    111. 

  111. !-- DT			time step (second)
    112. 

  112. !-- ANAL_INTERVAL       Interval of analyses used for soil moisture and temperature nudging
    113. 

  113. 

  114. !-- XLAND       	land mask (1 for land, 2 for water)
    115. 

  115. !-- XICE          Sea ice
    116. 

  116. !-- ALBBCK              Background Albedo
    117. 

  117. !-- ALBEDO      	surface albedo with snow cover effects
    118. 

  118. !-- SNOALB              Albedo of snow
    119. 

  119. 

  120. !-- SMOIS       	total soil moisture content (volumetric fraction)
    121. 

  121. !-- TSLB        	soil temp (k)
    122. 

  122. !-- MAVAIL              Moisture availibility of soil
    123. 

  123. !-- TA2                 2-m temperature
    124. 

  124. !-- QA2                 2-m mixing ratio
    125. 

  125. 

  126. !-- SVPT0       constant for saturation vapor pressure (K)
    127. 

  127. !-- SVP1        constant for saturation vapor pressure (kPa)
    128. 

  128. !-- SVP2        constant for saturation vapor pressure (dimensionless)
    129. 

  129. !-- SVP3        constant for saturation vapor pressure (K)
    130. 

  130. 

  131. !-- ZS          depths of centers of soil layers
    132. 

  132. !-- DZS         thicknesses of soil layers
    133. 

  133. !-- PSIH        similarity stability function for heat
    134. 

  134. 

  135. !-- LANDUSEF     Landuse fraction		
    136. 

  136. !-- SOILCTOP     Top soil fraction		
    137. 

  137. !-- SOILCBOT     Bottom soil fraction		
    138. 

  138. !-- VEGFRA       Vegetation fraction                   
    139. 

  139. !-- VEGF_PX      Veg fraction recomputed and used by PX LSM
    140. 

  140. !-- ISLTYP       Soil type
    141. 

  141. 

  142. !-- RA           Aerodynamic resistence			
    143. 

  143. !-- RS           Stomatal resistence			
    144. 

  144. !-- LAI          Leaf area index (weighted according to fractional landuse)
    145. 

  145. !-- NLCAT        Number of landuse categories		
    146. 

  146. !-- NSCAT        Number of soil categories		
    147. 

  147. 

  148. !-- HFX           net upward heat flux at the surface (W/m^2)
    149. 

  149. !-- QFX           net upward moisture flux at the surface (kg/m^2/s)
    150. 

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

  151. !-- TSK 	  surface skin temperature (K)
    152. 

  152. !-- SST     sea surface temperature
    153. 

  153. !-- ZNT           rougness length
    154. 

  154. !-- CANWAT        Canopy water (mm)
    155. 

  155. 

  156. !-- GRDFLX        Ground heat flux
    157. 

  157. !-- SFCEVP      Evaportation from surface
    158. 

  158. !-- SHDMIN       Minimum annual vegetation fraction for each grid cell
    159. 

  159. !-- SHDMAX       Maximum annual vegetation fraction for each grid cell
    160. 

  160. 

  161. !-- SNOWC         flag indicating snow coverage (1 for snow cover)
    162. 

  162. !-- PBLH          PBL height (m)
    163. 

  163. !-- RMOL          1/L Reciprocal of Monin-Obukhov length
    164. 

  164. !-- UST           u* in similarity theory (m/s)
    165. 

  165. !-- CAPG          heat capacity for soil (J/K/m^3)
    166. 

  166. !-- DTBL	  time step of boundary layer calls
    167. 

  167. 

  168. !-- T2_NDG_OLD    Analysis temperature prior to current time
    169. 

  169. !-- T2_NDG_NEW    Analysis temperature ahead of current time
    170. 

  170. !-- Q2_NDG_OLD    Analysis mixing ratio prior to current time
    171. 

  171. !-- Q2_NDG_NEW    Analysis mixing ratio ahead of current time
    172. 

  172. 

  173. !-- SN_NDG_OLD    Analysis snow water prior to current time
    174. 

  174. !-- SN_NDG_NEW    Analysis snow water ahead of current time
    175. 

  175. !-- SNOW          Snow water equivalent
    176. 

  176. !-- SNOWH         Physical snow depth
    177. 

  177. !-- SNOWNCV       Time step accumulated snow
    178. 

  178. 

  179. !-- T2OBS         Analysis temperature interpolated from prior and next in time analysese
    180. 

  180. !-- Q2OBS         Analysis moisture interpolated from prior and next in time analysese
    181. 

  181. !-- PXLSM_SMOIS_INIT  Flag to intialize deep soil moisture to a value derived from moisture availiability.
    182. 

  182. !-- PXLSM_SOIL_NUDGE  Flag to use soil moisture and temperature nudging in the PX LSM
    183. 

  183. !                     This is typically done for the first simulation.
    184. 

  184. 

  185. !-- ids             start index for i in domain
    186. 

  186. !-- ide             end index for i in domain
    187. 

  187. !-- jds             start index for j in domain
    188. 

  188. !-- jde             end index for j in domain
    189. 

  189. !-- kds             start index for k in domain
    190. 

  190. !-- kde             end index for k in domain
    191. 

  191. !-- ims             start index for i in memory
    192. 

  192. !-- ime             end index for i in memory
    193. 

  193. !-- jms             start index for j in memory
    194. 

  194. !-- jme             end index for j in memory
    195. 

  195. !-- kms             start index for k in memory
    196. 

  196. !-- kme             end index for k in memory
    197. 

  197. !-- jts             start index for j in tile
    198. 

  198. !-- jte             end index for j in tile
    199. 

  199. !-- kts             start index for k in tile
    200. 

  200. !-- kte             end index for k in tile
    201. 

  201. !... Outputs:
    202. 

  202. 

  203.      IMPLICIT NONE
    204. 

  204. 

  205. !.......Arguments
    206. 

  206. ! DECLARATIONS - INTEGER
    207. 

  207.     INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
    208. 

  208.                                    ims,ime, jms,jme, kms,kme, &
    209. 

  209.                                    its,ite, jts,jte, kts,kte 
    210. 

  210. 

  211.    INTEGER,   INTENT(IN)  ::      NSOIL, ITIMESTEP, NLCAT, NSCAT,             &
    212. 

  212.                                   ANAL_INTERVAL, PXLSM_SMOIS_INIT, PXLSM_SOIL_NUDGE
    213. 

  213. 

  214.    REAL,     INTENT(IN )  ::      DT,DTBL 
    215. 

  215. 

  216.    INTEGER,   DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::   ISLTYP
    217. 

  217. 

  218. ! DECLARATIONS - REAL
    219. 

  219.     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),  INTENT(IN) :: U3D, V3D, RHO, &
    220. 

  220.                                                                     T3D, TH3D, DZ8W, QV3D           
    221. 

  221. 

  222.    REAL,     DIMENSION(1:NSOIL), INTENT(IN)::ZS,DZS
    223. 

  223.    REAL,     DIMENSION( ims:ime , 1:NSOIL, jms:jme ),  INTENT(INOUT)   ::  SMOIS, TSLB     
    224. 

  224. 

  225.    REAL,     DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: RA, RS, LAI, ZNT
    226. 

  226.    REAL,     DIMENSION( ims:ime , jms:jme ), INTENT(OUT):: GRDFLX, TSK, TA2, QA2
    227. 

  227. 

  228.    REAL,     DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN):: LANDUSEF       
    229. 

  229.    REAL,     DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN):: SOILCTOP, SOILCBOT
    230. 

  230. 

  231. 

  232.    REAL,    DIMENSION( ims:ime, jms:jme ),                                     &
    233. 

  233.             INTENT(IN) ::                            PSFC, GSW, GLW, RAINBL,   &                
    234. 

  234.                                                      EMISS, SNOALB,            &
    235. 

  235.                                                      ALBBCK, SHDMIN, SHDMAX,   &
    236. 

  236.                                                      PBLH, RMOL, SNOWNCV,      &
    237. 

  237.                                                      UST, MAVAIL, SST
    238. 

  238.    REAL,    DIMENSION( ims:ime, jms:jme ),                                     &
    239. 

  239.             INTENT(IN) ::                            T2_NDG_OLD, T2_NDG_NEW,   &                
    240. 

  240.                                                      Q2_NDG_OLD, Q2_NDG_NEW,   & 
    241. 

  241.                                                      SN_NDG_OLD, SN_NDG_NEW   
    242. 

  242. 

  243.    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: T2OBS, Q2OBS                
    244. 

  244.                                                      
    245. 

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

  246.             INTENT(INOUT) ::           CAPG,CANWAT, QFX, HFX, LH,              &
    247. 

  247.                                        PSIH,VEGFRA, VEGF_PX, SNOW,             &
    248. 

  248.                                        SNOWH, SNOWC, ALBEDO, XLAND, XICE
    249. 

  249. 

  250. 

  251.    LOGICAL :: radiation
    252. 

  252. 

  253. !-------------------------------------------------------------------------
    254. 

  254. !     ---------- Local Variables --------------------------------
    255. 

  255. 

  256. 

  257.       !----  PARAMETERS
    258. 

  258.       INTEGER, PARAMETER  :: NSTPS  = 11     ! max. soil types
    259. 

  259.       REAL, PARAMETER     :: DTPBLX = 40.0   ! Max PX timestep = 40 sec
    260. 

  260. 

  261.       
    262. 

  262.       !----  INTEGERS
    263. 

  263.       INTEGER, DIMENSION( 1: NSTPS ) :: JP
    264. 

  264.       INTEGER:: J, I, NS, NUDGE, ISTI, WEIGHT
    265. 

  265.       INTEGER:: NTSPS, IT
    266. 

  266. 

  267.       !----  REALS
    268. 

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

  269.                                                  XVEG, XVEGMN, XSNUP, &
    270. 

  270.                                                  XALB
    271. 

  271.                                                   
    272. 

  272.       REAL,     DIMENSION( ims:ime, jms:jme ) :: RADNET, EG, ER, ETR, QST
    273. 

  273. 

  274.       REAL:: SFCPRS,TA1,DENS1,QV1,ZLVL,SOLDN,LWDN,    &
    275. 

  275.             EMISSI,PRECIP,THETA1,VAPPRS,QSBT,         &
    276. 

  276.             WG,W2,WR,TG,T2,USTAR,MOLX,Z0,             &
    277. 

  277.             RAIR,CPAIR,IFLAND,ISNOW,                  &
    278. 

  278.             ES,QSS,BETAP,                             &
    279. 

  279.             RH2_OLD, RH2_NEW, T2_OLD, T2_NEW,         &
    280. 

  280.             CORE, CORB, TIME_BETWEEN_ANALYSIS,        &
    281. 

  281.             G1000, ALN10,RH2OBS, HU, SNOBS,           &
    282. 

  282.             FWSAT,FWFC,FWWLT,FB,FCGSAT,FJP,FAS,       &
    283. 

  283.             FSEAS, T2I, HC_SNOW, SNOW_FRA,SNOWALB,    &
    284. 

  284.             QST12,ZFUNC,ZF1,ZA2,QV2, DT_FDDA,CURTIME, &
    285. 

  285.             FC2R,FC1SAT, DTPBL
    286. 

  286. 

  287.       CHARACTER (LEN = 5) :: LAND_USE_TYPE
    288. 

  288. 

  289. !-------------------------------------------------------------------------
    290. 

  290. !-------------------------------Executable starts here--------------------
    291. 

  291. !
    292. 

  292.       ! Determine Landuse Dataset by the number of categories
    293. 

  293.       IF (NLCAT == 50) THEN
    294. 

  294.          LAND_USE_TYPE = 'NLCD'
    295. 

  295.       ELSE IF (NLCAT == 20) THEN
    296. 

  296.          LAND_USE_TYPE = 'MODIS'
    297. 

  297.       ELSE IF (NLCAT == 24) THEN
    298. 

  298.          LAND_USE_TYPE = 'USGS'
    299. 

  299.       ELSE
    300. 

  300.          PRINT *, 'Error: Unknown Land Use Category'
    301. 

  301.          STOP
    302. 

  302.       END IF 
    303. 

  303. 

  304.       ALN10  = ALOG(10.0)
    305. 

  305.       G1000  = g*1.0E-3            ! G/1000
    306. 

  306.       WEIGHT = 0
    307. 

  307.       DT_FDDA = ANAL_INTERVAL * 1.0    ! Convert DT of Analysis to real 
    308. 

  308. 

  309.       !-----------------------------------------------------------------------------------
    310. 

  310.       ! Compute vegetation and land-use characteristics by land-use fraction weighting
    311. 

  311.       ! These parameters include LAI, VEGF, ZNT, ALBEDO, RS, etc.
    312. 

  312.       CALL VEGELAND(LANDUSEF, VEGFRA, SHDMIN, SHDMAX,         &
    313. 

  313.                     SOILCTOP, SOILCBOT, NLCAT, NSCAT,         &
    314. 

  314.                     ZNT,XLAI,XLAIMN,RSTMIN,XVEG,XVEGMN,XSNUP,	&
    315. 

  315.                     XLAND, XALB,                              &
    316. 

  316.                     ids,ide, jds,jde, kds,kde,			          &
    317. 

  317.                     ims,ime, jms,jme, kms,kme,			          &
    318. 

  318.                     its,ite, jts,jte, kts,kte, LAND_USE_TYPE   )
    319. 

  319.       !-----------------------------------------------------------------------------------
    320. 

  320.       !--- Compute time relatve to old and new analysis time for timestep interpolation
    321. 

  321.       CURTIME = (ITIMESTEP-1) * DT
    322. 

  322.       TIME_BETWEEN_ANALYSIS = MOD(CURTIME,DT_FDDA)
    323. 

  323.       IF (TIME_BETWEEN_ANALYSIS .EQ. 0.0) THEN   
    324. 

  324.         CORB = 1.0       
    325. 

  325.         CORE = 0.0    
    326. 

  326.       ELSE    
    327. 

  327.         CORE = TIME_BETWEEN_ANALYSIS  / DT_FDDA     
    328. 

  328.         CORB = 1.0 - CORE  
    329. 

  329.       ENDIF 
    330. 

  330.       !-----------------------------------------------------------------------------------
    331. 

  331. 

  332.       !-----------------------------------------------------------------------------------
    333. 

  333.       ! Main loop over individual grid cells
    334. 

  334.       DO J = jts,jte    !-- J LOOP
    335. 

  335.        DO I = its,ite   !-- I LOOP
    336. 

  336.           
    337. 

  337.           IFLAND = XLAND(I,J)
    338. 

  338. 

  339.           ! Compute soil properties via weighting of fractional components
    340. 

  340.           IF (IFLAND .LT. 1.5 ) THEN 
    341. 

  341. 

  342.           !---------------------------------------------------------
    343. 

  343.             CALL SOILPROP (SOILCBOT(I,:,J), WEIGHT,               &
    344. 

  344.                          ITIMESTEP, MAVAIL(I,J),                  &
    345. 

  345.                          PXLSM_SMOIS_INIT,                        &
    346. 

  346.                          FWSAT,FWFC,FWWLT,FB,FCGSAT,              &  
    347. 

  347.                          FJP,FAS,FC2R,FC1SAT,ISTI,SMOIS(I,2,J)    ) 
    348. 

  348.           !----------------------------------------------------------
    349. 

  349.             ISLTYP(I,J) = ISTI                      
    350. 

  350.           ELSE
    351. 

  351.             ISLTYP(I,J) = 14                    ! STATSGO type for water
    352. 

  352.           ENDIF
    353. 

  353. 

  354.           !--  Variables Sub. SURFPX needs
    355. 

  355.           SFCPRS = PSFC(i,j) / 1000.0                 ! surface pressure in cb
    356. 

  356.           TA1    = T3D(i,1,j)                         ! air temperature at first layer
    357. 

  357.           DENS1  = RHO(I,1,J)                         ! air density at first layer
    358. 

  358.           QV1    = QV3D(i,1,j)                        ! water vapor mixing ratio at first layer
    359. 

  359.           QV2    = QV3D(i,2,j)
    360. 

  360.           ZLVL   = 0.5 * DZ8W(i,1,j)                  ! thickness of lowest half level
    361. 

  361.           ZF1    = DZ8W(i,1,j)
    362. 

  362.           ZA2    = ZF1 + 0.5 * DZ8W(i,2,j)
    363. 

  363. 

  364.           LWDN   = GLW(I,J)                           ! longwave radiation
    365. 

  365.           EMISSI = EMISS(I,J)                         ! emissivity
    366. 

  366.           PRECIP = MAX ( 1.0E-3*RAINBL(i,j)/DTBL,0.0) ! accumulated precip. rate during DT (=dtpbl)
    367. 

  367.                                                       ! convert RAINBL from mm to m for PXLSM
    368. 

  368.           WR     = 1.0E-3*CANWAT(I,J)                 ! convert CANWAT from mm to m for PXLSM
    369. 

  369.           THETA1 = TH3D(i,1,j)			                  ! potential temp at first layer
    370. 

  370.           SNOBS  = SNOW(I,J)                          ! Set snow cover to existing model value
    371. 

  371.                                                       !   this is overwritten below if snow analysis is availiable
    372. 

  372.                                                       !   otherwise snow cover remains constant through simulation
    373. 

  373.                    
    374. 

  374.           IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN                   
    375. 

  375.             !-- 2 m Temp and RH for Nudging
    376. 

  376.             T2_OLD = T2_NDG_OLD(I,J)
    377. 

  377.             T2_NEW     = T2_NDG_NEW(I,J)
    378. 

  378.             VAPPRS     = SVP1 * EXP(SVP2 * (T2_OLD - SVPT0) / ( T2_OLD - SVP3))
    379. 

  379.             QSBT       = EP_2 * VAPPRS / (SFCPRS - VAPPRS)          
    380. 

  380.             RH2_OLD    = Q2_NDG_OLD(I,J) / QSBT
    381. 

  381.             VAPPRS     = SVP1 * EXP(SVP2 * (T2_NEW - SVPT0) / (T2_NEW - SVP3))
    382. 

  382.             QSBT       = EP_2 * VAPPRS / (SFCPRS - VAPPRS)          
    383. 

  383.             RH2_NEW    = Q2_NDG_NEW(I,J) / QSBT
    384. 

  384.             RH2OBS     = CORB * RH2_OLD +  CORE * RH2_NEW  
    385. 

  385.             T2OBS(I,J) = CORB * T2_OLD +  CORE * T2_NEW
    386. 

  386.             Q2OBS(I,J) = CORB * Q2_NDG_OLD(I,J) +  CORE * Q2_NDG_NEW(I,J)
    387. 

  387.             SNOBS = CORB * SN_NDG_OLD(I,J) +  CORE * SN_NDG_NEW(I,J)  
    388. 

  388.           ENDIF
    389. 

  389. 

  390.           USTAR  = MAX(UST(I,J),0.005)
    391. 

  391.           
    392. 

  392.           IF (IFLAND .GT. 1.5) THEN ! if over water                            
    393. 

  393.             ZNT(I,J) = CZO * UST(I,J) * UST(I,J) / G + OZO
    394. 

  394.           ENDIF                                                              
    395. 

  395.           
    396. 

  396.           Z0     = ZNT(I,J)
    397. 

  397.           CPAIR  = CPD * (1.0 + 0.84 * QV1)            ! J/(K KG)
    398. 

  398. 

  399.           !-------------------------------------------------------------
    400. 

  400.           ! Compute fractional snow area and snow albedo
    401. 

  401.           CALL PXSNOW (ITIMESTEP, SNOBS, SNOWNCV(I,J), SNOW(I,J),  &
    402. 

  402.                        SNOWH(I,J), XSNUP(I,J),  XALB(i,j),         &
    403. 

  403.                        SNOALB(I,J),VEGF_PX(I,J), SHDMIN(I,J),      &
    404. 

  404.                        HC_SNOW, SNOW_FRA, SNOWC(I,J),  ALBEDO(I,J) ) 
    405. 

  405.           !-------------------------------------------------------------
    406. 

  406.                                  
    407. 

  407.           !-------------------------------------------------------------
    408. 

  408.           ! Sea Ice from analysis and water cells that are very cold, but more than 50% water
    409. 

  409.           ! are converted to ice/snow for more reasonable treatment.
    410. 

  410.           IF( (XICE(I,J).GE.0.5) .OR.   &
    411. 

  411.               (SST(I,J).LE.270.0.AND.XLAND(I,J).GT.1.50) ) THEN
    412. 

  412.               XLAND(I,J) = 1.0
    413. 

  413.               IFLAND = 1.0
    414. 

  414.               ZNT(I,J) = 0.001  !  Ice
    415. 

  415.               SMOIS(I,1,J) = 1.0     ! FWSAT
    416. 

  416.               SMOIS(I,2,J) = 1.0     ! FWSAT
    417. 

  417.               XICE(I,J) = 1.0
    418. 

  418.               ALBEDO(I,J) = 0.7
    419. 

  419.               SNOWC(I,J) = 1.0
    420. 

  420.               SNOW_FRA = 1.0
    421. 

  421.               VEGF_PX(I,J) = 0.0
    422. 

  422.               LAI(I,J) = 0.0
    423. 

  423.           ENDIF
    424. 

  424.           !-------------------------------------------------------------
    425. 

  425. 

  426.           !-------------------------------------------------------------
    427. 

  427.           !-- Note that when IFGROW = 0 is selected in Vegeland then max and min           
    428. 

  428.           !-- LAI and Veg are the same                                                     
    429. 

  429.           T2I = TSLB(I,2,J)                                            
    430. 

  430. !         FSEAS = AMAX1(1.0 - 0.0016 * (298.0 - T2I) ** 2,0.0)  ! BATS            
    431. 

  431.           FSEAS = AMAX1(1.0 - 0.015625 * (290.0 - T2I) ** 2,0.0) ! JP97            
    432. 

  432.           IF (T2I .GE. 290.0) FSEAS = 1.0                                          
    433. 

  433.           
    434. 

  434.           LAI(I,J)       = XLAIMN(I,J) + FSEAS*(XLAI(I,J) - XLAIMN(I,J))                 
    435. 

  435.           VEGF_PX(I,J)   = XVEGMN(I,J) + FSEAS*(XVEG(I,J) - XVEGMN(I,J))                       
    436. 

  436.           
    437. 

  437.           ! Ensure veg algorithms not used for water 
    438. 

  438.           IF (IFLAND .GT. 1.5) THEN                        
    439. 

  439.              VEGF_PX(I,J) = 0.0         
    440. 

  440.           ENDIF                                                                   
    441. 

  441.           !-------------------------------------------------------------
    442. 

  442. 

  443. 

  444.           SOLDN  = GSW(I,J) / (1.0 - ALBEDO(I,J))     ! downward shortwave radiaton
    445. 

  445.           ISNOW = SNOWC(I,J)
    446. 

  446. 

  447. 

  448.           NUDGE=PXLSM_SOIL_NUDGE
    449. 

  449.           IF ( J .LE. 2 .OR. J .GE. (jde-1) ) NUDGE=0
    450. 

  450.           IF ( I .LE. 2 .OR. I .GE. (ide-1) ) NUDGE=0
    451. 

  451.           
    452. 

  452.           IF ( RMOL(I,J) .GT. 0.0 )  THEN
    453. 

  453.               MOLX = AMIN1(1/RMOL(I,J),1000.0)
    454. 

  454.           ELSE IF ( RMOL(I,J) .LT. 0.0 ) THEN
    455. 

  455.               MOLX = AMAX1(1/RMOL(I,J),-1000.0)
    456. 

  456.           ELSE
    457. 

  457.               MOLX = 1000
    458. 

  458.           ENDIF   
    459. 

  459.  
    460. 

  460.           ZFUNC = ZF1 * (1.0 - ZF1 / MAX(100.,PBLH(I,J))) ** 2
    461. 

  461.           QST12 = KARMAN * ZFUNC*(QV2-QV1) / (ZA2-ZLVL)
    462. 

  462.                
    463. 

  463. 

  464. !--- LSM sub-time loop too prevent dt > 40 sec       
    465. 

  465.           NTSPS = INT(DT / (DTPBLX + 0.000001) + 1.0)                                
    466. 

  466.           DTPBL = DT / NTSPS                                                       
    467. 

  467. 

  468.           DO IT=1,NTSPS                                                     
    469. 

  469.           
    470. 

  470.             !... SATURATION VAPOR PRESSURE (MB) 
    471. 

  471.             IF ( TSLB(I,1,J) .LE. SVPT0 ) THEN        ! For ground that is below freezing
    472. 

  472.               ES = SVP1 * EXP(22.514 - 6.15E3 / TSLB(I,1,J))      ! cb
    473. 

  473.             ELSE
    474. 

  474.               ES = SVP1 * EXP(SVP2 * (TSLB(I,1,J) - SVPT0) /  (TSLB(I,1,J) - SVP3))
    475. 

  475.             ENDIF
    476. 

  476.             QSS  = ES * 0.622 / (SFCPRS - ES)
    477. 

  477.           
    478. 

  478.             !... beta method, Lee & Pielke (JAM,May1992)
    479. 

  479.             BETAP = 1.0
    480. 

  480.             IF (IFLAND .LT. 1.5 .AND. ISNOW .LT. 0.5 .AND. SMOIS(I,1,J) .LE. FWFC) THEN       
    481. 

  481.               BETAP = 0.25 * (1.0 - COS(SMOIS(I,1,J) / FWFC * PI)) ** 2
    482. 

  482.             ENDIF
    483. 

  483.             
    484. 

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

  485.             CALL SURFPX (DTPBL, IFLAND, SNOWC(I,J),  NUDGE, XICE(I,J),          & !in
    486. 

  486.                       SOLDN,  GSW(I,J), LWDN,   EMISSI, ZLVL,                   & !in
    487. 

  487.                       MOLX,    Z0,    USTAR,                                    & !in
    488. 

  488.                       SFCPRS, DENS1,  QV1,    QSS,   TA1,                       & !in
    489. 

  489.                       THETA1,   PRECIP,                                         & !in
    490. 

  490.                       CPAIR, PSIH(I,J),                                         & !in
    491. 

  491.                       RH2OBS,T2OBS(I,J),                                        & !in
    492. 

  492.                       VEGF_PX(I,J), ISTI,   LAI(I,J),   BETAP,                  & !in
    493. 

  493.                       RSTMIN(I,J), HC_SNOW, SNOW_FRA,                           & !in          
    494. 

  494.                       FWWLT, FWFC, FCGSAT,  FWSAT, FB,                          & !in
    495. 

  495.                       FC1SAT,FC2R,FAS,FJP,DZS(1),DZS(2),QST12,                  & !in
    496. 

  496.                       RADNET(I,J), GRDFLX(I,J), HFX(I,J), QFX(I,J), LH(I,J),    & !out
    497. 

  497.                       EG(I,J), ER(I,J), ETR(I,J),                               & !out
    498. 

  498.                       QST(I,J), CAPG(I,J), RS(I,J), RA(I,J),                    & !out
    499. 

  499.                       TSLB(I,1,J), TSLB(I,2,J),                                 & !out
    500. 

  500.                       SMOIS(I,1,J), SMOIS(I,2,J), WR, TA2(I,J), QA2(I,J),       &
    501. 

  501.                       LAND_USE_TYPE )
    502. 

  502.             !-------------------------------------------------------------------------
    503. 

  503.           
    504. 

  504.           END DO                        ! Time internal PX time loop   
    505. 

  505. 

  506.           TSK(I,J)    = TSLB(I,1,J)     ! Skin temp set to 1 cm soil temperature in PX for now
    507. 

  507.           CANWAT(I,J) = WR * 1000.      ! convert WR back to mm for CANWAT
    508. 

  508.           
    509. 

  509.        ENDDO !  END MIAN I LOOP
    510. 

  510.       ENDDO  !  END MAIN J LOOP
    511. 

  511.       
    512. 

  512. !------------------------------------------------------
    513. 

  513.    END SUBROUTINE pxlsm
    514. 

  514. !-------------------------------------------------------------------------
    515. 

  515. !-------------------------------------------------------------------------
    516. 

  516. 

  517. 

  518. !-------------------------------------------------------------------------
    519. 

  519. !-------------------------------------------------------------------------
    520. 

  520.       SUBROUTINE VEGELAND( landusef, vegfra,    			                  & 
    521. 

  521.                             shdmin, shdmax,                             &                          
    522. 

  522.                             soilctop, soilcbot, nlcat, nscat,           &
    523. 

  523.                             znt,  xlai, xlaimn, rstmin, xveg, xvegmn,	  &
    524. 

  524.                             xsnup, xland, xalb,                         &
    525. 

  525.                             ids,ide, jds,jde, kds,kde,			            &
    526. 

  526.                             ims,ime, jms,jme, kms,kme,			            &
    527. 

  527.                             its,ite, jts,jte, kts,kte,                  &
    528. 

  528.                             LAND_USE_TYPE                                 )
    529. 

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

  530. !           
    531. 

  531. !   CALLED FROM Sub. bl_init in module_physics.init.F
    532. 

  532. !           
    533. 

  533. !   THIS PROGRAM PROCESSES THE USGS LANDUSE DATA
    534. 

  534. !   WHICH HAS BEEN GRIDDED BY THE WPS SYSTEM
    535. 

  535. !   AND PRODUCES 2-D FIELDS OF LU RELATED PARAMETERS
    536. 

  536. !   FOR USE IN THE PX SURFACE MODEL
    537. 

  537. !
    538. 

  538. !
    539. 

  539. ! REVISION HISTORY:
    540. 

  540. !  AX      Oct 2004 - developed WRF version based on VEGELAND in the MM5 PX LSM
    541. 

  541. !  RG      Dec 2006 - revised for WRFV2.1.2
    542. 

  542. !  JP      Dec 2007 - revised for WRFV3 - removed IFGROW options
    543. 

  543. !-------------------------------------------------------------------------
    544. 

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

  545.                               
    546. 

  546.       IMPLICIT NONE
    547. 

  547. !... 
    548. 

  548.       INTEGER,  INTENT(IN)        ::  ids,ide, jds,jde, kds,kde,  &
    549. 

  549.                                       ims,ime, jms,jme, kms,kme,  &
    550. 

  550.                                       its,ite, jts,jte, kts,kte
    551. 

  551.                                                                             
    552. 

  552.       INTEGER , INTENT(IN)        ::  NSCAT, NLCAT
    553. 

  553. 

  554.       REAL,    DIMENSION( ims:ime , 1:NLCAT, jms:jme ),  INTENT(IN) :: LANDUSEF                        
    555. 

  555.       REAL,    DIMENSION( ims:ime , 1:NSCAT, jms:jme ),  INTENT(IN) :: SOILCTOP, SOILCBOT
    556. 

  556.                                                                    
    557. 

  557.       REAL,    DIMENSION( ims:ime, jms:jme ),            INTENT(IN) ::  VEGFRA, SHDMIN, SHDMAX 
    558. 

  558. 

  559.       REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ZNT
    560. 

  560. 

  561.       REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: XLAI, XLAIMN, RSTMIN, XALB, &
    562. 

  562.                                                               XVEG, XVEGMN, XSNUP, XLAND 
    563. 

  563.       CHARACTER (LEN = 5), INTENT(IN) :: LAND_USE_TYPE
    564. 

  564. 

  565.                                      
    566. 

  566. !... local variables
    567. 

  567. 

  568.       INTEGER :: itf, jtf, k, j, i
    569. 

  569.       REAL    :: SUMLAI, SUMLMN, SUMRSI, SUMLZ0, SUMVEG, SUMVMN, ALAI, VEGF, SUMSNUP
    570. 

  570.       REAL    :: VFMX, VFMN, VSEAS, FAREA, FWAT, ZNOTC, SUMALB
    571. 

  571.       REAL,   DIMENSION( NLCAT )  :: LAIMX, LAIMN, Z0, VEG, VEGMN, SNUP
    572. 

  572.       REAL, PARAMETER      :: ZNOTCMN = 5.0  ! CM, MIN Zo FOR CROPS
    573. 

  573.       REAL, PARAMETER      :: ZNOTCMX = 15.0 ! CM, MAX Zo FOR CROPS
    574. 

  574. 

  575.       REAL, SAVE, DIMENSION(:), POINTER  :: RSMIN, Z00, VEG0, VEGMN0, LAI0, LAIMN0, SNUP0, ALBF
    576. 

  576. 

  577.       !---- INITIALIZE PARAMETERS
    578. 

  578. 

  579.       INTEGER, SAVE :: KWAT
    580. 

  580.       INTEGER, SAVE :: LIMIT1, LIMIT2
    581. 

  581. 

  582.      ! Initialize LU characteristics by LU Dataset
    583. 

  583.       IF (LAND_USE_TYPE == 'USGS') THEN
    584. 

  584.          KWAT = 16
    585. 

  585.          RSMIN  => RSMIN_USGS
    586. 

  586.          Z00    => Z00_USGS  
    587. 

  587.          VEG0   => VEG0_USGS 
    588. 

  588.          VEGMN0 => VEGMN0_USGS 
    589. 

  589.          LAI0   => LAI0_USGS 
    590. 

  590.          LAIMN0 => LAIMN0_USGS 
    591. 

  591.          SNUP0  => SNUP0_USGS 
    592. 

  592.          ALBF   => ALBF_USGS
    593. 

  593.          LIMIT1 = 2
    594. 

  594.          LIMIT1 = 6
    595. 

  595.       ELSE IF (LAND_USE_TYPE == 'NLCD') THEN
    596. 

  596.          KWAT = 1
    597. 

  597.          RSMIN  => RSMIN_NLCD
    598. 

  598.          Z00    => Z00_NLCD  
    599. 

  599.          VEG0   => VEG0_NLCD 
    600. 

  600.          VEGMN0 => VEGMN0_NLCD 
    601. 

  601.          LAI0   => LAI0_NLCD 
    602. 

  602.          LAIMN0 => LAIMN0_NLCD 
    603. 

  603.          SNUP0  => SNUP0_NLCD 
    604. 

  604.          ALBF   => ALBF_NLCD
    605. 

  605.          LIMIT1 = 20
    606. 

  606.          LIMIT1 = 43
    607. 

  607.       ELSE IF (LAND_USE_TYPE == 'MODIS') THEN
    608. 

  608.          KWAT = 17
    609. 

  609.          RSMIN  => RSMIN_MODIS
    610. 

  610.          Z00    => Z00_MODIS  
    611. 

  611.          VEG0   => VEG0_MODIS 
    612. 

  612.          VEGMN0 => VEGMN0_MODIS 
    613. 

  613.          LAI0   => LAI0_MODIS 
    614. 

  614.          LAIMN0 => LAIMN0_MODIS 
    615. 

  615.          SNUP0  => SNUP0_MODIS 
    616. 

  616.          ALBF   => ALBF_MODIS
    617. 

  617.          LIMIT1 = 12
    618. 

  618.          LIMIT1 = 14
    619. 

  619.       END IF
    620. 

  620. 

  621.       !--------------------------------------------------------------------
    622. 

  622.       DO J = jts,jte
    623. 

  623.         DO I = its,ite
    624. 

  624.           XLAI(I,J)   = 0.0                                                        
    625. 

  625.           XLAIMN(I,J) = 0.0                                                     
    626. 

  626.           RSTMIN(I,J) = 9999.0                                                    
    627. 

  627.           XVEG(I,J)   = 0.0                                                       
    628. 

  628.           XVEGMN(I,J) = 0.0
    629. 

  629.           XSNUP(I,J)  = 0.0                                                            
    630. 

  630.           XALB(I,J)   = 0.0                                                            
    631. 

  631.         ENDDO		! END LOOP THROUGH GRID CELLS
    632. 

  632.       ENDDO		! END LOOP THROUGH GRID CELLS
    633. 

  633.       !--------------------------------------------------------------------
    634. 

  634. 

  635.       DO J = jts,jte
    636. 

  636.         DO I = its,ite
    637. 

  637.       
    638. 

  638. 

  639.           !-- Initialize 2 and 3-D veg parameters to be caculated
    640. 

  640.           DO K=1,NLCAT                                                               
    641. 

  641.             LAIMX(K) = LAI0(K)                                                       
    642. 

  642.             LAIMN(K) = LAIMN0(K)                                                   
    643. 

  643.             Z0(K)    = Z00(K)                                                   
    644. 

  644.             VEG(K)   = VEG0(K)                                              
    645. 

  645.             VEGMN(K) = VEGMN0(K)
    646. 

  646.             SNUP(K)  = SNUP0(K)                                            
    647. 

  647.           ENDDO                                              
    648. 

  648. 

  649.           !--  INITIALIZE SUMS
    650. 

  650.           SUMLAI = 0.0
    651. 

  651.           SUMLMN = 0.0
    652. 

  652.           SUMRSI = 0.0
    653. 

  653.           SUMLZ0 = 0.0
    654. 

  654.           SUMVEG = 0.0
    655. 

  655.           SUMVMN = 0.0
    656. 

  656.           ALAI   = 0.0
    657. 

  657.           SUMSNUP= 0.0
    658. 

  658.           SUMALB = 0.0
    659. 

  659. 

  660.           !--   ESTIMATE CROP EMERGANCE DATE FROM VEGFRAC
    661. 

  661.           VFMX   = SHDMAX(I,J)
    662. 

  662.           VFMN   = SHDMIN(I,J)
    663. 

  663.           VEGF   = VEGFRA(I,J) 
    664. 

  664.                 
    665. 

  665.           !-- Computations for VEGETATION CELLS ONLY
    666. 

  666.           IF(VFMX.GT.0.0.AND.LANDUSEF(I,KWAT,J).LT.1.00) THEN
    667. 

  667.             VSEAS = VEGF/VFMX
    668. 

  668.             IF(VSEAS.GT.1.0.OR.VSEAS.LT.0.0) THEN
    669. 

  669.               VSEAS = MIN(VSEAS,1.0)
    670. 

  670.               VSEAS = MAX(VSEAS,0.0)
    671. 

  671.             ENDIF
    672. 

  672. 

  673.             ZNOTC = ZNOTCMN * (1-VSEAS) + ZNOTCMX * VSEAS  ! Zo FOR CROPS
    674. 

  674.             DO K = 1, NLCAT
    675. 

  675.               IF (LAND_USE_TYPE == 'MODIS') THEN
    676. 

  676. 

  677.                 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
    678. 

  678.                 IF (K.EQ.12 .OR. K.EQ.14) THEN
    679. 

  679.                    LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
    680. 

  680.                    LAIMN(K) = LAIMX(K)
    681. 

  681.                    VEG(K)   = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
    682. 

  682.                    VEGMN(K) = VEG(K)
    683. 

  683.                    !-- SEASONALLY VARY Zo FOR MODIS DryCrop (k=12)
    684. 

  684.                    IF (K .EQ. 12) THEN
    685. 

  685.                       Z0(K) = ZNOTC
    686. 

  686.                    !-- CrGrM (k=14) USE AVG WITH GRASS AND FOREST
    687. 

  687.                    ELSE IF (K .EQ.14) THEN
    688. 

  688.                       Z0(K)  = 0.5 * (ZNOTC + Z00(K))
    689. 

  689.                    ENDIF
    690. 

  690.                 ENDIF
    691. 

  691.               ELSE IF (LAND_USE_TYPE == 'NLCD') THEN
    692. 

  692.                 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
    693. 

  693.                 IF (K.EQ.20 .OR. K.EQ.43 .OR. K.EQ.45) THEN
    694. 

  694.                    LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
    695. 

  695.                    LAIMN(K) = LAIMX(K)
    696. 

  696.                    VEG(K)   = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
    697. 

  697.                    VEGMN(K) = VEG(K)
    698. 

  698.                    !-- SEASONALLY VARY Zo FOR DryCrop (k=20) OR Irigated Crop (k=43) OR  Mix Crop (k=4)
    699. 

  699.                    IF (K.EQ.20 .OR. K.EQ.43) THEN
    700. 

  700.                       Z0(K) = ZNOTC
    701. 

  701.                    !-- CrNatM (k=45) USE AVG WITH GRASS AND FOREST
    702. 

  702.                   ELSE IF (K.EQ.45) THEN
    703. 

  703.                       Z0(K)  = 0.5 * (ZNOTC + Z00(K))
    704. 

  704.                   ENDIF
    705. 

  705.                 ENDIF
    706. 

  706.               ELSE IF (LAND_USE_TYPE == 'USGS') THEN
    707. 

  707.                 !-- USE THE VEGFRAC DATA ONLY FOR CROPS 
    708. 

  708.                 IF (K .GE. 2 .AND. K .LE. 6) THEN
    709. 

  709.                    LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
    710. 

  710.                    LAIMN(K) = LAIMX(K)
    711. 

  711.                    VEG(K)   = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
    712. 

  712.                    VEGMN(K) = VEG(K)
    713. 

  713.                    !-- SEASONALLY VARY Zo FOR DryCrop (k=2) OR Irigated Crop (k=3) OR  Mix Crop (k=4)
    714. 

  714.                    IF (K .GE. 2 .AND. K .LE. 4) THEN
    715. 

  715.                       Z0(K) = ZNOTC
    716. 

  716.                    !-- CrGrM (k=5) or CrWdM (k=6) USE AVG WITH GRASS AND FOREST
    717. 

  717.                    ELSE IF (K .GE.5 .AND. K .LE. 6) THEN
    718. 

  718.                       Z0(K)  = 0.5 * (ZNOTC + Z00(K))
    719. 

  719.                    ENDIF
    720. 

  720.                 ENDIF
    721. 

  721.               END IF
    722. 

  722. 

  723.             ENDDO
    724. 

  724. 

  725.           ENDIF       !-- IF cell is vegetation
    726. 

  726.           !-------------------------------------
    727. 

  727.           !-- LOOP THROUGH LANDUSE Fraction and compute totals
    728. 

  728.           DO K = 1, NLCAT 
    729. 

  729.             FAREA  = LANDUSEF(I,K,J)
    730. 

  730.             SUMLAI = SUMLAI + LAIMX(K) * FAREA
    731. 

  731.             SUMLMN = SUMLMN + LAIMN(K) * FAREA
    732. 

  732.             ALAI   = ALAI + FAREA
    733. 

  733.             SUMRSI = SUMRSI + FAREA * LAIMX(K) / RSMIN(K)
    734. 

  734.             SUMLZ0 = SUMLZ0 + FAREA * ALOG(Z0(K))
    735. 

  735.             SUMVEG = SUMVEG + FAREA * VEG(K)
    736. 

  736.             SUMVMN = SUMVMN + FAREA * VEGMN(K)
    737. 

  737.             SUMSNUP= SUMSNUP+ FAREA * SNUP(K)
    738. 

  738.             SUMALB = SUMALB + FAREA * ALBF(K)
    739. 

  739.           ENDDO
    740. 

  740. 

  741.           !-- CHECK FOR WATER
    742. 

  742.           FWAT = LANDUSEF(I,KWAT,J)
    743. 

  743.           IF (FWAT .GT. 0.999) THEN
    744. 

  744.             XLAI(I,J)   = LAIMX(KWAT)
    745. 

  745.             XLAIMN(I,J) = LAIMN(KWAT)
    746. 

  746.             RSTMIN(I,J) = RSMIN(KWAT)
    747. 

  747.             ZNT(I,J)    = Z0(KWAT)
    748. 

  748.             XVEG(I,J)   = VEG(KWAT)
    749. 

  749.             XVEGMN(I,J) = VEGMN(KWAT)
    750. 

  750.             XSNUP(I,J)  = SNUP(KWAT)
    751. 

  751.             XALB(I,J)   = ALBF(KWAT)
    752. 

  752.           ELSE
    753. 

  753.             IF (FWAT .GT. 0.10) THEN
    754. 

  754.               ALAI   = ALAI - FWAT
    755. 

  755.               SUMLZ0 = SUMLZ0 - FWAT * ALOG(Z0(KWAT))
    756. 

  756.             ENDIF
    757. 

  757.             XLAI(I,J)   = SUMLAI / ALAI
    758. 

  758.             XLAIMN(I,J) = SUMLMN / ALAI
    759. 

  759.             RSTMIN(I,J) = SUMLAI / SUMRSI
    760. 

  760.             ZNT(I,J)    = EXP(SUMLZ0/ALAI)
    761. 

  761.             XVEG(I,J)   = SUMVEG / ALAI
    762. 

  762.             XVEGMN(I,J) = SUMVMN / ALAI
    763. 

  763.             XSNUP(I,J)  = SUMSNUP/ALAI
    764. 

  764.             XALB(I,J)   = SUMALB/ALAI
    765. 

  765.           ENDIF
    766. 

  766. 

  767.           IF (FWAT .GT. 0.50) THEN
    768. 

  768.             ZNT(I,J) = Z0(KWAT)
    769. 

  769.             XALB(I,J)= ALBF(KWAT)
    770. 

  770.           ENDIF
    771. 

  771. 

  772.           ZNT(I,J)    = ZNT(I,J) * 0.01           !CONVERT TO M
    773. 

  773.           XVEG(I,J)   = XVEG(I,J) * 0.01          !CONVERT TO FRAC
    774. 

  774.           XVEGMN(I,J) = XVEGMN(I,J) * 0.01
    775. 

  775.           XLAND(I,J)  = 1.0 + FWAT
    776. 

  776.           XALB(I,J)   = XALB(I,J) * 0.01
    777. 

  777.         
    778. 

  778.         ENDDO		! END LOOP THROUGH GRID CELLS
    779. 

  779.       ENDDO		  ! END LOOP THROUGH GRID CELLS
    780. 

  780.       !--------------------------------------------------------------------
    781. 

  781. 

  782.   END SUBROUTINE vegeland
    783. 

  783. 

  784. !------------------------------------------------------------------------------
    785. 

  785. !------------------------------------------------------------------------------
    786. 

  786.       SUBROUTINE SURFPX(DTPBL,  IFLAND, ISNOW,  NUDGEX, XICE1,               & !in
    787. 

  787.                         SOLDN,  GSW, LWDN,   EMISSI, Z1,                     & !in
    788. 

  788.                         MOL,    ZNT,    UST,                                 & !in 
    789. 

  789.                         PSURF,  DENS1,  QV1,    QSS,   TA1,                  & !in
    790. 

  790.                         THETA1, PRECIP,                                      & !in
    791. 

  791.                         CPAIR, PSIH,                                         & !in
    792. 

  792.                         RH2OBS, T2OBS,                                       & !in
    793. 

  793.                         VEGFRC, ISTI,   LAI,   BETAP,                        & !in
    794. 

  794.                         RSTMIN, HC_SNOW, SNOW_FRA,                           & !in
    795. 

  795.                         WWLT, WFC, CGSAT, WSAT, B,                           & !in
    796. 

  796.                         C1SAT, C2R, AS, JP, DS1, DS2,QST12,                  & !in
    797. 

  797.                         RADNET, GRDFLX, HFX, QFX, LH, EG, ER,  ETR,          & !out
    798. 

  798.                         QST,    CAPG,   RS, RA,  TG,     T2,                 & !out
    799. 

  799.                         WG,     W2,     WR, TA2, QA2, LAND_USE_TYPE   )        !out
    800. 

  800. 

  801. !------------------------------------------------------------------------------
    802. 

  802. ! 
    803. 

  803. !  FUNCTION:
    804. 

  804. !    THIS SUBROUTINE COMPUTES SOIL MOISTURE AND TEMPERATURE TENDENCIES
    805. 

  805. !    BY SOLVING THE PROGNOSTIC EQUATIONS IN PX95.
    806. 

  806. !
    807. 

  807. !  SUBROUTINES CALLED:
    808. 

  808. !    SUB. QFLUX   compute the soil and canopy evaporation, and transpiration
    809. 

  809. !    SUB. SMASS   compute nudging coefficients for soil moisture and temp nudging
    810. 

  810. !
    811. 

  811. !  ARGUMENTS:
    812. 

  812. !    DTPBL:   TIME STEP OF THE MINOR LOOP FOR THE LAND-SURFACE/PBL MODEL
    813. 

  813. !    IFLAND:  INDEX WHICH INDICATES THE TYPE OF SURFACE,=1,LAND;=2,SEA
    814. 

  814. !    ISNOW:   SNOW (=1) OR NOT (=0)
    815. 

  815. !    NUDGE:   SWITCH FOR SOIL MOISTURE NUDGING
    816. 

  816. !    SOLDN:   SHORT-WAVE RADIATION
    817. 

  817. !    LWDN:    LONG-WAVE RADIATION
    818. 

  818. !    EMISSI:  EMISSIVITY
    819. 

  819. !    Z1:      HEIGHT OF THE LOWEST HALF LAYER
    820. 

  820. !    MOL:     MONIN-OBUKOV LENGH (M)
    821. 

  821. !    ZNT:     ROUGHNESS LENGTH (M)
    822. 

  822. !    UST:     FRICTION VELOCITY (M/S)
    823. 

  823. !    TST:     Turbulent moisture scale
    824. 

  824. !    RA:      AERODYNAMIC RESISTENCE
    825. 

  825. !    PSURF:   P AT SURFACE (CB)
    826. 

  826. !    DENS1:   AIR DENSITY AT THE FIRST HALF LAYER
    827. 

  827. !    QV1:     Air humidity at first half layer
    828. 

  828. !    QSS:     Saturation mixing ratio at ground
    829. 

  829. !    TA1:     Air temperature at first half layer
    830. 

  830. !    THETA1:  Potential temperature at first half layer
    831. 

  831. !    PRECIP:  Precipitation rate in m/s
    832. 

  832. !    STBOLT:  STEFAN BOLTZMANN'S CONSTANT
    833. 

  833. !    KARMAN:  VON KARMAN CONSTANT
    834. 

  834. !    CPAIR:   Specific heat of moist air (M^2 S^-2 K^-1)
    835. 

  835. !    TAUINV:  1/1DAY(SEC)
    836. 

  836. !    VEGFRC:  Vegetation coverage
    837. 

  837. !    ISTI:    soil type
    838. 

  838. !    LAI:     Leaf area index
    839. 

  839. !    BETAP:   Coefficient for bare soil evaporation
    840. 

  840. !    THZ1OB:  Observed TEMP FROM ANAL OF OBS TEMPS AT SCREEN HT
    841. 

  841. !    RHOBS:   Observed relative humidity at SCREEN HT
    842. 

  842. !    RSTMIN   Minimum Stomatol resistence
    843. 

  843. !... Outputs from SURFPX
    844. 

  844. !    RADNET: Net radiation
    845. 

  845. !    HFX:    SENSIBLE HEAT FLUX (W / M^2)
    846. 

  846. !    QFX:    TOTAL EVAP FLUX (KG / M^2 S)
    847. 

  847. !    GRDFLX: Ground heat flux (W / M^2)
    848. 

  848. !    QST:    Turbulent moisture scale
    849. 

  849. !    CAPG:   THERMAL CAPACITY OF GROUND SLAB (J/M^2/K)
    850. 

  850. !    RS:     Surface resistence
    851. 

  851. !    RA:     Surface aerodynamic resistence
    852. 

  852. !    EG:     evaporation from ground (bare soil)
    853. 

  853. !    ER:     evaporation from canopy
    854. 

  854. !    ETR:    transpiration from vegetation
    855. 

  855. !    TA2:    diagnostic 2-m temperature from surface layer and lsm
    856. 

  856. !
    857. 

  857. !... Updated variables in this subroutine
    858. 

  858. !    TG:     Soil temperature at first soil layer
    859. 

  859. !    T2:     Soil temperature in root zone
    860. 

  860. !    WG:     Soil moisture at first soil layer
    861. 

  861. !    W2:     Soil moisture at root zone
    862. 

  862. !    WR:     Canopy water content in m
    863. 

  863. !
    864. 

  864. !  REVISION HISTORY:
    865. 

  865. !     AX     2/2005 - developed WRF version based on the MM5 PX LSM
    866. 

  866. !
    867. 

  867. !------------------------------------------------------------------------------
    868. 

  868. !------------------------------------------------------------------------------
    869. 

  869.       IMPLICIT NONE
    870. 

  870. 

  871. !.......Arguments
    872. 

  872. 

  873. !.. Integer
    874. 

  874.       INTEGER , INTENT(IN)  :: ISTI, NUDGEX
    875. 

  875. 

  876. !... Real
    877. 

  877.       REAL ,    INTENT(IN)  :: DTPBL, DS1, DS2
    878. 

  878.       REAL ,    INTENT(IN)  :: IFLAND, ISNOW, XICE1
    879. 

  879.       REAL ,    INTENT(IN)  :: SOLDN, GSW, LWDN, EMISSI, Z1
    880. 

  880.       REAL ,    INTENT(IN)  :: ZNT
    881. 

  881.       REAL ,    INTENT(IN)  :: PSURF, DENS1, QV1, QSS, TA1,  THETA1, PRECIP
    882. 

  882.       REAL ,    INTENT(IN)  :: CPAIR
    883. 

  883.       REAL ,    INTENT(IN)  :: VEGFRC, LAI
    884. 

  884.       REAL ,    INTENT(IN)  :: RSTMIN, HC_SNOW, SNOW_FRA 
    885. 

  885.       REAL ,    INTENT(IN)  :: WWLT, WFC, CGSAT, WSAT, B, C1SAT, C2R, AS, JP
    886. 

  886.       REAL ,    INTENT(IN)  :: RH2OBS,T2OBS
    887. 

  887.       REAL ,    INTENT(IN)  :: QST12
    888. 

  888. 

  889.       REAL ,    INTENT(OUT) :: RADNET, EG, ER, ETR
    890. 

  890.       REAL ,    INTENT(OUT) :: QST, CAPG, RS, TA2, QA2
    891. 

  891. 

  892.       REAL ,    INTENT(INOUT) :: TG, T2, WG, W2, WR, UST, RA, BETAP 
    893. 

  893.       REAL ,    INTENT(INOUT) :: GRDFLX, QFX, HFX, LH, PSIH, MOL
    894. 

  894. 

  895.       CHARACTER (LEN = 5), INTENT(IN) :: LAND_USE_TYPE
    896. 

  896. 

  897. !... Local Variables
    898. 

  898. 

  899. !... Real
    900. 

  900.       REAL        :: HF, LV, CQ4 
    901. 

  901.       REAL        :: RAH, RAW, ET, W2CG, CG, CT, SOILFLX, CPOT, THETAG
    902. 

  902.       REAL        :: ZOL, ZOBOL, ZNTOL, Y, Y0, PSIH15, YNT
    903. 

  903.       REAL        :: WGNUDG, W2NUDG, ALPH1, ALPH2, BET1, BET2, T1P5
    904. 

  904.       REAL        :: CQ1, CQ2, CQ3, COEFFNP1, COEFFN, TSNEW, TSHLF, T2NEW
    905. 

  905.       REAL        :: ROFF, WRMAX, PC, DWR, PNET, TENDWR, WRNEW
    906. 

  906.       REAL        :: COF1, CFNP1WR, CFNWR, PG, FASS
    907. 

  907.       REAL        :: TENDW2, W2NEW, W2HLF, W2REL, C1, C2, WEQ, CFNP1, CFN, WGNEW
    908. 

  908.       REAL        :: ALN10, TMP1, TMP2, TMP3, AA, AB, TST, RBH, CTVEG
    909. 

  909.       REAL        :: QST1,PHIH,PSIOB
    910. 

  910.       REAL        :: T2NUD, T2NUDF
    911. 

  911.       REAL        :: VAPPRS, QSBT, RH2MOD
    912. 

  912. 

  913. !... Parameters
    914. 

  914.       REAL        :: ZOBS, GAMAH, BETAH, SIGF, BH, CT_SNOW
    915. 

  915.       REAL, PARAMETER :: CV = 8.0E-6   ! K-M2/J
    916. 

  916.       PARAMETER (ZOBS  = 1.5)      ! height for observed screen temp., (m)
    917. 

  917.       PARAMETER (BH    = 15.7)
    918. 

  918.       PARAMETER (GAMAH = 16. )   !11.6)
    919. 

  919.       PARAMETER (BETAH = 5.0 )   !8.21)
    920. 

  920.       PARAMETER (SIGF  = 0.5)      ! rain interception see LSM (can be 0-1)
    921. 

  921.       !PARAMETER (CT_SNOW  = 5.54E-5)   
    922. 

  922. ! New value of CT_SNOW calibrated using multilayer soil model where csnow=6.9E5 J/(m3 K) 
    923. 

  923. ! from NCAR CSM                                         
    924. 

  924.       PARAMETER (CT_SNOW  = 2.0E-5)  
    925. 

  925.       ALN10  = ALOG(10.0)
    926. 

  926.         RADNET = SOLDN - (EMISSI *(STBOLT *TG **4 - LWDN))        ! NET RADIATION                             
    927. 

  927.         !--------------------------------------------------------------------
    928. 

  928.         CPOT= (100.0 / PSURF) ** ROVCP       ! rcp is global constant(module_model_constants)                            
    929. 

  929.         THETAG = TG * CPOT
    930. 

  930. 

  931.         ZOL   = Z1/MOL                                                       
    932. 

  932.         ZOBOL = ZOBS/MOL                                                      
    933. 

  933.         ZNTOL = ZNT/MOL                              
    934. 

  934. 

  935.         !-----------------------------------------------------------------------------------------
    936. 

  936.         IF (MOL .LT. 0.0) THEN
    937. 

  937.           Y      = ( 1.0 - GAMAH * ZOL ) ** 0.5
    938. 

  938.           Y0     = ( 1.0 - GAMAH * ZOBOL ) ** 0.5
    939. 

  939.           YNT    = ( 1.0 - GAMAH * ZNTOL ) ** 0.5
    940. 

  940.           PSIH15 =  2.0 * ALOG((Y + 1.0) / (Y0 + 1.0))
    941. 

  941.           PSIH   =  2.0 * ALOG((Y + 1.0) / (YNT + 1.0))
    942. 

  942.           PSIOB  =  2.0 * ALOG((Y0 + 1.0) / (YNT + 1.0))
    943. 

  943.           PHIH = 1.0 / Y
    944. 

  944.         ELSE
    945. 

  945.           IF((ZOL-ZNTOL).LE.1.0) THEN                                         
    946. 

  946.             PSIH = -BETAH*(ZOL-ZNTOL)                                        
    947. 

  947.           ELSE                                                               
    948. 

  948.             PSIH = 1.-BETAH-(ZOL-ZNTOL)                                      
    949. 

  949.           ENDIF             
    950. 

  950.           IF ((ZOBOL - ZNTOL) .LE. 1.0) THEN
    951. 

  951.             PSIOB = -BETAH * (ZOBOL - ZNTOL)
    952. 

  952.           ELSE
    953. 

  953.             PSIOB = 1.0 - BETAH - (ZOBOL - ZNTOL)
    954. 

  954.           ENDIF
    955. 

  955.           PSIH15 =  PSIH - PSIOB
    956. 

  956.           IF(ZOL.LE.1.0) THEN
    957. 

  957.             PHIH = 1.0 + BETAH * ZOL
    958. 

  958.           ELSE
    959. 

  959.             PHIH = BETAH + ZOL
    960. 

  960.           ENDIF
    961. 

  961.         ENDIF
    962. 

  962.         !-----------------------------------------------------------------------------------------
    963. 

  963.         !-- ADD RA AND RB FOR HEAT AND MOISTURE
    964. 

  964.         !... RB FOR HEAT = 5 /UST
    965. 

  965.         !... RB FOR WATER VAPOR =  5*(0.599/0.709)^2/3 /UST = 4.503/UST
    966. 

  966.         RA=PR0* ( ALOG(Z1/ZNT) - PSIH )/(KARMAN*UST)
    967. 

  967.         RAH = RA + 5.0 / UST
    968. 

  968.         RAW = RA + 4.503 / UST
    969. 

  969.         !--------------------------------------------------------------------
    970. 

  970.         !--  COMPUTE MOISTURE FLUX
    971. 

  971.         CALL QFLUX( DENS1,  QV1,    TA1,  SOLDN,  RAW, QSS,            &
    972. 

  972.                     VEGFRC, ISNOW,  ISTI, IFLAND, LAI, BETAP,          &
    973. 

  973.                     WG,     W2,     WR,                                &
    974. 

  974.                     RSTMIN, WWLT, WFC,                                 &
    975. 

  975.                     EG,     ER,     ETR,  CQ4,    RS,  FASS)
    976. 

  976.         !--------------------------------------------------------------------
    977. 

  977. 

  978.         !--------------------------------------------------------------------
    979. 

  979.         !..........Total evaporation (ET)
    980. 

  980.         ET  = EG + ER + ETR
    981. 

  981.         QST = -ET / (DENS1 * UST)
    982. 

  982.  
    983. 

  983.         LV  = 2.83E6                         !-- LATENT HEAT OF SUBLIMATION AT 0C FROM STULL(1988)
    984. 

  984.         IF (ISNOW .LT. 0.5.AND.TG.GT.273.15)                            &                                                                               
    985. 

  985.                     LV = (2.501 - 0.00237 * (TG - 273.15)) * 1.E6  !-- FROM STULL(1988) in J/KG
    986. 

  986.         
    987. 

  987.         IF (IFLAND .LT. 1.5 )  QFX = ET     !-- Recaculate QFX over land to account for P-X LSM EG, ER and ETR
    988. 

  988.         LH  = LV * QFX
    989. 

  989.         !-----------------------------------------------------------------------------------------
    990. 

  990. 

  991.         !-----------------------------------------------------------------------------------------
    992. 

  992.         ! Surface sensible heat flux
    993. 

  993.         TST = (THETA1 - THETAG ) / (UST*RAH)
    994. 

  994.         HF  = UST * TST           
    995. 

  995.         HFX = AMAX1(-DENS1 * CPAIR * HF, -250.0)  ! using -250. from MM5                                   
    996. 

  996.         !-----------------------------------------------------------------------------------------
    997. 

  997. 

  998.         !-----------------------------------------------------------------------------------------
    999. 

  999.         ! Compute the diagnosed 2m Q and T consistent with PX LSM
    1000. 

  1000.         QST1 = 0.5*(QST+QST12/PHIH)
    1001. 

  1001.         TA2  = (THETAG + TST * (PR0 / KARMAN * (ALOG(ZOBS / ZNT) - PSIOB)+5.))/CPOT
    1002. 

  1002.         QA2  = QV1 - QST1 * PR0/ KARMAN *  (ALOG(Z1 / ZOBS) - PSIH15)
    1003. 

  1003. 

  1004.         IF ((LAND_USE_TYPE == 'NLCD') .OR. (LAND_USE_TYPE == 'USGS')) THEN
    1005. 

  1005.            IF (QA2.LE.0.0) QA2 = QV1
    1006. 

  1006.         END IF
    1007. 

  1007. 

  1008.         !--  Relative humidity
    1009. 

  1009.         VAPPRS = SVP1 * EXP(SVP2 * (TA2 - SVPT0) / (TA2 - SVP3))
    1010. 

  1010.         QSBT   = EP_2 * VAPPRS / (PSURF - VAPPRS)
    1011. 

  1011.         RH2MOD = QA2 / QSBT
    1012. 

  1012.         !-----------------------------------------------------------------------------------------
    1013. 

  1013.         IF (IFLAND .LT. 1.5 ) THEN
    1014. 

  1014.           W2CG = AMAX1(W2,WWLT)
    1015. 

  1015.           CG   = CGSAT * 1.0E-6 * (WSAT/ W2CG) **    &  
    1016. 

  1016.                  (0.5 * B / ALN10)
    1017. 

  1017.           CT   = 1./((1-VEGFRC)/CG + VEGFRC/CV)
    1018. 

  1018.           CT   = 1./(SNOW_FRA/CT_SNOW + (1-SNOW_FRA)/CT)
    1019. 

  1019.           CAPG = 1.0/CT          
    1020. 

  1020. 

  1021.           SOILFLX = 2.0 * PI * TAUINV * (TG - T2)
    1022. 

  1022.           GRDFLX  = SOILFLX / CT
    1023. 

  1023.         ENDIF
    1024. 

  1024.         !-----------------------------------------------------------------------------------------
    1025. 

  1025. 

  1026.         !--------------------------------------------------------------------
    1027. 

  1027.         !-- ASSIMILATION --- COMPUTE SOIL MOISTURE NUDGING FROM TA2 and RH2               
    1028. 

  1028.         !-------COMPUTE ASSIMILATION COEFFICIENTS FOR ALL I                            
    1029. 

  1029.         IF (IFLAND .LT. 1.5) THEN
    1030. 

  1030.           IF (NUDGEX .EQ. 0) THEN                                          !-- NO NUDGING CASE                
    1031. 

  1031.             WGNUDG = 0.0                                                        
    1032. 

  1032.             W2NUDG = 0.0    
    1033. 

  1033.             T2NUD  = 0.0                                    
    1034. 

  1034.           ELSE                                                               !-- NUDGING CASE        
    1035. 

  1035.             
    1036. 

  1036.             CALL SMASS (ISTI,  FASS,  SOLDN,   VEGFRC, RA, WWLT, WFC,   &                       
    1037. 

  1037.                         ALPH1, ALPH2, BET1, BET2, T2NUDF)                             
    1038. 

  1038. 

  1039.             !--COMPUTE MODEL RH
    1040. 

  1040.             WGNUDG = ALPH1 * (T2OBS - TA2) + ALPH2 * (RH2OBS - RH2MOD) * 100  
    1041. 

  1041.             W2NUDG = BET1  * (T2OBS - TA2) + BET2  * (RH2OBS - RH2MOD) * 100
    1042. 

  1042.             IF (W2 .GE. WFC)  W2NUDG = AMIN1(W2NUDG,0.0)
    1043. 

  1043.             IF (W2 .LE. WWLT) W2NUDG = AMAX1(W2NUDG,0.0)
    1044. 

  1044.             T2NUD = T2NUDF * (T2OBS - TA2)
    1045. 

  1045.           !print *, 'T2NUD =',T2NUD,T2NUDF
    1046. 

  1046.           ENDIF
    1047. 

  1047.         ENDIF
    1048. 

  1048.         !-----------------------------------------------------------------------------------------
    1049. 

  1049. 

  1050.         !-----------------------------------------------------------------------------------------
    1051. 

  1051.         !-- Compute new values for TS,T2,WG,W2 and WR. No change over ice or water (XLAND > 1)
    1052. 

  1052.         IF (IFLAND .LT. 1.5) THEN
    1053. 

  1053.           !-- SOLVE BY CRANK-NIC --   TENDTS=CT*(RADNET-HFX-QFX)-SOILFLX 
    1054. 

  1054.           !-- Calculate the coefficients for implicit calculation of TG
    1055. 

  1055.           CQ1      = (1.0 - 0.622 * LV * CRANKP / (r_d * TG)) * QSS
    1056. 

  1056.           CQ2      = 0.622 * LV * QSS * CRANKP / (r_d * TG * TG)
    1057. 

  1057.           CQ3      = DENS1 * BETAP * (1.0 - VEGFRC) / RAW
    1058. 

  1058.           COEFFNP1 = 1.0 + DTPBL * CRANKP * (4.0 * EMISSI *  STBOLT * TG ** 3    &
    1059. 

  1059.                      * CT + DENS1 * CPAIR / RAH * CPOT * CT + 2.0 * PI           &
    1060. 

  1060.                      * TAUINV ) + DTPBL * (CT * LV * CQ2 * (CQ3 + CQ4))
    1061. 

  1061.           COEFFN   = CT * (GSW + EMISSI * (STBOLT * (4.0 * CRANKP - 1.0)         &
    1062. 

  1062.                      * TG*TG*TG*TG + LWDN)                                       & !NET RAD
    1063. 

  1063.                      + DENS1 * CPAIR / RAH * (THETA1 - (1.0 - CRANKP) * THETAG)  &
    1064. 

  1064.                      - LV * (CQ3 * (CQ1 - QV1) + CQ4 * (CQ1 - QV1)))             & !SFC HEAT FLUX 
    1065. 

  1065.                      - 2.0 * PI * TAUINV * ((1.0 - CRANKP) * TG - T2)              !SOIL FLUX
    1066. 

  1066.           TSNEW    = (TG + DTPBL * COEFFN) / COEFFNP1
    1067. 

  1067.           !-- FOR SNOW COVERED SURFACE TEMPERATURE IS NOT MORE THAN ZERO
    1068. 

  1068.           IF (XICE1 .GT. 0.5) TSNEW = AMIN1(TSNEW,273.15)                         ! Re-added Jan 2010 to keep ice surface at or below freezing (J. Pleim)
    1069. 

  1069.           TSHLF = 0.5 * ( TSNEW + TG)
    1070. 

  1070.           T2NEW = (T2 + DTPBL * TAUINV * T2TFAC * (TSHLF - (1 - CRANKP) * T2)     &
    1071. 

  1071.                 + DTPBL*T2NUD)  &              ! Added deep temperature nudging 
    1072. 

  1072.                   / (1.0 + DTPBL * TAUINV * T2TFAC * CRANKP)
    1073. 

  1073.           !-- REPLACE OLD with NEW Value
    1074. 

  1074.           TG = TSNEW
    1075. 

  1075.           T2 = T2NEW             
    1076. 

  1076.         ENDIF
    1077. 

  1077.         !-----------------------------------------------------------------------------------------
    1078. 

  1078. 

  1079.         !-----------------------------------------------------------------------------------------
    1080. 

  1080. 	      ! Compute new subsurface soil and canopy moisture values DENS1. No change required over ocean.
    1081. 

  1081.         IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN  
    1082. 

  1082.           !-- Compute WR
    1083. 

  1083.           ROFF  = 0.0
    1084. 

  1084.           WRMAX = 0.2E-3 * VEGFRC * LAI                     ! max. WRMAX IN m
    1085. 

  1085.           
    1086. 

  1086.           IF(WRMAX.GT.0.0) THEN
    1087. 

  1087.             !--  PC is precip. intercepted by veg.(M/S)
    1088. 

  1088.             PC    = VEGFRC * SIGF * PRECIP
    1089. 

  1089.             DWR   = (WRMAX - WR) / DTPBL                       !the tendency to reach max.
    1090. 

  1090.             PNET  = PC - ER/ DENW                              ! residual of precip. and evap.
    1091. 

  1091.             IF (PNET .GT. DWR) THEN
    1092. 

  1092.               ROFF = PNET - DWR
    1093. 

  1093.               PC   = PC - ROFF
    1094. 

  1094.             ENDIF
    1095. 

  1095.             IF (QSS .LT. QV1) THEN
    1096. 

  1096.               TENDWR = PC - ER / DENW
    1097. 

  1097.               WRNEW  = WR + DTPBL * TENDWR
    1098. 

  1098.             ELSE
    1099. 

  1099.               COF1    = DENS1 / DENW * VEGFRC * (QSS - QV1) / RAW
    1100. 

  1100.               !-- using delta=wr/wrmax
    1101. 

  1101.               CFNP1WR = 1.0 + DTPBL * COF1 * CRANKP / WRMAX  
    1102. 

  1102.               CFNWR   = PC - COF1 * (1.0 - CRANKP) * WR / WRMAX
    1103. 

  1103.               WRNEW   = (WR + DTPBL * CFNWR) / CFNP1WR
    1104. 

  1104.             ENDIF
    1105. 

  1105.           ELSE
    1106. 

  1106.             PC=0.0
    1107. 

  1107.             WRNEW=0.0  
    1108. 

  1108.           ENDIF
    1109. 

  1109.           !---------------------------------------------
    1110. 

  1110.           !-- Compute W2
    1111. 

  1111.           PG     = DENW * (PRECIP - PC)               ! PG is precip. reaching soil (PC already including ROFF)
    1112. 

  1112.           TENDW2 = 1.0 / (DENW * DS2) * (PG - EG - ETR) &
    1113. 

  1113.                    + (W2NUDG + WGNUDG) / DS2                    ! NUDGING
    1114. 

  1114.           W2NEW  = W2 + DTPBL * TENDW2
    1115. 

  1115.           W2NEW  = AMIN1(W2NEW,WSAT)
    1116. 

  1116.           W2NEW  = AMAX1(W2NEW,0.05)
    1117. 

  1117.           W2HLF  = 0.5 * (W2 + W2NEW)
    1118. 

  1118.           !.. new values
    1119. 

  1119.           W2     = W2NEW
    1120. 

  1120.           WR     = AMIN1(WRMAX,WRNEW)
    1121. 

  1121.         ENDIF
    1122. 

  1122.         !-----------------------------------------------------------------------------------------
    1123. 

  1123. 

  1124.         !-----------------------------------------------------------------------------------------
    1125. 

  1125. 	      ! Compute new surface soil moisture values (WR).
    1126. 

  1126.         IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN  ! over ocean no change to wg w2,wr
    1127. 

  1127.         !--    FOR SNOW COVERED SURFACE, ASSUME SURFACE IS SATURATED AND
    1128. 

  1128.         !      WG AND W2 ARE NOT CHANGED
    1129. 

  1129.           IF (ISNOW .GT.0.5) THEN
    1130. 

  1130.             WG = WSAT
    1131. 

  1131.           ELSE
    1132. 

  1132.             W2REL = W2HLF / WSAT
    1133. 

  1133.             IF (WG .GT. WWLT) THEN
    1134. 

  1134.               C1 = C1SAT * (WSAT / WG) ** (0.5 * B + 1.0) 
    1135. 

  1135.             ELSE   ! elimilate C1 for wg < wilting point
    1136. 

  1136.               C1 = C1SAT * (WSAT / WWLT) ** (0.5 * B + 1.0)
    1137. 

  1137.             ENDIF
    1138. 

  1138.             C2 = C2R * W2HLF / (WSAT - W2HLF + 1.E-11) 
    1139. 

  1139.             IF (W2HLF .EQ. WSAT) THEN
    1140. 

  1140.               WEQ = WSAT
    1141. 

  1141.             ELSE
    1142. 

  1142.               WEQ = W2HLF - AS * WSAT * W2REL ** JP *         &
    1143. 

  1143.                    (1.0 - W2REL ** (8.0 * JP))
    1144. 

  1144.             ENDIF
    1145. 

  1145. 

  1146.             !.... The beta method, Lee & Pielke (JAM, May 1992)
    1147. 

  1147.             CFNP1 = 1.0 + DTPBL * C2 * TAUINV * CRANKP
    1148. 

  1148.             CFN   = C1 / (DENW * DS1) * (PG - EG) - C2 * TAUINV *               &
    1149. 

  1149.                     ((1.0 - CRANKP) * WG - WEQ) + WGNUDG/ DS1                                    
    1150. 

  1150. 

  1151.             WGNEW = AMAX1((WG + DTPBL * CFN) / CFNP1,0.001)
    1152. 

  1152.             !-- NEW VALUES
    1153. 

  1153.             WG    = AMIN1(WGNEW,WSAT)
    1154. 

  1154.             
    1155. 

  1155.           ENDIF                  !endif for ISNOW
    1156. 

  1156.         ENDIF                    !endif for XLAND
    1157. 

  1157.             
    1158. 

  1158.       END SUBROUTINE surfpx
    1159. 

  1159. !-------------------------------------------------------------------------
    1160. 

  1160. !-------------------------------------------------------------------------
    1161. 

  1161. 

  1162. !-------------------------------------------------------------------------
    1163. 

  1163. !-------------------------------------------------------------------------
    1164. 

  1164.       SUBROUTINE QFLUX (DENS1,  QV1,   TA1,  RG,     RAW, QSS,           & ! in
    1165. 

  1165.                         VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP,         & ! in
    1166. 

  1166.                         WG,     W2,    WR,                               & ! in
    1167. 

  1167.                         RSTMIN, WWLT,  WFC,                              &
    1168. 

  1168.                         EG,     ER,    ETR,  CQ4,    RS,  FASS)            ! out
    1169. 

  1169. 

  1170. !-------------------------------------------------------------------------
    1171. 

  1171. !                                                              
    1172. 

  1172. !  FUNCTION:                                                          
    1173. 

  1173. !    THIS SUBROUTINE COMPUTES EVAPORATION FROM BARE SOIL (EG) AND FROM
    1174. 

  1174. !    THE WET PART OF CANOPY (ER) AND TRANSPIRATION FROM THE DRY PART OF
    1175. 

  1175. !    CANOPY (ETR).                       
    1176. 

  1176. !                                        
    1177. 

  1177. !  REVISION HISTORY: 
    1178. 

  1178. !    A. Xiu       Oct  2004  - adapted from the PX LSM in MM5 for the WRF system                            
    1179. 

  1179. !    R. Gilliam   Dec  2006  - Completed WRF V2.1.2 implementation                                 
    1180. 

  1180. !                                                               
    1181. 

  1181. !-------------------------------------------------------------------------
    1182. 

  1182. !  QFLUX ARGUMENT LIST:                                                        
    1183. 

  1183. !-------------------------------------------------------------------------
    1184. 

  1184. ! INPUTS:
    1185. 

  1185. !-- DENS1        air density at first layer
    1186. 

  1186. !-- QV1          air humidity at first layer
    1187. 

  1187. !-- TA1          air temperature at first layer
    1188. 

  1188. !-- RG           shortwave radition reaching the ground
    1189. 

  1189. !-- RAW          RA+RB for moisture
    1190. 

  1190. !-- QSS          saturation mixing ratio at ground
    1191. 

  1191. !-- VEGFRC       vegetation coverage
    1192. 

  1192. !-- ISNOW        if snow on the ground 
    1193. 

  1193. !-- ISTI         soil type
    1194. 

  1194. !-- IFLAND       if land (=1) or water (=2)
    1195. 

  1195. !-- LAI          leaf area index
    1196. 

  1196. !-- BETAP        
    1197. 

  1197. !-- WG           soil moisture at first soil layer
    1198. 

  1198. !-- W2           soil moisture at root zone
    1199. 

  1199. !-- WR           Canopy water
    1200. 

  1200. !
    1201. 

  1201. !  OUTPUTS FROM QFLUX:
    1202. 

  1202. !-- EG           evaporation from ground (bare soil)
    1203. 

  1203. !-- ER           evaporation from canopy
    1204. 

  1204. !-- ETR          transpiration from vegetation
    1205. 

  1205. !-- CQ4
    1206. 

  1206. !-- RS           surface resistence
    1207. 

  1207. !-- FASS         parameter for soil moisture nudging
    1208. 

  1208. !-------------------------------------------------------------------------
    1209. 

  1209. !-------------------------------------------------------------------------
    1210. 

  1210. 

  1211.       IMPLICIT NONE
    1212. 

  1212. 

  1213. ! DECLARATIONS - INTEGER
    1214. 

  1214.       INTEGER , INTENT(IN)  :: ISTI
    1215. 

  1215. 

  1216. ! DECLARATIONS - REAL
    1217. 

  1217.       REAL ,    INTENT(IN)  :: ISNOW, IFLAND
    1218. 

  1218.       REAL ,    INTENT(IN)  :: DENS1, QV1, TA1, RG, RAW, QSS,            &
    1219. 

  1219.                                VEGFRC, LAI,                              &
    1220. 

  1220.                                WG, W2, WR, RSTMIN   
    1221. 

  1221.       REAL ,    INTENT(INOUT)  :: BETAP
    1222. 

  1222.       REAL,     INTENT(IN)   :: WWLT, WFC
    1223. 

  1223. 

  1224.       REAL ,    INTENT(OUT) :: EG, ER, ETR, CQ4, RS, FASS
    1225. 

  1225. 

  1226. !... Local Variables
    1227. 

  1227. 

  1228. !... Real
    1229. 

  1229.       REAL    :: WRMAX, DELTA, SIGG, RADL, RADF, W2AVAIL, W2MXAV
    1230. 

  1230.       REAL    :: FTOT, F1, F2, F3,  F4
    1231. 

  1231.       REAL    :: FSHELT, GS, GA, FX
    1232. 

  1232. 

  1233. 

  1234. 

  1235. !... Parameters
    1236. 

  1236.       REAL, PARAMETER :: RSMAX = 5000.0                ! s/m
    1237. 

  1237.       REAL, PARAMETER :: FTMIN = 0.0000001             ! m/s
    1238. 

  1238.       REAL, PARAMETER :: F3MIN = 0.25
    1239. 

  1239. 

  1240. !
    1241. 

  1241. !... for water surface, no canopy evaporation and transpiration
    1242. 

  1242.       ER  = 0.0
    1243. 

  1243.       ETR = 0.0
    1244. 

  1244.       CQ4 = 0.0
    1245. 

  1245.       
    1246. 

  1246. !... GROUND EVAPORATION (DEPOSITION)
    1247. 

  1247.       IF (QSS .LT. QV1) BETAP = 1.0
    1248. 

  1248.       EG = DENS1 * (1.0 - VEGFRC) * BETAP * (QSS - QV1) / RAW
    1249. 

  1249. 

  1250.       !!---------------------------------------------------------------------
    1251. 

  1251. !... CANOPY
    1252. 

  1252.       IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0)  THEN
    1253. 

  1253.         WRMAX = 0.2E-3 * VEGFRC * LAI   ! in unit m
    1254. 

  1254.         IF (WR .LE. 0.0) THEN
    1255. 

  1255.           DELTA = 0.0
    1256. 

  1256.         ELSE
    1257. 

  1257. !         DELTA = (WR / WRMAX) ** 0.66667
    1258. 

  1258.           DELTA = WR / WRMAX           ! referred to SiB model
    1259. 

  1259.         ENDIF
    1260. 

  1260.         
    1261. 

  1261.         IF (QSS .GE. QV1) THEN
    1262. 

  1262.           SIGG = DELTA
    1263. 

  1263.         ELSE
    1264. 

  1264.           SIGG = 1.0
    1265. 

  1265.         ENDIF
    1266. 

  1266. 

  1267.         ER = DENS1 * VEGFRC * SIGG * (QSS - QV1) / RAW  
    1268. 

  1268.       ENDIF
    1269. 

  1269.       !!---------------------------------------------------------------------
    1270. 

  1270. 

  1271.       !-- TRANSPIRATION
    1272. 

  1272.       !!---------------------------------------------------------------------
    1273. 

  1273.       IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0)  THEN
    1274. 

  1274.         
    1275. 

  1275.         !!!-RADIATION
    1276. 

  1276.         IF (RSTMIN .GT. 130.0) THEN
    1277. 

  1277.           RADL = 30.0                                            ! W/M2
    1278. 

  1278.         ELSE
    1279. 

  1279.           RADL = 100.0                                           ! W/M2
    1280. 

  1280.         ENDIF
    1281. 

  1281.         RADF = 1.1 * RG / (RADL * LAI)                  ! NP89 - EQN34
    1282. 

  1282.         F1   = (RSTMIN / RSMAX + RADF) / (1.0 + RADF)
    1283. 

  1283. 

  1284.         !!!-SOIL MOISTURE
    1285. 

  1285.         W2AVAIL = W2 - WWLT
    1286. 

  1286.         W2MXAV  = WFC - WWLT
    1287. 

  1287.         F2      = 1.0 / (1.0 + EXP(-5.0 * ( W2AVAIL / W2MXAV -               &
    1288. 

  1288.                   (W2MXAV / 3.0 + WWLT))))    ! according JP, 9/94
    1289. 

  1289. 

  1290.         !-AIR TEMP
    1291. 

  1291.         !... according to Avissar (1985) and AX 7/95
    1292. 

  1292.         IF (TA1 .LE. 302.15) THEN
    1293. 

  1293.           F4 = 1.0 / (1.0 + EXP(-0.41 * (TA1 - 282.05)))
    1294. 

  1294.         ELSE
    1295. 

  1295.           F4 = 1.0 / (1.0 + EXP(0.5 * (TA1 - 314.0)))
    1296. 

  1296.         ENDIF
    1297. 

  1297. 

  1298.         FTOT = LAI * F1 * F2 * F4
    1299. 

  1299.       ENDIF
    1300. 

  1300. 

  1301.       !!---------------------------------------------------------------------
    1302. 

  1302.       IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0)  THEN
    1303. 

  1303.         FSHELT = 1.0   ! go back to NP89
    1304. 

  1304.         GS     = FTOT / (RSTMIN * FSHELT)
    1305. 

  1305.         GA     = 1.0 / RAW
    1306. 

  1306.         !-- Compute humidity effect according to RH at leaf surf
    1307. 

  1307.         F3 = 0.5 * (GS - GA + SQRT(GA * GA + GA * GS *                       &
    1308. 

  1308.              (4.0 * QV1 / QSS - 2.0) + GS * GS)) / GS
    1309. 

  1309.         F3 = AMIN1(AMAX1(F3,F3MIN),1.0)
    1310. 

  1310.         RS = 1.0 / (GS * F3)
    1311. 

  1311.         
    1312. 

  1312.         !--- Compute Assimilation factor for soil moisture nudging - jp 12/94
    1313. 

  1313.         !--  Note that the 30 coef is to result in order 1 factor for max
    1314. 

  1314.         IF (RG .LT. 0.00001) THEN   ! do not nudge during night
    1315. 

  1315.           FX = 0.0
    1316. 

  1316.         ELSE
    1317. 

  1317.           FX = 30.0 * F1 * F4 * LAI / (RSTMIN * FSHELT)
    1318. 

  1318.         ENDIF
    1319. 

  1319. 

  1320.         FASS = FX
    1321. 

  1321.         ETR = DENS1 * VEGFRC * (1.0 - SIGG) * (QSS - QV1) / (RAW + RS)
    1322. 

  1322.         !..... CQ4 is used for the implicit calculation of TG in SURFACE
    1323. 

  1323.         CQ4 = DENS1 * VEGFRC * ((1.0 - SIGG) / (RAW + RS) + SIGG / RAW)
    1324. 

  1324.       ENDIF         
    1325. 

  1325. 

  1326.     END SUBROUTINE qflux
    1327. 

  1327. !------------------------------------------------------------------------------------------
    1328. 

  1328. !------------------------------------------------------------------------------------------
    1329. 

  1329. 

  1330. !------------------------------------------------------------------------------------------
    1331. 

  1331. !------------------------------------------------------------------------------------------
    1332. 

  1332.       SUBROUTINE SMASS (ISTI,  FASS,  RG,   VEGFRC, RA,              & !in                      
    1333. 

  1333.                         WWLT, WFC,                                   & !in
    1334. 

  1334.                         ALPH1, ALPH2, BET1, BET2, T2NUDF )             !out  
    1335. 

  1335. !------------------------------------------------------------------------------------------
    1336. 

  1336. !------------------------------------------------------------------------------------------
    1337. 

  1337.       IMPLICIT NONE
    1338. 

  1338. !------------------------------------------------------------------------------------------
    1339. 

  1339. !     SMASS COMPUTES SOIL MOISTURE NUDGING COEFFICIENTS
    1340. 

  1340. !------------------------------------------------------------------------------------------
    1341. 

  1341. !
    1342. 

  1342. !.........Arguments
    1343. 

  1343.       INTEGER, PARAMETER             :: NSCAT    = 16     ! max. soil types
    1344. 

  1344.       
    1345. 

  1345.       INTEGER,                   INTENT(IN)   :: ISTI 
    1346. 

  1346.       REAL,                      INTENT(IN)   :: FASS, RG, VEGFRC, RA
    1347. 

  1347.       REAL,                      INTENT(IN)   :: WWLT, WFC
    1348. 

  1348.       REAL,                      INTENT(OUT)  :: ALPH1, ALPH2, BET1, BET2, T2NUDF
    1349. 

  1349. 

  1350. 

  1351. !........Local variables
    1352. 

  1352. 

  1353. !... Real
    1354. 

  1354.       REAL    :: FBET, FALPH, FRA, FTEXT
    1355. 

  1355.       REAL,    DIMENSION( 1: NSCAT ) :: WFCX, WWLTX
    1356. 

  1356.       
    1357. 

  1357. !... Parameters
    1358. 

  1358.       REAL, PARAMETER  :: A1MAX = -10.E-5, A2MAX = 1.E-5  ! m/K, m for 6hr period
    1359. 

  1359.       REAL, PARAMETER  :: B1MAX = -10.E-3, B2MAX = 1.E-3  ! m/K, m (Bouttier et al 1993)
    1360. 

  1360.       REAL, PARAMETER  :: TASSI = 4.6296E-5               ! 1/6hr in 1/sec
    1361. 

  1361.       REAL, PARAMETER  :: RAMIN = 10.0                    ! 0.1 s/cm
    1362. 

  1362. !
    1363. 

  1363. !-- WFC is field capacity (M^3/M^3) (JN90)
    1364. 

  1364.       DATA WFCX   /  0.135, 0.150, 0.195, 0.255, 0.240, 0.255, 0.322,    &
    1365. 

  1365.                     0.325, 0.310, 0.370, 0.367, 0.367, 0.367, 0.367, 0.367, 0.367 /
    1366. 

  1366. !
    1367. 

  1367. !-- WWLT is wilting point (M^3/M^3) (JN90)
    1368. 

  1368.       DATA WWLTX  /  0.068, 0.075, 0.114, 0.179, 0.155, 0.175, 0.218,    &
    1369. 

  1369.                     0.250, 0.219, 0.283, 0.286, 0.286, 0.286, 0.286, 0.286, 0.286 /
    1370. 

  1370. 

  1371. !
    1372. 

  1372.       FBET  = FASS
    1373. 

  1373.       FALPH = RG / 1370.0
    1374. 

  1374. !--TEXTURE FACTOR NORMALIZED BY LOAM (IST=5)
    1375. 

  1375.       FRA   = RAMIN / RA            ! scale by aerodynamic resistance
    1376. 

  1376.       FTEXT = TASSI * (WWLT + WFC) / (WWLTX(5) + WFCX(5)) * FRA
    1377. 

  1377. !        write(6,*) ' ftot, fbet=',ftot, fbet,' ftext=',ftext/tassi
    1378. 

  1378. !
    1379. 

  1379.       ALPH1 = A1MAX * FALPH * (1.0 - VEGFRC) * FTEXT
    1380. 

  1380.       ALPH2 = A2MAX * FALPH * (1.0 - VEGFRC) * FTEXT
    1381. 

  1381.       BET1  = B1MAX * FBET  *        VEGFRC  * FTEXT
    1382. 

  1382.       BET2  = B2MAX * FBET  *        VEGFRC  * FTEXT
    1383. 

  1383.       T2NUDF = 1.0E-5 * MAX((1.0 - 5.0 * FALPH),0.0)   ! T2 Nudging at night
    1384. 

  1384. 

  1385.     END SUBROUTINE smass
    1386. 

  1386. !------------------------------------------------------------------------------------------
    1387. 

  1387. !------------------------------------------------------------------------------------------
    1388. 

  1388. 

  1389. !------------------------------------------------------------------------------------------
    1390. 

  1390. !------------------------------------------------------------------------------------------
    1391. 

  1391.       SUBROUTINE SOILPROP (SOILCBOT,WEIGHT, ITIMESTEP, MAVAIL, &  ! IN
    1392. 

  1392.                            PXLSM_SMOIS_INIT,                   &  ! IN
    1393. 

  1393.                            FWSAT,FWFC,FWWLT,FB,FCGSAT,         &  ! OUT
    1394. 

  1394.                            FJP,FAS,FC2R,FC1SAT,ISTI, W2        )  ! OUT
    1395. 

  1395.      !------------------------------------------------------------------------
    1396. 

  1396.      !     SOILPROP COMPUTES SOIL PARAMETERS FOR BOTH BOTTOM AND TOP LAYERS
    1397. 

  1397.      !              USING FRACTIONAL SOIL TYPE. A HARD CODED OPTION IS AVAILIABLE
    1398. 

  1398.      !              TO COMPUTE THE SOIL PARAMETERS USING FRACTIONAL INFORMATION, OR
    1399. 

  1399.      !              TO JUST USE SOIL PARAMETERS OF THE DOMINANT SOIL TYPE
    1400. 

  1400.      !------------------------------------------------------------------------
    1401. 

  1401.      !-- SOIL PARAMETERS ARE SPECIFIED BY SOIL TYPE: 
    1402. 

  1402.      !   #  SOIL TYPE  WSAT  WFC  WWLT    B   CGSAT   JP   AS   C2R  C1SAT
    1403. 

  1403.      !   _  _________  ____  ___  ____  ____  _____   ___  ___  ___  _____
    1404. 

  1404.      !   1  SAND       .395 .135  .068  4.05  3.222    4  .387  3.9  .082
    1405. 

  1405.      !   2  LOAMY SAND .410 .150  .075  4.38  3.057    4  .404  3.7  .098
    1406. 

  1406.      !   3  SANDY LOAM .435 .195  .114  4.90  3.560    4  .219  1.8  .132
    1407. 

  1407.      !   4  SILT LOAM  .485 .255  .179  5.30  4.418    6  .105  0.8  .153
    1408. 

  1408.      !   4  SILT       .485 .255  .179  5.30  4.418    6  .105  0.8  .153 NP89 does not have Silt so mapped to Silt Loam
    1409. 

  1409.      !   5  LOAM       .451 .240  .155  5.39  4.111    6  .148  0.8  .191
    1410. 

  1410.      !   6  SND CLY LM .420 .255  .175  7.12  3.670    6  .135  0.8  .213
    1411. 

  1411.      !   7  SLT CLY LM .477 .322  .218  7.75  3.593    8  .127  0.4  .385
    1412. 

  1412.      !   8  CLAY LOAM  .476 .325  .250  8.52  3.995   10  .084  0.6  .227
    1413. 

  1413.      !   9  SANDY CLAY .426 .310  .219 10.40  3.058    8  .139  0.3  .421
    1414. 

  1414.      !  10  SILTY CLAY .482 .370  .283 10.40  3.729   10  .075  0.3  .375
    1415. 

  1415.      !  11  CLAY       .482 .367  .286 11.40  3.600   12  .083  0.3  .342
    1416. 

  1416.      !------------------------------------------------------------------------
    1417. 

  1417. !------------------------------------------------------------------------------------------
    1418. 

  1418. !------------------------------------------------------------------------------------------
    1419. 

  1419.       IMPLICIT NONE
    1420. 

  1420. 

  1421.      !.........Arguments
    1422. 

  1422.       INTEGER, PARAMETER             :: NSCAT    = 16     ! max. soil types
    1423. 

  1423.       INTEGER, PARAMETER             :: NSCATMIN = 16     ! min soil types
    1424. 

  1424. 

  1425.       INTEGER,                      INTENT(IN)   :: WEIGHT, ITIMESTEP, PXLSM_SMOIS_INIT
    1426. 

  1426.       REAL,                         INTENT(IN)   :: MAVAIL
    1427. 

  1427.       REAL,     DIMENSION(1:NSCAT), INTENT(IN)   :: SOILCBOT
    1428. 

  1428.       REAL,                         INTENT(OUT)  :: FWSAT,FWFC,FWWLT,FB,FCGSAT,          &
    1429. 

  1429.                                                     FJP,FAS,FC2R,FC1SAT
    1430. 

  1430.       REAL,                         INTENT(INOUT)  :: W2
    1431. 

  1431. 

  1432. 

  1433.       INTEGER,                      INTENT(OUT)  :: ISTI
    1434. 

  1434. !........Local variables
    1435. 

  1435.       CHARACTER*4, AVCLASS
    1436. 

  1436.       CHARACTER*4, DIMENSION( 1: NSCAT ) ::  TEXID 
    1437. 

  1437. !... Integer
    1438. 

  1438.       INTEGER:: S
    1439. 

  1439. !... Real
    1440. 

  1440.       REAL:: TFRACBOT, CFRAC, SUMSND, SUMCLY, AVS, AVC, AVSLT
    1441. 

  1441.       REAL,    DIMENSION( 1: NSCAT ) :: WSAT, WFC, WWLT, B, CGSAT, AS,  &
    1442. 

  1442.                                         JP, C2R, C1SAT
    1443. 

  1443. 

  1444.       REAL,    DIMENSION( 1: NSCATMIN ) ::  SAND, CLAY
    1445. 

  1445. !.......... DATA statement for SOIL PARAMETERS for the 11 soil types
    1446. 

  1446.       DATA SAND /92.5,80.5,61.1,19.6,4.0,40.0,57.1,11.3,26.8,       &
    1447. 

  1447.                  52.0,6.5,10.2,1.0,1.0,1.0,1.0/                                                  
    1448. 

  1448.       DATA CLAY/2.1,4.1,10.9,19.1,7.3,18.8,23.3,32.2,36.6,          &                     
    1449. 

  1449.                 43.0,46.2,58.8,1.0,1.0,1.0,1.0/                                                 
    1450. 

  1450.       DATA TEXID/'Sand','Lsan','Sloa','Sill','Silt','Loam','Sclo',  &             
    1451. 

  1451.                  'Sicl','Cllo','Sacl','Sicy','Clay','Ormt','Wate',  &           
    1452. 

  1452.                  'Bedr','Othe'/                                         
    1453. 

  1453. 

  1454. !
    1455. 

  1455. !-- WSAT is saturated soil moisture (M^3/M^3) (JN90)
    1456. 

  1456.       DATA WSAT  /  0.395, 0.410, 0.435, 0.485, 0.451, 0.420, 0.477,    &
    1457. 

  1457.                     0.476, 0.426, 0.482, 0.482, 0.482, 0.482, 0.482, 0.482, 0.482 /
    1458. 

  1458. !
    1459. 

  1459. !-- WFC is field capacity (M^3/M^3) (JN90)
    1460. 

  1460.       DATA WFC   /  0.135, 0.150, 0.195, 0.255, 0.240, 0.255, 0.322,    &
    1461. 

  1461.                     0.325, 0.310, 0.370, 0.367, 0.367, 0.367, 0.367, 0.367, 0.367 /
    1462. 

  1462. !
    1463. 

  1463. !-- WWLT is wilting point (M^3/M^3) (JN90)
    1464. 

  1464.       DATA WWLT  /  0.068, 0.075, 0.114, 0.179, 0.155, 0.175, 0.218,    &
    1465. 

  1465.                     0.250, 0.219, 0.283, 0.286, 0.286, 0.286, 0.286, 0.286, 0.286 /
    1466. 

  1466. !
    1467. 

  1467. !-- B is slop of the retention curve (NP89)
    1468. 

  1468.       DATA B     /  4.05,  4.38,  4.90,  5.30,  5.39,  7.12,  7.75,     &
    1469. 

  1469.                     8.52, 10.40, 10.40, 11.40, 11.40, 11.40, 11.40, 11.40, 11.40  /
    1470. 

  1470. !
    1471. 

  1471. !-- CGSAT is soil thermal coef. at saturation (10^-6 K M^2 J^-1) (NP89)
    1472. 

  1472.       DATA CGSAT /  3.222, 3.057, 3.560, 4.418, 4.111, 3.670, 3.593,    &
    1473. 

  1473.                     3.995, 3.058, 3.729, 3.600, 3.600, 3.600, 3.600, 3.600, 3.600 /
    1474. 

  1474. !
    1475. 

  1475. !-- JP is coefficient of WGEQ formulation (NP89)
    1476. 

  1476.       DATA JP    /  4,     4,     4,     6,     6,     6,     8,        &
    1477. 

  1477.                    10,     8,    10,    12,    12,    12,    12,    12,    12     /
    1478. 

  1478. !
    1479. 

  1479. !-- AS is coefficient of WGEQ formulation (NP89)
    1480. 

  1480.       DATA AS    /  0.387, 0.404, 0.219, 0.105, 0.148, 0.135, 0.127,    &
    1481. 

  1481.                     0.084, 0.139, 0.075, 0.083, 0.083, 0.083, 0.083, 0.083, 0.083 /
    1482. 

  1482. !
    1483. 

  1483. !-- C2R is the value of C2 for W2=0.5WSAT (NP89)
    1484. 

  1484.       DATA C2R   /  3.9,   3.7,   1.8,   0.8,   0.8,   0.8,   0.4,      &
    1485. 

  1485.                     0.6,   0.3,   0.3,   0.3,   0.3,   0.3,   0.3,   0.3,   0.3   /
    1486. 

  1486. !
    1487. 

  1487. !-- C1SAT is the value of C1 at saturation (NP89)
    1488. 

  1488.       DATA C1SAT /  0.082, 0.098, 0.132, 0.153, 0.191, 0.213, 0.385,    &
    1489. 

  1489.                     0.227, 0.421, 0.375, 0.342, 0.342, 0.342, 0.342, 0.342, 0.342 /
    1490. 

  1490. !   
    1491. 

  1491. !-------------------------------Exicutable starts here--------------------
    1492. 

  1492. 

  1493.     IF(WEIGHT.GE.1.0) THEN   !Compute soil characteristics using weighting determined by fractional soil content
    1494. 

  1494.      FWSAT =0
    1495. 

  1495.      FWFC  =0
    1496. 

  1496.      FWWLT =0
    1497. 

  1497.      FB    =0
    1498. 

  1498.      FCGSAT=0
    1499. 

  1499.      FJP   =0
    1500. 

  1500.      FAS   =0
    1501. 

  1501.      FC2R  =0
    1502. 

  1502.      FC1SAT=0
    1503. 

  1503.      TFRACBOT =0
    1504. 

  1504.      CFRAC=0
    1505. 

  1505.      
    1506. 

  1506.      DO S=1,NSCAT
    1507. 

  1507. 

  1508.        IF(SOILCBOT(S).GE.CFRAC) THEN
    1509. 

  1509.          ISTI=S
    1510. 

  1510.          CFRAC=SOILCBOT(S)
    1511. 

  1511.        ENDIF
    1512. 

  1512.       
    1513. 

  1513.        TFRACBOT=TFRACBOT+SOILCBOT(S)
    1514. 

  1514.        FWSAT =FWSAT  +  WSAT(S)  *SOILCBOT(S)
    1515. 

  1515.        FWFC  =FWFC   +  WFC(S)   *SOILCBOT(S)
    1516. 

  1516.        FWWLT =FWWLT  +  WWLT(S)  *SOILCBOT(S)
    1517. 

  1517.        FB    =FB     +  B(S)     *SOILCBOT(S)
    1518. 

  1518.        FCGSAT=FCGSAT +  CGSAT(S) *SOILCBOT(S)
    1519. 

  1519.        FJP   =FJP    +  JP(S)    *SOILCBOT(S)
    1520. 

  1520.        FAS   =FAS    +  AS(S)    *SOILCBOT(S)
    1521. 

  1521.        FC2R  =FC2R   +  C2R(S)   *SOILCBOT(S)
    1522. 

  1522.        FC1SAT=FC1SAT +  C1SAT(S) *SOILCBOT(S)
    1523. 

  1523.              
    1524. 

  1524.      ENDDO
    1525. 

  1525. 

  1526.      TFRACBOT = 1/TFRACBOT
    1527. 

  1527.      FWSAT =FWSAT   *  TFRACBOT
    1528. 

  1528.      FWFC  =FWFC    *  TFRACBOT
    1529. 

  1529.      FWWLT =FWWLT   *  TFRACBOT
    1530. 

  1530.      FB    =FB      *  TFRACBOT
    1531. 

  1531.      FCGSAT=FCGSAT  *  TFRACBOT
    1532. 

  1532.      FJP   =FJP     *  TFRACBOT
    1533. 

  1533.      FAS   =FAS     *  TFRACBOT
    1534. 

  1534.      FC2R  =FC2R    *  TFRACBOT
    1535. 

  1535.      FC1SAT=FC1SAT  *  TFRACBOT
    1536. 

  1536.     
    1537. 

  1537.     ELSE                !Compute soil characteristics by sand and clay fraction
    1538. 

  1538.     
    1539. 

  1539.       CFRAC      = 0.0                                                               
    1540. 

  1540.       SUMSND     = 0.0                                                               
    1541. 

  1541.       SUMCLY     = 0.0                                                               
    1542. 

  1542.       TFRACBOT   = 0.0                                                              
    1543. 

  1543.                                                                                 
    1544. 

  1544.       DO S = 1,12                                                         
    1545. 

  1545. 

  1546.         TFRACBOT  = TFRACBOT  + SOILCBOT(S)
    1547. 

  1547.         SUMSND    = SUMSND    + SAND(S) * SOILCBOT(S)
    1548. 

  1548.         SUMCLY    = SUMCLY    + CLAY(S) * SOILCBOT(S)
    1549. 

  1549.          
    1550. 

  1550.         IF(SOILCBOT(S).GE.CFRAC) THEN   ! Find Dominant Category and fraction
    1551. 

  1551.           ISTI=S
    1552. 

  1552.           CFRAC=SOILCBOT(S)
    1553. 

  1553.         ENDIF
    1554. 

  1554. 

  1555.       ENDDO
    1556. 

  1556.          
    1557. 

  1557.       IF(TFRACBOT.GT.0.001) THEN                                                  
    1558. 

  1558.         AVS = SUMSND / TFRACBOT                                               
    1559. 

  1559.         AVC = SUMCLY / TFRACBOT                                           
    1560. 

  1560.         AVSLT = 100 - AVS - AVC                                   
    1561. 

  1561.         
    1562. 

  1562.         IF(AVS.GT.(85.+ 0.5*AVC)) THEN
    1563. 

  1563.           AVCLASS= 'Sand'                                                   
    1564. 

  1564.           ISTI = 1                                                          
    1565. 

  1565.         ELSE IF(AVS.GT.(70.+ AVC)) THEN
    1566. 

  1566.           AVCLASS= 'Lsan'                                             
    1567. 

  1567.           ISTI = 2                                                  
    1568. 

  1568.         ELSE IF((AVC.LT.20..AND.AVS.GT.52.) &
    1569. 

  1569.                 .OR.(AVC.LE.7.5.AND.AVSLT.LT.50.)) THEN                
    1570. 

  1570.           AVCLASS= 'Sloa'                                             
    1571. 

  1571.           ISTI = 3                                                     
    1572. 

  1572.         ELSE IF(AVC.LT.35..AND.AVS.GT.45..AND.AVSLT.LT.28.) THEN      
    1573. 

  1573.           AVCLASS= 'Sclo'                                                  
    1574. 

  1574.           ISTI = 6                                                     
    1575. 

  1575.         ELSE IF(AVC.GE.35..AND.AVS.GT.45.) THEN                           
    1576. 

  1576.           AVCLASS = 'Sacl'                                                  
    1577. 

  1577.           ISTI = 9                                                         
    1578. 

  1578.         ELSE IF(AVC.LT.27.0.AND.AVSLT.LT.50.) THEN                        
    1579. 

  1579.           AVCLASS= 'Loam'                                                
    1580. 

  1580.           ISTI = 5                                                        
    1581. 

  1581.         ELSE IF(AVC.LT.12..AND.AVSLT.GT.80.) THEN                
    1582. 

  1582.           AVCLASS = 'Silt'                                         
    1583. 

  1583.           ISTI = 4                                                           
    1584. 

  1584.         ELSE IF(AVC.LT.27.) THEN                                                 
    1585. 

  1585.           AVCLASS = 'Sill'                                                       
    1586. 

  1586.           ISTI = 4                                                              
    1587. 

  1587.         ELSE IF(AVC.LT.40..AND.AVS.GT.20.) THEN                                  
    1588. 

  1588.           AVCLASS = 'Cllo'                                                       
    1589. 

  1589.           ISTI = 8                                                               
    1590. 

  1590.         ELSE IF(AVC.LT.40.) THEN                                                 
    1591. 

  1591.           AVCLASS = 'Sicl'                                                       
    1592. 

  1592.           ISTI = 7                                                               
    1593. 

  1593.         ELSE IF(AVSLT.GE.40.) THEN                                               
    1594. 

  1594.           AVCLASS = 'Sicy'                                                       
    1595. 

  1595.           ISTI = 10                                                              
    1596. 

  1596.         ELSE                                                                     
    1597. 

  1597.           AVCLASS = 'Clay'                                                       
    1598. 

  1598.           ISTI = 11                                                              
    1599. 

  1599.         ENDIF                                                                          
    1600. 

  1600.       ELSE
    1601. 

  1601.         ISTI=5
    1602. 

  1602.         AVCLASS = TEXID(ISTI)  
    1603. 

  1603.       ENDIF
    1604. 

  1604.   
    1605. 

  1605.       FWSAT =WSAT(ISTI)   
    1606. 

  1606.       FWFC  =WFC(ISTI)    
    1607. 

  1607.       FWWLT =WWLT(ISTI)   
    1608. 

  1608.       FB    =B(ISTI)      
    1609. 

  1609.       FCGSAT=CGSAT(ISTI)  
    1610. 

  1610.       FJP   =JP(ISTI)     
    1611. 

  1611.       FAS   =AS(ISTI)     
    1612. 

  1612.       FC2R  =C2R(ISTI)   
    1613. 

  1613.       FC1SAT=C1SAT(ISTI)  
    1614. 

  1614.                             
    1615. 

  1615.     ENDIF
    1616. 

  1616.     
    1617. 

  1617.     ! Compute W2 using soil moisture availiability if pxlsm_smois_init (in namelist) is not zero
    1618. 

  1618.     IF (ITIMESTEP .EQ. 1 .AND. PXLSM_SMOIS_INIT .GT. 0) THEN
    1619. 

  1619.        W2 = FWWLT + (0.5*(FWSAT+FWFC) - FWWLT)  * MAVAIL 
    1620. 

  1620.     ENDIF                                                                    
    1621. 

  1621.     
    1622. 

  1622.     END SUBROUTINE soilprop
    1623. 

  1623. !------------------------------------------------------------------------------------------
    1624. 

  1624. !------------------------------------------------------------------------------------------
    1625. 

  1625. 

  1626. 

  1627. !------------------------------------------------------------------------------------------
    1628. 

  1628. !------------------------------------------------------------------------------------------
    1629. 

  1629.       SUBROUTINE PXSNOW (ITIMESTEP, ASNOW, CSNOW, SNOW,       &
    1630. 

  1630.                          SNOWH, SNUP,                         &
    1631. 

  1631.                          ALB, SNOALB, SHDFAC, SHDMIN,         & 
    1632. 

  1632.                          HC_SNOW, SNOW_FRA, SNOWC, SNOWALB)                   
    1633. 

  1633. !------------------------------------------------------------------------------------------
    1634. 

  1634. !------------------------------------------------------------------------------------------
    1635. 

  1635. !     Pleim-Xiu LSM Simple Snow model
    1636. 

  1636. !---------------------------------------------------------------------------------------------------
    1637. 

  1637. !     ITIMESTEP -- Model time step index
    1638. 

  1638. !     ASNOW --  Analyzed snow water equivalent (mm)
    1639. 

  1639. !     CSNOW --  Current snow water equivalent (mm)
    1640. 

  1640. !     SNOW  --  Snow water equivalent in model (mm)
    1641. 

  1641. !     SNOWH --  Physical snow depth (m)
    1642. 

  1642. !     SNUP --   Physical snow depth (landuse dependent) where when below, snow cover is fractional
    1643. 

  1643. !
    1644. 

  1644. !     HC_SNOW -- Heat capacity of snow as a function of depth
    1645. 

  1645. !     SNOW_FRA-- Factional snow area 
    1646. 

  1646. !     IFSNOW  -- Snow flag
    1647. 

  1647. !     SNOWALB -- Snow albedo
    1648. 

  1648. !---------------------------------------------------------------------------------------------------
    1649. 

  1649. 

  1650.       IMPLICIT NONE
    1651. 

  1651. 

  1652. !---  Arguments
    1653. 

  1653.       REAL, PARAMETER  :: W2SN_CONV   =   10.0
    1654. 

  1654.       REAL, PARAMETER  :: CS_SNOWPACK = 2092.0
    1655. 

  1655.       REAL, PARAMETER  :: SALP        =    2.6
    1656. 

  1656. !-----------------------------------------------
    1657. 

  1657.       INTEGER,       INTENT(IN)    :: ITIMESTEP
    1658. 

  1658.       REAL,          INTENT(IN)    :: ASNOW, CSNOW, SNUP, ALB, SNOALB, SHDFAC, SHDMIN
    1659. 

  1659.       REAL,          INTENT(INOUT) :: SNOW, SNOWH, SNOWC
    1660. 

  1660.       REAL,          INTENT(OUT)   :: HC_SNOW, SNOW_FRA, SNOWALB
    1661. 

  1661. !------------------------------------------------------------------------------------
    1662. 

  1662. 

  1663.                                                    
    1664. 

  1664. !-----------------------------------------------
    1665. 

  1665. !    Local variables
    1666. 

  1666.      !... Real
    1667. 

  1667.      REAL:: CONV_WAT2SNOW, CSNOWW, RHO_SNOWPACK,   &
    1668. 

  1668.             LIQSN_RATIO, SNEQV, RSNOW
    1669. 

  1669. !-----------------------------------------------
    1670. 

  1670.              
    1671. 

  1671.      SNEQV=ASNOW*0.001			              ! Snow water in meters
    1672. 

  1672.      RHO_SNOWPACK = 450                   ! Snowpack density (kg/m3), this should be computed in the future
    1673. 

  1673.      LIQSN_RATIO  = DENW/RHO_SNOWPACK     ! Ratio of water density and snowpack density
    1674. 

  1674.      
    1675. 

  1675.      CONV_WAT2SNOW = LIQSN_RATIO/1000     ! Conversion factor for snow liquid equiv. (mm) to snowpack (m)
    1676. 

  1676.       
    1677. 

  1677.      SNOW = ASNOW                         ! Set snow water to analysis value for now, implement a nudging scheme later
    1678. 

  1678.      SNOWH= SNOW * CONV_WAT2SNOW          ! Conversion of snow water (mm) to snow depth (m)     
    1679. 

  1679.      
    1680. 

  1680.      
    1681. 

  1681.      ! Is snow cover now present. The limit is 0.45 mm of water eqivalent or about 2 inches snow depth 
    1682. 

  1682.      SNOWC = 0.0
    1683. 

  1683.      IF (SNOWH .GT. 0.005) SNOWC = 1.0  
    1684. 

  1684.      
    1685. 

  1685.      HC_SNOW = RHO_SNOWPACK * CS_SNOWPACK * SNOWH
    1686. 

  1686.      
    1687. 

  1687.       IF (SNEQV < SNUP) THEN                                                  
    1688. 

  1688.          RSNOW = SNEQV / SNUP                                                    
    1689. 

  1689.          SNOW_FRA = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))          
    1690. 

  1690.       ELSE                                                                       
    1691. 

  1691.          SNOW_FRA = 1.0                                                            
    1692. 

  1692.       END IF   
    1693. 

  1693.       
    1694. 

  1694.       SNOWC  = SNOW_FRA
    1695. 

  1695. 

  1696.       SNOWALB = ALB + SNOWC*(SNOALB-ALB)
    1697. 

  1697. 

  1698.     END SUBROUTINE pxsnow
    1699. 

  1699. !------------------------------------------------------------------------------------------
    1700. 

  1700. !------------------------------------------------------------------------------------------
    1701. 

  1701. 

  1702. END MODULE module_sf_pxlsm
    1703. 

  1703. 

  1704.