| blob | a4bbcb77bd3cb5a56294b0b5f2966fdf6f97a3d4 |
1 MODULE module_sf_pxlsm
3 USE module_model_constants
4 USE module_sf_pxlsm_data
6 INTEGER, PARAMETER :: NSOLD=20
7 REAL, PARAMETER :: RD = 287.04, CPD = 1004.67, &
8 CPH2O = 4.218E+3, CPICE = 2.106E+3, &
9 LSUBF = 3.335E+5, SIGMA = 5.67E-8, &
10 ROVCP = RD / CPD
12 REAL, PARAMETER :: CRANKP = 0.5 ! CRANK-NIC PARAMETER
13 REAL, PARAMETER :: RIC = 0.25 ! critical Richardson number
14 REAL, PARAMETER :: DENW = 1000.0 ! water density in KG/M3
15 REAL, PARAMETER :: TAUINV = 1.0 / 86400.0 ! 1/1DAY(SEC)
16 REAL, PARAMETER :: T2TFAC = 1.0 / 10.0 ! Bottom soil temp response factor
17 REAL, PARAMETER :: PI = 3.1415926
18 REAL, PARAMETER :: PR0 = 0.95
19 REAL, PARAMETER :: CZO = 0.032
20 REAL, PARAMETER :: OZO = 1.E-4
22 CONTAINS
23 !
25 !-------------------------------------------------------------------------
26 !-------------------------------------------------------------------------
27 SUBROUTINE pxlsm(U3D, V3D, DZ8W, QV3D,T3D,TH3D, RHO, &
28 PSFC, GSW, GLW, RAINBL, EMISS, &
29 ITIMESTEP,CURR_SECS,NSOIL,DT,ANAL_INTERVAL, &
30 XLAND, XICE, ALBBCK, ALBEDO, &
31 SNOALB, SMOIS, TSLB, MAVAIL, TA2, &
32 QA2, QSFC, ZS,DZS, PSIH, &
33 LANDUSEF,SOILCTOP,SOILCBOT,VEGFRA,VEGF_PX, &
34 ISLTYP,RA,RS,LAI,IMPERV,CANFRA,NLCAT,NSCAT, &
35 HFX,QFX,LH,TSK,SST,ZNT,CANWAT, &
36 GRDFLX,SHDMIN,SHDMAX, &
37 SNOWC,PBLH,RMOL,UST,CAPG,DTBL, &
38 T2_NDG_OLD, T2_NDG_NEW, &
39 Q2_NDG_OLD, Q2_NDG_NEW, &
40 SN_NDG_OLD, SN_NDG_NEW, SNOW, SNOWH,SNOWNCV,&
41 T2OBS, Q2OBS, PXLSM_SMOIS_INIT, &
42 PXLSM_SOIL_NUDGE, &
43 pxlsm_modis_veg, &
44 LAI_PX, &
45 WWLT_PX, WFC_PX, WSAT_PX, &
46 CLAY_PX, CSAND_PX, FMSAND_PX, &
47 ids,ide, jds,jde, kds,kde, &
48 ims,ime, jms,jme, kms,kme, &
49 its,ite, jts,jte, kts,kte )
51 !-------------------------------------------------------------------------
52 ! THIS MODULE CONTAINS THE PLEIM-XIU LAND-SURFACE MODEL (PX-LSM).
53 ! IT IS DESIGNED TO SIMULATE CHARACTERISTICS OF THE LAND SURFACE AND
54 ! VEGETATION AND EXCHANGE WITH THE PLANETARY BOUNDARY LAYER (PBL). THE
55 ! SOIL MOISTURE MODEL IS BASED ON THE ISBA SCHEME DEVELOPED BY NOILHAN
56 ! AND PLANTON (1989) AND JACQUEMIN AND NOILHAN (1990) AND INCLUDES
57 ! PROGNOSTIC EQUATIONS FOR SOIL MOISTURE AND SOIL TEMPERATURE IN TWO
58 ! LAYERS (1 CM AND 1 M) AS WELL AS CANOPY WATER CONTENT. SURFACE
59 ! MOISTURE FLUXES ARE MODELED BY 3 PATHWAYS: SOIL EVAPORATION, CANOPY
60 ! EVAPORATION, AND VEGETATIVE EVAPOTRANSPIRRATION.
61 ! EVAPOTRANSPIRATION DIRECTLY FROM THE ROOT ZONE SOIL LAYER IS MODELED
62 ! VIA A CANOPY RESISTANCE ANALOG ALGORITHM WHERE STOMATAL CONDUCTANCE
63 ! IS CONTROLLED BY SOLAR RADIATION, AIR TEMPERATURE, AIR HUMIDITY, AND
64 ! ROOT ZONE SOIL MOISTURE. REQUIRED VEGETATION CHARACTERISTICS DERIVED
65 ! FROM THE USGS LANDUSE DATA INCLUDE: LEAF AREA INDEX, FRACTIONAL VEGETATION
66 ! COVERAGE, ROUGHNESS LENGTH, AND MINIMUM STOMATAL RESISTANCE. AN INDIRECT
67 ! NUDGING SCHEME ADJUSTS SOIL MOISTURE ACCORDING TO DIFFERENCES BETWEEN
68 ! MODELED TEMPERATURE AND HUMIDITY AND ANALYSED SURFACE FIELDS.
69 !
70 ! References:
71 ! Pleim and Xiu, 1995: Development and testing of a surface flux and planetary
72 ! boundary layer model for application in mesoscale models.
73 ! J. Appl. Meteoro., Vol. 34, 16-32.
74 ! Xiu and Pleim, 2001: Development of a land surface model. Part I: Application
75 ! in a mesoscale meteorological model. J. Appl. Meteoro.,
76 ! Vol. 40, 192-209.
77 ! Pleim and Xiu, 2003: Development of a land surface model. Part II: Data
78 ! assimilation. J. Appl. Meteoro., Vol. 42, 1811-1822.
79 !
80 ! Pleim and Gilliam, 2009: An Indirect Data Assimilation Scheme for Deep Soil Temperature in the
81 ! Pleim-Xiu Land Surface Model. J. Appl. Meteor. Climatol., 48, 1362-1376.
82 !
83 ! Gilliam and Pleim, 2010: Performance assessment of new land-surface and planetary boundary layer
84 ! physics in the WRF-ARW. Journal of Applied Meteorology and Climatology, 49, 760-774.
85 ! REVISION HISTORY:
86 ! AX 4/2005 - developed the initial WRF version based on the MM5 PX LSM
87 ! RG 2/2008 - Completed testing of the intial working version of PX LSM, released in WRFV3.0 early 2008
88 ! DW 8/2011 - Landuse specific versions of PX (USGS/NLCD/MODIS) were unified into a single code with
89 ! landuse characteristics defined in module_sf_pxsfclay.F.
90 ! RG 12/2011 - Basic code clean, removed commented out debug statements, lined up columns, etc.
91 ! RG 01/2012 - Removed FIRSTIME Logic that computes PX Landuse characteristics at first time step only. This resulted
92 ! in different solutions when OpenMP was used and would not work with moving domains.
93 ! RG 08/2012 - Added CURR_SECS variable in argument list as replacement for PX internal CURRTIME internally comp var.
94 ! This is neccessary for PX to correctly interpolate analyses for soil nudging. In this same calculation
95 ! logic was added for cases where user does not specify the analysis interval, or no analysis interval is
96 ! relevant as in the no PX soil nudging via namelist (pxlsm_soil_nudge = 0). Prior to this fix the default
97 ! analysis interval was zero, so if not speficied a divide by zero was the result. Also, changes were made to
98 ! ensure PX LSM will work with not only MODIS and USGS, but also both the 40 and 50 class NLCD-MODIS data.
99 ! Also, coupled module_sf_pxlsm_data.F was updated so landuse characteristics across datasets are more
100 ! consistent. Albedo for NLCD two grassland categories were lowered from 23 and 25 to 18 and 19.
101 ! For the NLCD40 and NLCD50 roughness and leaf area were made consistent between the US NLCD and
102 ! outside US MODIS datasets. Prior, US boundaries created boundaries of roughness and LAI.
103 !
104 ! RG 10/2014 - Wetlands soil moisture treatment. Grid cell soil moisture cannot fall less than fraction of a grid
105 ! cells wetland area * soil saturation (e.g., SMOIS of cell with 50% wetlands cannot fall below 50% of WSAT)
106 ! - Both soil levels are initialized using MVAVAIL (Soil moisture availability) instead of just layer 2.
107 ! - Veg Cv (heat capacity) changed from 8x10-6 to 1.2x10-5 (K-M2/J)
108 ! - Alternate empirical stomatal function of PAR (F1) to better replicate photosynthesis-conductance models.
109 ! The main effect is to reduce stomatal conductance for low PAR.
110 ! - Snow albedo is now computed using fractional land-use weighting. Values for each land-use class
111 ! are defined like other PX landuse parameters in module_sf_pxlsm_data.F. These are based on values
112 ! used by NOAH LSM MODIS in VEGPARM.TBL (MAXALB), but tuned to better match satellite values in maxsnowalb
113 ! dataset. Tuning reduced the MAXALB for all forest classes from values in the 50-60% range to 30-40% range.
114 ! These static values are more representative of albedo after snow has melted of fallen from trees. These
115 ! values were also cross verified with http://www.globalbedo.org/global.php
116 ! - USGS 28 category added as an option
117 ! - Impervious surface and canopy fraction data can be used if processed (otherwise 0% so no impact)
118 ! to alter surface heat capacity (See SURFPX subroutine for details) in urban areas and refine
119 ! LAI and VEGF_PX estimations (see VEGLAND subroutine).
120 !
121 ! JP 12/2015 - Surface water vapor mixing ratio calculation added for land surface, which is passed to PX-SFCLAY
122 ! for use over all non-water and non-frozen surfaces.
123 ! - PAR function and impact on transpiration modified according to Echer et al.(2015). See P-X LSM documentation
124 ! for full reference. These act to reduce moisture bias near surface during PBL transition.
125 !
126 ! JP 11/2017 - Updated vegetation table for different land cover types. Added in WRFv4.0.
127 !
128 ! LR 11/2017 - Update for MODIS vegetation: many changes in soil properties. Added in WRFv4.1.
129 ! (Ran et al., 2016 JGR-atmosphere, Ran et al. 2017 in preparation)
130 ! JP 12/2018 - revised soil type categories (ISTI) to conform to WRF soil type input data
131 ! soil types Sand through Clay are now 1-12 rather than 1-11
132 ! JP 12/2021 - Added new pathway for evaporation from the ground in the vegetated fraction of the grid cell.
133 !
134 !--------------------------------------------------------------------------------------------------------------
135 !--------------------------------------------------------------------------------------------------------------
136 ! ARGUMENT LIST:
137 !
138 !... Inputs:
139 !-- U3D 3D u-velocity interpolated to theta points (m/s)
140 !-- V3D 3D v-velocity interpolated to theta points (m/s)
141 !-- DZ8W dz between full levels (m)
142 !-- QV3D 3D mixing ratio
143 !-- T3D Temperature (K)
144 !-- TH3D Theta (K)
145 !-- RHO 3D dry air density (kg/m^3)
147 !-- PSFC surface pressure (Pa)
148 !-- GSW downward short wave flux at ground surface (W/m^2)
149 !-- GLW downward long wave flux at ground surface (W/m^2)
150 !-- RAINBL Timestep rainfall
151 !-- EMISS surface emissivity (between 0 and 1)
153 !-- ITIMESTEP time step number
154 !-- NSOIL number of soil layers
155 !-- DT time step (second)
156 !-- CURR_SECS time on model domain in seconds, universal WRF variable
157 !-- ANAL_INTERVAL Interval of analyses used for soil moisture and temperature nudging
159 !-- XLAND land mask (1 for land, 2 for water)
160 !-- XICE Sea ice
161 !-- ALBBCK Background Albedo
162 !-- ALBEDO surface albedo with snow cover effects
163 !-- SNOALB Albedo of snow
165 !-- SMOIS total soil moisture content (volumetric fraction)
166 !-- TSLB soil temp (k)
167 !-- MAVAIL Moisture availibility of soil
168 !-- TA2 2-m temperature
169 !-- QA2 2-m mixing ratio
171 !-- SVPT0 constant for saturation vapor pressure (K)
172 !-- SVP1 constant for saturation vapor pressure (kPa)
173 !-- SVP2 constant for saturation vapor pressure (dimensionless)
174 !-- SVP3 constant for saturation vapor pressure (K)
176 !-- ZS depths of centers of soil layers
177 !-- DZS thicknesses of soil layers
178 !-- PSIH similarity stability function for heat
180 !-- LANDUSEF Landuse fraction
181 !-- SOILCTOP Top soil fraction
182 !-- SOILCBOT Bottom soil fraction
183 !-- VEGFRA Vegetation fraction
184 !-- VEGF_PX Veg fraction recomputed and used by PX LSM
185 !-- ISLTYP Soil type
187 !-- RA Aerodynamic resistence
188 !-- RS Stomatal resistence
189 !-- LAI read in Leaf area index (weighted according to fractional landuse)
190 !-- ZNT rougness length
191 !-- QSFC Sat. water vapor mixing ratio at the surface interface
193 !-- IMPERV Fraction (percent) of grid cell that is impervious surface (concrete/road/non-veg)
194 !-- CANFRA Fraction (percent) of grid cell that is covered with tree canopy
196 !-- NLCAT Number of landuse categories
197 !-- NSCAT Number of soil categories
199 !-- HFX net upward heat flux at the surface (W/m^2)
200 !-- QFX net upward moisture flux at the surface (kg/m^2/s)
201 !-- LH net upward latent heat flux at surface (W/m^2)
202 !-- TSK surface skin temperature (K)
203 !-- SST sea surface temperature
204 !-- CANWAT Canopy water (mm)
206 !-- GRDFLX Ground heat flux
207 !-- SFCEVP Evaportation from surface
208 !-- SHDMIN Minimum annual vegetation fraction for each grid cell
209 !-- SHDMAX Maximum annual vegetation fraction for each grid cell
211 !-- SNOWC flag indicating snow coverage (1 for snow cover)
212 !-- PBLH PBL height (m)
213 !-- RMOL 1/L Reciprocal of Monin-Obukhov length
214 !-- UST u* in similarity theory (m/s)
215 !-- CAPG heat capacity for soil (J/K/m^3)
216 !-- DTBL time step of boundary layer calls
218 !-- T2_NDG_OLD Analysis temperature prior to current time
219 !-- T2_NDG_NEW Analysis temperature ahead of current time
220 !-- Q2_NDG_OLD Analysis mixing ratio prior to current time
221 !-- Q2_NDG_NEW Analysis mixing ratio ahead of current time
223 !-- SN_NDG_OLD Analysis snow water prior to current time
224 !-- SN_NDG_NEW Analysis snow water ahead of current time
225 !-- SNOW Snow water equivalent
226 !-- SNOWH Physical snow depth
227 !-- SNOWNCV Time step accumulated snow
229 !-- T2OBS Analysis temperature interpolated from prior and next in time analysese
230 !-- Q2OBS Analysis moisture interpolated from prior and next in time analysese
231 !-- PXLSM_SMOIS_INIT Flag to intialize deep soil moisture to a value derived from moisture availiability.
232 !-- PXLSM_SOIL_NUDGE Flag to use soil moisture and temperature nudging in the PX LSM
233 ! This is typically done for the first simulation.
234 !-- pxlsm_modis_veg Use MODIS vegetation option: 1 yes, 0 no
235 !-- LAI_PX LAI used for PX (m^2/m^2)
236 !-- WWLT_PX Computed soil wilting point for PX (m^3/m^3)
237 !-- WFC_PX Computed soil field capacity for PX (m^3/m^3)
238 !-- WSAT_PX Computed soil saturation for PX (m^3/m^3)
239 !-- CLAY_PX Aggregated soil clay fraction for PX (%)
240 !-- CSAND_PX Aggregated soil coarse sand fraction for PX (%)
241 !-- FMSAND_PX Aggregated soil fine-medium sand fraction for PX (%)
242 !
243 !-- ids start index for i in domain
244 !-- ide end index for i in domain
245 !-- jds start index for j in domain
246 !-- jde end index for j in domain
247 !-- kds start index for k in domain
248 !-- kde end index for k in domain
249 !-- ims start index for i in memory
250 !-- ime end index for i in memory
251 !-- jms start index for j in memory
252 !-- jme end index for j in memory
253 !-- kms start index for k in memory
254 !-- kme end index for k in memory
255 !-- jts start index for j in tile
256 !-- jte end index for j in tile
257 !-- kts start index for k in tile
258 !-- kte end index for k in tile
259 !... Outputs:
261 IMPLICIT NONE
263 !.......Arguments
264 ! DECLARATIONS - INTEGER
265 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
266 ims,ime, jms,jme, kms,kme, &
267 its,ite, jts,jte, kts,kte
269 INTEGER, INTENT(IN) :: NSOIL, ITIMESTEP, NLCAT, NSCAT, &
270 ANAL_INTERVAL, PXLSM_SMOIS_INIT, PXLSM_SOIL_NUDGE
272 REAL, INTENT(IN ),OPTIONAL :: curr_secs
274 REAL, INTENT(IN ) :: DT,DTBL
276 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ISLTYP
278 ! DECLARATIONS - REAL
279 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U3D, V3D, RHO, &
280 T3D, TH3D, DZ8W, QV3D
282 REAL, DIMENSION(1:NSOIL), INTENT(IN)::ZS,DZS
283 REAL, DIMENSION( ims:ime , 1:NSOIL, jms:jme ), INTENT(INOUT) :: SMOIS, TSLB
285 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: RA, RS, LAI, ZNT, QSFC
286 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT):: GRDFLX, TSK, TA2, QA2
288 REAL, DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN):: LANDUSEF
289 REAL, DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN):: SOILCTOP, SOILCBOT
292 REAL, DIMENSION( ims:ime, jms:jme ), &
293 INTENT(IN) :: PSFC, GSW, GLW, RAINBL, &
294 ALBBCK, SHDMIN, SHDMAX, &
295 PBLH, RMOL, SNOWNCV, &
296 UST, MAVAIL, SST, EMISS
298 REAL, DIMENSION( ims:ime, jms:jme ), &
299 INTENT(IN) :: T2_NDG_OLD, T2_NDG_NEW, &
300 Q2_NDG_OLD, Q2_NDG_NEW, &
301 SN_NDG_OLD, SN_NDG_NEW
303 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: T2OBS, Q2OBS
305 REAL, DIMENSION( ims:ime, jms:jme ), &
306 INTENT(INOUT) :: CAPG,CANWAT, QFX, HFX, LH, &
307 PSIH,VEGFRA, VEGF_PX, SNOW, SNOALB, &
308 SNOWH, SNOWC, ALBEDO, XLAND, XICE, &
309 IMPERV, CANFRA
311 INTEGER, OPTIONAL, INTENT(IN) :: pxlsm_modis_veg
313 REAL, DIMENSION( ims:ime, jms:jme ), &
314 OPTIONAL, INTENT(OUT) :: LAI_PX, WWLT_PX, WFC_PX, WSAT_PX, &
315 CLAY_PX, CSAND_PX, FMSAND_PX
316 INTEGER :: KWAT
318 LOGICAL :: radiation
320 !-------------------------------------------------------------------------
321 ! ---------- Local Variables --------------------------------
324 !---- PARAMETERS
325 INTEGER, PARAMETER :: NSTPS = 11 ! max. soil types
326 REAL, PARAMETER :: DTPBLX = 40.0 ! Max PX timestep = 40 sec
329 !---- INTEGERS
330 INTEGER, DIMENSION( 1: NSTPS ) :: JP
331 INTEGER:: J, I, NS, NUDGE, ISTI, WEIGHT
332 INTEGER:: NTSPS, IT
334 !---- REALS
335 REAL, DIMENSION( ims:ime, jms:jme ) :: XLAI, XLAIMN, RSTMIN, &
336 XVEG, XVEGMN, XSNUP, &
337 XALB, XSNOALB, WETFRA
339 REAL, DIMENSION( ims:ime, jms:jme ) :: RADNET, EG, ER, ETR, QST
341 REAL:: SFCPRS,TA1,DENS1,QV1,ZLVL,SOLDN,LWDN, &
342 EMISSI,PRECIP,THETA1,VAPPRS,QSBT, &
343 WG,W2,WR,TG,T2,USTAR,MOLX,Z0, &
344 RAIR,CPAIR,IFLAND,ISNOW, &
345 ES,QSS,BETAP, &
346 RH2_OLD, RH2_NEW, T2_OLD, T2_NEW, &
347 CORE, CORB, TIME_BETWEEN_ANALYSIS, &
348 G1000, ALN10,RH2OBS, HU, SNOBS, &
349 FWSAT,FWFC,FWWLT,FB,FCGSAT,FJP,FAS, &
350 FWRES, FC3, FCLAY, FCSAND, FFMSAND, &
351 FSEAS, T2I, HC_SNOW, SNOW_FRA,SNOWALB, &
352 QST12,ZFUNC,ZF1,ZA2,QV2, DT_FDDA, &
353 FC2R,FC1SAT, DTPBL, RAW
356 CHARACTER (LEN = 6) :: LAND_USE_TYPE
358 !-------------------------------------------------------------------------
359 !-------------------------------Executable starts here--------------------
360 !
361 ALN10 = ALOG(10.0)
362 G1000 = g*1.0E-3 ! G/1000
363 WEIGHT = 0
364 ! Determine Landuse Dataset by the number of categories
365 IF (NLCAT == 50) THEN
366 LAND_USE_TYPE = 'NLCD50'
367 ELSE IF (NLCAT == 40) THEN
368 LAND_USE_TYPE = 'NLCD40'
369 ELSE IF (NLCAT == 20) THEN
370 LAND_USE_TYPE = 'MODIS'
371 ELSE IF (NLCAT == 21) THEN
372 LAND_USE_TYPE = 'MODIS'
373 ELSE IF (NLCAT == 24) THEN
374 LAND_USE_TYPE = 'USGS'
375 ELSE IF (NLCAT == 28) THEN
376 LAND_USE_TYPE = 'USGS28'
377 ELSE
378 call wrf_error_fatal("Error: Unknown Land Use Category")
379 END IF
381 IF (ITIMESTEP .EQ. 1) THEN
382 CALL wrf_message( 'PX LSM will use the ' // TRIM(LAND_USE_TYPE) // ' landuse tables' )
383 PRINT *, 'The analysis interval for surface soil and temp nudging = ',ANAL_INTERVAL,'sec.'
384 ENDIF
386 !-----------------------------------------------------------------------------------
387 ! Kill WRF if user specifies soil nudging but provides no analysis interval, then provide helpful message.
388 IF (ANAL_INTERVAL .LE. 0.0 .AND. PXLSM_SOIL_NUDGE .EQ. 1) THEN
389 CALL wrf_message('PX LSM Error: The User specified analysis interval is zero or negative.')
390 CALL wrf_message('If the PX LSM is used with soil nudging (pxlsm_soil_nudge=1) a wrfsfdda_d0* file is required.')
391 CALL wrf_message('Make sure these files are present and')
392 CALL wrf_message('Check the namelist to ensure sgfdda_interval_m is set to proper sfc analysis interval')
393 STOP
394 ENDIF
395 !-----------------------------------------------------------------------------------
396 !--- Compute time relatve to old and new analysis time for timestep interpolation
397 IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN
398 DT_FDDA = ANAL_INTERVAL * 1.0 ! Convert DT of Analysis to real
399 TIME_BETWEEN_ANALYSIS = MOD(CURR_SECS,DT_FDDA)
400 IF (TIME_BETWEEN_ANALYSIS .EQ. 0.0) THEN
401 CORB = 1.0
402 CORE = 0.0
403 ELSE
404 CORE = TIME_BETWEEN_ANALYSIS / DT_FDDA
405 CORB = 1.0 - CORE
406 ENDIF
407 ENDIF
408 !-----------------------------------------------------------------------------------
409 ! Compute vegetation and land-use characteristics by land-use fraction weighting
410 ! These parameters include LAI, VEGF, ZNT, ALBEDO, SNOALB, RS, etc.
411 CALL VEGELAND(LANDUSEF, VEGFRA, SHDMIN, SHDMAX, &
412 SOILCTOP, SOILCBOT, NLCAT, NSCAT, &
413 ZNT,XLAI,XLAIMN,RSTMIN,XVEG,XVEGMN,XSNUP, &
414 XLAND, XALB,XSNOALB,WETFRA,IMPERV,CANFRA, &
415 ids,ide, jds,jde, kds,kde, &
416 ims,ime, jms,jme, kms,kme, &
417 its,ite, jts,jte, kts,kte, LAND_USE_TYPE, &
418 KWAT )
419 !-----------------------------------------------------------------------------------
422 !-----------------------------------------------------------------------------------
423 ! Main loop over individual grid cells
424 DO J = jts,jte !-- J LOOP
425 DO I = its,ite !-- I LOOP
427 IFLAND = XLAND(I,J)
429 ! Compute soil properties via weighting of fractional components
430 IF (IFLAND .LT. 1.5 ) THEN
432 !---------------------------------------------------------
433 CALL SOILPROP (SOILCBOT(I,:,J), WEIGHT, &
434 ITIMESTEP, MAVAIL(I,J), &
435 PXLSM_SMOIS_INIT, &
436 FWSAT,FWFC,FWWLT,FCLAY,FCSAND,FFMSAND, &
437 FB,FCGSAT, &
438 FJP,FAS,FC2R,FC1SAT,FWRES, FC3, ISTI, &
439 SMOIS(I,1,J), SMOIS(I,2,J) )
440 !----------------------------------------------------------
441 !----------------------------------------------------------
442 ISLTYP(I,J) = ISTI
443 ELSE
444 ISLTYP(I,J) = 14 ! STATSGO type for water
446 !-- aded for MODIS model
447 FWWLT = 0.1
448 FWFC = 1.0
449 FWSAT = 1.0
451 FCLAY = 0.0
452 FCSAND = 0.0
453 FFMSAND = 0.0
454 !-- end
456 ENDIF
458 !-- aded for MODIS model
459 WWLT_PX(I,J) = FWWLT
460 WFC_PX(I,J) = FWFC
461 WSAT_PX(I,J) = FWSAT
463 CLAY_PX(I,J) = FCLAY * 0.01 ! percent to fraction
464 CSAND_PX(I,J) = FCSAND * 0.01
465 FMSAND_PX(I,J) = FFMSAND * 0.01
466 !-- end
468 !-- Variables Sub. SURFPX needs
469 SFCPRS = PSFC(i,j) / 1000.0 ! surface pressure in cb
470 TA1 = T3D(i,1,j) ! air temperature at first layer
471 DENS1 = RHO(I,1,J) ! air density at first layer
472 QV1 = QV3D(i,1,j) ! water vapor mixing ratio at first layer
473 QV2 = QV3D(i,2,j)
474 ZLVL = 0.5 * DZ8W(i,1,j) ! thickness of lowest half level
475 ZF1 = DZ8W(i,1,j)
476 ZA2 = ZF1 + 0.5 * DZ8W(i,2,j)
478 LWDN = GLW(I,J) ! longwave radiation
479 EMISSI = EMISS(I,J) ! emissivity
480 PRECIP = MAX ( 1.0E-3*RAINBL(i,j)/DTBL,0.0) ! accumulated precip. rate during DT (=dtpbl)
481 ! convert RAINBL from mm to m for PXLSM
482 WR = 1.0E-3*CANWAT(I,J) ! convert CANWAT from mm to m for PXLSM
483 THETA1 = TH3D(i,1,j) ! potential temp at first layer
484 SNOBS = SNOW(I,J) ! Set snow cover to existing model value
485 ! this is overwritten below if snow analysis is availiable
486 ! otherwise snow cover remains constant through simulation
488 IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN
489 !-- 2 m Temp and RH for Nudging
490 T2_OLD = T2_NDG_OLD(I,J)
491 T2_NEW = T2_NDG_NEW(I,J)
492 VAPPRS = SVP1 * EXP(SVP2 * (T2_OLD - SVPT0) / ( T2_OLD - SVP3))
493 QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS)
494 RH2_OLD = Q2_NDG_OLD(I,J) / QSBT
495 VAPPRS = SVP1 * EXP(SVP2 * (T2_NEW - SVPT0) / (T2_NEW - SVP3))
496 QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS)
497 RH2_NEW = Q2_NDG_NEW(I,J) / QSBT
498 RH2OBS = CORB * RH2_OLD + CORE * RH2_NEW
499 T2OBS(I,J) = CORB * T2_OLD + CORE * T2_NEW
500 Q2OBS(I,J) = CORB * Q2_NDG_OLD(I,J) + CORE * Q2_NDG_NEW(I,J)
501 SNOBS = CORB * SN_NDG_OLD(I,J) + CORE * SN_NDG_NEW(I,J)
502 ENDIF
504 USTAR = MAX(UST(I,J),0.005)
506 IF (IFLAND .GE. 1.5) THEN ! if over water
507 ZNT(I,J) = CZO * UST(I,J) * UST(I,J) / G + OZO
508 ENDIF
510 Z0 = ZNT(I,J)
511 CPAIR = CPD * (1.0 + 0.84 * QV1) ! J/(K KG)
513 ! Set WRF Snow albedo to PX snow albedo based on fractional landuse
514 ! Snow albedo for each landuse class is defined in module_sf_pxlsm_data.F
515 SNOALB(I,J) = XSNOALB(I,J)
516 !-------------------------------------------------------------
517 ! Compute fractional snow area and snow albedo
518 CALL PXSNOW (ITIMESTEP, SNOBS, SNOWNCV(I,J), SNOW(I,J), &
519 SNOWH(I,J), XSNUP(I,J), XALB(i,j), &
520 SNOALB(I,J),VEGF_PX(I,J), SHDMIN(I,J), &
521 HC_SNOW, SNOW_FRA, SNOWC(I,J), ALBEDO(I,J) )
522 !-------------------------------------------------------------
524 !-------------------------------------------------------------
525 ! Sea Ice from analysis and water cells that are very cold, but more than 50% water
526 ! are converted to ice/snow for more reasonable treatment.
527 IF( (XICE(I,J).GE.0.5) .OR. &
528 (SST(I,J).LE.270.0.AND.XLAND(I,J).GE.1.50) ) THEN
529 XLAND(I,J) = 1.0
530 IFLAND = 1.0
531 ZNT(I,J) = 0.001 ! Ice
532 SMOIS(I,1,J) = 1.0 ! FWSAT
533 SMOIS(I,2,J) = 1.0 ! FWSAT
534 XICE(I,J) = 1.0
535 ALBEDO(I,J) = 0.7
536 SNOWC(I,J) = 1.0
537 SNOW_FRA = 1.0
538 VEGF_PX(I,J) = 0.0
539 LAI(I,J) = 0.0
540 ENDIF
541 !-------------------------------------------------------------
543 !-------------------------------------------------------------
544 !-- Note that when IFGROW = 0 is selected in Vegeland then max and min
545 !-- LAI and Veg are the same
546 T2I = TSLB(I,2,J)
547 FSEAS = AMAX1(1.0 - 0.015625 * (290.0 - T2I) ** 2,0.0)
548 IF (T2I .GE. 290.0) FSEAS = 1.0
550 !get PX table vegetation
551 LAI_PX(I,J) = XLAIMN(I,J) + FSEAS*(XLAI(I,J) - XLAIMN(I,J))
552 VEGF_PX(I,J) = XVEGMN(I,J) + FSEAS*(XVEG(I,J) - XVEGMN(I,J))
554 !... use MODIS LAI and VEGFRA from wrflowinp
555 IF ( pxlsm_modis_veg .EQ. 1 ) THEN
557 ! get LAI for vegetated area
558 IF ( VEGFRA(I,J) .GT. 0.0 ) THEN
559 LAI_PX(I,J) = LAI(I,J) / ( VEGFRA(I,J) / 100.0)
560 ELSE
561 LAI_PX(I,J) = 0.0
562 ENDIF
564 VEGF_PX(I,J) = VEGFRA(I,J) / 100.0
566 !vegF is just for the land
567 IF ( LANDUSEF(I,KWAT,J) .LT. 1.0 ) THEN
568 VEGF_PX(I,J) = VEGF_PX(I,J) / (1.0 - LANDUSEF(I,KWAT,J))
569 ELSE
570 VEGF_PX(I,J) = 0.0
571 ENDIF
573 ENDIF
575 LAI_PX(I,J) = MIN(LAI_PX(I,J), 8.0)
576 LAI_PX(I,J) = MAX(LAI_PX(I,J), 0.0001)
578 VEGF_PX(I,J) = MIN(VEGF_PX(I,J), 1.0)
579 VEGF_PX(I,J) = MAX(VEGF_PX(I,J), 0.0001)
581 !... END OF MODIS LAI and FPAR
584 ! Ensure veg algorithms not used for water
585 IF (IFLAND .GE. 1.5) THEN
586 VEGF_PX(I,J) = 0.0
587 LAI_PX(I,J) = 0.0
588 ENDIF
589 !-------------------------------------------------------------
592 SOLDN = GSW(I,J) / (1.0 - ALBEDO(I,J)) ! downward shortwave radiaton
593 ISNOW = SNOWC(I,J)
596 NUDGE=PXLSM_SOIL_NUDGE
597 IF ( J .LE. 2 .OR. J .GE. (jde-1) ) NUDGE=0
598 IF ( I .LE. 2 .OR. I .GE. (ide-1) ) NUDGE=0
600 IF ( RMOL(I,J) .GT. 0.0 ) THEN
601 MOLX = AMIN1(1/RMOL(I,J),1000.0)
602 ELSE IF ( RMOL(I,J) .LT. 0.0 ) THEN
603 MOLX = AMAX1(1/RMOL(I,J),-1000.0)
604 ELSE
605 MOLX = 1000
606 ENDIF
608 ZFUNC = ZF1 * (1.0 - ZF1 / MAX(100.,PBLH(I,J))) ** 2
609 QST12 = KARMAN * ZFUNC*(QV2-QV1) / (ZA2-ZLVL)
612 !-------------------------------------------------------------
613 !-- LSM sub-time loop too prevent dt > 40 sec
614 NTSPS = INT(DT / (DTPBLX + 0.000001) + 1.0)
615 DTPBL = DT / NTSPS
617 DO IT=1,NTSPS
619 !... SATURATION VAPOR PRESSURE (MB)
620 IF ( TSLB(I,1,J) .LE. SVPT0 ) THEN ! For ground that is below freezing
621 ES = SVP1 * EXP(22.514 - 6.15E3 / TSLB(I,1,J)) ! cb
622 ELSE
623 ES = SVP1 * EXP(SVP2 * (TSLB(I,1,J) - SVPT0) / (TSLB(I,1,J) - SVP3))
624 ENDIF
625 QSS = ES * 0.622 / (SFCPRS - ES)
627 !... beta method, Lee & Pielke (JAM,May1992)
628 BETAP = 1.0
629 IF (IFLAND .LT. 1.5 .AND. ISNOW .LT. 0.5 .AND. SMOIS(I,1,J) .LE. FWFC) THEN
630 BETAP = 0.25 * (1.0 - COS(SMOIS(I,1,J) / FWFC * PI)) ** 2
631 ENDIF
633 !-------------------------------------------------------------------------
634 CALL SURFPX (DTPBL, IFLAND, SNOWC(I,J), NUDGE, XICE(I,J), & !in
635 SOLDN, GSW(I,J), LWDN, EMISSI, ZLVL, & !in
636 MOLX, Z0, USTAR, & !in
637 SFCPRS, DENS1, QV1, QSS, TA1, & !in
638 THETA1, PRECIP, & !in
639 CPAIR, PSIH(I,J), & !in
640 RH2OBS,T2OBS(I,J), & !in
641 VEGF_PX(I,J), ISTI, LAI_PX(I,J), IMPERV(I,J), CANFRA(I,J), & !in
642 BETAP, RSTMIN(I,J), HC_SNOW, SNOW_FRA, WETFRA(I,J), & !in
643 FWWLT, FWFC, FWRES, FCGSAT, FWSAT, FB, & !in ! Soil model updates - JEP 12/14
644 FC1SAT,FC2R,FAS,FJP,FC3,DZS(1),DZS(2),QST12, & !in
645 RADNET(I,J), GRDFLX(I,J), HFX(I,J), QFX(I,J), LH(I,J), & !out
646 EG(I,J), ER(I,J), ETR(I,J), & !out
647 QST(I,J), CAPG(I,J), RS(I,J), RA(I,J), & !out
648 TSLB(I,1,J), TSLB(I,2,J), & !out
649 SMOIS(I,1,J), SMOIS(I,2,J), WR, &
650 TA2(I,J), QA2(I,J), LAND_USE_TYPE,I,J )
651 !-------------------------------------------------------------------------
653 END DO ! Time internal PX time loop
655 IF (IFLAND .LT. 1.5) TSK(I,J) = TSLB(I,1,J) ! ATV 02/20: changed for compatibility with sst_skin = 1.
656 CANWAT(I,J)= WR * 1000. ! convert WR back to mm for CANWAT
657 RAW = RA(I,J) + 4.503 / USTAR
658 QSFC(I,J) = QFX(I,J) * RAW / DENS1 + QV1
660 ENDDO ! END MIAN I LOOP
661 ENDDO ! END MAIN J LOOP
663 !------------------------------------------------------
664 END SUBROUTINE pxlsm
665 !-------------------------------------------------------------------------
666 !-------------------------------------------------------------------------
669 !-------------------------------------------------------------------------
670 !-------------------------------------------------------------------------
671 SUBROUTINE VEGELAND( landusef, vegfra, &
672 shdmin, shdmax, &
673 soilctop, soilcbot, nlcat, nscat,znt, xlai, &
674 xlaimn, rstmin, xveg, xvegmn, xsnup, xland, &
675 xalb, xsnoalb, wetfra, imperv, canfra, &
676 ids,ide, jds,jde, kds,kde, &
677 ims,ime, jms,jme, kms,kme, &
678 its,ite, jts,jte, kts,kte, &
679 LAND_USE_TYPE, KWAT_OUT )
680 !-------------------------------------------------------------------------
681 !
682 ! CALLED FROM Sub. bl_init in module_physics.init.F
683 !
684 ! THIS PROGRAM PROCESSES THE USGS LANDUSE DATA
685 ! WHICH HAS BEEN GRIDDED BY THE WPS SYSTEM
686 ! AND PRODUCES 2-D FIELDS OF LU RELATED PARAMETERS
687 ! FOR USE IN THE PX SURFACE MODEL
688 !
689 !
690 ! REVISION HISTORY:
691 ! AX Oct 2004 - developed WRF version based on VEGELAND in the MM5 PX LSM
692 ! RG Dec 2006 - revised for WRFV2.1.2
693 ! JP Dec 2007 - revised for WRFV3 - removed IFGROW options
694 !-------------------------------------------------------------------------
695 !-------------------------------------------------------------------------
697 IMPLICIT NONE
698 !...
699 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
700 ims,ime, jms,jme, kms,kme, &
701 its,ite, jts,jte, kts,kte
703 INTEGER , INTENT(IN) :: NSCAT, NLCAT
705 REAL, DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN) :: LANDUSEF
706 REAL, DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN) :: SOILCTOP, SOILCBOT
708 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: VEGFRA, SHDMIN, SHDMAX
710 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ZNT, IMPERV, CANFRA
712 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: XLAI, XLAIMN, RSTMIN, XALB, &
713 XVEG, XVEGMN, XSNUP, XLAND, &
714 WETFRA, XSNOALB
716 INTEGER, INTENT(OUT) :: KWAT_OUT
718 CHARACTER (LEN = 6), INTENT(IN) :: LAND_USE_TYPE
721 !... local variables
723 INTEGER :: ITF, JTF, K, J, I
724 REAL :: SUMLAI, SUMLMN, SUMRSI, SUMLZ0, SUMVEG, SUMVMN, &
725 ALAI, VEGF, SUMSNUP, SUMALB, SUMSNOALB
727 REAL :: VFMX, VFMN, VSEAS, FAREA, FWAT, ZNOTC, FCAN, FIMP, FORFRA, EXTFOR
728 REAL, DIMENSION( NLCAT ) :: LAIMX, LAIMN, Z0, VEG, VEGMN, SNUP, ALB, SNOALB
730 REAL, PARAMETER :: ZNOTCMN = 5.0 ! CM, MIN Zo FOR CROPS
731 REAL, PARAMETER :: ZNOTCMX = 15.0 ! CM, MAX Zo FOR CROPS
733 REAL, SAVE, DIMENSION(:), POINTER :: RSMIN, Z00, VEG0, VEGMN0, LAI0, &
734 LAIMN0, SNUP0, ALBF, SNOALBF
736 !---- INITIALIZE PARAMETERS
737 INTEGER, SAVE :: KWAT, LIMIT1, LIMIT2
739 ! Initialize LU characteristics by LU Dataset
740 IF (LAND_USE_TYPE == 'USGS') THEN
741 KWAT = 16
742 RSMIN => RSMIN_USGS
743 Z00 => Z00_USGS
744 VEG0 => VEG0_USGS
745 VEGMN0 => VEGMN0_USGS
746 LAI0 => LAI0_USGS
747 LAIMN0 => LAIMN0_USGS
748 SNUP0 => SNUP0_USGS
749 ALBF => ALBF_USGS
750 SNOALBF=> SNOALB_USGS
751 LIMIT1 = 2
752 LIMIT1 = 6
753 ELSE IF (LAND_USE_TYPE == 'USGS28') THEN
754 KWAT = 16
755 RSMIN => RSMIN_USGS28
756 Z00 => Z00_USGS28
757 VEG0 => VEG0_USGS28
758 VEGMN0 => VEGMN0_USGS28
759 LAI0 => LAI0_USGS28
760 LAIMN0 => LAIMN0_USGS28
761 SNUP0 => SNUP0_USGS28
762 ALBF => ALBF_USGS28
763 SNOALBF=> SNOALB_USGS28
764 LIMIT1 = 2
765 LIMIT1 = 6
766 ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN
767 KWAT = 1
768 RSMIN => RSMIN_NLCD50
769 Z00 => Z00_NLCD50
770 VEG0 => VEG0_NLCD50
771 VEGMN0 => VEGMN0_NLCD50
772 LAI0 => LAI0_NLCD50
773 LAIMN0 => LAIMN0_NLCD50
774 SNUP0 => SNUP0_NLCD50
775 ALBF => ALBF_NLCD50
776 SNOALBF=> SNOALB_NLCD50
777 LIMIT1 = 20
778 LIMIT1 = 43
779 ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN
780 KWAT = 17
781 RSMIN => RSMIN_NLCD40
782 Z00 => Z00_NLCD40
783 VEG0 => VEG0_NLCD40
784 VEGMN0 => VEGMN0_NLCD40
785 LAI0 => LAI0_NLCD40
786 LAIMN0 => LAIMN0_NLCD40
787 SNUP0 => SNUP0_NLCD40
788 ALBF => ALBF_NLCD40
789 SNOALBF=> SNOALB_NLCD40
790 LIMIT1 = 20
791 LIMIT1 = 43
792 ELSE IF (LAND_USE_TYPE == 'MODIS') THEN
793 KWAT = 17
794 RSMIN => RSMIN_MODIS
795 Z00 => Z00_MODIS
796 VEG0 => VEG0_MODIS
797 VEGMN0 => VEGMN0_MODIS
798 LAI0 => LAI0_MODIS
799 LAIMN0 => LAIMN0_MODIS
800 SNUP0 => SNUP0_MODIS
801 ALBF => ALBF_MODIS
802 SNOALBF=> SNOALB_MODIS
803 LIMIT1 = 12
804 LIMIT1 = 14
805 END IF
807 KWAT_OUT = KWAT
808 !--------------------------------------------------------------------
809 DO J = jts,jte
810 DO I = its,ite
811 XLAI(I,J) = 0.0
812 XLAIMN(I,J) = 0.0
813 RSTMIN(I,J) = 9999.0
814 XVEG(I,J) = 0.0
815 XVEGMN(I,J) = 0.0
816 XSNUP(I,J) = 0.0
817 XALB(I,J) = 0.0
818 XSNOALB(I,J)= 0.0
820 ! Code that may be needed in case these arrays are not intialized
821 ! with zero by real.exe when not defined by GEOGRID or present
822 ! in met_em* files
823 !IMPERV(I,J) = AMAX1(0.0001,IMPERV(I,J))
824 !CANFRA(I,J) = AMAX1(0.0001,CANFRA(I,J))
826 ENDDO ! END LOOP THROUGH GRID CELLS
827 ENDDO ! END LOOP THROUGH GRID CELLS
828 !--------------------------------------------------------------------
830 DO J = jts,jte
831 DO I = its,ite
832 !-- Initialize 2 and 3-D veg parameters to be caculated
833 DO K=1,NLCAT
834 LAIMX(K) = LAI0(K)
835 LAIMN(K) = LAIMN0(K)
836 Z0(K) = Z00(K)
837 VEG(K) = VEG0(K)
838 VEGMN(K) = VEGMN0(K)
839 SNUP(K) = SNUP0(K)
840 ALB(K) = ALBF(K)
841 SNOALB(K)= SNOALBF(K)
842 ENDDO
844 !-- INITIALIZE SUMS
845 SUMLAI = 0.0
846 SUMLMN = 0.0
847 SUMRSI = 0.0
848 SUMLZ0 = 0.0
849 SUMVEG = 0.0
850 SUMVMN = 0.0
851 ALAI = 0.0
852 SUMSNUP = 0.0
853 SUMALB = 0.0
854 SUMSNOALB = 0.0
856 !-- ESTIMATE CROP EMERGANCE DATE FROM VEGFRAC
857 VFMX = SHDMAX(I,J)
858 VFMN = SHDMIN(I,J)
859 VEGF = VEGFRA(I,J)
861 !-- Computations for VEGETATION CELLS ONLY
862 IF(VFMX.GT.0.0.AND.LANDUSEF(I,KWAT,J).LT.1.00) THEN
863 VSEAS = VEGF/VFMX
864 IF(VSEAS.GT.1.0.OR.VSEAS.LT.0.0) THEN
865 VSEAS = MIN(VSEAS,1.0)
866 VSEAS = MAX(VSEAS,0.0)
867 ENDIF
869 ZNOTC = ZNOTCMN * (1-VSEAS) + ZNOTCMX * VSEAS ! Zo FOR CROPS
870 DO K = 1, NLCAT
871 IF (LAND_USE_TYPE == 'MODIS') THEN
872 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
873 IF (K.EQ.12 .OR. K.EQ.14) THEN
874 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
875 LAIMN(K) = LAIMX(K)
876 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
877 VEGMN(K) = VEG(K)
878 !-- SEASONALLY VARY Zo FOR MODIS DryCrop (k=12)
879 IF (K .EQ. 12) THEN
880 Z0(K) = ZNOTC
881 !-- CrGrM (k=14) USE AVG WITH GRASS AND FOREST
882 ELSE IF (K .EQ.14) THEN
883 Z0(K) = 0.5 * (ZNOTC + Z00(K))
884 ENDIF
885 ENDIF
886 ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN
887 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
888 IF (K.EQ.20 .OR. K.EQ.43 .OR. K.EQ.45) THEN
889 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
890 LAIMN(K) = LAIMX(K)
891 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
892 VEGMN(K) = VEG(K)
893 !-- SEASONALLY VARY Zo FOR DryCrop (k=20) OR Irigated Crop (k=43) OR Mix Crop (k=4)
894 IF (K.EQ.20 .OR. K.EQ.43) THEN
895 Z0(K) = ZNOTC
896 !-- CrNatM (k=45) USE AVG WITH GRASS AND FOREST
897 ELSE IF (K.EQ.45) THEN
898 Z0(K) = 0.5 * (ZNOTC + Z00(K))
899 ENDIF
900 ENDIF
901 ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN
902 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
903 IF (K.EQ.12 .OR. K.EQ.14 .OR. K.EQ.38) THEN
904 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
905 LAIMN(K) = LAIMX(K)
906 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
907 VEGMN(K) = VEG(K)
908 !-- SEASONALLY VARY Zo FOR Crop (k=12 for MODIS or 38 for NLCD)
909 IF (K.EQ.12 .OR. K.EQ.38) THEN
910 Z0(K) = ZNOTC
911 !-- CrNatM (k=14) USE AVG WITH GRASS AND FOREST
912 ELSE IF (K.EQ.14) THEN
913 Z0(K) = 0.5 * (ZNOTC + Z00(K))
914 ENDIF
915 ENDIF
916 ELSE IF (LAND_USE_TYPE == 'USGS') THEN
917 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
918 IF (K .GE. 2 .AND. K .LE. 6) THEN
919 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
920 LAIMN(K) = LAIMX(K)
921 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
922 VEGMN(K) = VEG(K)
923 !-- SEASONALLY VARY Zo FOR DryCrop (k=2) OR Irigated Crop (k=3) OR Mix Crop (k=4)
924 IF (K .GE. 2 .AND. K .LE. 4) THEN
925 Z0(K) = ZNOTC
926 !-- CrGrM (k=5) or CrWdM (k=6) USE AVG WITH GRASS AND FOREST
927 ELSE IF (K .GE.5 .AND. K .LE. 6) THEN
928 Z0(K) = 0.5 * (ZNOTC + Z00(K))
929 ENDIF
930 ENDIF
931 ELSE IF (LAND_USE_TYPE == 'USGS28') THEN
932 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
933 IF (K .GE. 2 .AND. K .LE. 6) THEN
934 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
935 LAIMN(K) = LAIMX(K)
936 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
937 VEGMN(K) = VEG(K)
938 !-- SEASONALLY VARY Zo FOR DryCrop (k=2) OR Irigated Crop (k=3) OR Mix Crop (k=4)
939 IF (K .GE. 2 .AND. K .LE. 4) THEN
940 Z0(K) = ZNOTC
941 !-- CrGrM (k=5) or CrWdM (k=6) USE AVG WITH GRASS AND FOREST
942 ELSE IF (K .GE.5 .AND. K .LE. 6) THEN
943 Z0(K) = 0.5 * (ZNOTC + Z00(K))
944 ENDIF
945 ENDIF
947 END IF
949 ENDDO
951 ENDIF !-- IF cell is vegetation
953 !-------------------------------------
954 !-- LOOP THROUGH LANDUSE Fraction and compute totals
955 DO K = 1, NLCAT
956 FAREA = LANDUSEF(I,K,J)
957 SUMLAI = SUMLAI + LAIMX(K) * FAREA
958 SUMLMN = SUMLMN + LAIMN(K) * FAREA
959 ALAI = ALAI + FAREA
960 SUMRSI = SUMRSI + FAREA * LAIMX(K) / RSMIN(K)
961 SUMLZ0 = SUMLZ0 + FAREA * ALOG(Z0(K))
962 SUMVEG = SUMVEG + FAREA * VEG(K)
963 SUMVMN = SUMVMN + FAREA * VEGMN(K)
964 SUMSNUP = SUMSNUP+ FAREA * SNUP(K)
965 SUMALB = SUMALB + FAREA * ALB(K)
966 SUMSNOALB= SUMSNOALB + FAREA * SNOALB(K)
967 ENDDO
969 FWAT = LANDUSEF(I,KWAT,J)
970 !-- CHECK FOR WATER
971 IF (FWAT .GE. 0.50) THEN ! Changed WRFV3.7
972 XLAI(I,J) = LAIMX(KWAT)
973 XLAIMN(I,J) = LAIMN(KWAT)
974 RSTMIN(I,J) = RSMIN(KWAT)
975 ZNT(I,J) = Z0(KWAT)
976 XVEG(I,J) = VEG(KWAT)
977 XVEGMN(I,J) = VEGMN(KWAT)
978 XSNUP(I,J) = SNUP(KWAT)
979 XALB(I,J) = ALB(KWAT)
980 XSNOALB(I,J)= SNOALB(KWAT)
981 ELSE
982 IF (FWAT .GT. 0.10) THEN
983 ALAI = ALAI - FWAT
984 SUMLZ0 = SUMLZ0 - FWAT * ALOG(Z0(KWAT))
985 ENDIF
986 XLAI(I,J) = SUMLAI / ALAI
987 XLAIMN(I,J) = SUMLMN / ALAI
988 RSTMIN(I,J) = SUMLAI / SUMRSI
989 ZNT(I,J) = EXP(SUMLZ0/ALAI)
990 XVEG(I,J) = SUMVEG / ALAI
991 XVEGMN(I,J) = SUMVMN / ALAI
992 XSNUP(I,J) = SUMSNUP
993 XALB(I,J) = SUMALB
994 XSNOALB(I,J)= SUMSNOALB
995 ENDIF
997 IF (FWAT .GT. 0.50) THEN
998 ZNT(I,J) = Z0(KWAT)
999 XALB(I,J) = ALB(KWAT)
1000 XSNOALB(I,J)= SNOALB(KWAT)
1001 ENDIF
1003 !-- Compute wetlands fraction for proper MMLUIN data set
1004 !-- Note: if LU categories change, these hard coded indicies must be updated
1005 IF (LAND_USE_TYPE == 'USGS') THEN
1006 WETFRA(I,J) = LANDUSEF(I,17,J)+LANDUSEF(I,18,J)
1007 ELSE IF (LAND_USE_TYPE == 'USGS28') THEN
1008 WETFRA(I,J) = LANDUSEF(I,17,J)+LANDUSEF(I,18,J)
1009 ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN
1010 WETFRA(I,J) = LANDUSEF(I,22,J)+LANDUSEF(I,23,J)+LANDUSEF(I,27,J)+LANDUSEF(I,28,J)+LANDUSEF(I,42,J)
1011 ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN
1012 WETFRA(I,J) = LANDUSEF(I,39,J)+LANDUSEF(I,40,J)+LANDUSEF(I,11,J)
1013 ELSE IF (LAND_USE_TYPE == 'MODIS') THEN
1014 WETFRA(I,J) = LANDUSEF(I,11,J)
1015 END IF
1017 ZNT(I,J) = ZNT(I,J) * 0.01 !CONVERT TO M
1018 XVEG(I,J) = XVEG(I,J) * 0.01 !CONVERT TO FRAC
1019 XVEGMN(I,J) = XVEGMN(I,J) * 0.01
1020 XLAND(I,J) = 1.0 + FWAT
1021 XALB(I,J) = XALB(I,J) * 0.01
1022 XSNOALB(I,J)= XSNOALB(I,J) * 0.01
1024 !-------Adjustment according to CANFRA and IMPERV fo NLCD40 only -----------
1025 FIMP = IMPERV(I,J) * 0.01
1026 FCAN = CANFRA(I,J) * 0.01
1027 IF (LAND_USE_TYPE == 'NLCD40') THEN
1028 XVEG(I,J) = MIN(XVEG(I,J),1.0-FIMP)
1029 XVEGMN(I,J) = MIN(XVEGMN(I,J),1.0-FIMP)
1030 XVEG(I,J) = MAX(XVEG(I,J),FCAN)
1031 XVEGMN(I,J) = MAX(XVEGMN(I,J),FCAN)
1033 FORFRA = LANDUSEF(I,39,J)+LANDUSEF(I,30,J)+LANDUSEF(I,29,J)+LANDUSEF(I,28,J)
1034 EXTFOR = FCAN - FORFRA
1035 IF (EXTFOR.GE.0.01) THEN
1036 XLAI(I,J) = LAIMX(30) * EXTFOR + XLAI(I,J) * (1-EXTFOR)
1037 XLAIMN(I,J) = LAIMN(30) * EXTFOR + XLAIMN(I,J) * (1-EXTFOR)
1038 ENDIF
1039 ENDIF
1040 !--------------------------------------------------------------------------
1042 ENDDO ! END LOOP THROUGH GRID CELLS
1043 ENDDO ! END LOOP THROUGH GRID CELLS
1044 !--------------------------------------------------------------------
1046 END SUBROUTINE vegeland
1048 !------------------------------------------------------------------------------
1049 !------------------------------------------------------------------------------
1050 SUBROUTINE SURFPX(DTPBL, IFLAND, ISNOW, NUDGEX, XICE1, SOLDN, GSW, & !in
1051 LWDN, EMISSI, Z1, MOL, ZNT, UST, PSURF, DENS1, & !in
1052 QV1, QSS, TA1, THETA1, PRECIP, CPAIR, PSIH, & !in
1053 RH2OBS, T2OBS, VEGFRC, ISTI,LAI,IMPERV,CANFRA,BETAP, & !in
1054 RSTMIN, HC_SNOW, SNOW_FRA, WETFRA, WWLT, WFC, & !in
1055 WRES, CGSAT, WSAT, B, C1SAT, C2R, AS, JP, C3, DS1, & !in
1056 DS2, QST12, & !in
1057 RADNET, GRDFLX, HFX, QFX, LH, EG, ER, ETR, & !out
1058 QST, CAPG, RS, RA, TG, T2, WG, W2, WR, & !out
1059 TA2, QA2, LAND_USE_TYPE, I, J ) !out
1061 !------------------------------------------------------------------------------
1062 !
1063 ! FUNCTION:
1064 ! THIS SUBROUTINE COMPUTES SOIL MOISTURE AND TEMPERATURE TENDENCIES
1065 ! BY SOLVING THE PROGNOSTIC EQUATIONS IN PX95.
1066 !
1067 ! SUBROUTINES CALLED:
1068 ! SUB. QFLUX compute the soil and canopy evaporation, and transpiration
1069 ! SUB. SMASS compute nudging coefficients for soil moisture and temp nudging
1070 !
1071 ! ARGUMENTS:
1072 ! DTPBL: TIME STEP OF THE MINOR LOOP FOR THE LAND-SURFACE/PBL MODEL
1073 ! IFLAND: INDEX WHICH INDICATES THE TYPE OF SURFACE,=1,LAND;=2,SEA
1074 ! ISNOW: SNOW (=1) OR NOT (=0)
1075 ! NUDGE: SWITCH FOR SOIL MOISTURE NUDGING
1076 ! SOLDN: SHORT-WAVE RADIATION
1077 ! LWDN: LONG-WAVE RADIATION
1078 ! EMISSI: EMISSIVITY
1079 ! Z1: HEIGHT OF THE LOWEST HALF LAYER
1080 ! MOL: MONIN-OBUKOV LENGH (M)
1081 ! ZNT: ROUGHNESS LENGTH (M)
1082 ! UST: FRICTION VELOCITY (M/S)
1083 ! TST: Turbulent moisture scale
1084 ! RA: AERODYNAMIC RESISTENCE
1085 ! PSURF: P AT SURFACE (CB)
1086 ! DENS1: AIR DENSITY AT THE FIRST HALF LAYER
1087 ! QV1: Air humidity at first half layer
1088 ! QSS: Saturation mixing ratio at ground
1089 ! TA1: Air temperature at first half layer
1090 ! THETA1: Potential temperature at first half layer
1091 ! PRECIP: Precipitation rate in m/s
1092 ! STBOLT: STEFAN BOLTZMANN'S CONSTANT
1093 ! KARMAN: VON KARMAN CONSTANT
1094 ! CPAIR: Specific heat of moist air (M^2 S^-2 K^-1)
1095 ! TAUINV: 1/1DAY(SEC)
1096 ! VEGFRC: Vegetation coverage
1097 ! ISTI: soil type
1098 ! LAI: Leaf area index
1099 ! IMPERV: Percentage of IMPERVIOUS Fraction
1100 ! CANFRA: Percentage of Canopy/Tree Fraction
1101 ! BETAP: Coefficient for bare soil evaporation
1102 ! WETFRA: Fraction of Wetlands area
1103 ! THZ1OB: Observed TEMP FROM ANAL OF OBS TEMPS AT SCREEN HT
1104 ! RHOBS: Observed relative humidity at SCREEN HT
1105 ! RSTMIN Minimum Stomatol resistence
1106 !... Outputs from SURFPX
1107 ! RADNET: Net radiation
1108 ! HFX: SENSIBLE HEAT FLUX (W / M^2)
1109 ! QFX: TOTAL EVAP FLUX (KG / M^2 S)
1110 ! GRDFLX: Ground heat flux (W / M^2)
1111 ! QST: Turbulent moisture scale
1112 ! CAPG: THERMAL CAPACITY OF GROUND SLAB (J/M^2/K)
1113 ! RS: Surface resistence
1114 ! RA: Surface aerodynamic resistence
1115 ! EG: evaporation from ground (bare soil)
1116 ! ER: evaporation from canopy
1117 ! ETR: transpiration from vegetation
1118 ! TA2: diagnostic 2-m temperature from surface layer and lsm
1119 !
1120 !... Updated variables in this subroutine
1121 ! TG: Soil temperature at first soil layer
1122 ! T2: Soil temperature in root zone
1123 ! WG: Soil moisture at first soil layer
1124 ! W2: Soil moisture at root zone
1125 ! WR: Canopy water content in m
1126 !
1127 ! REVISION HISTORY:
1128 ! AX 2/2005 - developed WRF version based on the MM5 PX LSM
1129 !
1130 !------------------------------------------------------------------------------
1131 !------------------------------------------------------------------------------
1132 IMPLICIT NONE
1134 !.......Arguments
1136 !.. Integer
1137 INTEGER , INTENT(IN) :: ISTI, NUDGEX, I, J
1139 !... Real
1140 REAL , INTENT(IN) :: DTPBL, DS1, DS2
1141 REAL , INTENT(IN) :: IFLAND, ISNOW, XICE1
1142 REAL , INTENT(IN) :: SOLDN, GSW, LWDN, EMISSI, Z1
1143 REAL , INTENT(IN) :: ZNT
1144 REAL , INTENT(IN) :: PSURF, DENS1, QV1, QSS, TA1, THETA1, PRECIP
1145 REAL , INTENT(IN) :: CPAIR
1146 REAL , INTENT(IN) :: VEGFRC, LAI, IMPERV, CANFRA
1147 REAL , INTENT(IN) :: RSTMIN, HC_SNOW, SNOW_FRA, WETFRA
1148 REAL , INTENT(IN) :: WWLT, WFC, WRES, CGSAT, WSAT, B, C1SAT, C2R, AS, JP, C3
1149 REAL , INTENT(IN) :: RH2OBS,T2OBS
1150 REAL , INTENT(IN) :: QST12
1152 REAL , INTENT(OUT) :: RADNET, EG, ER, ETR
1153 REAL , INTENT(OUT) :: QST, CAPG, RS, TA2, QA2
1155 REAL , INTENT(INOUT) :: TG, T2, WG, W2, WR, UST, RA, BETAP
1156 REAL , INTENT(INOUT) :: GRDFLX, QFX, HFX, LH, PSIH, MOL
1158 CHARACTER (LEN = 6), INTENT(IN) :: LAND_USE_TYPE ! Check Limei Ran, wrf38 set to 5
1160 !... Local Variables
1162 !... Real
1163 REAL :: HF, LV, CQ4, WETSAT, SM2
1164 REAL :: RAH, RAW, ET, W2CG, CG, CT, SOILFLX, CPOT, THETAG
1165 REAL :: ZOL, ZOBOL, ZNTOL, Y, Y0, PSIH15, YNT
1166 REAL :: WGNUDG, W2NUDG, ALPH1, ALPH2, BET1, BET2, T1P5
1167 REAL :: CQ1, CQ2, CQ3, COEFFNP1, COEFFN, TSNEW, TSHLF, T2NEW
1168 REAL :: ROFF, WRMAX, PC, DWR, PNET, TENDWR, WRNEW
1169 REAL :: COF1, CFNP1WR, CFNWR, PG, FASS
1170 REAL :: TENDW2, W2NEW, W2HLF, W2REL, C1, C2, WEQ, CFNP1, CFN, WGNEW
1171 REAL :: ALN10, TMP1, TMP2, TMP3, AA, AB, TST, RBH, CTVEG
1172 REAL :: QST1,PHIH,PSIOB
1173 REAL :: T2NUD, T2NUDF
1174 REAL :: VAPPRS, QSBT, RH2MOD, IMF, VEGF, SOILF
1175 REAL :: RSOIL, LDRY, DP ! Soil model updates - JEP 12/14
1176 REAL :: C1MAX,ZZA,ZZB,ZDEL,ZLY,ZA,ZB,ZY2
1177 REAL :: Rinc, Hcan ! JEP 2020
1179 !... Parameters
1180 REAL :: ZOBS, GAMAH, BETAH, SIGF, BH, CT_SNOW, CT_IMPERV
1182 REAL, PARAMETER :: CV = 1.2E-5 ! K-M2/J Note: Update from 8E-6 10/14 Jon Pleim
1184 PARAMETER (ZOBS = 1.5) ! height for observed screen temp., (m)
1185 PARAMETER (BH = 15.7)
1186 PARAMETER (GAMAH = 16. ) !11.6)
1187 PARAMETER (BETAH = 5.0 ) !8.21)
1188 PARAMETER (SIGF = 0.5) ! rain interception see LSM (can be 0-1)
1189 REAL, PARAMETER :: DWAT = 0.2178 ! [cm^2 / s] at 273.15K
1190 !--------------------------------------------------------------------
1191 ! OLD PX legacy value from MM5 ... unknown origin PARAMETER (CT_SNOW = 5.54E-5)
1192 ! New value of CT_SNOW calibrated using multilayer soil model where csnow=6.9E5 J/(m3 K)
1193 ! from NCAR (WRFV3.2 -WRFV3.6.1) CSM PARAMETER (CT_SNOW = 2.0E-5)
1194 PARAMETER (CT_SNOW = 2.0E-5)
1196 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1197 ! CS for concrete/asphalt/road material (http://www.engineeringtoolbox.com/specific-heat-capacity-d_391.html)
1198 ! cs_imperv = 920 J kg-1 K-1
1199 ! CS for asphalt and concrete from 0.1785 to 0.20 cal g-1 K-1 (http://pages.towson.edu/morgan/files/Impervious.pdf)
1200 ! the values above translate to ~750 to 850 J kg-1 K-1
1201 !
1202 ! CAPG used for WRF urban physics ranges from 1.0E6 for roof and building walls to 1.4E6 J m-3 K-1 for roads/urban ground
1203 ! Using these values to back out CS along with 0.15 m thickness 1.4E6 J m-3 K-1 * 0.15 m = 2.1E5 J m-2 K-1
1204 ! inverse of the value above gives CT_IMPERV value of 1/0.000021 = 4.762E-6 K m2 J-1
1206 ! The middle value will be used for now. 850 J kg-1 K-1. This needs to be converted from J/K per kg to area using
1207 ! the approxiate concrete/asphalt density and layer thickness or represenative thickness. For starters (12/2011)
1208 ! well not use the PX first layer thickness, but representative thickness of most roads/parkinglots/buildings.
1209 ! for now welll use 6 inches or about 15 cm or 0.15 m. Density of concrete (normal) from
1210 ! Dorf, Richard. Engineering Handbook. New York: CRC Press, 1996. is ~2400 kg m-3.
1211 ! Using these values 850 J kg-1 K-1 * 2400 kg m-3 * 0.15 m = 3.06E5 J m-2 K-1 or in CT form (inverse) 3.268E-6 K m2 J-1
1212 ! If you look at the range of possible values considering density differences of concrete/asphalt/etc
1213 ! Values can range from 2.5 to 6.0 E-6
1214 PARAMETER (CT_IMPERV = 3.268E-6)
1215 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1217 ALN10 = ALOG(10.0)
1218 RADNET = SOLDN - (EMISSI *(STBOLT *TG **4 - LWDN)) ! NET RADIATION
1219 !--------------------------------------------------------------------
1220 CPOT= (100.0 / PSURF) ** ROVCP ! rcp is global constant(module_model_constants)
1221 THETAG = TG * CPOT
1223 ZOL = Z1/MOL
1224 ZOBOL = ZOBS/MOL
1225 ZNTOL = ZNT/MOL
1227 !-----------------------------------------------------------------------------------------
1228 IF (MOL .LT. 0.0) THEN
1229 Y = ( 1.0 - GAMAH * ZOL ) ** 0.5
1230 Y0 = ( 1.0 - GAMAH * ZOBOL ) ** 0.5
1231 YNT = ( 1.0 - GAMAH * ZNTOL ) ** 0.5
1232 PSIH15 = 2.0 * ALOG((Y + 1.0) / (Y0 + 1.0))
1233 PSIH = 2.0 * ALOG((Y + 1.0) / (YNT + 1.0))
1234 PSIOB = 2.0 * ALOG((Y0 + 1.0) / (YNT + 1.0))
1235 PHIH = 1.0 / Y
1236 ELSE
1237 IF((ZOL-ZNTOL).LE.1.0) THEN
1238 PSIH = -BETAH*(ZOL-ZNTOL)
1239 ELSE
1240 PSIH = 1.-BETAH-(ZOL-ZNTOL)
1241 ENDIF
1242 IF ((ZOBOL - ZNTOL) .LE. 1.0) THEN
1243 PSIOB = -BETAH * (ZOBOL - ZNTOL)
1244 ELSE
1245 PSIOB = 1.0 - BETAH - (ZOBOL - ZNTOL)
1246 ENDIF
1247 PSIH15 = PSIH - PSIOB
1248 IF(ZOL.LE.1.0) THEN
1249 PHIH = 1.0 + BETAH * ZOL
1250 ELSE
1251 PHIH = BETAH + ZOL
1252 ENDIF
1253 ENDIF
1254 !-----------------------------------------------------------------------------------------
1255 !-- ADD RA AND RB FOR HEAT AND MOISTURE
1256 !... RB FOR HEAT = 5 /UST
1257 !... RB FOR WATER VAPOR = 5*(0.599/0.709)^2/3 /UST = 4.503/UST
1258 RA=PR0* ( ALOG(Z1/ZNT) - PSIH )/(KARMAN*UST)
1259 RAH = RA + 5.0 / UST
1260 RAW = RA + 4.503 / UST
1261 IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN
1262 LDRY = 1.75*DS1*(EXP((1.-WG/WSAT)**5)-1.)/1.718 ! 1.75 cm is the layer thickness used by S&Z09
1263 DP = DWAT*1.E-4 * WSAT**2 * (1.-WRES/WSAT)**(2.+3./B)
1264 RSOIL=LDRY/DP
1266 ELSE
1267 RSOIL = 0.0
1268 ENDIF
1269 !--------------------------------------------------------------------
1270 ! Compute soil moisture layer 2 that considers fraction of saturated
1271 ! wetlands. If 100% of cell is wetland, soil moisture can be no lower
1272 ! than full soil saturation. If half wetland, no less than half saturated
1273 IF (IFLAND .LT. 1.5 ) THEN
1274 WETSAT = 1.00 * WSAT ! Wetlands soil moisture
1275 SM2 = (WETFRA * WETSAT) !
1276 W2 = AMAX1(SM2, W2) ! In case that W2 > Field capacity (heavy precip), use wetter W2
1277 ENDIF
1279 !--------------------------------------------------------------------
1280 ! In-canopy resistance - Erisman et al (1994)
1281 Hcan = ZNT * 10.0
1282 Rinc = 14.0 * LAI* Hcan / UST
1284 !--------------------------------------------------------------------
1285 !-- COMPUTE MOISTURE FLUX
1286 CALL QFLUX( DENS1, QV1, TA1, SOLDN, RAW, QSS, &
1287 VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, &
1288 WG, W2, WR, &
1289 RSTMIN, WWLT, WFC, RSOIL, RINC, &
1290 EG, ER, ETR, CQ4, RS, FASS)
1291 !--------------------------------------------------------------------
1293 !--------------------------------------------------------------------
1294 ! Compute Total evaporation (ET) from various modes of moisture flux
1295 ET = EG + ER + ETR
1296 QST = -ET / (DENS1 * UST)
1298 LV = 2.83E6 !-- LATENT HEAT OF SUBLIMATION AT 0C FROM STULL(1988)
1299 IF (ISNOW .LT. 0.5.AND.TG.GT.273.15) &
1300 LV = (2.501 - 0.00237 * (TG - 273.15)) * 1.E6 !-- FROM STULL(1988) in J/KG
1302 QFX = ET
1303 LH = LV * QFX
1304 !-----------------------------------------------------------------------------------------
1306 !-----------------------------------------------------------------------------------------
1307 ! Surface sensible heat flux
1308 TST = (THETA1 - THETAG ) / (UST*RAH)
1309 HF = UST * TST
1310 HFX = AMAX1(-DENS1 * CPAIR * HF, -250.0) ! using -250. from MM5
1311 !-----------------------------------------------------------------------------------------
1313 !-----------------------------------------------------------------------------------------
1314 ! Compute the diagnosed 2m Q and T consistent with PX LSM
1315 QST1 = 0.5*(QST+QST12/PHIH)
1316 TA2 = (THETAG + TST * (PR0 / KARMAN * (ALOG(ZOBS / ZNT) - PSIOB)+5.))/CPOT
1317 QA2 = QV1 - QST1 * PR0/ KARMAN * (ALOG(Z1 / ZOBS) - PSIH15)
1319 IF (QA2.LE.0.0) QA2 = QV1
1321 !-- Relative humidity
1322 VAPPRS = SVP1 * EXP(SVP2 * (TA2 - SVPT0) / (TA2 - SVP3))
1323 QSBT = EP_2 * VAPPRS / (PSURF - VAPPRS)
1324 RH2MOD = QA2 / QSBT
1325 !-----------------------------------------------------------------------------------------
1326 IF (IFLAND .LT. 1.5 ) THEN
1327 W2CG = AMAX1(W2,WWLT)
1328 CG = CGSAT * 1.0E-6 * (WSAT/ W2CG) ** &
1329 (0.5 * B / ALN10)
1330 ! IMPERVIOUS weighting scheme -- Subtract highly accurate impervious fraction from cell
1331 ! remainder is split between ground and vegetation. CT is a weighted fractional average.
1332 ! Snow CT is then applied for final heat capacity
1333 IMF = AMAX1(0.0,IMPERV/100.0)
1334 VEGF = (1.0 - IMF) * VEGFRC
1335 SOILF= (1.0 - IMF) * (1.0 - VEGFRC)
1336 CT = 1./( IMF/CT_IMPERV + VEGF/CV + SOILF/CG)
1337 CT = 1./(SNOW_FRA/CT_SNOW + (1-SNOW_FRA)/CT)
1338 CAPG = 1.0/CT
1340 SOILFLX = 2.0 * PI * TAUINV * (TG - T2)
1341 GRDFLX = SOILFLX / CT
1342 ENDIF
1343 !-----------------------------------------------------------------------------------------
1345 !--------------------------------------------------------------------
1346 !-- ASSIMILATION --- COMPUTE SOIL MOISTURE NUDGING FROM TA2 and RH2
1347 !-------COMPUTE ASSIMILATION COEFFICIENTS FOR ALL I
1348 IF (IFLAND .LT. 1.5) THEN
1349 IF (NUDGEX .EQ. 0) THEN !-- NO NUDGING CASE
1350 WGNUDG = 0.0
1351 W2NUDG = 0.0
1352 T2NUD = 0.0
1353 ELSE !-- NUDGING CASE
1355 CALL SMASS (ISTI, FASS, SOLDN, VEGFRC, RA, WWLT, WFC, &
1356 ALPH1, ALPH2, BET1, BET2, T2NUDF)
1358 !--COMPUTE SOIL MOISTURE NUDGING
1359 WGNUDG = ALPH1 * (T2OBS - TA2) + ALPH2 * (RH2OBS - RH2MOD) * 100 !NUDGING W2 FOR NON-VEG
1360 W2NUDG = BET1 * (T2OBS - TA2) + BET2 * (RH2OBS - RH2MOD) * 100
1361 IF (W2 .GE. WFC) W2NUDG = AMIN1(W2NUDG,0.0)
1362 IF (W2 .LE. WWLT) W2NUDG = AMAX1(W2NUDG,0.0)
1363 IF (W2 .GE. WFC) WGNUDG = AMIN1(WGNUDG,0.0)
1364 IF (W2 .LE. WWLT) WGNUDG = AMAX1(WGNUDG,0.0)
1365 T2NUD = T2NUDF * (T2OBS - TA2)
1366 ENDIF
1367 ENDIF
1368 !-----------------------------------------------------------------------------------------
1370 !-----------------------------------------------------------------------------------------
1371 !-- Compute new values for TS,T2,WG,W2 and WR. No change over ice or water (XLAND > 1)
1372 IF (IFLAND .LT. 1.5) THEN
1373 !-- SOLVE BY CRANK-NIC -- TENDTS=CT*(RADNET-HFX-QFX)-SOILFLX
1374 !-- Calculate the coefficients for implicit calculation of TG
1375 CQ1 = (1.0 - 0.622 * LV * CRANKP / (r_d * TG)) * QSS
1376 CQ2 = 0.622 * LV * QSS * CRANKP / (r_d * TG * TG)
1377 CQ3 = DENS1 * (1.0 - VEGFRC) / (RAW + RSOIL)
1378 COEFFNP1 = 1.0 + DTPBL * CRANKP * (4.0 * EMISSI * STBOLT * TG ** 3 &
1379 * CT + DENS1 * CPAIR / RAH * CPOT * CT + 2.0 * PI &
1380 * TAUINV ) + DTPBL * (CT * LV * CQ2 * (CQ3 + CQ4))
1381 COEFFN = CT * (GSW + EMISSI * (STBOLT * (4.0 * CRANKP - 1.0) &
1382 * TG*TG*TG*TG + LWDN) & !NET RAD
1383 + DENS1 * CPAIR / RAH * (THETA1 - (1.0 - CRANKP) * THETAG) &
1384 - LV * (CQ3 * (CQ1 - QV1) + CQ4 * (CQ1 - QV1))) & !SFC HEAT FLUX
1385 - 2.0 * PI * TAUINV * ((1.0 - CRANKP) * TG - T2) !SOIL FLUX
1386 TSNEW = (TG + DTPBL * COEFFN) / COEFFNP1
1387 !-- FOR SNOW COVERED SURFACE TEMPERATURE IS NOT MORE THAN ZERO
1388 IF (XICE1 .GT. 0.5) TSNEW = AMIN1(TSNEW,273.15) ! Re-added Jan 2010 to keep ice surface at or below freezing (J. Pleim)
1389 TSHLF = 0.5 * ( TSNEW + TG)
1390 T2NEW = (T2 + DTPBL * TAUINV * T2TFAC * (TSHLF - (1 - CRANKP) * T2) &
1391 + DTPBL*T2NUD) & ! Added deep temperature nudging
1392 / (1.0 + DTPBL * TAUINV * T2TFAC * CRANKP)
1393 !-- REPLACE OLD with NEW Value
1394 TG = TSNEW
1395 T2 = T2NEW
1396 ENDIF
1397 !-----------------------------------------------------------------------------------------
1399 !-----------------------------------------------------------------------------------------
1400 ! Compute new subsurface soil and canopy moisture values DENS1. No change required over ocean.
1401 IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN
1402 !-- Compute WR
1403 ROFF = 0.0
1404 WRMAX = 0.2E-3 * VEGFRC * LAI ! max. WRMAX IN m
1406 IF(WRMAX.GT.0.0) THEN
1407 !-- PC is precip. intercepted by veg.(M/S)
1408 PC = VEGFRC * SIGF * PRECIP
1409 DWR = (WRMAX - WR) / DTPBL !the tendency to reach max.
1410 PNET = PC - ER/ DENW ! residual of precip. and evap.
1411 IF (PNET .GT. DWR) THEN
1412 ROFF = PNET - DWR
1413 PC = PC - ROFF
1414 ENDIF
1415 IF (QSS .LT. QV1) THEN
1416 TENDWR = PC - ER / DENW
1417 WRNEW = WR + DTPBL * TENDWR
1418 ELSE
1419 COF1 = DENS1 / DENW * VEGFRC * (QSS - QV1) / RAW
1420 !-- using delta=wr/wrmax
1421 CFNP1WR = 1.0 + DTPBL * COF1 * CRANKP / WRMAX
1422 CFNWR = PC - COF1 * (1.0 - CRANKP) * WR / WRMAX
1423 WRNEW = (WR + DTPBL * CFNWR) / CFNP1WR
1424 ENDIF
1425 ELSE
1426 PC=0.0
1427 WRNEW=0.0
1428 ENDIF
1429 !---------------------------------------------
1430 !-- Compute W2
1431 PG = DENW * (PRECIP - PC) ! PG is precip. reaching soil (PC already including ROFF)
1432 TENDW2 = 1.0 / (DENW * DS2) * (PG - EG - ETR) &
1433 - C3/DS2 * TAUINV * AMAX1(0.0,(W2 - WFC)) &
1434 + (W2NUDG + WGNUDG) / DS2 ! NUDGING
1435 W2NEW = W2 + DTPBL * TENDW2
1436 W2NEW = AMIN1(W2NEW,WSAT)
1437 W2NEW = AMAX1(W2NEW,WRES) !0.05) !Limei 08/02/2017
1438 W2HLF = 0.5 * (W2 + W2NEW)
1439 !.. new values
1440 W2 = W2NEW
1441 WR = AMIN1(WRMAX,WRNEW)
1442 ENDIF
1443 !-----------------------------------------------------------------------------------------
1445 !-----------------------------------------------------------------------------------------
1446 ! Compute new surface soil moisture values (WR).
1447 IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN ! over ocean no change to wg w2,wr
1448 !-- FOR SNOW COVERED SURFACE, ASSUME SURFACE IS SATURATED AND
1449 ! WG AND W2 ARE NOT CHANGED
1450 IF (ISNOW .GT.0.5) THEN
1451 WG = WSAT
1452 ELSE
1453 W2REL = W2HLF / WSAT
1454 IF (WG .GT. WWLT) THEN
1455 C1 = DS1*C1SAT * (WSAT / WG) ** (0.5 * B + 1.0)
1456 ELSE ! revise C1 for wg < wilting point Noilhan & Mahfouf (1996)
1457 ZY2 = C1SAT * (WSAT / WWLT) ** (0.5 * B + 1.0)
1458 C1MAX = (1.19*WWLT - 5.09)*TG - 146.*WWLT + 1786.
1459 C1MAX = MAX(MAX(C1MAX,ZY2),10.)
1460 !* Giard-Bazile formulation (resolution of a second order equation)
1461 !
1462 ZLY = LOG( C1MAX/10.)
1463 ZZA = - LOG( ZY2 /10.)
1464 ZZB = 2. * WWLT * ZLY
1465 ZDEL = 4. * (ZLY+ZZA) * ZLY * WWLT**2
1466 ZA = (-ZZB+SQRT(ZDEL)) / (2.*ZZA)
1467 ZB = ZA**2 / ZLY
1468 C1 = DS1*C1MAX * EXP(-(WG-ZA)**2/ZB)
1469 ENDIF
1471 C2 = C2R * W2HLF / (WSAT - W2HLF + 1.E-11)
1472 IF (W2HLF .GE. WSAT) THEN
1473 WEQ = WSAT
1474 ELSE
1475 WEQ = W2HLF - AS * WSAT * W2REL ** JP * &
1476 (1.0 - W2REL ** (8.0 * JP))
1477 ENDIF
1479 !.... The diffusion method in Sakaguchi and Zeng [2009] (JGR-Atmos.)
1480 CFNP1 = 1.0 + DTPBL * C2 * TAUINV * CRANKP
1481 CFN = C1 / (DENW * DS1) * (PG - EG) - C2 * TAUINV * &
1482 ((1.0 - CRANKP) * WG - WEQ) + WGNUDG/ DS1
1484 WGNEW = AMAX1((WG + DTPBL * CFN) / CFNP1, WRES ) !0.001) ! Limei 08/02/2017
1485 !-- NEW VALUES
1486 WG = AMIN1(WGNEW,WSAT)
1488 ENDIF !endif for ISNOW
1489 ENDIF !endif for XLAND
1491 END SUBROUTINE surfpx
1492 !-------------------------------------------------------------------------
1493 !-------------------------------------------------------------------------
1495 !-------------------------------------------------------------------------
1496 !-------------------------------------------------------------------------
1497 SUBROUTINE QFLUX (DENS1, QV1, TA1, RG, RAW, QSS, & ! in
1498 VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, & ! in
1499 WG, W2, WR, & ! in
1500 RSTMIN, WWLT, WFC, RSOIL, RINC, & ! in !Soil model updates - JEP 12/14
1501 EG, ER, ETR, CQ4, RS, FASS) ! out
1503 !-------------------------------------------------------------------------
1504 !
1505 ! FUNCTION:
1506 ! THIS SUBROUTINE COMPUTES EVAPORATION FROM BARE SOIL (EG) AND FROM
1507 ! THE WET PART OF CANOPY (ER) AND TRANSPIRATION FROM THE DRY PART OF
1508 ! CANOPY (ETR).
1509 !
1510 ! REVISION HISTORY:
1511 ! A. Xiu Oct 2004 - adapted from the PX LSM in MM5 for the WRF system
1512 ! R. Gilliam Dec 2006 - Completed WRF V2.1.2 implementation
1513 !
1514 !-------------------------------------------------------------------------
1515 ! QFLUX ARGUMENT LIST:
1516 !-------------------------------------------------------------------------
1517 ! INPUTS:
1518 !-- DENS1 air density at first layer
1519 !-- QV1 air humidity at first layer
1520 !-- TA1 air temperature at first layer
1521 !-- RG shortwave radition reaching the ground
1522 !-- RAW RA+RB for moisture
1523 !-- QSS saturation mixing ratio at ground
1524 !-- VEGFRC vegetation coverage
1525 !-- ISNOW if snow on the ground
1526 !-- ISTI soil type
1527 !-- IFLAND if land (=1) or water (=2)
1528 !-- LAI leaf area index
1529 !-- BETAP
1530 !-- WG soil moisture at first soil layer
1531 !-- W2 soil moisture at root zone
1532 !-- WR Canopy water
1533 !
1534 ! OUTPUTS FROM QFLUX:
1535 !-- EG evaporation from ground (bare soil)
1536 !-- ER evaporation from canopy
1537 !-- ETR transpiration from vegetation
1538 !-- CQ4
1539 !-- RS surface resistence
1540 !-- FASS parameter for soil moisture nudging
1541 !-------------------------------------------------------------------------
1542 !-------------------------------------------------------------------------
1544 IMPLICIT NONE
1546 ! DECLARATIONS - INTEGER
1547 INTEGER , INTENT(IN) :: ISTI
1549 ! DECLARATIONS - REAL
1550 REAL , INTENT(IN) :: ISNOW, IFLAND
1551 REAL , INTENT(IN) :: DENS1, QV1, TA1, RG, RAW, QSS, &
1552 VEGFRC, LAI, &
1553 WG, W2, WR, RSTMIN
1554 REAL , INTENT(INOUT) :: BETAP, RSOIL
1555 REAL, INTENT(IN) :: WWLT, WFC, RINC
1557 REAL , INTENT(OUT) :: EG, ER, ETR, CQ4, RS, FASS
1559 !... Local Variables
1561 !... Real
1562 REAL :: WRMAX, DELTA, SIGG, RADL, RADF, W2AVAIL, W2MXAV
1563 REAL :: FTOT, F1, F2, F3, F4
1564 REAL :: FSHELT, GS, GA, FX
1565 REAL :: PAR, F1MAX
1568 !... Parameters
1569 REAL, PARAMETER :: RSMAX = 5000.0 ! s/m
1570 REAL, PARAMETER :: FTMIN = 0.0000001 ! m/s
1571 REAL, PARAMETER :: F3MIN = 0.25
1573 !
1574 !... for water surface, no canopy evaporation and transpiration
1575 ER = 0.0
1576 ETR = 0.0
1577 CQ4 = 0.0
1579 !... GROUND EVAPORATION (DEPOSITION)
1580 IF (QSS .LT. QV1) RSOIL = 0.0
1581 EG = DENS1 * (QSS - QV1) * ((1.0 - VEGFRC) / (RAW + RSOIL) + VEGFRC/(RAW + RINC +RSOIL))
1583 !... CANOPY
1584 IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1585 WRMAX = 0.2E-3 * VEGFRC * LAI ! in unit m
1586 IF (WR .LE. 0.0) THEN
1587 DELTA = 0.0
1588 ELSE
1589 DELTA = WR / WRMAX ! referred to SiB model
1590 ENDIF
1592 IF (QSS .GE. QV1) THEN
1593 SIGG = DELTA
1594 ELSE
1595 SIGG = 1.0
1596 ENDIF
1598 ER = DENS1 * VEGFRC * SIGG * (QSS - QV1) / RAW
1599 ENDIF
1600 !!---------------------------------------------------------------------
1602 !-- TRANSPIRATION
1603 !!---------------------------------------------------------------------
1604 IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1606 !!!-RADIATION
1607 IF (RSTMIN .GT. 130.0) THEN
1608 F1MAX = 1.-0.02*LAI !Echer2015 Trees
1609 ELSE
1610 F1MAX = 1.-0.07*LAI !Echer2015 crops/grass
1611 ENDIF
1612 PAR = 0.45 * RG * 4.566 ! converted from W/m2 to umoles/m2/s (1/.219) Echer2015
1613 F1 = F1MAX*(1.0-exp(-0.0017*PAR)) !Echer2015
1614 F1 = AMAX1(F1,RSTMIN / RSMAX)
1616 !!!-SOIL MOISTURE
1617 W2AVAIL = W2 - WWLT
1618 W2MXAV = WFC - WWLT
1619 F2 = 1.0 / (1.0 + EXP(-5.0 * ( W2AVAIL / W2MXAV - &
1620 (W2MXAV / 3.0 + WWLT)))) ! according JP, 9/94
1622 !-AIR TEMP
1623 !... according to Avissar (1985) and AX 7/95
1624 IF (TA1 .LE. 302.15) THEN
1625 F4 = 1.0 / (1.0 + EXP(-0.41 * (TA1 - 282.05)))
1626 ELSE
1627 F4 = 1.0 / (1.0 + EXP(0.5 * (TA1 - 314.0)))
1628 ENDIF
1630 FTOT = LAI * F1 * F2 * F4
1631 ENDIF
1633 !!---------------------------------------------------------------------
1634 IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1635 FSHELT = 1.0 ! go back to NP89
1636 GS = FTOT / (RSTMIN * FSHELT)
1637 GA = 1.0 / RAW
1638 !-- Compute humidity effect according to RH at leaf surf
1639 F3 = 0.5 * (GS - GA + SQRT(GA * GA + GA * GS * &
1640 (4.0 * QV1 / QSS - 2.0) + GS * GS)) / GS
1641 F3 = AMIN1(AMAX1(F3,F3MIN),1.0)
1642 RS = 1.0 / (GS * F3)
1644 !--- Compute Assimilation factor for soil moisture nudging - jp 12/94
1645 !-- Note that the 30 coef is to result in order 1 factor for max
1646 IF (RG .LT. 0.00001) THEN ! do not nudge during night
1647 FX = 0.0
1648 ELSE
1649 FX = 30.0 * F1 * F4 * LAI / (RSTMIN * FSHELT)
1650 ENDIF
1652 FASS = FX
1653 ETR = DENS1 * VEGFRC * (1.0 - SIGG) * (QSS - QV1) / (RAW + RS)
1654 !..... CQ4 is used for the implicit calculation of TG in SURFACE
1655 CQ4 = DENS1 * VEGFRC * ((1.0 - SIGG) / (RAW + RS) + SIGG / RAW)
1656 ENDIF
1658 END SUBROUTINE qflux
1659 !------------------------------------------------------------------------------------------
1660 !------------------------------------------------------------------------------------------
1662 !------------------------------------------------------------------------------------------
1663 !------------------------------------------------------------------------------------------
1664 SUBROUTINE SMASS (ISTI, FASS, RG, VEGFRC, RA, & !in
1665 WWLT, WFC, & !in
1666 ALPH1, ALPH2, BET1, BET2, T2NUDF ) !out
1667 !------------------------------------------------------------------------------------------
1668 !------------------------------------------------------------------------------------------
1669 IMPLICIT NONE
1670 !------------------------------------------------------------------------------------------
1671 ! SMASS COMPUTES SOIL MOISTURE NUDGING COEFFICIENTS
1672 !------------------------------------------------------------------------------------------
1673 !
1674 !.........Arguments
1675 INTEGER, PARAMETER :: NSCAT = 16 ! max. soil types
1677 INTEGER, INTENT(IN) :: ISTI
1678 REAL, INTENT(IN) :: FASS, RG, VEGFRC, RA
1679 REAL, INTENT(IN) :: WWLT, WFC
1680 REAL, INTENT(OUT) :: ALPH1, ALPH2, BET1, BET2, T2NUDF
1683 !........Local variables
1685 !... Real
1686 REAL :: FBET, FALPH, FRA, FTEXT
1688 !... Parameters
1689 REAL, PARAMETER :: A1MAX = -10.E-5, A2MAX = 1.E-5 ! m/K, m for 6hr period
1690 REAL, PARAMETER :: B1MAX = -10.E-3, B2MAX = 1.E-3 ! m/K, m (Bouttier et al 1993)
1691 REAL, PARAMETER :: TASSI = 4.6296E-5 ! 1/6hr in 1/sec
1692 REAL, PARAMETER :: RAMIN = 10.0 ! 0.1 s/cm
1693 REAL, PARAMETER :: WFCX = 0.243 ! middle of WFC range
1694 REAL, PARAMETER :: WWLTX = 0.169 ! middle of WWLT range
1696 !
1697 FBET = FASS
1698 FALPH = RG / 1370.0
1699 !--TEXTURE FACTOR NORMALIZED BY LOAM (IST=5)
1700 FRA = RAMIN / RA ! scale by aerodynamic resistance
1701 FTEXT = TASSI * (WWLT + WFC) / (WWLTX + WFCX) * FRA
1702 ! write(6,*) ' ftot, fbet=',ftot, fbet,' ftext=',ftext/tassi
1703 !
1704 ALPH1 = A1MAX * FALPH * (1.0 - VEGFRC) * FTEXT
1705 ALPH2 = A2MAX * FALPH * (1.0 - VEGFRC) * FTEXT
1706 BET1 = B1MAX * FBET * VEGFRC * FTEXT
1707 BET2 = B2MAX * FBET * VEGFRC * FTEXT
1708 T2NUDF = 1.0E-5 * ( VEGFRC*MAX((1.0 - 5.0 * FALPH),0.0) + (1-VEGFRC) ) ! T2 Nudging at night and day for non-veg frac - jp 10/30/14
1710 END SUBROUTINE smass
1711 !------------------------------------------------------------------------------------------
1712 !------------------------------------------------------------------------------------------
1714 !------------------------------------------------------------------------------------------
1715 !------------------------------------------------------------------------------------------
1717 SUBROUTINE SOILPROP (SOILCBOT,WEIGHT, ITIMESTEP, MAVAIL, & ! IN
1718 PXLSM_SMOIS_INIT, & ! IN
1719 FWSAT,FWFC,FWWLT,FCLAY,FCSAND, & ! OUT
1720 FFMSAND,FB,FCGSAT, & ! OUT
1721 FJP,FAS,FC2R,FC1SAT,FWRES,FC3,ISTI, & ! OUT
1722 WG, W2 ) ! OUT
1723 !------------------------------------------------------------------------
1724 ! SOILPROP COMPUTES SOIL PARAMETERS FOR BOTH BOTTOM AND TOP LAYERS
1725 ! USING FRACTIONAL SOIL TYPE. A HARD CODED OPTION IS AVAILIABLE
1726 ! TO COMPUTE THE SOIL PARAMETERS USING FRACTIONAL INFORMATION, OR
1727 ! TO JUST USE SOIL PARAMETERS OF THE DOMINANT SOIL TYPE
1728 !------------------------------------------------------------------------
1729 !-- SOIL PARAMETERS ARE SPECIFIED BY SOIL TYPE:
1730 ! # SOIL TYPE WSAT WFC WWLT B CGSAT JP AS C2R C1SAT WRES
1731 ! _ _________ ____ ___ ____ ____ _____ ___ ___ ___ _____ ____
1732 ! 1 SAND .395 .135 .068 4.05 3.222 4 .387 3.9 .082 0.020
1733 ! 2 LOAMY SAND .410 .150 .075 4.38 3.057 4 .404 3.7 .098 0.035
1734 ! 3 SANDY LOAM .435 .195 .114 4.90 3.560 4 .219 1.8 .132 0.041
1735 ! 4 SILT LOAM .485 .255 .179 5.30 4.418 6 .105 0.8 .153 0.015
1736 ! 5 SILT .485 .255 .179 5.30 4.418 6 .105 0.8 .153 0.015
1737 ! 6 LOAM .451 .240 .155 5.39 4.111 6 .148 0.8 .191 0.027
1738 ! 7 SND CLY LM .420 .255 .175 7.12 3.670 6 .135 0.8 .213 0.068
1739 ! 8 SLT CLY LM .477 .322 .218 7.75 3.593 8 .127 0.4 .385 0.040
1740 ! 9 CLAY LOAM .476 .325 .250 8.52 3.995 10 .084 0.6 .227 0.075
1741 ! 10 SANDY CLAY .426 .310 .219 10.40 3.058 8 .139 0.3 .421 0.109
1742 ! 11 SILTY CLAY .482 .370 .283 10.40 3.729 10 .075 0.3 .375 0.056
1743 ! 12 CLAY .482 .367 .286 11.40 3.600 12 .083 0.3 .342 0.090
1744 !------------------------------------------------------------------------
1745 !------------------------------------------------------------------------------------------
1746 !------------------------------------------------------------------------------------------
1747 IMPLICIT NONE
1749 !.........Arguments
1750 INTEGER, PARAMETER :: NSCAT = 16 ! max. soil types
1751 INTEGER, PARAMETER :: NSCATMIN = 16 ! min soil types
1753 INTEGER, INTENT(IN) :: WEIGHT, ITIMESTEP, PXLSM_SMOIS_INIT
1754 REAL, INTENT(IN) :: MAVAIL
1755 REAL, DIMENSION(1:NSCAT), INTENT(IN) :: SOILCBOT
1756 REAL, INTENT(OUT) :: FWSAT,FWFC,FWWLT,FCLAY, &
1757 FCSAND,FFMSAND,FB,FCGSAT, &
1758 FJP,FAS,FC2R,FC1SAT,FWRES,FC3
1759 REAL, INTENT(INOUT) :: W2, WG
1762 INTEGER, INTENT(OUT) :: ISTI
1763 !........Local variables
1764 CHARACTER*4, AVCLASS
1765 CHARACTER*4, DIMENSION( 1: NSCAT ) :: TEXID
1766 !... Integer
1767 INTEGER:: S
1768 !... Real
1769 REAL:: TFRACBOT, CFRAC, SUMCSND, SUMFMSND, SUMCLY, &
1770 AVS, AVCS, AVFMS, AVC, AVSLT, &
1771 SSMPOT, DSMPOT ! saturated and air-dry soil matric potential
1773 REAL, DIMENSION( 1: NSCAT ) :: WSAT, WFC, WWLT, B, CGSAT, AS, &
1774 JP, C2R, C1SAT, WRES
1776 REAL, DIMENSION( 1: NSCATMIN ) :: CSAND,FMSAND, CLAY
1777 !.......... DATA statement for SOIL PARAMETERS for the 11 soil types
1778 !...........Follow Menut et al., 2013 JGR
1779 DATA CSAND /46.0,40.0,29.0, 0.0, 0.0, &
1780 0.0,29.0, 0.0, 0.0, 0.0, &
1781 0.0, 0.0, 0.0, 0.0,46.0, &
1782 0.0/
1783 DATA FMSAND /46.0,40.0,29.0,17.0,10.0, &
1784 43.0,29.0,10.0,32.0,52.0, &
1785 6.0,22.0,43.0,43.0,46.0, &
1786 32.0/
1788 DATA CLAY / 3.0, 4.0,10.0,13.0, 5.0, &
1789 18.0,27.0,34.0,34.0,42.0, &
1790 47.0,58.0,18.0,18.0, 3.0, &
1791 34.0/
1793 DATA TEXID/'Sand','Lsan','Sloa','Sill','Silt', &
1794 'Loam','Sclo','Sicl','Cllo','Sacl', &
1795 'Sicy','Clay','Ormt','Wate','Bedr', &
1796 'Othe'/
1797 !-- Removed soil parameter lookup table data since these parameters are now computed analytically (Noilhan and Mahfouf 1996)
1799 DSMPOT = -1.0E7 ! mm air-dry water matric potential
1800 !
1801 !-------------------------------Exicutable starts here--------------------
1803 ! Compute soil characteristics by sand (coarse and fine-medium) and clay fraction
1805 CFRAC = 0.0
1806 SUMCSND = 0.0
1807 SUMFMSND = 0.0
1808 SUMCLY = 0.0
1809 TFRACBOT = 0.0
1811 DO S = 1,NSCAT
1813 TFRACBOT = TFRACBOT + SOILCBOT(S)
1814 SUMCSND = SUMCSND + CSAND(S) * SOILCBOT(S)
1815 SUMFMSND = SUMFMSND + FMSAND(S) * SOILCBOT(S)
1816 SUMCLY = SUMCLY + CLAY(S) * SOILCBOT(S)
1818 IF(SOILCBOT(S).GE.CFRAC) THEN ! Find Dominant Category and fraction
1819 ISTI=S
1820 CFRAC=SOILCBOT(S)
1821 ENDIF
1823 ENDDO
1826 IF(TFRACBOT.GT.0.001) THEN
1827 AVCS = SUMCSND / TFRACBOT
1828 AVFMS = SUMFMSND / TFRACBOT
1829 AVS = AVCS + AVFMS
1831 AVC = SUMCLY / TFRACBOT
1832 AVSLT = 100.0 - AVS - AVC
1834 IF(AVS.GT.(85.+ 0.5*AVC)) THEN
1835 AVCLASS= 'Sand'
1836 ISTI = 1
1837 ELSE IF(AVS.GT.(70.+ AVC)) THEN
1838 AVCLASS= 'Lsan'
1839 ISTI = 2
1840 ELSE IF((AVC.LT.20..AND.AVS.GT.52.) &
1841 .OR.(AVC.LE.7.5.AND.AVSLT.LT.50.)) THEN
1842 AVCLASS= 'Sloa'
1843 ISTI = 3
1844 ELSE IF(AVC.LT.35..AND.AVS.GT.45..AND.AVSLT.LT.28.) THEN
1845 AVCLASS= 'Sclo'
1846 ISTI = 7
1847 ELSE IF(AVC.GE.35..AND.AVS.GT.45.) THEN
1848 AVCLASS = 'Sacl'
1849 ISTI = 10
1850 ELSE IF(AVC.LT.27.0.AND.AVSLT.LT.50.) THEN
1851 AVCLASS= 'Loam'
1852 ISTI = 6
1853 ELSE IF(AVC.LT.12..AND.AVSLT.GT.80.) THEN
1854 AVCLASS = 'Silt'
1855 ISTI = 5
1856 ELSE IF(AVC.LT.27.) THEN
1857 AVCLASS = 'Sill'
1858 ISTI = 4
1859 ELSE IF(AVC.LT.40..AND.AVS.GT.20.) THEN
1860 AVCLASS = 'Cllo'
1861 ISTI = 9
1862 ELSE IF(AVC.LT.40.) THEN
1863 AVCLASS = 'Sicl'
1864 ISTI = 8
1865 ELSE IF(AVSLT.GE.40.) THEN
1866 AVCLASS = 'Sicy'
1867 ISTI = 11
1868 ELSE
1869 AVCLASS = 'Clay'
1870 ISTI = 12
1871 ENDIF
1872 ELSE
1873 ! set no soil to 9 - clay loam
1874 ISTI = 9
1875 AVCLASS = TEXID(ISTI)
1877 AVCS = CSAND(ISTI)
1878 AVFMS = FMSAND(ISTI)
1879 AVS = AVCS + AVFMS
1881 AVC = CLAY(ISTI)
1882 AVSLT = 100.0 - AVS - AVC
1884 ENDIF
1886 FCSAND = AVCS
1887 FFMSAND = AVFMS
1888 FCLAY = AVC
1890 ! Continous formulation of secondary soil parmeters (Noilhan and Mahfouf 1996)
1891 FWSAT = (-1.08 * AVS + 494.305) * 1.0E-3
1892 FWWLT = 37.1342E-3 * SQRT(AVC)
1893 FWFC = 89.0467E-3 * AVC**0.3496
1894 FB = 0.137 * AVC + 3.501
1895 FCGSAT= -1.557E-2 * AVS - 1.441E-2 * AVC + 4.7021
1896 FC1SAT= (5.58 * AVC + 84.88) * 1.0E-2
1897 FC2R = 13.815 * AVC**(-0.954)
1898 FC3 = 5.327 * AVC **(-1.043)
1899 FAS = 732.42E-3 * AVC **(-0.539)
1900 FJP = 0.134 * AVC + 3.4
1901 FWRES = 0.00123 * AVC - 0.00066 * AVSLT + 0.0405 !J. Pleim fitted function
1902 FWRES = AMAX1(FWRES, 0.01) !L. Ran set minimum
1904 ! Compute W2 using soil moisture availiability if pxlsm_smois_init (in namelist) is not zero
1905 IF (ITIMESTEP .EQ. 1 .AND. PXLSM_SMOIS_INIT .GT. 0) THEN
1906 WG = FWWLT + (0.5*(FWSAT+FWFC) - FWWLT) * MAVAIL
1907 W2 = FWWLT + (0.5*(FWSAT+FWFC) - FWWLT) * MAVAIL
1908 ENDIF
1910 END SUBROUTINE soilprop
1911 !------------------------------------------------------------------------------------------
1912 !------------------------------------------------------------------------------------------
1914 !------------------------------------------------------------------------------------------
1915 !------------------------------------------------------------------------------------------
1916 SUBROUTINE PXSNOW (ITIMESTEP, ASNOW, CSNOW, SNOW, &
1917 SNOWH, SNUP, &
1918 ALB, SNOALB, SHDFAC, SHDMIN, &
1919 HC_SNOW, SNOW_FRA, SNOWC, SNOWALB)
1920 !------------------------------------------------------------------------------------------
1921 !------------------------------------------------------------------------------------------
1922 ! Pleim-Xiu LSM Simple Snow model
1923 !---------------------------------------------------------------------------------------------------
1924 ! ITIMESTEP -- Model time step index
1925 ! ASNOW -- Analyzed snow water equivalent (mm)
1926 ! CSNOW -- Current snow water equivalent (mm)
1927 ! SNOW -- Snow water equivalent in model (mm)
1928 ! SNOWH -- Physical snow depth (m)
1929 ! SNUP -- Physical snow depth (landuse dependent) where when below, snow cover is fractional
1930 !
1931 ! HC_SNOW -- Heat capacity of snow as a function of depth
1932 ! SNOW_FRA-- Factional snow area
1933 ! IFSNOW -- Snow flag
1934 ! SNOWALB -- Snow albedo
1935 !---------------------------------------------------------------------------------------------------
1937 IMPLICIT NONE
1939 !--- Arguments
1940 REAL, PARAMETER :: W2SN_CONV = 10.0
1941 REAL, PARAMETER :: CS_SNOWPACK = 2092.0
1942 REAL, PARAMETER :: SALP = 2.6
1943 !-----------------------------------------------
1944 INTEGER, INTENT(IN) :: ITIMESTEP
1945 REAL, INTENT(IN) :: ASNOW, CSNOW, SNUP, ALB, SNOALB, SHDFAC, SHDMIN
1946 REAL, INTENT(INOUT) :: SNOW, SNOWH, SNOWC
1947 REAL, INTENT(OUT) :: HC_SNOW, SNOW_FRA, SNOWALB
1948 !------------------------------------------------------------------------------------
1951 !-----------------------------------------------
1952 ! Local variables
1953 !... Real
1954 REAL:: CONV_WAT2SNOW, CSNOWW, RHO_SNOWPACK, &
1955 LIQSN_RATIO, SNEQV, RSNOW
1956 !-----------------------------------------------
1958 SNEQV=ASNOW*0.001 ! Snow water in meters
1959 RHO_SNOWPACK = 450 ! Snowpack density (kg/m3), this should be computed in the future
1960 LIQSN_RATIO = DENW/RHO_SNOWPACK ! Ratio of water density and snowpack density
1962 CONV_WAT2SNOW = LIQSN_RATIO/1000 ! Conversion factor for snow liquid equiv. (mm) to snowpack (m)
1964 SNOW = ASNOW ! Set snow water to analysis value for now, implement a nudging scheme later
1965 SNOWH= SNOW * CONV_WAT2SNOW ! Conversion of snow water (mm) to snow depth (m)
1968 ! Is snow cover now present. The limit is 0.45 mm of water eqivalent or about 2 inches snow depth
1969 SNOWC = 0.0
1970 IF (SNOWH .GT. 0.005) SNOWC = 1.0
1972 HC_SNOW = RHO_SNOWPACK * CS_SNOWPACK * SNOWH
1974 IF (SNEQV < SNUP) THEN
1975 RSNOW = SNEQV / SNUP
1976 SNOW_FRA = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
1977 ELSE
1978 SNOW_FRA = 1.0
1979 END IF
1981 SNOWC = SNOW_FRA
1983 SNOWALB = ALB + SNOWC*(SNOALB-ALB)
1985 END SUBROUTINE pxsnow
1986 !------------------------------------------------------------------------------------------
1987 !------------------------------------------------------------------------------------------
1989 END MODULE module_sf_pxlsm