| blob | b2c0ce1aad21b81fb997826276052367a0551911 |
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 == 61) THEN
374 LAND_USE_TYPE = 'MODIS'
375 ELSE IF (NLCAT == 24) THEN
376 LAND_USE_TYPE = 'USGS'
377 ELSE IF (NLCAT == 28) THEN
378 LAND_USE_TYPE = 'USGS28'
379 ELSE
380 call wrf_error_fatal("Error: Unknown Land Use Category")
381 END IF
383 IF (ITIMESTEP .EQ. 1) THEN
384 CALL wrf_message( 'PX LSM will use the ' // TRIM(LAND_USE_TYPE) // ' landuse tables' )
385 PRINT *, 'The analysis interval for surface soil and temp nudging = ',ANAL_INTERVAL,'sec.'
386 ENDIF
388 !-----------------------------------------------------------------------------------
389 ! Kill WRF if user specifies soil nudging but provides no analysis interval, then provide helpful message.
390 IF (ANAL_INTERVAL .LE. 0.0 .AND. PXLSM_SOIL_NUDGE .EQ. 1) THEN
391 CALL wrf_message('PX LSM Error: The User specified analysis interval is zero or negative.')
392 CALL wrf_message('If the PX LSM is used with soil nudging (pxlsm_soil_nudge=1) a wrfsfdda_d0* file is required.')
393 CALL wrf_message('Make sure these files are present and')
394 CALL wrf_message('Check the namelist to ensure sgfdda_interval_m is set to proper sfc analysis interval')
395 STOP
396 ENDIF
397 !-----------------------------------------------------------------------------------
398 !--- Compute time relatve to old and new analysis time for timestep interpolation
399 IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN
400 DT_FDDA = ANAL_INTERVAL * 1.0 ! Convert DT of Analysis to real
401 TIME_BETWEEN_ANALYSIS = MOD(CURR_SECS,DT_FDDA)
402 IF (TIME_BETWEEN_ANALYSIS .EQ. 0.0) THEN
403 CORB = 1.0
404 CORE = 0.0
405 ELSE
406 CORE = TIME_BETWEEN_ANALYSIS / DT_FDDA
407 CORB = 1.0 - CORE
408 ENDIF
409 ENDIF
410 !-----------------------------------------------------------------------------------
411 ! Compute vegetation and land-use characteristics by land-use fraction weighting
412 ! These parameters include LAI, VEGF, ZNT, ALBEDO, SNOALB, RS, etc.
413 CALL VEGELAND(LANDUSEF, VEGFRA, SHDMIN, SHDMAX, &
414 SOILCTOP, SOILCBOT, NLCAT, NSCAT, &
415 ZNT,XLAI,XLAIMN,RSTMIN,XVEG,XVEGMN,XSNUP, &
416 XLAND, XALB,XSNOALB,WETFRA,IMPERV,CANFRA, &
417 ids,ide, jds,jde, kds,kde, &
418 ims,ime, jms,jme, kms,kme, &
419 its,ite, jts,jte, kts,kte, LAND_USE_TYPE, &
420 KWAT )
421 !-----------------------------------------------------------------------------------
424 !-----------------------------------------------------------------------------------
425 ! Main loop over individual grid cells
426 DO J = jts,jte !-- J LOOP
427 DO I = its,ite !-- I LOOP
429 IFLAND = XLAND(I,J)
431 ! Compute soil properties via weighting of fractional components
432 IF (IFLAND .LT. 1.5 ) THEN
434 !---------------------------------------------------------
435 CALL SOILPROP (SOILCBOT(I,:,J), WEIGHT, &
436 ITIMESTEP, MAVAIL(I,J), &
437 PXLSM_SMOIS_INIT, &
438 FWSAT,FWFC,FWWLT,FCLAY,FCSAND,FFMSAND, &
439 FB,FCGSAT, &
440 FJP,FAS,FC2R,FC1SAT,FWRES, FC3, ISTI, &
441 SMOIS(I,1,J), SMOIS(I,2,J) )
442 !----------------------------------------------------------
443 !----------------------------------------------------------
444 ISLTYP(I,J) = ISTI
445 ELSE
446 ISLTYP(I,J) = 14 ! STATSGO type for water
448 !-- aded for MODIS model
449 FWWLT = 0.1
450 FWFC = 1.0
451 FWSAT = 1.0
453 FCLAY = 0.0
454 FCSAND = 0.0
455 FFMSAND = 0.0
456 !-- end
458 ENDIF
460 !-- aded for MODIS model
461 WWLT_PX(I,J) = FWWLT
462 WFC_PX(I,J) = FWFC
463 WSAT_PX(I,J) = FWSAT
465 CLAY_PX(I,J) = FCLAY * 0.01 ! percent to fraction
466 CSAND_PX(I,J) = FCSAND * 0.01
467 FMSAND_PX(I,J) = FFMSAND * 0.01
468 !-- end
470 !-- Variables Sub. SURFPX needs
471 SFCPRS = PSFC(i,j) / 1000.0 ! surface pressure in cb
472 TA1 = T3D(i,1,j) ! air temperature at first layer
473 DENS1 = RHO(I,1,J) ! air density at first layer
474 QV1 = QV3D(i,1,j) ! water vapor mixing ratio at first layer
475 QV2 = QV3D(i,2,j)
476 ZLVL = 0.5 * DZ8W(i,1,j) ! thickness of lowest half level
477 ZF1 = DZ8W(i,1,j)
478 ZA2 = ZF1 + 0.5 * DZ8W(i,2,j)
480 LWDN = GLW(I,J) ! longwave radiation
481 EMISSI = EMISS(I,J) ! emissivity
482 PRECIP = MAX ( 1.0E-3*RAINBL(i,j)/DTBL,0.0) ! accumulated precip. rate during DT (=dtpbl)
483 ! convert RAINBL from mm to m for PXLSM
484 WR = 1.0E-3*CANWAT(I,J) ! convert CANWAT from mm to m for PXLSM
485 THETA1 = TH3D(i,1,j) ! potential temp at first layer
486 SNOBS = SNOW(I,J) ! Set snow cover to existing model value
487 ! this is overwritten below if snow analysis is availiable
488 ! otherwise snow cover remains constant through simulation
490 IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN
491 !-- 2 m Temp and RH for Nudging
492 T2_OLD = T2_NDG_OLD(I,J)
493 T2_NEW = T2_NDG_NEW(I,J)
494 VAPPRS = SVP1 * EXP(SVP2 * (T2_OLD - SVPT0) / ( T2_OLD - SVP3))
495 QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS)
496 RH2_OLD = Q2_NDG_OLD(I,J) / QSBT
497 VAPPRS = SVP1 * EXP(SVP2 * (T2_NEW - SVPT0) / (T2_NEW - SVP3))
498 QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS)
499 RH2_NEW = Q2_NDG_NEW(I,J) / QSBT
500 RH2OBS = CORB * RH2_OLD + CORE * RH2_NEW
501 T2OBS(I,J) = CORB * T2_OLD + CORE * T2_NEW
502 Q2OBS(I,J) = CORB * Q2_NDG_OLD(I,J) + CORE * Q2_NDG_NEW(I,J)
503 SNOBS = CORB * SN_NDG_OLD(I,J) + CORE * SN_NDG_NEW(I,J)
504 ENDIF
506 USTAR = MAX(UST(I,J),0.005)
508 IF (IFLAND .GE. 1.5) THEN ! if over water
509 ZNT(I,J) = CZO * UST(I,J) * UST(I,J) / G + OZO
510 ENDIF
512 Z0 = ZNT(I,J)
513 CPAIR = CPD * (1.0 + 0.84 * QV1) ! J/(K KG)
515 ! Set WRF Snow albedo to PX snow albedo based on fractional landuse
516 ! Snow albedo for each landuse class is defined in module_sf_pxlsm_data.F
517 SNOALB(I,J) = XSNOALB(I,J)
518 !-------------------------------------------------------------
519 ! Compute fractional snow area and snow albedo
520 CALL PXSNOW (ITIMESTEP, SNOBS, SNOWNCV(I,J), SNOW(I,J), &
521 SNOWH(I,J), XSNUP(I,J), XALB(i,j), &
522 SNOALB(I,J),VEGF_PX(I,J), SHDMIN(I,J), &
523 HC_SNOW, SNOW_FRA, SNOWC(I,J), ALBEDO(I,J) )
524 !-------------------------------------------------------------
526 !-------------------------------------------------------------
527 ! Sea Ice from analysis and water cells that are very cold, but more than 50% water
528 ! are converted to ice/snow for more reasonable treatment.
529 IF( (XICE(I,J).GE.0.5) .OR. &
530 (SST(I,J).LE.270.0.AND.XLAND(I,J).GE.1.50) ) THEN
531 XLAND(I,J) = 1.0
532 IFLAND = 1.0
533 ZNT(I,J) = 0.001
534 SMOIS(I,1,J) = 1.0
535 SMOIS(I,2,J) = 1.0
536 XICE(I,J) = 1.0
537 ALBEDO(I,J) = 0.7
538 SNOWC(I,J) = 1.0
539 SNOW_FRA = 1.0
540 VEGF_PX(I,J) = 0.0
541 LAI_PX(I,J) = 0.0
542 LAI(I,J) = 0.0
543 FCGSAT = 3.670
544 ENDIF
545 !-------------------------------------------------------------
547 !-------------------------------------------------------------
548 !-- Note that when IFGROW = 0 is selected in Vegeland then max and min
549 !-- LAI and Veg are the same
550 T2I = TSLB(I,2,J)
551 FSEAS = AMAX1(1.0 - 0.015625 * (290.0 - T2I) ** 2,0.0)
552 IF (T2I .GE. 290.0) FSEAS = 1.0
554 !get PX table vegetation
555 LAI_PX(I,J) = XLAIMN(I,J) + FSEAS*(XLAI(I,J) - XLAIMN(I,J))
556 VEGF_PX(I,J) = XVEGMN(I,J) + FSEAS*(XVEG(I,J) - XVEGMN(I,J))
558 !... use MODIS LAI and VEGFRA from wrflowinp
559 IF ( pxlsm_modis_veg .EQ. 1 ) THEN
561 ! get LAI for vegetated area
562 IF ( VEGFRA(I,J) .GT. 0.0 ) THEN
563 LAI_PX(I,J) = LAI(I,J) / ( VEGFRA(I,J) / 100.0)
564 ELSE
565 LAI_PX(I,J) = 0.0
566 ENDIF
568 VEGF_PX(I,J) = VEGFRA(I,J) / 100.0
570 !vegF is just for the land
571 IF ( LANDUSEF(I,KWAT,J) .LT. 1.0 ) THEN
572 VEGF_PX(I,J) = VEGF_PX(I,J) / (1.0 - LANDUSEF(I,KWAT,J))
573 ELSE
574 VEGF_PX(I,J) = 0.0
575 ENDIF
577 ENDIF
579 LAI_PX(I,J) = MIN(LAI_PX(I,J), 8.0)
580 LAI_PX(I,J) = MAX(LAI_PX(I,J), 0.0001)
582 VEGF_PX(I,J) = MIN(VEGF_PX(I,J), 1.0)
583 VEGF_PX(I,J) = MAX(VEGF_PX(I,J), 0.0001)
585 !... END OF MODIS LAI and FPAR
588 ! Ensure veg algorithms not used for water
589 IF (IFLAND .GE. 1.5) THEN
590 VEGF_PX(I,J) = 0.0
591 LAI_PX(I,J) = 0.0
592 ENDIF
593 !-------------------------------------------------------------
596 SOLDN = GSW(I,J) / (1.0 - ALBEDO(I,J)) ! downward shortwave radiaton
597 ISNOW = SNOWC(I,J)
600 NUDGE=PXLSM_SOIL_NUDGE
601 IF ( J .LE. 2 .OR. J .GE. (jde-1) ) NUDGE=0
602 IF ( I .LE. 2 .OR. I .GE. (ide-1) ) NUDGE=0
604 IF ( RMOL(I,J) .GT. 0.0 ) THEN
605 MOLX = AMIN1(1/RMOL(I,J),1000.0)
606 ELSE IF ( RMOL(I,J) .LT. 0.0 ) THEN
607 MOLX = AMAX1(1/RMOL(I,J),-1000.0)
608 ELSE
609 MOLX = 1000
610 ENDIF
612 ZFUNC = ZF1 * (1.0 - ZF1 / MAX(100.,PBLH(I,J))) ** 2
613 QST12 = KARMAN * ZFUNC*(QV2-QV1) / (ZA2-ZLVL)
616 !-------------------------------------------------------------
617 !-- LSM sub-time loop too prevent dt > 40 sec
618 NTSPS = INT(DT / (DTPBLX + 0.000001) + 1.0)
619 DTPBL = DT / NTSPS
621 DO IT=1,NTSPS
623 !... SATURATION VAPOR PRESSURE (MB)
624 IF ( TSLB(I,1,J) .LE. SVPT0 ) THEN ! For ground that is below freezing
625 ES = SVP1 * EXP(22.514 - 6.15E3 / TSLB(I,1,J)) ! cb
626 ELSE
627 ES = SVP1 * EXP(SVP2 * (TSLB(I,1,J) - SVPT0) / (TSLB(I,1,J) - SVP3))
628 ENDIF
629 QSS = ES * 0.622 / (SFCPRS - ES)
631 !... beta method, Lee & Pielke (JAM,May1992)
632 BETAP = 1.0
633 IF (IFLAND .LT. 1.5 .AND. ISNOW .LT. 0.5 .AND. SMOIS(I,1,J) .LE. FWFC) THEN
634 BETAP = 0.25 * (1.0 - COS(SMOIS(I,1,J) / FWFC * PI)) ** 2
635 ENDIF
637 !-------------------------------------------------------------------------
638 CALL SURFPX (DTPBL, IFLAND, SNOWC(I,J), NUDGE, XICE(I,J), & !in
639 SOLDN, GSW(I,J), LWDN, EMISSI, ZLVL, & !in
640 MOLX, Z0, USTAR, & !in
641 SFCPRS, DENS1, QV1, QSS, TA1, & !in
642 THETA1, PRECIP, & !in
643 CPAIR, PSIH(I,J), & !in
644 RH2OBS,T2OBS(I,J), & !in
645 VEGF_PX(I,J), ISTI, LAI_PX(I,J), IMPERV(I,J), CANFRA(I,J), & !in
646 BETAP, RSTMIN(I,J), HC_SNOW, SNOW_FRA, WETFRA(I,J), & !in
647 FWWLT, FWFC, FWRES, FCGSAT, FWSAT, FB, & !in ! Soil model updates - JEP 12/14
648 FC1SAT,FC2R,FAS,FJP,FC3,DZS(1),DZS(2),QST12, & !in
649 RADNET(I,J), GRDFLX(I,J), HFX(I,J), QFX(I,J), LH(I,J), & !out
650 EG(I,J), ER(I,J), ETR(I,J), & !out
651 QST(I,J), CAPG(I,J), RS(I,J), RA(I,J), & !out
652 TSLB(I,1,J), TSLB(I,2,J), & !out
653 SMOIS(I,1,J), SMOIS(I,2,J), WR, &
654 TA2(I,J), QA2(I,J), LAND_USE_TYPE,I,J )
655 !-------------------------------------------------------------------------
657 END DO ! Time internal PX time loop
659 IF (IFLAND .GE. 1.5) THEN
660 TSK(I,J) = SST(I,J) ! Skin temp set to sea surface temperature for open water
661 GRDFLX(I,J) = 0.0
662 ELSE
663 TSK(I,J) = TSLB(I,1,J) ! Skin temp set to 1 cm soil temperature in PX for now
664 ENDIF
665 CANWAT(I,J) = WR * 1000. ! convert WR back to mm for CANWAT
666 RAW = RA(I,J) + 4.503 / USTAR
667 QSFC(I,J) = QFX(I,J) * RAW / DENS1 + QV1
669 ENDDO ! END MIAN I LOOP
670 ENDDO ! END MAIN J LOOP
672 !------------------------------------------------------
673 END SUBROUTINE pxlsm
674 !-------------------------------------------------------------------------
675 !-------------------------------------------------------------------------
678 !-------------------------------------------------------------------------
679 !-------------------------------------------------------------------------
680 SUBROUTINE VEGELAND( landusef, vegfra, &
681 shdmin, shdmax, &
682 soilctop, soilcbot, nlcat, nscat,znt, xlai, &
683 xlaimn, rstmin, xveg, xvegmn, xsnup, xland, &
684 xalb, xsnoalb, wetfra, imperv, canfra, &
685 ids,ide, jds,jde, kds,kde, &
686 ims,ime, jms,jme, kms,kme, &
687 its,ite, jts,jte, kts,kte, &
688 LAND_USE_TYPE, KWAT_OUT )
689 !-------------------------------------------------------------------------
690 !
691 ! CALLED FROM Sub. bl_init in module_physics.init.F
692 !
693 ! THIS PROGRAM PROCESSES THE USGS LANDUSE DATA
694 ! WHICH HAS BEEN GRIDDED BY THE WPS SYSTEM
695 ! AND PRODUCES 2-D FIELDS OF LU RELATED PARAMETERS
696 ! FOR USE IN THE PX SURFACE MODEL
697 !
698 !
699 ! REVISION HISTORY:
700 ! AX Oct 2004 - developed WRF version based on VEGELAND in the MM5 PX LSM
701 ! RG Dec 2006 - revised for WRFV2.1.2
702 ! JP Dec 2007 - revised for WRFV3 - removed IFGROW options
703 !-------------------------------------------------------------------------
704 !-------------------------------------------------------------------------
706 IMPLICIT NONE
707 !...
708 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
709 ims,ime, jms,jme, kms,kme, &
710 its,ite, jts,jte, kts,kte
712 INTEGER , INTENT(IN) :: NSCAT, NLCAT
714 REAL, DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN) :: LANDUSEF
715 REAL, DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN) :: SOILCTOP, SOILCBOT
717 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: VEGFRA, SHDMIN, SHDMAX
719 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ZNT, IMPERV, CANFRA
721 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: XLAI, XLAIMN, RSTMIN, XALB, &
722 XVEG, XVEGMN, XSNUP, XLAND, &
723 WETFRA, XSNOALB
725 INTEGER, INTENT(OUT) :: KWAT_OUT
727 CHARACTER (LEN = 6), INTENT(IN) :: LAND_USE_TYPE
730 !... local variables
732 INTEGER :: ITF, JTF, K, J, I
733 REAL :: SUMLAI, SUMLMN, SUMRSI, SUMLZ0, SUMVEG, SUMVMN, &
734 ALAI, VEGF, SUMSNUP, SUMALB, SUMSNOALB
736 REAL :: VFMX, VFMN, VSEAS, FAREA, FWAT, ZNOTC, FCAN, FIMP, FORFRA, EXTFOR
737 REAL, DIMENSION( NLCAT ) :: LAIMX, LAIMN, Z0, VEG, VEGMN, SNUP, ALB, SNOALB
739 REAL, PARAMETER :: ZNOTCMN = 5.0 ! CM, MIN Zo FOR CROPS
740 REAL, PARAMETER :: ZNOTCMX = 15.0 ! CM, MAX Zo FOR CROPS
742 REAL, SAVE, DIMENSION(:), POINTER :: RSMIN, Z00, VEG0, VEGMN0, LAI0, &
743 LAIMN0, SNUP0, ALBF, SNOALBF
745 !---- INITIALIZE PARAMETERS
746 INTEGER, SAVE :: KWAT, LIMIT1, LIMIT2
748 ! Initialize LU characteristics by LU Dataset
749 IF (LAND_USE_TYPE == 'USGS') THEN
750 KWAT = 16
751 RSMIN => RSMIN_USGS
752 Z00 => Z00_USGS
753 VEG0 => VEG0_USGS
754 VEGMN0 => VEGMN0_USGS
755 LAI0 => LAI0_USGS
756 LAIMN0 => LAIMN0_USGS
757 SNUP0 => SNUP0_USGS
758 ALBF => ALBF_USGS
759 SNOALBF=> SNOALB_USGS
760 LIMIT1 = 2
761 LIMIT1 = 6
762 ELSE IF (LAND_USE_TYPE == 'USGS28') THEN
763 KWAT = 16
764 RSMIN => RSMIN_USGS28
765 Z00 => Z00_USGS28
766 VEG0 => VEG0_USGS28
767 VEGMN0 => VEGMN0_USGS28
768 LAI0 => LAI0_USGS28
769 LAIMN0 => LAIMN0_USGS28
770 SNUP0 => SNUP0_USGS28
771 ALBF => ALBF_USGS28
772 SNOALBF=> SNOALB_USGS28
773 LIMIT1 = 2
774 LIMIT1 = 6
775 ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN
776 KWAT = 1
777 RSMIN => RSMIN_NLCD50
778 Z00 => Z00_NLCD50
779 VEG0 => VEG0_NLCD50
780 VEGMN0 => VEGMN0_NLCD50
781 LAI0 => LAI0_NLCD50
782 LAIMN0 => LAIMN0_NLCD50
783 SNUP0 => SNUP0_NLCD50
784 ALBF => ALBF_NLCD50
785 SNOALBF=> SNOALB_NLCD50
786 LIMIT1 = 20
787 LIMIT1 = 43
788 ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN
789 KWAT = 17
790 RSMIN => RSMIN_NLCD40
791 Z00 => Z00_NLCD40
792 VEG0 => VEG0_NLCD40
793 VEGMN0 => VEGMN0_NLCD40
794 LAI0 => LAI0_NLCD40
795 LAIMN0 => LAIMN0_NLCD40
796 SNUP0 => SNUP0_NLCD40
797 ALBF => ALBF_NLCD40
798 SNOALBF=> SNOALB_NLCD40
799 LIMIT1 = 20
800 LIMIT1 = 43
801 ELSE IF (LAND_USE_TYPE == 'MODIS') THEN
802 KWAT = 17
803 RSMIN => RSMIN_MODIS
804 Z00 => Z00_MODIS
805 VEG0 => VEG0_MODIS
806 VEGMN0 => VEGMN0_MODIS
807 LAI0 => LAI0_MODIS
808 LAIMN0 => LAIMN0_MODIS
809 SNUP0 => SNUP0_MODIS
810 ALBF => ALBF_MODIS
811 SNOALBF=> SNOALB_MODIS
812 LIMIT1 = 12
813 LIMIT1 = 14
814 END IF
816 KWAT_OUT = KWAT
817 !--------------------------------------------------------------------
818 DO J = jts,jte
819 DO I = its,ite
820 XLAI(I,J) = 0.0
821 XLAIMN(I,J) = 0.0
822 RSTMIN(I,J) = 9999.0
823 XVEG(I,J) = 0.0
824 XVEGMN(I,J) = 0.0
825 XSNUP(I,J) = 0.0
826 XALB(I,J) = 0.0
827 XSNOALB(I,J)= 0.0
829 ! Code that may be needed in case these arrays are not intialized
830 ! with zero by real.exe when not defined by GEOGRID or present
831 ! in met_em* files
832 !IMPERV(I,J) = AMAX1(0.0001,IMPERV(I,J))
833 !CANFRA(I,J) = AMAX1(0.0001,CANFRA(I,J))
835 ENDDO ! END LOOP THROUGH GRID CELLS
836 ENDDO ! END LOOP THROUGH GRID CELLS
837 !--------------------------------------------------------------------
839 DO J = jts,jte
840 DO I = its,ite
841 !-- Initialize 2 and 3-D veg parameters to be caculated
842 DO K=1,NLCAT
843 LAIMX(K) = LAI0(K)
844 LAIMN(K) = LAIMN0(K)
845 Z0(K) = Z00(K)
846 VEG(K) = VEG0(K)
847 VEGMN(K) = VEGMN0(K)
848 SNUP(K) = SNUP0(K)
849 ALB(K) = ALBF(K)
850 SNOALB(K)= SNOALBF(K)
851 ENDDO
853 !-- INITIALIZE SUMS
854 SUMLAI = 0.0
855 SUMLMN = 0.0
856 SUMRSI = 0.0
857 SUMLZ0 = 0.0
858 SUMVEG = 0.0
859 SUMVMN = 0.0
860 ALAI = 0.0
861 SUMSNUP = 0.0
862 SUMALB = 0.0
863 SUMSNOALB = 0.0
865 !-- ESTIMATE CROP EMERGANCE DATE FROM VEGFRAC
866 VFMX = SHDMAX(I,J)
867 VFMN = SHDMIN(I,J)
868 VEGF = VEGFRA(I,J)
870 !-- Computations for VEGETATION CELLS ONLY
871 IF(VFMX.GT.0.0.AND.LANDUSEF(I,KWAT,J).LT.1.00) THEN
872 VSEAS = VEGF/VFMX
873 IF(VSEAS.GT.1.0.OR.VSEAS.LT.0.0) THEN
874 VSEAS = MIN(VSEAS,1.0)
875 VSEAS = MAX(VSEAS,0.0)
876 ENDIF
878 ZNOTC = ZNOTCMN * (1-VSEAS) + ZNOTCMX * VSEAS ! Zo FOR CROPS
879 DO K = 1, NLCAT
880 IF (LAND_USE_TYPE == 'MODIS') THEN
881 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
882 IF (K.EQ.12 .OR. K.EQ.14) THEN
883 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
884 LAIMN(K) = LAIMX(K)
885 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
886 VEGMN(K) = VEG(K)
887 !-- SEASONALLY VARY Zo FOR MODIS DryCrop (k=12)
888 IF (K .EQ. 12) THEN
889 Z0(K) = ZNOTC
890 !-- CrGrM (k=14) USE AVG WITH GRASS AND FOREST
891 ELSE IF (K .EQ.14) THEN
892 Z0(K) = 0.5 * (ZNOTC + Z00(K))
893 ENDIF
894 ENDIF
895 ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN
896 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
897 IF (K.EQ.20 .OR. K.EQ.43 .OR. K.EQ.45) THEN
898 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
899 LAIMN(K) = LAIMX(K)
900 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
901 VEGMN(K) = VEG(K)
902 !-- SEASONALLY VARY Zo FOR DryCrop (k=20) OR Irigated Crop (k=43) OR Mix Crop (k=4)
903 IF (K.EQ.20 .OR. K.EQ.43) THEN
904 Z0(K) = ZNOTC
905 !-- CrNatM (k=45) USE AVG WITH GRASS AND FOREST
906 ELSE IF (K.EQ.45) THEN
907 Z0(K) = 0.5 * (ZNOTC + Z00(K))
908 ENDIF
909 ENDIF
910 ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN
911 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
912 IF (K.EQ.12 .OR. K.EQ.14 .OR. K.EQ.38) THEN
913 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
914 LAIMN(K) = LAIMX(K)
915 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
916 VEGMN(K) = VEG(K)
917 !-- SEASONALLY VARY Zo FOR Crop (k=12 for MODIS or 38 for NLCD)
918 IF (K.EQ.12 .OR. K.EQ.38) THEN
919 Z0(K) = ZNOTC
920 !-- CrNatM (k=14) USE AVG WITH GRASS AND FOREST
921 ELSE IF (K.EQ.14) THEN
922 Z0(K) = 0.5 * (ZNOTC + Z00(K))
923 ENDIF
924 ENDIF
925 ELSE IF (LAND_USE_TYPE == 'USGS') THEN
926 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
927 IF (K .GE. 2 .AND. K .LE. 6) THEN
928 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
929 LAIMN(K) = LAIMX(K)
930 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
931 VEGMN(K) = VEG(K)
932 !-- SEASONALLY VARY Zo FOR DryCrop (k=2) OR Irigated Crop (k=3) OR Mix Crop (k=4)
933 IF (K .GE. 2 .AND. K .LE. 4) THEN
934 Z0(K) = ZNOTC
935 !-- CrGrM (k=5) or CrWdM (k=6) USE AVG WITH GRASS AND FOREST
936 ELSE IF (K .GE.5 .AND. K .LE. 6) THEN
937 Z0(K) = 0.5 * (ZNOTC + Z00(K))
938 ENDIF
939 ENDIF
940 ELSE IF (LAND_USE_TYPE == 'USGS28') THEN
941 !-- USE THE VEGFRAC DATA ONLY FOR CROPS
942 IF (K .GE. 2 .AND. K .LE. 6) THEN
943 LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
944 LAIMN(K) = LAIMX(K)
945 VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
946 VEGMN(K) = VEG(K)
947 !-- SEASONALLY VARY Zo FOR DryCrop (k=2) OR Irigated Crop (k=3) OR Mix Crop (k=4)
948 IF (K .GE. 2 .AND. K .LE. 4) THEN
949 Z0(K) = ZNOTC
950 !-- CrGrM (k=5) or CrWdM (k=6) USE AVG WITH GRASS AND FOREST
951 ELSE IF (K .GE.5 .AND. K .LE. 6) THEN
952 Z0(K) = 0.5 * (ZNOTC + Z00(K))
953 ENDIF
954 ENDIF
956 END IF
958 ENDDO
960 ENDIF !-- IF cell is vegetation
962 !-------------------------------------
963 !-- LOOP THROUGH LANDUSE Fraction and compute totals
964 DO K = 1, NLCAT
965 FAREA = LANDUSEF(I,K,J)
966 SUMLAI = SUMLAI + LAIMX(K) * FAREA
967 SUMLMN = SUMLMN + LAIMN(K) * FAREA
968 ALAI = ALAI + FAREA
969 SUMRSI = SUMRSI + FAREA * LAIMX(K) / RSMIN(K)
970 SUMLZ0 = SUMLZ0 + FAREA * ALOG(Z0(K))
971 SUMVEG = SUMVEG + FAREA * VEG(K)
972 SUMVMN = SUMVMN + FAREA * VEGMN(K)
973 SUMSNUP = SUMSNUP+ FAREA * SNUP(K)
974 SUMALB = SUMALB + FAREA * ALB(K)
975 SUMSNOALB= SUMSNOALB + FAREA * SNOALB(K)
976 ENDDO
978 FWAT = LANDUSEF(I,KWAT,J)
979 !-- CHECK FOR WATER
980 IF (FWAT .GE. 0.50) THEN ! Changed WRFV3.7
981 XLAI(I,J) = LAIMX(KWAT)
982 XLAIMN(I,J) = LAIMN(KWAT)
983 RSTMIN(I,J) = RSMIN(KWAT)
984 ZNT(I,J) = Z0(KWAT)
985 XVEG(I,J) = VEG(KWAT)
986 XVEGMN(I,J) = VEGMN(KWAT)
987 XSNUP(I,J) = SNUP(KWAT)
988 XALB(I,J) = ALB(KWAT)
989 XSNOALB(I,J)= SNOALB(KWAT)
990 ELSE
991 IF (FWAT .GT. 0.10) THEN
992 ALAI = ALAI - FWAT
993 SUMLZ0 = SUMLZ0 - FWAT * ALOG(Z0(KWAT))
994 ENDIF
995 XLAI(I,J) = SUMLAI / ALAI
996 XLAIMN(I,J) = SUMLMN / ALAI
997 RSTMIN(I,J) = SUMLAI / SUMRSI
998 ZNT(I,J) = EXP(SUMLZ0/ALAI)
999 XVEG(I,J) = SUMVEG / ALAI
1000 XVEGMN(I,J) = SUMVMN / ALAI
1001 XSNUP(I,J) = SUMSNUP
1002 XALB(I,J) = SUMALB
1003 XSNOALB(I,J)= SUMSNOALB
1004 ENDIF
1006 IF (FWAT .GT. 0.50) THEN
1007 ZNT(I,J) = Z0(KWAT)
1008 XALB(I,J) = ALB(KWAT)
1009 XSNOALB(I,J)= SNOALB(KWAT)
1010 ENDIF
1012 !-- Compute wetlands fraction for proper MMLUIN data set
1013 !-- Note: if LU categories change, these hard coded indicies must be updated
1014 IF (LAND_USE_TYPE == 'USGS') THEN
1015 WETFRA(I,J) = LANDUSEF(I,17,J)+LANDUSEF(I,18,J)
1016 ELSE IF (LAND_USE_TYPE == 'USGS28') THEN
1017 WETFRA(I,J) = LANDUSEF(I,17,J)+LANDUSEF(I,18,J)
1018 ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN
1019 WETFRA(I,J) = LANDUSEF(I,22,J)+LANDUSEF(I,23,J)+LANDUSEF(I,27,J)+LANDUSEF(I,28,J)+LANDUSEF(I,42,J)
1020 ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN
1021 WETFRA(I,J) = LANDUSEF(I,39,J)+LANDUSEF(I,40,J)+LANDUSEF(I,11,J)
1022 ELSE IF (LAND_USE_TYPE == 'MODIS') THEN
1023 WETFRA(I,J) = LANDUSEF(I,11,J)
1024 END IF
1026 ZNT(I,J) = ZNT(I,J) * 0.01 !CONVERT TO M
1027 XVEG(I,J) = XVEG(I,J) * 0.01 !CONVERT TO FRAC
1028 XVEGMN(I,J) = XVEGMN(I,J) * 0.01
1029 XLAND(I,J) = 1.0 + FWAT
1030 XALB(I,J) = XALB(I,J) * 0.01
1031 XSNOALB(I,J)= XSNOALB(I,J) * 0.01
1033 !-------Adjustment according to CANFRA and IMPERV fo NLCD40 only -----------
1034 FIMP = IMPERV(I,J) * 0.01
1035 FCAN = CANFRA(I,J) * 0.01
1036 IF (LAND_USE_TYPE == 'NLCD40') THEN
1037 XVEG(I,J) = MIN(XVEG(I,J),1.0-FIMP)
1038 XVEGMN(I,J) = MIN(XVEGMN(I,J),1.0-FIMP)
1039 XVEG(I,J) = MAX(XVEG(I,J),FCAN)
1040 XVEGMN(I,J) = MAX(XVEGMN(I,J),FCAN)
1042 FORFRA = LANDUSEF(I,39,J)+LANDUSEF(I,30,J)+LANDUSEF(I,29,J)+LANDUSEF(I,28,J)
1043 EXTFOR = FCAN - FORFRA
1044 IF (EXTFOR.GE.0.01) THEN
1045 XLAI(I,J) = LAIMX(30) * EXTFOR + XLAI(I,J) * (1-EXTFOR)
1046 XLAIMN(I,J) = LAIMN(30) * EXTFOR + XLAIMN(I,J) * (1-EXTFOR)
1047 ENDIF
1048 ENDIF
1049 !--------------------------------------------------------------------------
1051 ENDDO ! END LOOP THROUGH GRID CELLS
1052 ENDDO ! END LOOP THROUGH GRID CELLS
1053 !--------------------------------------------------------------------
1055 END SUBROUTINE vegeland
1057 !------------------------------------------------------------------------------
1058 !------------------------------------------------------------------------------
1059 SUBROUTINE SURFPX(DTPBL, IFLAND, ISNOW, NUDGEX, XICE1, SOLDN, GSW, & !in
1060 LWDN, EMISSI, Z1, MOL, ZNT, UST, PSURF, DENS1, & !in
1061 QV1, QSS, TA1, THETA1, PRECIP, CPAIR, PSIH, & !in
1062 RH2OBS, T2OBS, VEGFRC, ISTI,LAI,IMPERV,CANFRA,BETAP, & !in
1063 RSTMIN, HC_SNOW, SNOW_FRA, WETFRA, WWLT, WFC, & !in
1064 WRES, CGSAT, WSAT, B, C1SAT, C2R, AS, JP, C3, DS1, & !in
1065 DS2, QST12, & !in
1066 RADNET, GRDFLX, HFX, QFX, LH, EG, ER, ETR, & !out
1067 QST, CAPG, RS, RA, TG, T2, WG, W2, WR, & !out
1068 TA2, QA2, LAND_USE_TYPE, I, J ) !out
1070 !------------------------------------------------------------------------------
1071 !
1072 ! FUNCTION:
1073 ! THIS SUBROUTINE COMPUTES SOIL MOISTURE AND TEMPERATURE TENDENCIES
1074 ! BY SOLVING THE PROGNOSTIC EQUATIONS IN PX95.
1075 !
1076 ! SUBROUTINES CALLED:
1077 ! SUB. QFLUX compute the soil and canopy evaporation, and transpiration
1078 ! SUB. SMASS compute nudging coefficients for soil moisture and temp nudging
1079 !
1080 ! ARGUMENTS:
1081 ! DTPBL: TIME STEP OF THE MINOR LOOP FOR THE LAND-SURFACE/PBL MODEL
1082 ! IFLAND: INDEX WHICH INDICATES THE TYPE OF SURFACE,=1,LAND;=2,SEA
1083 ! ISNOW: SNOW (=1) OR NOT (=0)
1084 ! NUDGE: SWITCH FOR SOIL MOISTURE NUDGING
1085 ! SOLDN: SHORT-WAVE RADIATION
1086 ! LWDN: LONG-WAVE RADIATION
1087 ! EMISSI: EMISSIVITY
1088 ! Z1: HEIGHT OF THE LOWEST HALF LAYER
1089 ! MOL: MONIN-OBUKOV LENGH (M)
1090 ! ZNT: ROUGHNESS LENGTH (M)
1091 ! UST: FRICTION VELOCITY (M/S)
1092 ! TST: Turbulent moisture scale
1093 ! RA: AERODYNAMIC RESISTENCE
1094 ! PSURF: P AT SURFACE (CB)
1095 ! DENS1: AIR DENSITY AT THE FIRST HALF LAYER
1096 ! QV1: Air humidity at first half layer
1097 ! QSS: Saturation mixing ratio at ground
1098 ! TA1: Air temperature at first half layer
1099 ! THETA1: Potential temperature at first half layer
1100 ! PRECIP: Precipitation rate in m/s
1101 ! STBOLT: STEFAN BOLTZMANN'S CONSTANT
1102 ! KARMAN: VON KARMAN CONSTANT
1103 ! CPAIR: Specific heat of moist air (M^2 S^-2 K^-1)
1104 ! TAUINV: 1/1DAY(SEC)
1105 ! VEGFRC: Vegetation coverage
1106 ! ISTI: soil type
1107 ! LAI: Leaf area index
1108 ! IMPERV: Percentage of IMPERVIOUS Fraction
1109 ! CANFRA: Percentage of Canopy/Tree Fraction
1110 ! BETAP: Coefficient for bare soil evaporation
1111 ! WETFRA: Fraction of Wetlands area
1112 ! THZ1OB: Observed TEMP FROM ANAL OF OBS TEMPS AT SCREEN HT
1113 ! RHOBS: Observed relative humidity at SCREEN HT
1114 ! RSTMIN Minimum Stomatol resistence
1115 !... Outputs from SURFPX
1116 ! RADNET: Net radiation
1117 ! HFX: SENSIBLE HEAT FLUX (W / M^2)
1118 ! QFX: TOTAL EVAP FLUX (KG / M^2 S)
1119 ! GRDFLX: Ground heat flux (W / M^2)
1120 ! QST: Turbulent moisture scale
1121 ! CAPG: THERMAL CAPACITY OF GROUND SLAB (J/M^2/K)
1122 ! RS: Surface resistence
1123 ! RA: Surface aerodynamic resistence
1124 ! EG: evaporation from ground (bare soil)
1125 ! ER: evaporation from canopy
1126 ! ETR: transpiration from vegetation
1127 ! TA2: diagnostic 2-m temperature from surface layer and lsm
1128 !
1129 !... Updated variables in this subroutine
1130 ! TG: Soil temperature at first soil layer
1131 ! T2: Soil temperature in root zone
1132 ! WG: Soil moisture at first soil layer
1133 ! W2: Soil moisture at root zone
1134 ! WR: Canopy water content in m
1135 !
1136 ! REVISION HISTORY:
1137 ! AX 2/2005 - developed WRF version based on the MM5 PX LSM
1138 !
1139 !------------------------------------------------------------------------------
1140 !------------------------------------------------------------------------------
1141 IMPLICIT NONE
1143 !.......Arguments
1145 !.. Integer
1146 INTEGER , INTENT(IN) :: ISTI, NUDGEX, I, J
1148 !... Real
1149 REAL , INTENT(IN) :: DTPBL, DS1, DS2
1150 REAL , INTENT(IN) :: IFLAND, ISNOW, XICE1
1151 REAL , INTENT(IN) :: SOLDN, GSW, LWDN, EMISSI, Z1
1152 REAL , INTENT(IN) :: ZNT
1153 REAL , INTENT(IN) :: PSURF, DENS1, QV1, QSS, TA1, THETA1, PRECIP
1154 REAL , INTENT(IN) :: CPAIR
1155 REAL , INTENT(IN) :: VEGFRC, LAI, IMPERV, CANFRA
1156 REAL , INTENT(IN) :: RSTMIN, HC_SNOW, SNOW_FRA, WETFRA
1157 REAL , INTENT(IN) :: WWLT, WFC, WRES, CGSAT, WSAT, B, C1SAT, C2R, AS, JP, C3
1158 REAL , INTENT(IN) :: RH2OBS,T2OBS
1159 REAL , INTENT(IN) :: QST12
1161 REAL , INTENT(OUT) :: RADNET, EG, ER, ETR
1162 REAL , INTENT(OUT) :: QST, CAPG, RS, TA2, QA2
1164 REAL , INTENT(INOUT) :: TG, T2, WG, W2, WR, UST, RA, BETAP
1165 REAL , INTENT(INOUT) :: GRDFLX, QFX, HFX, LH, PSIH, MOL
1167 CHARACTER (LEN = 6), INTENT(IN) :: LAND_USE_TYPE ! Check Limei Ran, wrf38 set to 5
1169 !... Local Variables
1171 !... Real
1172 REAL :: HF, LV, CQ4, WETSAT, SM2
1173 REAL :: RAH, RAW, ET, W2CG, CG, CT, SOILFLX, CPOT, THETAG
1174 REAL :: ZOL, ZOBOL, ZNTOL, Y, Y0, PSIH15, YNT
1175 REAL :: WGNUDG, W2NUDG, ALPH1, ALPH2, BET1, BET2, T1P5
1176 REAL :: CQ1, CQ2, CQ3, COEFFNP1, COEFFN, TSNEW, TSHLF, T2NEW
1177 REAL :: ROFF, WRMAX, PC, DWR, PNET, TENDWR, WRNEW
1178 REAL :: COF1, CFNP1WR, CFNWR, PG, FASS
1179 REAL :: TENDW2, W2NEW, W2HLF, W2REL, C1, C2, WEQ, CFNP1, CFN, WGNEW
1180 REAL :: ALN10, TMP1, TMP2, TMP3, AA, AB, TST, RBH, CTVEG
1181 REAL :: QST1,PHIH,PSIOB
1182 REAL :: T2NUD, T2NUDF
1183 REAL :: VAPPRS, QSBT, RH2MOD, IMF, VEGF, SOILF
1184 REAL :: RSOIL, LDRY, DP ! Soil model updates - JEP 12/14
1185 REAL :: C1MAX,ZZA,ZZB,ZDEL,ZLY,ZA,ZB,ZY2
1186 REAL :: Rinc, Hcan ! JEP 2020
1187 REAL :: CQ3BG,CQ3VW,CQ3VG,SIGG
1189 !... Parameters
1190 REAL :: ZOBS, GAMAH, BETAH, SIGF, BH, CT_SNOW, CT_IMPERV
1192 REAL, PARAMETER :: CV = 1.2E-5 ! K-M2/J Note: Update from 8E-6 10/14 Jon Pleim
1194 PARAMETER (ZOBS = 1.5) ! height for observed screen temp., (m)
1195 PARAMETER (BH = 15.7)
1196 PARAMETER (GAMAH = 16. ) !11.6)
1197 PARAMETER (BETAH = 5.0 ) !8.21)
1198 PARAMETER (SIGF = 0.5) ! rain interception see LSM (can be 0-1)
1199 REAL, PARAMETER :: DWAT = 0.2178 ! [cm^2 / s] at 273.15K
1200 !--------------------------------------------------------------------
1201 ! OLD PX legacy value from MM5 ... unknown origin PARAMETER (CT_SNOW = 5.54E-5)
1202 ! New value of CT_SNOW calibrated using multilayer soil model where csnow=6.9E5 J/(m3 K)
1203 ! from NCAR (WRFV3.2 -WRFV3.6.1) CSM PARAMETER (CT_SNOW = 2.0E-5)
1204 PARAMETER (CT_SNOW = 2.0E-5)
1206 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1207 ! CS for concrete/asphalt/road material (http://www.engineeringtoolbox.com/specific-heat-capacity-d_391.html)
1208 ! cs_imperv = 920 J kg-1 K-1
1209 ! CS for asphalt and concrete from 0.1785 to 0.20 cal g-1 K-1 (http://pages.towson.edu/morgan/files/Impervious.pdf)
1210 ! the values above translate to ~750 to 850 J kg-1 K-1
1211 !
1212 ! 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
1213 ! 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
1214 ! inverse of the value above gives CT_IMPERV value of 1/0.000021 = 4.762E-6 K m2 J-1
1216 ! 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
1217 ! the approxiate concrete/asphalt density and layer thickness or represenative thickness. For starters (12/2011)
1218 ! well not use the PX first layer thickness, but representative thickness of most roads/parkinglots/buildings.
1219 ! for now welll use 6 inches or about 15 cm or 0.15 m. Density of concrete (normal) from
1220 ! Dorf, Richard. Engineering Handbook. New York: CRC Press, 1996. is ~2400 kg m-3.
1221 ! 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
1222 ! If you look at the range of possible values considering density differences of concrete/asphalt/etc
1223 ! Values can range from 2.5 to 6.0 E-6
1224 PARAMETER (CT_IMPERV = 3.268E-6)
1225 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1227 ALN10 = ALOG(10.0)
1228 RADNET = SOLDN - (EMISSI *(STBOLT *TG **4 - LWDN)) ! NET RADIATION
1229 !--------------------------------------------------------------------
1230 CPOT= (100.0 / PSURF) ** ROVCP ! rcp is global constant(module_model_constants)
1231 THETAG = TG * CPOT
1233 ZOL = Z1/MOL
1234 ZOBOL = ZOBS/MOL
1235 ZNTOL = ZNT/MOL
1237 !-----------------------------------------------------------------------------------------
1238 IF (MOL .LT. 0.0) THEN
1239 Y = ( 1.0 - GAMAH * ZOL ) ** 0.5
1240 Y0 = ( 1.0 - GAMAH * ZOBOL ) ** 0.5
1241 YNT = ( 1.0 - GAMAH * ZNTOL ) ** 0.5
1242 PSIH15 = 2.0 * ALOG((Y + 1.0) / (Y0 + 1.0))
1243 PSIH = 2.0 * ALOG((Y + 1.0) / (YNT + 1.0))
1244 PSIOB = 2.0 * ALOG((Y0 + 1.0) / (YNT + 1.0))
1245 PHIH = 1.0 / Y
1246 ELSE
1247 IF((ZOL-ZNTOL).LE.1.0) THEN
1248 PSIH = -BETAH*(ZOL-ZNTOL)
1249 ELSE
1250 PSIH = 1.-BETAH-(ZOL-ZNTOL)
1251 ENDIF
1252 IF ((ZOBOL - ZNTOL) .LE. 1.0) THEN
1253 PSIOB = -BETAH * (ZOBOL - ZNTOL)
1254 ELSE
1255 PSIOB = 1.0 - BETAH - (ZOBOL - ZNTOL)
1256 ENDIF
1257 PSIH15 = PSIH - PSIOB
1258 IF(ZOL.LE.1.0) THEN
1259 PHIH = 1.0 + BETAH * ZOL
1260 ELSE
1261 PHIH = BETAH + ZOL
1262 ENDIF
1263 ENDIF
1264 !-----------------------------------------------------------------------------------------
1265 !-- ADD RA AND RB FOR HEAT AND MOISTURE
1266 !... RB FOR HEAT = 5 /UST
1267 !... RB FOR WATER VAPOR = 5*(0.599/0.709)^2/3 /UST = 4.503/UST
1268 RA=PR0* ( ALOG(Z1/ZNT) - PSIH )/(KARMAN*UST)
1269 RAH = RA + 5.0 / UST
1270 RAW = RA + 4.503 / UST
1271 IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN
1272 LDRY = 1.75*DS1*(EXP((1.-WG/WSAT)**5)-1.)/1.718 ! 1.75 cm is the layer thickness used by S&Z09
1273 DP = DWAT*1.E-4 * WSAT**2 * (1.-WRES/WSAT)**(2.+3./B)
1274 RSOIL=LDRY/DP
1276 ELSE
1277 RSOIL = 0.0
1278 ENDIF
1279 !--------------------------------------------------------------------
1280 ! Compute soil moisture layer 2 that considers fraction of saturated
1281 ! wetlands. If 100% of cell is wetland, soil moisture can be no lower
1282 ! than full soil saturation. If half wetland, no less than half saturated
1283 IF (IFLAND .LT. 1.5 ) THEN
1284 WETSAT = 1.00 * WSAT ! Wetlands soil moisture
1285 SM2 = (WETFRA * WETSAT) !
1286 W2 = AMAX1(SM2, W2) ! In case that W2 > Field capacity (heavy precip), use wetter W2
1287 ENDIF
1289 !--------------------------------------------------------------------
1290 ! In-canopy resistance - Erisman et al (1994)
1291 Hcan = ZNT * 10.0
1292 Rinc = 14.0 * LAI* Hcan / UST
1294 !--------------------------------------------------------------------
1295 !-- COMPUTE MOISTURE FLUX
1296 CALL QFLUX( DENS1, QV1, TA1, SOLDN, RAW, QSS, &
1297 VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, &
1298 WG, W2, WR, &
1299 RSTMIN, WWLT, WFC, RSOIL, RINC, &
1300 EG, ER, ETR, CQ4, RS, FASS, SIGG)
1301 !--------------------------------------------------------------------
1303 !--------------------------------------------------------------------
1304 ! Compute Total evaporation (ET) from various modes of moisture flux
1305 ET = EG + ER + ETR
1306 QST = -ET / (DENS1 * UST)
1308 LV = 2.83E6 !-- LATENT HEAT OF SUBLIMATION AT 0C FROM STULL(1988)
1309 IF (ISNOW .LT. 0.5.AND.TG.GT.273.15) &
1310 LV = (2.501 - 0.00237 * (TG - 273.15)) * 1.E6 !-- FROM STULL(1988) in J/KG
1312 QFX = ET
1313 LH = LV * QFX
1314 !-----------------------------------------------------------------------------------------
1316 !-----------------------------------------------------------------------------------------
1317 ! Surface sensible heat flux
1318 TST = (THETA1 - THETAG ) / (UST*RAH)
1319 HF = UST * TST
1320 HFX = AMAX1(-DENS1 * CPAIR * HF, -250.0) ! using -250. from MM5
1321 !-----------------------------------------------------------------------------------------
1323 !-----------------------------------------------------------------------------------------
1324 ! Compute the diagnosed 2m Q and T consistent with PX LSM
1325 QST1 = 0.5*(QST+QST12/PHIH)
1326 TA2 = (THETAG + TST * (PR0 / KARMAN * (ALOG(ZOBS / ZNT) - PSIOB)+5.))/CPOT
1327 QA2 = QV1 - QST1 * PR0/ KARMAN * (ALOG(Z1 / ZOBS) - PSIH15)
1329 IF (QA2.LE.0.0) QA2 = QV1
1331 !-- Relative humidity
1332 VAPPRS = SVP1 * EXP(SVP2 * (TA2 - SVPT0) / (TA2 - SVP3))
1333 QSBT = EP_2 * VAPPRS / (PSURF - VAPPRS)
1334 RH2MOD = QA2 / QSBT
1335 !-----------------------------------------------------------------------------------------
1336 IF (IFLAND .LT. 1.5 ) THEN
1337 W2CG = AMAX1(W2,WWLT)
1338 CG = CGSAT * 1.0E-6 * (WSAT/ W2CG) ** &
1339 (0.5 * B / ALN10)
1340 ! IMPERVIOUS weighting scheme -- Subtract highly accurate impervious fraction from cell
1341 ! remainder is split between ground and vegetation. CT is a weighted fractional average.
1342 ! Snow CT is then applied for final heat capacity
1343 IMF = AMAX1(0.0,IMPERV/100.0)
1344 VEGF = (1.0 - IMF) * VEGFRC
1345 SOILF = (1.0 - IMF) * (1.0 - VEGFRC)
1346 CT = 1./( IMF/CT_IMPERV + VEGF/CV + SOILF/CG)
1347 CT = 1./(SNOW_FRA/CT_SNOW + (1-SNOW_FRA)/CT)
1348 CAPG = 1.0/CT
1349 SOILFLX = 2.0 * PI * TAUINV * (TG - T2)
1350 GRDFLX = SOILFLX / CT
1351 ENDIF
1352 !-----------------------------------------------------------------------------------------
1354 !--------------------------------------------------------------------
1355 !-- ASSIMILATION --- COMPUTE SOIL MOISTURE NUDGING FROM TA2 and RH2
1356 !-------COMPUTE ASSIMILATION COEFFICIENTS FOR ALL I
1357 IF (IFLAND .LT. 1.5) THEN
1358 IF (NUDGEX .EQ. 0) THEN !-- NO NUDGING CASE
1359 WGNUDG = 0.0
1360 W2NUDG = 0.0
1361 T2NUD = 0.0
1362 ELSE !-- NUDGING CASE
1364 CALL SMASS (ISTI, FASS, SOLDN, VEGFRC, RA, WWLT, WFC, &
1365 ALPH1, ALPH2, BET1, BET2, T2NUDF)
1367 !--COMPUTE SOIL MOISTURE NUDGING
1368 WGNUDG = ALPH1 * (T2OBS - TA2) + ALPH2 * (RH2OBS - RH2MOD) * 100 !NUDGING W2 FOR NON-VEG
1369 W2NUDG = BET1 * (T2OBS - TA2) + BET2 * (RH2OBS - RH2MOD) * 100
1370 IF (W2 .GE. WFC) W2NUDG = AMIN1(W2NUDG,0.0)
1371 IF (W2 .LE. WWLT) W2NUDG = AMAX1(W2NUDG,0.0)
1372 IF (W2 .GE. WFC) WGNUDG = AMIN1(WGNUDG,0.0)
1373 IF (W2 .LE. WWLT) WGNUDG = AMAX1(WGNUDG,0.0)
1374 T2NUD = T2NUDF * (T2OBS - TA2)
1375 ENDIF
1376 ENDIF
1377 !-----------------------------------------------------------------------------------------
1379 !-----------------------------------------------------------------------------------------
1380 !-- Compute new values for TS,T2,WG,W2 and WR. No change over ice or water (XLAND > 1)
1381 IF (IFLAND .LT. 1.5) THEN
1382 !-- SOLVE BY CRANK-NIC -- TENDTS=CT*(RADNET-HFX-QFX)-SOILFLX
1383 !-- Calculate the coefficients for implicit calculation of TG
1384 CQ1 = (1.0 - 0.622 * LV * CRANKP / (r_d * TG)) * QSS
1385 CQ2 = 0.622 * LV * QSS * CRANKP / (r_d * TG * TG)
1386 CQ3BG = DENS1 * (1.0 - VEGFRC) / (RAW + RSOIL)
1387 CQ3VW = DENS1 * VEGFRC*SIGG / RAW
1388 CQ3VG = DENS1 * VEGFRC / (RAW + RSOIL + RINC)
1389 CQ3 = CQ3BG + CQ3VW + CQ3VG
1390 COEFFNP1 = 1.0 + DTPBL * CRANKP * (4.0 * EMISSI * STBOLT * TG ** 3 &
1391 * CT + DENS1 * CPAIR / RAH * CPOT * CT + 2.0 * PI &
1392 * TAUINV ) + DTPBL * (CT * LV * CQ2 * (CQ3 + CQ4))
1393 COEFFN = CT * (GSW + EMISSI * (STBOLT * (4.0 * CRANKP - 1.0) &
1394 * TG*TG*TG*TG + LWDN) & !NET RAD
1395 + DENS1 * CPAIR / RAH * (THETA1 - (1.0 - CRANKP) * THETAG) &
1396 - LV * (CQ3 * (CQ1 - QV1) + CQ4 * (CQ1 - QV1))) & !SFC HEAT FLUX
1397 - 2.0 * PI * TAUINV * ((1.0 - CRANKP) * TG - T2) !SOIL FLUX
1398 TSNEW = (TG + DTPBL * COEFFN) / COEFFNP1
1399 !-- FOR SNOW COVERED SURFACE TEMPERATURE IS NOT MORE THAN ZERO
1400 IF (XICE1 .GT. 0.5) TSNEW = AMIN1(TSNEW,273.15) ! Re-added Jan 2010 to keep ice surface at or below freezing (J. Pleim)
1401 TSHLF = 0.5 * ( TSNEW + TG)
1402 T2NEW = (T2 + DTPBL * TAUINV * T2TFAC * (TSHLF - (1 - CRANKP) * T2) &
1403 + DTPBL*T2NUD) & ! Added deep temperature nudging
1404 / (1.0 + DTPBL * TAUINV * T2TFAC * CRANKP)
1405 !-- REPLACE OLD with NEW Value
1406 TG = TSNEW
1407 T2 = T2NEW
1408 ENDIF
1409 !-----------------------------------------------------------------------------------------
1411 !-----------------------------------------------------------------------------------------
1412 ! Compute new subsurface soil and canopy moisture values DENS1. No change required over ocean.
1413 IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN
1414 !-- Compute WR
1415 ROFF = 0.0
1416 WRMAX = 0.2E-3 * VEGFRC * LAI ! max. WRMAX IN m
1418 IF(WRMAX.GT.0.0) THEN
1419 !-- PC is precip. intercepted by veg.(M/S)
1420 PC = VEGFRC * SIGF * PRECIP
1421 DWR = (WRMAX - WR) / DTPBL !the tendency to reach max.
1422 PNET = PC - ER/ DENW ! residual of precip. and evap.
1423 IF (PNET .GT. DWR) THEN
1424 ROFF = PNET - DWR
1425 PC = PC - ROFF
1426 ENDIF
1427 IF (QSS .LT. QV1) THEN
1428 TENDWR = PC - ER / DENW
1429 WRNEW = WR + DTPBL * TENDWR
1430 ELSE
1431 COF1 = DENS1 / DENW * VEGFRC * (QSS - QV1) / RAW
1432 !-- using delta=wr/wrmax
1433 CFNP1WR = 1.0 + DTPBL * COF1 * CRANKP / WRMAX
1434 CFNWR = PC - COF1 * (1.0 - CRANKP) * WR / WRMAX
1435 WRNEW = (WR + DTPBL * CFNWR) / CFNP1WR
1436 ENDIF
1437 ELSE
1438 PC=0.0
1439 WRNEW=0.0
1440 ENDIF
1441 !---------------------------------------------
1442 !-- Compute W2
1443 PG = DENW * (PRECIP - PC) ! PG is precip. reaching soil (PC already including ROFF)
1444 TENDW2 = 1.0 / (DENW * DS2) * (PG - EG - ETR) &
1445 - C3/DS2 * TAUINV * AMAX1(0.0,(W2 - WFC)) &
1446 + (W2NUDG + WGNUDG) / DS2 ! NUDGING
1447 W2NEW = W2 + DTPBL * TENDW2
1448 W2NEW = AMIN1(W2NEW,WSAT)
1449 W2NEW = AMAX1(W2NEW,WRES) !0.05) !Limei 08/02/2017
1450 W2HLF = 0.5 * (W2 + W2NEW)
1451 !.. new values
1452 W2 = W2NEW
1453 WR = AMIN1(WRMAX,WRNEW)
1454 ENDIF
1455 !-----------------------------------------------------------------------------------------
1457 !-----------------------------------------------------------------------------------------
1458 ! Compute new surface soil moisture values (WR).
1459 IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN ! over ocean no change to wg w2,wr
1460 !-- FOR SNOW COVERED SURFACE, ASSUME SURFACE IS SATURATED AND
1461 ! WG AND W2 ARE NOT CHANGED
1462 IF (ISNOW .GT.0.5) THEN
1463 WG = WSAT
1464 ELSE
1465 W2REL = W2HLF / WSAT
1466 IF (WG .GT. WWLT) THEN
1467 C1 = DS1*C1SAT * (WSAT / WG) ** (0.5 * B + 1.0)
1468 ELSE ! revise C1 for wg < wilting point Noilhan & Mahfouf (1996)
1469 ZY2 = C1SAT * (WSAT / WWLT) ** (0.5 * B + 1.0)
1470 C1MAX = (1.19*WWLT - 5.09)*TG - 146.*WWLT + 1786.
1471 C1MAX = MAX(MAX(C1MAX,ZY2),10.)
1472 !* Giard-Bazile formulation (resolution of a second order equation)
1473 !
1474 ZLY = LOG( C1MAX/10.)
1475 ZZA = - LOG( ZY2 /10.)
1476 ZZB = 2. * WWLT * ZLY
1477 ZDEL = 4. * (ZLY+ZZA) * ZLY * WWLT**2
1478 ZA = (-ZZB+SQRT(ZDEL)) / (2.*ZZA)
1479 ZB = ZA**2 / ZLY
1480 C1 = DS1*C1MAX * EXP(-(WG-ZA)**2/ZB)
1481 ENDIF
1483 C2 = C2R * W2HLF / (WSAT - W2HLF + 1.E-11)
1484 IF (W2HLF .GE. WSAT) THEN
1485 WEQ = WSAT
1486 ELSE
1487 WEQ = W2HLF - AS * WSAT * W2REL ** JP * &
1488 (1.0 - W2REL ** (8.0 * JP))
1489 ENDIF
1491 !.... The diffusion method in Sakaguchi and Zeng [2009] (JGR-Atmos.)
1492 CFNP1 = 1.0 + DTPBL * C2 * TAUINV * CRANKP
1493 CFN = C1 / (DENW * DS1) * (PG - EG) - C2 * TAUINV * &
1494 ((1.0 - CRANKP) * WG - WEQ) + WGNUDG/ DS1
1496 WGNEW = AMAX1((WG + DTPBL * CFN) / CFNP1, WRES ) !0.001) ! Limei 08/02/2017
1497 !-- NEW VALUES
1498 WG = AMIN1(WGNEW,WSAT)
1500 ENDIF !endif for ISNOW
1501 ENDIF !endif for XLAND
1503 END SUBROUTINE surfpx
1504 !-------------------------------------------------------------------------
1505 !-------------------------------------------------------------------------
1507 !-------------------------------------------------------------------------
1508 !-------------------------------------------------------------------------
1509 SUBROUTINE QFLUX (DENS1, QV1, TA1, RG, RAW, QSS, & ! in
1510 VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, & ! in
1511 WG, W2, WR, & ! in
1512 RSTMIN, WWLT, WFC, RSOIL, RINC, & ! in
1513 EG, ER, ETR, CQ4, RS, FASS, SIGG ) ! out
1515 !-------------------------------------------------------------------------
1516 !
1517 ! FUNCTION:
1518 ! THIS SUBROUTINE COMPUTES EVAPORATION FROM BARE SOIL (EG) AND FROM
1519 ! THE WET PART OF CANOPY (ER) AND TRANSPIRATION FROM THE DRY PART OF
1520 ! CANOPY (ETR).
1521 !
1522 ! REVISION HISTORY:
1523 ! A. Xiu Oct 2004 - adapted from the PX LSM in MM5 for the WRF system
1524 ! R. Gilliam Dec 2006 - Completed WRF V2.1.2 implementation
1525 !
1526 !-------------------------------------------------------------------------
1527 ! QFLUX ARGUMENT LIST:
1528 !-------------------------------------------------------------------------
1529 ! INPUTS:
1530 !-- DENS1 air density at first layer
1531 !-- QV1 air humidity at first layer
1532 !-- TA1 air temperature at first layer
1533 !-- RG shortwave radition reaching the ground
1534 !-- RAW RA+RB for moisture
1535 !-- QSS saturation mixing ratio at ground
1536 !-- VEGFRC vegetation coverage
1537 !-- ISNOW if snow on the ground
1538 !-- ISTI soil type
1539 !-- IFLAND if land (=1) or water (=2)
1540 !-- LAI leaf area index
1541 !-- BETAP
1542 !-- WG soil moisture at first soil layer
1543 !-- W2 soil moisture at root zone
1544 !-- WR Canopy water
1545 !
1546 ! OUTPUTS FROM QFLUX:
1547 !-- EG evaporation from ground (bare soil)
1548 !-- ER evaporation from canopy
1549 !-- ETR transpiration from vegetation
1550 !-- CQ4 CQ4 is used for the implicit calculation of TG in SURFACE
1551 !-- RS surface resistence
1552 !-- FASS parameter for soil moisture nudging
1553 !-- SIGG near ground-first layer moisture difference for evap or dew (1)
1554 !-------------------------------------------------------------------------
1555 !-------------------------------------------------------------------------
1557 IMPLICIT NONE
1559 ! DECLARATIONS - INTEGER
1560 INTEGER , INTENT(IN) :: ISTI
1562 ! DECLARATIONS - REAL
1563 REAL , INTENT(IN) :: ISNOW, IFLAND
1564 REAL , INTENT(IN) :: DENS1, QV1, TA1, RG, RAW, QSS, &
1565 VEGFRC, LAI, &
1566 WG, W2, WR, RSTMIN
1567 REAL , INTENT(INOUT) :: BETAP, RSOIL
1568 REAL, INTENT(IN) :: WWLT, WFC, RINC
1570 REAL , INTENT(OUT) :: EG, ER, ETR, CQ4, RS, FASS, SIGG
1572 !... Local Variables
1574 !... Real
1575 REAL :: WRMAX, DELTA, RADL, RADF, W2AVAIL, W2MXAV
1576 REAL :: FTOT, F1, F2, F3, F4
1577 REAL :: FSHELT, GS, GA, FX
1578 REAL :: PAR, F1MAX
1581 !... Parameters
1582 REAL, PARAMETER :: RSMAX = 5000.0 ! s/m
1583 REAL, PARAMETER :: FTMIN = 0.0000001 ! m/s
1584 REAL, PARAMETER :: F3MIN = 0.25
1586 !
1587 !... for water surface, no canopy evaporation and transpiration
1588 ER = 0.0
1589 ETR = 0.0
1590 CQ4 = 0.0
1592 !... GROUND EVAPORATION (DEPOSITION)
1593 IF (QSS .LT. QV1) RSOIL = 0.0
1594 EG = DENS1 * (QSS - QV1) * ((1.0 - VEGFRC) / (RAW + RSOIL) + VEGFRC/(RAW + RINC +RSOIL))
1596 !... CANOPY
1597 IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1598 WRMAX = 0.2E-3 * VEGFRC * LAI ! in unit m
1599 IF (WR .LE. 0.0) THEN
1600 DELTA = 0.0
1601 ELSE
1602 DELTA = WR / WRMAX ! referred to SiB model
1603 ENDIF
1605 IF (QSS .GE. QV1) THEN
1606 SIGG = DELTA
1607 ELSE
1608 SIGG = 1.0
1609 ENDIF
1611 ER = DENS1 * VEGFRC * SIGG * (QSS - QV1) / RAW
1612 ENDIF
1613 !!---------------------------------------------------------------------
1615 !-- TRANSPIRATION
1616 !!---------------------------------------------------------------------
1617 IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1619 !!!-RADIATION
1620 IF (RSTMIN .GT. 130.0) THEN
1621 F1MAX = 1.-0.02*LAI !Echer2015 Trees
1622 ELSE
1623 F1MAX = 1.-0.07*LAI !Echer2015 crops/grass
1624 ENDIF
1625 PAR = 0.45 * RG * 4.566 ! converted from W/m2 to umoles/m2/s (1/.219) Echer2015
1626 F1 = F1MAX*(1.0-exp(-0.0017*PAR)) !Echer2015
1627 F1 = AMAX1(F1,RSTMIN / RSMAX)
1629 !!!-SOIL MOISTURE
1630 W2AVAIL = W2 - WWLT
1631 W2MXAV = WFC - WWLT
1632 F2 = 1.0 / (1.0 + EXP(-5.0 * ( W2AVAIL / W2MXAV - &
1633 (W2MXAV / 3.0 + WWLT)))) ! according JP, 9/94
1635 !-AIR TEMP
1636 !... according to Avissar (1985) and AX 7/95
1637 IF (TA1 .LE. 302.15) THEN
1638 F4 = 1.0 / (1.0 + EXP(-0.41 * (TA1 - 282.05)))
1639 ELSE
1640 F4 = 1.0 / (1.0 + EXP(0.5 * (TA1 - 314.0)))
1641 ENDIF
1643 FTOT = LAI * F1 * F2 * F4
1644 ENDIF
1646 !!---------------------------------------------------------------------
1647 IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1648 FSHELT = 1.0 ! go back to NP89
1649 GS = FTOT / (RSTMIN * FSHELT)
1650 GA = 1.0 / RAW
1651 !-- Compute humidity effect according to RH at leaf surf
1652 F3 = 0.5 * (GS - GA + SQRT(GA * GA + GA * GS * &
1653 (4.0 * QV1 / QSS - 2.0) + GS * GS)) / GS
1654 F3 = AMIN1(AMAX1(F3,F3MIN),1.0)
1655 RS = 1.0 / (GS * F3)
1657 !--- Compute Assimilation factor for soil moisture nudging - jp 12/94
1658 !-- Note that the 30 coef is to result in order 1 factor for max
1659 IF (RG .LT. 0.00001) THEN ! do not nudge during night
1660 FX = 0.0
1661 ELSE
1662 FX = 30.0 * F1 * F4 * LAI / (RSTMIN * FSHELT)
1663 ENDIF
1665 FASS = FX
1666 ETR = DENS1 * VEGFRC * (1.0 - SIGG) * (QSS - QV1) / (RAW + RS)
1667 !..... CQ4 is used for the implicit calculation of TG in SURFACE
1668 CQ4 = DENS1 * VEGFRC * ((1.0 - SIGG) / (RAW + RS) + SIGG / RAW)
1669 ENDIF
1671 END SUBROUTINE qflux
1672 !------------------------------------------------------------------------------------------
1673 !------------------------------------------------------------------------------------------
1675 !------------------------------------------------------------------------------------------
1676 !------------------------------------------------------------------------------------------
1677 SUBROUTINE SMASS (ISTI, FASS, RG, VEGFRC, RA, & !in
1678 WWLT, WFC, & !in
1679 ALPH1, ALPH2, BET1, BET2, T2NUDF ) !out
1680 !------------------------------------------------------------------------------------------
1681 !------------------------------------------------------------------------------------------
1682 IMPLICIT NONE
1683 !------------------------------------------------------------------------------------------
1684 ! SMASS COMPUTES SOIL MOISTURE NUDGING COEFFICIENTS
1685 !------------------------------------------------------------------------------------------
1686 !
1687 !.........Arguments
1688 INTEGER, PARAMETER :: NSCAT = 16 ! max. soil types
1690 INTEGER, INTENT(IN) :: ISTI
1691 REAL, INTENT(IN) :: FASS, RG, VEGFRC, RA
1692 REAL, INTENT(IN) :: WWLT, WFC
1693 REAL, INTENT(OUT) :: ALPH1, ALPH2, BET1, BET2, T2NUDF
1696 !........Local variables
1698 !... Real
1699 REAL :: FBET, FALPH, FRA, FTEXT
1701 !... Parameters
1702 REAL, PARAMETER :: A1MAX = -10.E-5, A2MAX = 1.E-5 ! m/K, m for 6hr period
1703 REAL, PARAMETER :: B1MAX = -10.E-3, B2MAX = 1.E-3 ! m/K, m (Bouttier et al 1993)
1704 REAL, PARAMETER :: TASSI = 4.6296E-5 ! 1/6hr in 1/sec
1705 REAL, PARAMETER :: RAMIN = 10.0 ! 0.1 s/cm
1706 REAL, PARAMETER :: WFCX = 0.243 ! middle of WFC range
1707 REAL, PARAMETER :: WWLTX = 0.169 ! middle of WWLT range
1709 !
1710 FBET = FASS
1711 FALPH = RG / 1370.0
1712 !--TEXTURE FACTOR NORMALIZED BY LOAM (IST=5)
1713 FRA = RAMIN / RA ! scale by aerodynamic resistance
1714 FTEXT = TASSI * (WWLT + WFC) / (WWLTX + WFCX) * FRA
1715 ! write(6,*) ' ftot, fbet=',ftot, fbet,' ftext=',ftext/tassi
1716 !
1717 ALPH1 = A1MAX * FALPH * (1.0 - VEGFRC) * FTEXT
1718 ALPH2 = A2MAX * FALPH * (1.0 - VEGFRC) * FTEXT
1719 BET1 = B1MAX * FBET * VEGFRC * FTEXT
1720 BET2 = B2MAX * FBET * VEGFRC * FTEXT
1721 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
1723 END SUBROUTINE smass
1724 !------------------------------------------------------------------------------------------
1725 !------------------------------------------------------------------------------------------
1727 !------------------------------------------------------------------------------------------
1728 !------------------------------------------------------------------------------------------
1730 SUBROUTINE SOILPROP (SOILCBOT,WEIGHT, ITIMESTEP, MAVAIL, & ! IN
1731 PXLSM_SMOIS_INIT, & ! IN
1732 FWSAT,FWFC,FWWLT,FCLAY,FCSAND, & ! OUT
1733 FFMSAND,FB,FCGSAT, & ! OUT
1734 FJP,FAS,FC2R,FC1SAT,FWRES,FC3,ISTI, & ! OUT
1735 WG, W2 ) ! OUT
1736 !------------------------------------------------------------------------
1737 ! SOILPROP COMPUTES SOIL PARAMETERS FOR BOTH BOTTOM AND TOP LAYERS
1738 ! USING FRACTIONAL SOIL TYPE. A HARD CODED OPTION IS AVAILIABLE
1739 ! TO COMPUTE THE SOIL PARAMETERS USING FRACTIONAL INFORMATION, OR
1740 ! TO JUST USE SOIL PARAMETERS OF THE DOMINANT SOIL TYPE
1741 !------------------------------------------------------------------------
1742 !-- SOIL PARAMETERS ARE SPECIFIED BY SOIL TYPE:
1743 ! # SOIL TYPE WSAT WFC WWLT B CGSAT JP AS C2R C1SAT WRES
1744 ! _ _________ ____ ___ ____ ____ _____ ___ ___ ___ _____ ____
1745 ! 1 SAND .395 .135 .068 4.05 3.222 4 .387 3.9 .082 0.020
1746 ! 2 LOAMY SAND .410 .150 .075 4.38 3.057 4 .404 3.7 .098 0.035
1747 ! 3 SANDY LOAM .435 .195 .114 4.90 3.560 4 .219 1.8 .132 0.041
1748 ! 4 SILT LOAM .485 .255 .179 5.30 4.418 6 .105 0.8 .153 0.015
1749 ! 5 SILT .485 .255 .179 5.30 4.418 6 .105 0.8 .153 0.015
1750 ! 6 LOAM .451 .240 .155 5.39 4.111 6 .148 0.8 .191 0.027
1751 ! 7 SND CLY LM .420 .255 .175 7.12 3.670 6 .135 0.8 .213 0.068
1752 ! 8 SLT CLY LM .477 .322 .218 7.75 3.593 8 .127 0.4 .385 0.040
1753 ! 9 CLAY LOAM .476 .325 .250 8.52 3.995 10 .084 0.6 .227 0.075
1754 ! 10 SANDY CLAY .426 .310 .219 10.40 3.058 8 .139 0.3 .421 0.109
1755 ! 11 SILTY CLAY .482 .370 .283 10.40 3.729 10 .075 0.3 .375 0.056
1756 ! 12 CLAY .482 .367 .286 11.40 3.600 12 .083 0.3 .342 0.090
1757 !------------------------------------------------------------------------
1758 !------------------------------------------------------------------------------------------
1759 !------------------------------------------------------------------------------------------
1760 IMPLICIT NONE
1762 !.........Arguments
1763 INTEGER, PARAMETER :: NSCAT = 16 ! max. soil types
1764 INTEGER, PARAMETER :: NSCATMIN = 16 ! min soil types
1766 INTEGER, INTENT(IN) :: WEIGHT, ITIMESTEP, PXLSM_SMOIS_INIT
1767 REAL, INTENT(IN) :: MAVAIL
1768 REAL, DIMENSION(1:NSCAT), INTENT(IN) :: SOILCBOT
1769 REAL, INTENT(OUT) :: FWSAT,FWFC,FWWLT,FCLAY, &
1770 FCSAND,FFMSAND,FB,FCGSAT, &
1771 FJP,FAS,FC2R,FC1SAT,FWRES,FC3
1772 REAL, INTENT(INOUT) :: W2, WG
1775 INTEGER, INTENT(OUT) :: ISTI
1776 !........Local variables
1777 CHARACTER*4, AVCLASS
1778 CHARACTER*4, DIMENSION( 1: NSCAT ) :: TEXID
1779 !... Integer
1780 INTEGER:: S
1781 !... Real
1782 REAL:: TFRACBOT, CFRAC, SUMCSND, SUMFMSND, SUMCLY, &
1783 AVS, AVCS, AVFMS, AVC, AVSLT, &
1784 SSMPOT, DSMPOT ! saturated and air-dry soil matric potential
1786 REAL, DIMENSION( 1: NSCAT ) :: WSAT, WFC, WWLT, B, CGSAT, AS, &
1787 JP, C2R, C1SAT, WRES
1789 REAL, DIMENSION( 1: NSCATMIN ) :: CSAND,FMSAND, CLAY
1790 !.......... DATA statement for SOIL PARAMETERS for the 11 soil types
1791 !...........Follow Menut et al., 2013 JGR
1792 DATA CSAND /46.0,40.0,29.0, 0.0, 0.0, &
1793 0.0,29.0, 0.0, 0.0, 0.0, &
1794 0.0, 0.0, 0.0, 0.0,46.0, &
1795 0.0/
1796 DATA FMSAND /46.0,40.0,29.0,17.0,10.0, &
1797 43.0,29.0,10.0,32.0,52.0, &
1798 6.0,22.0,43.0,43.0,46.0, &
1799 32.0/
1801 DATA CLAY / 3.0, 4.0,10.0,13.0, 5.0, &
1802 18.0,27.0,34.0,34.0,42.0, &
1803 47.0,58.0,18.0,18.0, 3.0, &
1804 34.0/
1806 DATA TEXID/'Sand','Lsan','Sloa','Sill','Silt', &
1807 'Loam','Sclo','Sicl','Cllo','Sacl', &
1808 'Sicy','Clay','Ormt','Wate','Bedr', &
1809 'Othe'/
1810 !-- Removed soil parameter lookup table data since these parameters are now computed analytically (Noilhan and Mahfouf 1996)
1812 DSMPOT = -1.0E7 ! mm air-dry water matric potential
1813 !
1814 !-------------------------------Exicutable starts here--------------------
1816 ! Compute soil characteristics by sand (coarse and fine-medium) and clay fraction
1818 CFRAC = 0.0
1819 SUMCSND = 0.0
1820 SUMFMSND = 0.0
1821 SUMCLY = 0.0
1822 TFRACBOT = 0.0
1824 DO S = 1,NSCAT
1826 TFRACBOT = TFRACBOT + SOILCBOT(S)
1827 SUMCSND = SUMCSND + CSAND(S) * SOILCBOT(S)
1828 SUMFMSND = SUMFMSND + FMSAND(S) * SOILCBOT(S)
1829 SUMCLY = SUMCLY + CLAY(S) * SOILCBOT(S)
1831 IF(SOILCBOT(S).GE.CFRAC) THEN ! Find Dominant Category and fraction
1832 ISTI=S
1833 CFRAC=SOILCBOT(S)
1834 ENDIF
1836 ENDDO
1839 IF(TFRACBOT.GT.0.001) THEN
1840 AVCS = SUMCSND / TFRACBOT
1841 AVFMS = SUMFMSND / TFRACBOT
1842 AVS = AVCS + AVFMS
1844 AVC = SUMCLY / TFRACBOT
1845 AVSLT = 100.0 - AVS - AVC
1847 IF(AVS.GT.(85.+ 0.5*AVC)) THEN
1848 AVCLASS= 'Sand'
1849 ISTI = 1
1850 ELSE IF(AVS.GT.(70.+ AVC)) THEN
1851 AVCLASS= 'Lsan'
1852 ISTI = 2
1853 ELSE IF((AVC.LT.20..AND.AVS.GT.52.) &
1854 .OR.(AVC.LE.7.5.AND.AVSLT.LT.50.)) THEN
1855 AVCLASS= 'Sloa'
1856 ISTI = 3
1857 ELSE IF(AVC.LT.35..AND.AVS.GT.45..AND.AVSLT.LT.28.) THEN
1858 AVCLASS= 'Sclo'
1859 ISTI = 7
1860 ELSE IF(AVC.GE.35..AND.AVS.GT.45.) THEN
1861 AVCLASS = 'Sacl'
1862 ISTI = 10
1863 ELSE IF(AVC.LT.27.0.AND.AVSLT.LT.50.) THEN
1864 AVCLASS= 'Loam'
1865 ISTI = 6
1866 ELSE IF(AVC.LT.12..AND.AVSLT.GT.80.) THEN
1867 AVCLASS = 'Silt'
1868 ISTI = 5
1869 ELSE IF(AVC.LT.27.) THEN
1870 AVCLASS = 'Sill'
1871 ISTI = 4
1872 ELSE IF(AVC.LT.40..AND.AVS.GT.20.) THEN
1873 AVCLASS = 'Cllo'
1874 ISTI = 9
1875 ELSE IF(AVC.LT.40.) THEN
1876 AVCLASS = 'Sicl'
1877 ISTI = 8
1878 ELSE IF(AVSLT.GE.40.) THEN
1879 AVCLASS = 'Sicy'
1880 ISTI = 11
1881 ELSE
1882 AVCLASS = 'Clay'
1883 ISTI = 12
1884 ENDIF
1885 ELSE
1886 ! set no soil to 9 - clay loam
1887 ISTI = 9
1888 AVCLASS = TEXID(ISTI)
1890 AVCS = CSAND(ISTI)
1891 AVFMS = FMSAND(ISTI)
1892 AVS = AVCS + AVFMS
1894 AVC = CLAY(ISTI)
1895 AVSLT = 100.0 - AVS - AVC
1897 ENDIF
1899 FCSAND = AVCS
1900 FFMSAND = AVFMS
1901 FCLAY = AVC
1903 ! Continous formulation of secondary soil parmeters (Noilhan and Mahfouf 1996)
1904 FWSAT = (-1.08 * AVS + 494.305) * 1.0E-3
1905 FWWLT = 37.1342E-3 * SQRT(AVC)
1906 FWFC = 89.0467E-3 * AVC**0.3496
1907 FB = 0.137 * AVC + 3.501
1908 FCGSAT= -1.557E-2 * AVS - 1.441E-2 * AVC + 4.7021
1909 FC1SAT= (5.58 * AVC + 84.88) * 1.0E-2
1910 FC2R = 13.815 * AVC**(-0.954)
1911 FC3 = 5.327 * AVC **(-1.043)
1912 FAS = 732.42E-3 * AVC **(-0.539)
1913 FJP = 0.134 * AVC + 3.4
1914 FWRES = 0.00123 * AVC - 0.00066 * AVSLT + 0.0405 !J. Pleim fitted function
1915 FWRES = AMAX1(FWRES, 0.01) !L. Ran set minimum
1917 ! Compute W2 using soil moisture availiability if pxlsm_smois_init (in namelist) is not zero
1918 IF (ITIMESTEP .EQ. 1 .AND. PXLSM_SMOIS_INIT .GT. 0) THEN
1919 WG = FWWLT + (0.5*(FWSAT+FWFC) - FWWLT) * MAVAIL
1920 W2 = FWWLT + (0.5*(FWSAT+FWFC) - FWWLT) * MAVAIL
1921 ENDIF
1923 END SUBROUTINE soilprop
1924 !------------------------------------------------------------------------------------------
1925 !------------------------------------------------------------------------------------------
1927 !------------------------------------------------------------------------------------------
1928 !------------------------------------------------------------------------------------------
1929 SUBROUTINE PXSNOW (ITIMESTEP, ASNOW, CSNOW, SNOW, &
1930 SNOWH, SNUP, &
1931 ALB, SNOALB, SHDFAC, SHDMIN, &
1932 HC_SNOW, SNOW_FRA, SNOWC, SNOWALB)
1933 !------------------------------------------------------------------------------------------
1934 !------------------------------------------------------------------------------------------
1935 ! Pleim-Xiu LSM Simple Snow model
1936 !---------------------------------------------------------------------------------------------------
1937 ! ITIMESTEP -- Model time step index
1938 ! ASNOW -- Analyzed snow water equivalent (mm)
1939 ! CSNOW -- Current snow water equivalent (mm)
1940 ! SNOW -- Snow water equivalent in model (mm)
1941 ! SNOWH -- Physical snow depth (m)
1942 ! SNUP -- Physical snow depth (landuse dependent) where when below, snow cover is fractional
1943 !
1944 ! HC_SNOW -- Heat capacity of snow as a function of depth
1945 ! SNOW_FRA-- Factional snow area
1946 ! IFSNOW -- Snow flag
1947 ! SNOWALB -- Snow albedo
1948 !---------------------------------------------------------------------------------------------------
1950 IMPLICIT NONE
1952 !--- Arguments
1953 REAL, PARAMETER :: W2SN_CONV = 10.0
1954 REAL, PARAMETER :: CS_SNOWPACK = 2092.0
1955 REAL, PARAMETER :: SALP = 2.6
1956 !-----------------------------------------------
1957 INTEGER, INTENT(IN) :: ITIMESTEP
1958 REAL, INTENT(IN) :: ASNOW, CSNOW, SNUP, ALB, SNOALB, SHDFAC, SHDMIN
1959 REAL, INTENT(INOUT) :: SNOW, SNOWH, SNOWC
1960 REAL, INTENT(OUT) :: HC_SNOW, SNOW_FRA, SNOWALB
1961 !------------------------------------------------------------------------------------
1964 !-----------------------------------------------
1965 ! Local variables
1966 !... Real
1967 REAL:: CONV_WAT2SNOW, CSNOWW, RHO_SNOWPACK, &
1968 LIQSN_RATIO, SNEQV, RSNOW
1969 !-----------------------------------------------
1971 SNEQV=ASNOW*0.001 ! Snow water in meters
1972 RHO_SNOWPACK = 450 ! Snowpack density (kg/m3), this should be computed in the future
1973 LIQSN_RATIO = DENW/RHO_SNOWPACK ! Ratio of water density and snowpack density
1975 CONV_WAT2SNOW = LIQSN_RATIO/1000 ! Conversion factor for snow liquid equiv. (mm) to snowpack (m)
1977 SNOW = ASNOW ! Set snow water to analysis value for now, implement a nudging scheme later
1978 SNOWH= SNOW * CONV_WAT2SNOW ! Conversion of snow water (mm) to snow depth (m)
1981 ! Is snow cover now present. The limit is 0.45 mm of water eqivalent or about 2 inches snow depth
1982 SNOWC = 0.0
1983 IF (SNOWH .GT. 0.005) SNOWC = 1.0
1985 HC_SNOW = RHO_SNOWPACK * CS_SNOWPACK * SNOWH
1987 IF (SNEQV < SNUP) THEN
1988 RSNOW = SNEQV / SNUP
1989 SNOW_FRA = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
1990 ELSE
1991 SNOW_FRA = 1.0
1992 END IF
1994 SNOWC = SNOW_FRA
1996 SNOWALB = ALB + SNOWC*(SNOALB-ALB)
1998 END SUBROUTINE pxsnow
1999 !------------------------------------------------------------------------------------------
2000 !------------------------------------------------------------------------------------------
2002 END MODULE module_sf_pxlsm