Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_pxlsm.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 77.5 KB

Line
1	MODULE module_sf_pxlsm
2
3	USE module_model_constants
4
5	INTEGER, PARAMETER :: NSOLD=20
6	REAL, PARAMETER :: RD = 287.04, CPD = 1004.67, &
7	CPH2O = 4.218E+3, CPICE = 2.106E+3, &
8	LSUBF = 3.335E+5, SIGMA = 5.67E-8, &
9	ROVCP = RD / CPD
10
11	REAL, PARAMETER :: CRANKP = 0.5 ! CRANK-NIC PARAMETER
12	REAL, PARAMETER :: RIC = 0.25 ! critical Richardson number
13	REAL, PARAMETER :: DENW = 1000.0 ! water density in KG/M3
14	REAL, PARAMETER :: TAUINV = 1.0 / 86400.0 ! 1/1DAY(SEC)
15	REAL, PARAMETER :: T2TFAC = 1.0 / 10.0 ! Bottom soil temp response factor
16	REAL, PARAMETER :: PI = 3.1415926
17	REAL, PARAMETER :: PR0 = 0.95
18	REAL, PARAMETER :: CZO = 0.032
19	REAL, PARAMETER :: OZO = 1.E-4
20
21	CONTAINS
22	!
23
24	!-------------------------------------------------------------------------
25	SUBROUTINE pxlsm(U3D, V3D, DZ8W, QV3D,T3D,TH3D, RHO, &
26	PSFC, GSW, GLW, RAINBL, EMISS, &
27	ITIMESTEP, NSOIL, DT, ANAL_INTERVAL, &
28	XLAND, XICE, ALBBCK, ALBEDO, SNOALB, &
29	SMOIS, TSLB, MAVAIL, TA2, QA2, &
30	ZS,DZS, PSIH, &
31	LANDUSEF,SOILCTOP,SOILCBOT,VEGFRA,VEGF_PX, &
32	ISLTYP,RA,RS,LAI,NLCAT,NSCAT, &
33	HFX,QFX,LH,TSK,SST,ZNT,CANWAT, &
34	GRDFLX,SHDMIN,SHDMAX, &
35	SNOWC,PBLH,RMOL,UST,CAPG,DTBL, &
36	T2_NDG_OLD, T2_NDG_NEW, &
37	Q2_NDG_OLD, Q2_NDG_NEW, &
38	SN_NDG_OLD, SN_NDG_NEW, SNOW, SNOWH,SNOWNCV,&
39	T2OBS, Q2OBS, PXLSM_SMOIS_INIT, &
40	PXLSM_SOIL_NUDGE, &
41	ids,ide, jds,jde, kds,kde, &
42	ims,ime, jms,jme, kms,kme, &
43	its,ite, jts,jte, kts,kte )
44
45	!*************************************************************************
46	! THIS MODULE CONTAINS THE PLEIM-XIU LAND-SURFACE MODEL (PX-LSM).
47	! IT IS DESIGNED TO SIMULATE CHARACTERISTICS OF THE LAND SURFACE AND
48	! VEGETATION AND EXCHANGE WITH THE PLANETARY BOUNDARY LAYER (PBL). THE
49	! SOIL MOISTURE MODEL IS BASED ON THE ISBA SCHEME DEVELOPED BY NOILHAN
50	! AND PLANTON (1989) AND JACQUEMIN AND NOILHAN (1990) AND INCLUDES
51	! PROGNOSTIC EQUATIONS FOR SOIL MOISTURE AND SOIL TEMPERATURE IN TWO
52	! LAYERS (1 CM AND 1 M) AS WELL AS CANOPY WATER CONTENT. SURFACE
53	! MOISTURE FLUXES ARE MODELED BY 3 PATHWAYS: SOIL EVAPORATION, CANOPY
54	! EVAPORATION, AND VEGETATIVE EVAPOTRANSPIRRATION.
55	! EVAPOTRANSPIRATION DIRECTLY FROM THE ROOT ZONE SOIL LAYER IS MODELED
56	! VIA A CANOPY RESISTANCE ANALOG ALGORITHM WHERE STOMATAL CONDUCTANCE
57	! IS CONTROLLED BY SOLAR RADIATION, AIR TEMPERATURE, AIR HUMIDITY, AND
58	! ROOT ZONE SOIL MOISTURE. REQUIRED VEGETATION CHARACTERISTICS DERIVED
59	! FROM THE USGS LANDUSE DATA INCLUDE: LEAF AREA INDEX, FRACTIONAL VEGETATION
60	! COVERAGE, ROUGHNESS LENGTH, AND MINIMUM STOMATAL RESISTANCE. AN INDIRECT
61	! NUDGING SCHEME ADJUSTS SOIL MOISTURE ACCORDING TO DIFFERENCES BETWEEN
62	! MODELED TEMPERATURE AND HUMIDITY AND ANALYSED SURFACE FIELDS.
63	!
64	! References:
65	! Pleim and Xiu, 1995: Development and testing of a surface flux and planetary
66	! boundary layer model for application in mesoscale models.
67	! J. Appl. Meteoro., Vol. 34, 16-32.
68	! Xiu and Pleim, 2001: Development of a land surface model. Part I: Application
69	! in a mesoscale meteorological model. J. Appl. Meteoro.,
70	! Vol. 40, 192-209.
71	! Pleim and Xiu, 2003: Development of a land surface model. Part II: Data
72	! assimilation. J. Appl. Meteoro., Vol. 42, 1811-1822.
73	!
74	! REVISION HISTORY:
75	! AX 4/2005 - developed WRF version based on the MM5 PX LSM
76	!************************************************************************
77	! ARGUMENT LIST:
78	!
79	!... Inputs:
80	!-- U3D 3D u-velocity interpolated to theta points (m/s)
81	!-- V3D 3D v-velocity interpolated to theta points (m/s)
82	!-- DZ8W dz between full levels (m)
83	!-- QV3D 3D mixing ratio
84	!-- T3D Temperature (K)
85	!-- TH3D Theta (K)
86	!-- RHO 3D dry air density (kg/m^3)
87
88	!-- PSFC surface pressure (Pa)
89	!-- GSW downward short wave flux at ground surface (W/m^2)
90	!-- GLW downward long wave flux at ground surface (W/m^2)
91	!-- RAINBL Timestep rainfall
92	!-- EMISS surface emissivity (between 0 and 1)
93
94	!-- ITIMESTEP time step number
95	!-- NSOIL number of soil layers
96	!-- DT time step (second)
97	!-- ANAL_INTERVAL Interval of analyses used for soil moisture and temperature nudging
98
99	!-- XLAND land mask (1 for land, 2 for water)
100	!-- XICE Sea ice
101	!-- ALBBCK Background Albedo
102	!-- ALBEDO surface albedo with snow cover effects
103	!-- SNOALB Albedo of snow
104
105	!-- SMOIS total soil moisture content (volumetric fraction)
106	!-- TSLB soil temp (k)
107	!-- MAVAIL Moisture availibility of soil
108	!-- TA2 2-m temperature
109	!-- QA2 2-m mixing ratio
110
111	!-- SVPT0 constant for saturation vapor pressure (K)
112	!-- SVP1 constant for saturation vapor pressure (kPa)
113	!-- SVP2 constant for saturation vapor pressure (dimensionless)
114	!-- SVP3 constant for saturation vapor pressure (K)
115
116	!-- ZS depths of centers of soil layers
117	!-- DZS thicknesses of soil layers
118	!-- PSIH similarity stability function for heat
119
120	!-- LANDUSEF Landuse fraction
121	!-- SOILCTOP Top soil fraction
122	!-- SOILCBOT Bottom soil fraction
123	!-- VEGFRA Vegetation fraction
124	!-- VEGF_PX Veg fraction recomputed and used by PX LSM
125	!-- ISLTYP Soil type
126
127	!-- RA Aerodynamic resistence
128	!-- RS Stomatal resistence
129	!-- LAI Leaf area index (weighted according to fractional landuse)
130	!-- NLCAT Number of landuse categories
131	!-- NSCAT Number of soil categories
132
133	!-- HFX net upward heat flux at the surface (W/m^2)
134	!-- QFX net upward moisture flux at the surface (kg/m^2/s)
135	!-- LH net upward latent heat flux at surface (W/m^2)
136	!-- TSK surface skin temperature (K)
137	!-- SST sea surface temperature
138	!-- ZNT rougness length
139	!-- CANWAT Canopy water (mm)
140
141	!-- GRDFLX Ground heat flux
142	!-- SFCEVP Evaportation from surface
143	!-- SHDMIN Minimum annual vegetation fraction for each grid cell
144	!-- SHDMAX Maximum annual vegetation fraction for each grid cell
145
146	!-- SNOWC flag indicating snow coverage (1 for snow cover)
147	!-- PBLH PBL height (m)
148	!-- RMOL 1/L Reciprocal of Monin-Obukhov length
149	!-- UST u* in similarity theory (m/s)
150	!-- CAPG heat capacity for soil (J/K/m^3)
151	!-- DTBL time step of boundary layer calls
152
153	!-- T2_NDG_OLD Analysis temperature prior to current time
154	!-- T2_NDG_NEW Analysis temperature ahead of current time
155	!-- Q2_NDG_OLD Analysis mixing ratio prior to current time
156	!-- Q2_NDG_NEW Analysis mixing ratio ahead of current time
157
158	!-- SN_NDG_OLD Analysis snow water prior to current time
159	!-- SN_NDG_NEW Analysis snow water ahead of current time
160	!-- SNOW Snow water equivalent
161	!-- SNOWH Physical snow depth
162	!-- SNOWNCV Time step accumulated snow
163
164	!-- T2OBS Analysis temperature interpolated from prior and next in time analysese
165	!-- Q2OBS Analysis moisture interpolated from prior and next in time analysese
166	!-- PXLSM_SMOIS_INIT Flag to intialize deep soil moisture to a value derived from moisture availiability.
167	!-- PXLSM_SOIL_NUDGE Flag to use soil moisture and temperature nudging in the PX LSM
168	! This is typically done for the first simulation.
169
170	!-- ids start index for i in domain
171	!-- ide end index for i in domain
172	!-- jds start index for j in domain
173	!-- jde end index for j in domain
174	!-- kds start index for k in domain
175	!-- kde end index for k in domain
176	!-- ims start index for i in memory
177	!-- ime end index for i in memory
178	!-- jms start index for j in memory
179	!-- jme end index for j in memory
180	!-- kms start index for k in memory
181	!-- kme end index for k in memory
182	!-- jts start index for j in tile
183	!-- jte end index for j in tile
184	!-- kts start index for k in tile
185	!-- kte end index for k in tile
186	!... Outputs:
187
188	IMPLICIT NONE
189
190	!.......Arguments
191	! DECLARATIONS - INTEGER
192	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
193	ims,ime, jms,jme, kms,kme, &
194	its,ite, jts,jte, kts,kte
195
196	INTEGER, INTENT(IN) :: NSOIL, ITIMESTEP, NLCAT, NSCAT, &
197	ANAL_INTERVAL, PXLSM_SMOIS_INIT, PXLSM_SOIL_NUDGE
198
199	REAL, INTENT(IN ) :: DT,DTBL
200
201	INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ISLTYP
202
203	! DECLARATIONS - REAL
204	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U3D, V3D, RHO, &
205	T3D, TH3D, DZ8W, QV3D
206
207	REAL, DIMENSION(1:NSOIL), INTENT(IN)::ZS,DZS
208	REAL, DIMENSION( ims:ime , 1:NSOIL, jms:jme ), INTENT(INOUT) :: SMOIS, TSLB
209
210	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: RA, RS, LAI, ZNT
211	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT):: GRDFLX, TSK, TA2, QA2
212
213	REAL, DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN):: LANDUSEF
214	REAL, DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN):: SOILCTOP, SOILCBOT
215
216
217	REAL, DIMENSION( ims:ime, jms:jme ), &
218	INTENT(IN) :: PSFC, GSW, GLW, RAINBL, &
219	EMISS, SNOALB, &
220	ALBBCK, SHDMIN, SHDMAX, &
221	PBLH, RMOL, SNOWNCV, &
222	UST, MAVAIL, SST
223	REAL, DIMENSION( ims:ime, jms:jme ), &
224	INTENT(IN) :: T2_NDG_OLD, T2_NDG_NEW, &
225	Q2_NDG_OLD, Q2_NDG_NEW, &
226	SN_NDG_OLD, SN_NDG_NEW
227
228	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: T2OBS, Q2OBS
229
230	REAL, DIMENSION( ims:ime, jms:jme ), &
231	INTENT(INOUT) :: CAPG,CANWAT, QFX, HFX, LH, &
232	PSIH,VEGFRA, VEGF_PX, SNOW, &
233	SNOWH, SNOWC, ALBEDO, XLAND, XICE
234
235
236	LOGICAL :: radiation
237
238	!-------------------------------------------------------------------------
239	! ---------- Local Variables --------------------------------
240
241
242	!---- PARAMETERS
243	INTEGER, PARAMETER :: NSTPS = 11 ! max. soil types
244	REAL, PARAMETER :: DTPBLX = 40.0 ! Max PX timestep = 40 sec
245
246
247	!---- INTEGERS
248	INTEGER, DIMENSION( 1: NSTPS ) :: JP
249	INTEGER:: J, I, NS, NUDGE, ISTI, WEIGHT
250	INTEGER:: NTSPS, IT
251
252	!---- REALS
253	REAL, DIMENSION( ims:ime, jms:jme ) :: XLAI, XLAIMN, RSTMIN, &
254	XVEG, XVEGMN, XSNUP, &
255	XALB
256
257	REAL, DIMENSION( ims:ime, jms:jme ) :: RADNET, EG, ER, ETR, QST
258
259	REAL:: SFCPRS,TA1,DENS1,QV1,ZLVL,SOLDN,LWDN, &
260	EMISSI,PRECIP,THETA1,VAPPRS,QSBT, &
261	WG,W2,WR,TG,T2,USTAR,MOLX,Z0, &
262	RAIR,CPAIR,IFLAND,ISNOW, &
263	ES,QSS,BETAP, &
264	RH2_OLD, RH2_NEW, T2_OLD, T2_NEW, &
265	CORE, CORB, TIME_BETWEEN_ANALYSIS, &
266	G1000, ALN10,RH2OBS, HU, SNOBS, &
267	FWSAT,FWFC,FWWLT,FB,FCGSAT,FJP,FAS, &
268	FSEAS, T2I, HC_SNOW, SNOW_FRA,SNOWALB, &
269	QST12,ZFUNC,ZF1,ZA2,QV2, DT_FDDA,CURTIME, &
270	FC2R,FC1SAT, DTPBL
271
272	!-------------------------------------------------------------------------
273	!-------------------------------Executable starts here--------------------
274	!
275	ALN10 = ALOG(10.0)
276	G1000 = g*1.0E-3 ! G/1000
277	WEIGHT = 0
278	DT_FDDA = ANAL_INTERVAL * 1.0 ! Convert DT of Analysis to real
279
280	! Compute vegetation and land-use characteristics by land-use fraction weighting
281	! These parameters include LAI, VEGF, ZNT, ALBEDO, RS, etc.
282	CALL VEGELAND(LANDUSEF, VEGFRA, SHDMIN, SHDMAX, &
283	SOILCTOP, SOILCBOT, NLCAT, NSCAT, &
284	ZNT,XLAI,XLAIMN,RSTMIN,XVEG,XVEGMN,XSNUP, &
285	XLAND, XALB, &
286	ids,ide, jds,jde, kds,kde, &
287	ims,ime, jms,jme, kms,kme, &
288	its,ite, jts,jte, kts,kte)
289	!-----------------------------------------------------------------------------------
290	!--- Compute time relatve to old and new analysis time for timestep interpolation
291	CURTIME = (ITIMESTEP-1) * DT
292	TIME_BETWEEN_ANALYSIS = MOD(CURTIME,DT_FDDA)
293	IF (TIME_BETWEEN_ANALYSIS .EQ. 0.0) THEN
294	CORB = 1.0
295	CORE = 0.0
296	ELSE
297	CORE = TIME_BETWEEN_ANALYSIS / DT_FDDA
298	CORB = 1.0 - CORE
299	ENDIF
300	!-----------------------------------------------------------------------------------
301
302	!-----------------------------------------------------------------------------------
303	! Main loop over individual grid cells
304	DO J = jts,jte !-- J LOOP
305	DO I = its,ite !-- I LOOP
306
307	IFLAND = XLAND(I,J)
308
309	! Compute soil properties via weighting of fractional components
310	IF (IFLAND .LT. 1.5 ) THEN
311	CALL SOILPROP (SOILCBOT(I,:,J), WEIGHT, &
312	ITIMESTEP, MAVAIL(I,J), &
313	PXLSM_SMOIS_INIT, &
314	FWSAT,FWFC,FWWLT,FB,FCGSAT, &
315	FJP,FAS,FC2R,FC1SAT,ISTI,SMOIS(I,2,J) )
316	ISLTYP(I,J) = ISTI
317	ENDIF
318
319	!-- Variables Sub. SURFPX needs
320	SFCPRS = PSFC(i,j) / 1000.0 ! surface pressure in cb
321	TA1 = T3D(i,1,j) ! air temperature at first layer
322	DENS1 = RHO(I,1,J) ! air density at first layer
323	QV1 = QV3D(i,1,j) ! water vapor mixing ratio at first layer
324	QV2 = QV3D(i,2,j)
325	ZLVL = 0.5 * DZ8W(i,1,j) ! thickness of lowest half level
326	ZF1 = DZ8W(i,1,j)
327	ZA2 = ZF1 + 0.5 * DZ8W(i,2,j)
328
329	LWDN = GLW(I,J) ! longwave radiation
330	EMISSI = EMISS(I,J) ! emissivity
331	PRECIP = MAX ( 1.0E-3*RAINBL(i,j)/DTBL,0.0) ! accumulated precip. rate during DT (=dtpbl)
332	! convert RAINBL from mm to m for PXLSM
333	WR = 1.0E-3*CANWAT(I,J) ! convert CANWAT from mm to m for PXLSM
334	THETA1 = TH3D(i,1,j) ! potential temp at first layer
335	SNOBS = SNOW(I,J) ! Set snow cover to existing model value
336	! this is overwritten below if snow analysis is availiable
337	! otherwise snow cover remains constant through simulation
338
339	IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN
340	!-- 2 m Temp and RH for Nudging
341	T2_OLD = T2_NDG_OLD(I,J)
342	T2_NEW = T2_NDG_NEW(I,J)
343	VAPPRS = SVP1 * EXP(SVP2 * (T2_OLD - SVPT0) / ( T2_OLD - SVP3))
344	QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS)
345	RH2_OLD = Q2_NDG_OLD(I,J) / QSBT
346	VAPPRS = SVP1 * EXP(SVP2 * (T2_NEW - SVPT0) / (T2_NEW - SVP3))
347	QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS)
348	RH2_NEW = Q2_NDG_NEW(I,J) / QSBT
349	RH2OBS = CORB * RH2_OLD + CORE * RH2_NEW
350	T2OBS(I,J) = CORB * T2_OLD + CORE * T2_NEW
351	Q2OBS(I,J) = CORB * Q2_NDG_OLD(I,J) + CORE * Q2_NDG_NEW(I,J)
352	SNOBS = CORB * SN_NDG_OLD(I,J) + CORE * SN_NDG_NEW(I,J)
353	ENDIF
354
355	USTAR = MAX(UST(I,J),0.005)
356	IF (IFLAND .GT. 1.5) THEN ! if over water
357	ZNT(I,J) = CZO * UST(I,J) * UST(I,J) / G + OZO
358	ENDIF
359	Z0 = ZNT(I,J)
360	CPAIR = CPD * (1.0 + 0.84 * QV1) ! J/(K KG)
361
362	! Compute fractional snow area and snow albedo
363	CALL PXSNOW (ITIMESTEP, SNOBS, SNOWNCV(I,J), SNOW(I,J), &
364	SNOWH(I,J), XSNUP(I,J), XALB(i,j), &
365	SNOALB(I,J),VEGF_PX(I,J), SHDMIN(I,J), &
366	HC_SNOW, SNOW_FRA, SNOWC(I,J), ALBEDO(I,J) )
367
368	!-------------------------------------------------------------
369	! Sea Ice from analysis and water cells that are very cold, but more than 50% water
370	! are converted to ice/snow for more reasonable treatment.
371	IF( (XICE(I,J).GE.0.5) .OR. &
372	(SST(I,J).LE.270.0.AND.XLAND(I,J).GT.1.50) ) THEN
373	XLAND(I,J) = 1.0
374	IFLAND = 1.0
375	ZNT(I,J) = 0.001 ! Ice
376	SMOIS(I,1,J) = 1.0 ! FWSAT
377	SMOIS(I,2,J) = 1.0 ! FWSAT
378	XICE(I,J) = 1.0
379	ALBEDO(I,J) = 0.7
380	SNOWC(I,J) = 1.0
381	SNOW_FRA = 1.0
382	VEGF_PX(I,J) = 0.0
383	LAI(I,J) = 0.0
384	ENDIF
385	!-------------------------------------------------------------
386
387	!-------------------------------------------------------------
388	!-- Note that when IFGROW = 0 is selected in Vegeland then max and min
389	!-- LAI and Veg are the same
390	T2I = TSLB(I,2,J)
391	! FSEAS = AMAX1(1.0 - 0.0016 * (298.0 - T2I) ** 2,0.0) ! BATS
392	FSEAS = AMAX1(1.0 - 0.015625 * (290.0 - T2I) ** 2,0.0) ! JP97
393	IF (T2I .GE. 290.0) FSEAS = 1.0
394
395	LAI(I,J) = XLAIMN(I,J) + FSEAS*(XLAI(I,J) - XLAIMN(I,J))
396	VEGF_PX(I,J) = XVEGMN(I,J) + FSEAS*(XVEG(I,J) - XVEGMN(I,J))
397
398	! Ensure veg algorithms not used for water
399	IF (IFLAND .GT. 1.5) THEN
400	VEGF_PX(I,J) = 0.0
401	ENDIF
402	!-------------------------------------------------------------
403
404
405	SOLDN = GSW(I,J) / (1.0 - ALBEDO(I,J)) ! downward shortwave radiaton
406	ISNOW = SNOWC(I,J)
407
408
409	NUDGE=PXLSM_SOIL_NUDGE
410	IF ( J .LE. 2 .OR. J .GE. (jde-1) ) NUDGE=0
411	IF ( I .LE. 2 .OR. I .GE. (ide-1) ) NUDGE=0
412
413	IF ( RMOL(I,J) .GT. 0.0 ) THEN
414	MOLX = AMIN1(1/RMOL(I,J),1000.0)
415	ELSE IF ( RMOL(I,J) .LT. 0.0 ) THEN
416	MOLX = AMAX1(1/RMOL(I,J),-1000.0)
417	ELSE
418	MOLX = 1000
419	ENDIF
420
421	ZFUNC = ZF1 * (1.0 - ZF1 / MAX(100.,PBLH(I,J))) ** 2
422	QST12 = KARMAN * ZFUNC*(QV2-QV1) &
423	/ (ZA2-ZLVL)
424
425	!--- LSM sub-time loop too prevent dt > 40 sec
426	NTSPS = INT(DT / (DTPBLX + 0.000001) + 1.0)
427	DTPBL = DT / NTSPS
428
429	DO IT=1,NTSPS
430
431	!... SATURATION VAPOR PRESSURE (MB)
432	IF ( TSLB(I,1,J) .LE. SVPT0 ) THEN ! For ground that is below freezing
433	ES = SVP1 * EXP(22.514 - 6.15E3 / TSLB(I,1,J)) ! cb
434	ELSE
435	ES = SVP1 * EXP(SVP2 * (TSLB(I,1,J) - SVPT0) / (TSLB(I,1,J) - SVP3))
436	ENDIF
437	QSS = ES * 0.622 / (SFCPRS - ES)
438
439	!... beta method, Lee & Pielke (JAM,May1992)
440	BETAP = 1.0
441	IF (IFLAND .LT. 1.5 .AND. ISNOW .LT. 0.5 .AND. SMOIS(I,1,J) .LE. FWFC) THEN
442	BETAP = 0.25 * (1.0 - COS(SMOIS(I,1,J) / FWFC * PI)) ** 2
443	ENDIF
444
445	CALL SURFPX (DTPBL, IFLAND, SNOWC(I,J), NUDGE, XICE(I,J), & !in
446	SOLDN, GSW(I,J), LWDN, EMISSI, ZLVL, & !in
447	MOLX, Z0, USTAR, & !in
448	SFCPRS, DENS1, QV1, QSS, TA1, & !in
449	THETA1, PRECIP, & !in
450	CPAIR, PSIH(I,J), & !in
451	RH2OBS,T2OBS(I,J), & !in
452	VEGF_PX(I,J), ISTI, LAI(I,J), BETAP, & !in
453	RSTMIN(I,J), HC_SNOW, SNOW_FRA, & !in
454	FWWLT, FWFC, FCGSAT, FWSAT, FB, & !in
455	FC1SAT,FC2R,FAS,FJP,DZS(1),DZS(2),QST12, & !in
456	RADNET(I,J), GRDFLX(I,J), HFX(I,J), QFX(I,J), LH(I,J), & !out
457	EG(I,J), ER(I,J), ETR(I,J), & !out
458	QST(I,J), CAPG(I,J), RS(I,J), RA(I,J), & !out
459	TSLB(I,1,J), TSLB(I,2,J), & !out
460	SMOIS(I,1,J), SMOIS(I,2,J), WR, TA2(I,J), QA2(I,J) )
461
462	END DO ! Time internal PX time loop
463
464	TSK(I,J)= TSLB(I,1,J) ! Skin temp set to 1 cm soil temperature in PX for now
465	CANWAT(I,J) = WR * 1000. ! convert WR back to mm for CANWAT
466
467	ENDDO ! END MIAN I LOOP
468	ENDDO ! END MAIN J LOOP
469
470	!------------------------------------------------------
471	END SUBROUTINE pxlsm
472	!------------------------------------------------------
473	!***************************************************************************
474	!***************************************************************************
475	SUBROUTINE VEGELAND( landusef, vegfra, & ! in
476	shdmin, shdmax, &
477	soilctop, soilcbot, nlcat, nscat, &
478	znt, xlai, xlaimn, rstmin, xveg, xvegmn, &
479	xsnup, xland, xalb, &
480	ids,ide, jds,jde, kds,kde, &
481	ims,ime, jms,jme, kms,kme, &
482	its,ite, jts,jte, kts,kte )
483	!***************************************************************************
484	!
485	! CALLED FROM Sub. bl_init in module_physics.init.F
486	!
487	! THIS PROGRAM PROCESSES THE USGS LANDUSE DATA
488	! WHICH HAS BEEN GRIDDED BY THE WPS SYSTEM
489	! AND PRODUCES 2-D FIELDS OF LU RELATED PARAMETERS
490	! FOR USE IN THE PX SURFACE MODEL
491	!
492	!
493	! REVISION HISTORY:
494	! AX Oct 2004 - developed WRF version based on VEGELAND in the MM5 PX LSM
495	! RG Dec 2006 - revised for WRFV2.1.2
496	! JP Dec 2007 - revised for WRFV3 - removed IFGROW options
497	!******************************************************************************
498
499
500	IMPLICIT NONE
501	!...
502	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
503	ims,ime, jms,jme, kms,kme, &
504	its,ite, jts,jte, kts,kte
505
506	INTEGER , INTENT(IN) :: NSCAT, NLCAT
507
508	REAL, DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN) :: LANDUSEF
509	REAL, DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN) :: SOILCTOP, SOILCBOT
510
511	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: VEGFRA, SHDMIN, SHDMAX
512
513	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ZNT
514
515	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: XLAI, XLAIMN, RSTMIN, XALB, &
516	XVEG, XVEGMN, XSNUP, XLAND
517
518	!... local variables
519
520	INTEGER :: itf, jtf, k, j, i
521	REAL :: SUMLAI, SUMLMN, SUMRSI, SUMLZ0, SUMVEG, SUMVMN, ALAI, VEGF, SUMSNUP
522	REAL :: VFMX, VFMN, VSEAS, FAREA, FWAT, ZNOTC, SUMALB
523	REAL, DIMENSION( NLCAT ) :: LAIMX, LAIMN, Z0, VEG, VEGMN, SNUP
524	REAL, PARAMETER :: ZNOTCMN = 5.0 ! CM, MIN Zo FOR CROPS
525	REAL, PARAMETER :: ZNOTCMX = 15.0 ! CM, MAX Zo FOR CROPS
526
527	!*****************************************************************************
528	! USGS LU characterization
529	!---------------------------
530	! Name Rstmin Zo Mxfr MnFr MxLA MnLA
531	! 1 Urban 150. 50. 40. 20. 2.0 0.5 Urban or Built-up Land
532	! 2 DrCrp 70. 10. 95. 15. 3.0 0.5 Dryland Cropland and Pasture
533	! 3 IrCrp 60. 10. 95. 10. 3.0 0.5 Irrigated Cropland and Pasture
534	! 4 MixCp 70. 10. 95. 15. 3.0 0.5 Mixed Dry/Irr Crop and Past
535	! 5 CrGrM 80. 10. 95. 35. 2.5 1.0 Grassland/Cropland Mosaic
536	! 6 CrWdM 180. 40. 95. 40. 4.0 1.5 Woodland/Cropland Mosaic
537	! 7 GrsLd 100. 7. 95. 70. 2.5 1.0 Grassland
538	! 8 ShrLd 200. 20. 70. 50. 3.0 1.0 Shrubland
539	! 9 ShrGr 150. 20. 85. 60. 3.0 1.0 Mixed Shrubland/Grassland
540	! 10 Savan 120. 20. 80. 60. 2.0 1.0 Savanna
541	! 11 DBFst 200. 50. 95. 50. 5.0 1.0 Broadleaf Deciduous Forest
542	! 12 DNFst 175. 50. 95. 50. 5.0 1.0 Deciduous Coniferous Forest
543	! 13 EBFst 120. 40. 95. 85. 5.0 4.0 Evergreen Broadleaf Forest (Palm?)
544	! 14 ENFst 175. 50. 90. 80. 4.0 3.0 Evergreen Coniferous Forest
545	! 15 MxFst 200. 50. 95. 60. 5.0 2.0 Mixed forest
546	! 16 Water 9999. 0.1 00. 00. 0.0 0.0 Water
547	! 17 HWtld 164. 15. 60. 40. 2.0 1.0 Herbaceous Wetland (none in east)
548	! 18 WWtld 200. 45. 90. 80. 5.0 3.0 Forested Wetlands (e.g. Everglades)
549	! 19 BarSp 100. 5. 10. 05. 0.5 0.2 Barren or Sparsely Vegetated
550	! 20 HrTun 150. 10. 20. 10. 1.0 0.5 Herbaceous Tundra
551	! 21 WdTun 200. 10. 30. 10. 1.0 0.5 Shrub and Brush Tundra
552	! 22 MxTun 150. 5. 20. 05. 1.0 0.5 Mixed Tundra
553	! 23 BGTun 100. 5. 5. 02. 0.1 0.1 Bare Ground Tundra
554	! 24 SnwIc 300. 5. 5. 02. 0.1 0.1 Perennial Snowfields or Glaciers
555	!-----------------------------------------------------------------------------
556
557	REAL, DIMENSION(24) :: RSMIN, Z00, VEG0, VEGMN0, &
558	LAI0, LAIMN0, SNUP0, ALBF
559
560	DATA RSMIN /150.0, 70.0, 60.0, 70.0, 80.0, 180.0, 100.0, 200.0, &
561	150.0, 120.0, 200.0, 175.0, 120.0, 175.0, 200.0,9999.0, &
562	164.0, 200.0, 100.0, 150.0, 200.0, 150.0, 100.0, 300.0/
563	DATA Z00 / 50.0, 10.0, 10.0, 10.0, 10.0, 40.0, 7.0, 20.0, &
564	20.0, 20.0, 50.0, 50.0, 40.0, 50.0, 50.0, 0.1, &
565	15.0, 45.0, 5.0, 10.0, 10.0, 5.0, 5.0, 5.0/
566	DATA VEG0 / 40.0, 95.0, 95.0, 95.0, 95.0, 95.0, 95.0, 70.0, &
567	85.0, 80.0, 95.0, 95.0, 95.0, 90.0, 95.0, 0.00, &
568	60.0, 90.0, 10.0, 20.0, 30.0, 20.0, 5.0, 5.0/
569	DATA VEGMN0/ 20.0, 15.0, 10.0, 15.0, 35.0, 40.0, 70.0, 50.0, &
570	60.0, 60.0, 50.0, 50.0, 85.0, 80.0, 60.0, 0.0, &
571	40.0, 80.0, 5.0, 10.0, 10.0, 5.0, 2.0, 2.0/
572	DATA LAI0 / 2.0, 3.0, 3.0, 3.0, 2.5, 4.0, 2.5, 3.0, &
573	3.0, 2.0, 5.0, 5.0, 5.0, 4.0, 5.0, 0.0, &
574	2.0, 5.0, 0.50, 1.0, 1.0, 1.0, 0.1, 0.1/
575	DATA LAIMN0/ 0.50, 0.50, 0.50, 0.50, 1.0, 1.5, 1.0, 1.0, &
576	1.0, 1.0, 1.0, 1.0, 4.0, 3.0, 2.0, 0.0, &
577	1.0, 3.0, 0.20, 0.50, 0.50, 0.50, 0.10, 0.10/
578	DATA SNUP0/ 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.03, &
579	0.035, 0.04, 0.08, 0.08, 0.08, 0.08, 0.08, 0.01, &
580	0.01, 0.01, 0.02, 0.02, 0.025, 0.025, 0.025, 0.02/
581	DATA ALBF/ 15.0, 17.0, 18.0, 18.0, 18.0, 16.0, 19.0, 22.0, &
582	20.0, 20.0, 16.0, 14.0, 12.0, 12.0, 13.0, 8.0, &
583	14.0, 14.0, 25.0, 15.0, 15.0, 15.0, 25.0, 55.0/
584
585	!---- INITIALIZE PARAMETERS
586
587	INTEGER :: KWAT
588	DATA KWAT / 16/
589
590	!--------------------------------------------------------------------
591	DO J = jts,jte
592	DO I = its,ite
593	XLAI(I,J) = 0.0
594	XLAIMN(I,J) = 0.0
595	RSTMIN(I,J) = 9999.0
596	XVEG(I,J) = 0.0
597	XVEGMN(I,J) = 0.0
598	XSNUP(I,J) = 0.0
599	XALB(I,J) = 0.0
600	ENDDO ! END LOOP THROUGH GRID CELLS
601	ENDDO ! END LOOP THROUGH GRID CELLS
602	!--------------------------------------------------------------------
603
604	DO J = jts,jte
605	DO I = its,ite
606
607
608	!-- Initialize 2 and 3-D veg parameters to be caculated
609	DO K=1,NLCAT
610	LAIMX(K) = LAI0(K)
611	LAIMN(K) = LAIMN0(K)
612	Z0(K) = Z00(K)
613	VEG(K) = VEG0(K)
614	VEGMN(K) = VEGMN0(K)
615	SNUP(K) = SNUP0(K)
616	ENDDO
617
618	!-- INITIALIZE SUMS
619	SUMLAI = 0.0
620	SUMLMN = 0.0
621	SUMRSI = 0.0
622	SUMLZ0 = 0.0
623	SUMVEG = 0.0
624	SUMVMN = 0.0
625	ALAI = 0.0
626	SUMSNUP= 0.0
627	SUMALB = 0.0
628
629	!-- ESTIMATE CROP EMERGANCE DATE FROM VEGFRAC
630	VFMX = SHDMAX(I,J)
631	VFMN = SHDMIN(I,J)
632	VEGF = VEGFRA(I,J)
633
634	!-- Computations for VEGETATION CELLS ONLY
635	IF(VFMX.GT.0.0.AND.LANDUSEF(I,KWAT,J).LT.1.00) THEN
636	VSEAS = VEGF/VFMX
637	IF(VSEAS.GT.1.0.OR.VSEAS.LT.0.0) THEN
638	VSEAS = MIN(VSEAS,1.0)
639	VSEAS = MAX(VSEAS,0.0)
640	ENDIF
641
642	ZNOTC = ZNOTCMN * (1-VSEAS) + ZNOTCMX * VSEAS ! Zo FOR CROPS
643	DO K = 1, NLCAT
644	!-- USE THE VEGFRAC DATA ONLY FOR CROPS
645	IF (K .GE. 2 .AND. K .LE. 6) THEN
646	LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS
647	LAIMN(K) = LAIMX(K)
648	VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS
649	VEGMN(K) = VEG(K)
650	!-- SEASONALLY VARY Zo FOR DryCrop (k=2) OR Irigated Crop (k=3) OR Mix Crop (k=4)
651	IF (K .GE. 2 .AND. K .LE. 4) THEN
652	Z0(K) = ZNOTC
653	!-- CrGrM (k=5) or CrWdM (k=6) USE AVG WITH GRASS AND FOREST
654	ELSE IF (K .GE.5 .AND. K .LE. 6) THEN
655	Z0(K) = 0.5 * (ZNOTC + Z00(K))
656	ENDIF
657	ENDIF
658	ENDDO
659	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
660
661	ENDIF !-- IF cell is vegetation
662	!-------------------------------------
663	!-- LOOP THROUGH LANDUSE Fraction and compute totals
664	DO K = 1, NLCAT
665	FAREA = LANDUSEF(I,K,J)
666	SUMLAI = SUMLAI + LAIMX(K) * FAREA
667	SUMLMN = SUMLMN + LAIMN(K) * FAREA
668	ALAI = ALAI + FAREA
669	SUMRSI = SUMRSI + FAREA * LAIMX(K) / RSMIN(K)
670	SUMLZ0 = SUMLZ0 + FAREA * ALOG(Z0(K))
671	SUMVEG = SUMVEG + FAREA * VEG(K)
672	SUMVMN = SUMVMN + FAREA * VEGMN(K)
673	SUMSNUP= SUMSNUP+ FAREA * SNUP(K)
674	SUMALB = SUMALB + FAREA * ALBF(K)
675	ENDDO
676
677	!-- CHECK FOR WATER
678	FWAT = LANDUSEF(I,KWAT,J)
679	IF (FWAT .GT. 0.999) THEN
680	XLAI(I,J) = LAIMX(KWAT)
681	XLAIMN(I,J) = LAIMN(KWAT)
682	RSTMIN(I,J) = RSMIN(KWAT)
683	ZNT(I,J) = Z0(KWAT)
684	XVEG(I,J) = VEG(KWAT)
685	XVEGMN(I,J) = VEGMN(KWAT)
686	XSNUP(I,J) = SNUP(KWAT)
687	XALB(I,J) = ALBF(KWAT)
688	ELSE
689	IF (FWAT .GT. 0.10) THEN
690	ALAI = ALAI - FWAT
691	SUMLZ0 = SUMLZ0 - FWAT * ALOG(Z0(KWAT))
692	ENDIF
693	XLAI(I,J) = SUMLAI / ALAI
694	XLAIMN(I,J) = SUMLMN / ALAI
695	RSTMIN(I,J) = SUMLAI / SUMRSI
696	ZNT(I,J) = EXP(SUMLZ0/ALAI)
697	XVEG(I,J) = SUMVEG / ALAI
698	XVEGMN(I,J) = SUMVMN / ALAI
699	XSNUP(I,J) = SUMSNUP/ALAI
700	XALB(I,J) = SUMALB/ALAI
701	ENDIF
702
703	IF (FWAT .GT. 0.50) THEN
704	ZNT(I,J) = Z0(KWAT)
705	XALB(I,J) = ALBF(KWAT)
706	ENDIF
707
708	ZNT(I,J) = ZNT(I,J) * 0.01 !CONVERT TO M
709	XVEG(I,J) = XVEG(I,J) * 0.01 !CONVERT TO FRAC
710	XVEGMN(I,J) = XVEGMN(I,J) * 0.01
711	XLAND(I,J) = 1.0 + FWAT
712	XALB(I,J) = XALB(I,J) * 0.01
713
714	ENDDO ! END LOOP THROUGH GRID CELLS
715	ENDDO ! END LOOP THROUGH GRID CELLS
716	!--------------------------------------------------------------------
717
718	!------------------------------------------------------------------------------
719	END SUBROUTINE vegeland
720	!------------------------------------------------------------------------------
721	!**********************************************************************
722	SUBROUTINE SURFPX(DTPBL, IFLAND, ISNOW, NUDGEX, XICE1, & !in
723	SOLDN, GSW, LWDN, EMISSI, Z1, & !in
724	MOL, ZNT, UST, & !in
725	PSURF, DENS1, QV1, QSS, TA1, & !in
726	THETA1, PRECIP, & !in
727	CPAIR, PSIH, & !in
728	RH2OBS, T2OBS, & !in
729	VEGFRC, ISTI, LAI, BETAP, & !in
730	RSTMIN, HC_SNOW, SNOW_FRA, & !in
731	WWLT, WFC, CGSAT, WSAT, B, & !in
732	C1SAT, C2R, AS, JP, DS1, DS2,QST12, & !in
733	RADNET, GRDFLX, HFX, QFX, LH, EG, ER, ETR, & !out
734	QST, CAPG, RS, RA, TG, T2, & !out
735	WG, W2, WR, TA2, QA2 ) !out
736
737	IMPLICIT NONE
738	!**********************************************************************
739	!
740	! FUNCTION:
741	! THIS SUBROUTINE COMPUTES SOIL MOISTURE AND TEMPERATURE TENDENCIES
742	! BY SOLVING THE PROGNOSTIC EQUATIONS IN PX95.
743	!
744	! SUBROUTINES CALLED:
745	! SUB. QFLUX compute the soil and canopy evaporation, and transpiration
746	! SUB. SMASS compute nudging coefficients for soil moisture and temp nudging
747	!
748	! ARGUMENTS:
749	! DTPBL: TIME STEP OF THE MINOR LOOP FOR THE LAND-SURFACE/PBL MODEL
750	! IFLAND: INDEX WHICH INDICATES THE TYPE OF SURFACE,=1,LAND;=2,SEA
751	! ISNOW: SNOW (=1) OR NOT (=0)
752	! NUDGE: SWITCH FOR SOIL MOISTURE NUDGING
753	! SOLDN: SHORT-WAVE RADIATION
754	! LWDN: LONG-WAVE RADIATION
755	! EMISSI: EMISSIVITY
756	! Z1: HEIGHT OF THE LOWEST HALF LAYER
757	! MOL: MONIN-OBUKOV LENGH (M)
758	! ZNT: ROUGHNESS LENGTH (M)
759	! UST: FRICTION VELOCITY (M/S)
760	! TST: Turbulent moisture scale
761	! RA: AERODYNAMIC RESISTENCE
762	! PSURF: P AT SURFACE (CB)
763	! DENS1: AIR DENSITY AT THE FIRST HALF LAYER
764	! QV1: Air humidity at first half layer
765	! QSS: Saturation mixing ratio at ground
766	! TA1: Air temperature at first half layer
767	! THETA1: Potential temperature at first half layer
768	! PRECIP: Precipitation rate in m/s
769	! STBOLT: STEFAN BOLTZMANN'S CONSTANT
770	! KARMAN: VON KARMAN CONSTANT
771	! CPAIR: Specific heat of moist air (M^2 S^-2 K^-1)
772	! TAUINV: 1/1DAY(SEC)
773	! VEGFRC: Vegetation coverage
774	! ISTI: soil type
775	! LAI: Leaf area index
776	! BETAP: Coefficient for bare soil evaporation
777	! THZ1OB: Observed TEMP FROM ANAL OF OBS TEMPS AT SCREEN HT
778	! RHOBS: Observed relative humidity at SCREEN HT
779	! RSTMIN Minimum Stomatol resistence
780	!... Outputs from SURFPX
781	! RADNET: Net radiation
782	! HFX: SENSIBLE HEAT FLUX (W / M^2)
783	! QFX: TOTAL EVAP FLUX (KG / M^2 S)
784	! GRDFLX: Ground heat flux (W / M^2)
785	! QST: Turbulent moisture scale
786	! CAPG: THERMAL CAPACITY OF GROUND SLAB (J/M^2/K)
787	! RS: Surface resistence
788	! RA: Surface aerodynamic resistence
789	! EG: evaporation from ground (bare soil)
790	! ER: evaporation from canopy
791	! ETR: transpiration from vegetation
792	! TA2: diagnostic 2-m temperature from surface layer and lsm
793	!
794	!... Updated variables in this subroutine
795	! TG: Soil temperature at first soil layer
796	! T2: Soil temperature in root zone
797	! WG: Soil moisture at first soil layer
798	! W2: Soil moisture at root zone
799	! WR: Canopy water content in m
800	!
801	! REVISION HISTORY:
802	! AX 2/2005 - developed WRF version based on the MM5 PX LSM
803	!
804	!**********************************************************************
805	!
806	!----------------------------------------------------------------------
807	!
808	!.......Arguments
809
810	!.. Integer
811	INTEGER , INTENT(IN) :: ISTI, NUDGEX
812
813	!... Real
814	REAL , INTENT(IN) :: DTPBL, DS1, DS2
815	REAL , INTENT(IN) :: IFLAND, ISNOW, XICE1
816	REAL , INTENT(IN) :: SOLDN, GSW, LWDN, EMISSI, Z1
817	REAL , INTENT(IN) :: ZNT
818	REAL , INTENT(IN) :: PSURF, DENS1, QV1, QSS, TA1, THETA1, PRECIP
819	REAL , INTENT(IN) :: CPAIR
820	REAL , INTENT(IN) :: VEGFRC, LAI
821	REAL , INTENT(IN) :: RSTMIN, HC_SNOW, SNOW_FRA
822	REAL , INTENT(IN) :: WWLT, WFC, CGSAT, WSAT, B, C1SAT, C2R, AS, JP
823	REAL , INTENT(IN) :: RH2OBS,T2OBS
824	REAL , INTENT(IN) :: QST12
825
826	REAL , INTENT(OUT) :: RADNET, EG, ER, ETR
827	REAL , INTENT(OUT) :: QST, CAPG, RS, TA2, QA2
828
829	REAL , INTENT(INOUT) :: TG, T2, WG, W2, WR, UST, RA, BETAP
830	REAL , INTENT(INOUT) :: GRDFLX, QFX, HFX, LH, PSIH, MOL
831
832	!... Local Variables
833
834	!... Real
835	REAL :: HF, LV, CQ4
836	REAL :: RAH, RAW, ET, W2CG, CG, CT, SOILFLX, CPOT, THETAG
837	REAL :: ZOL, ZOBOL, ZNTOL, Y, Y0, PSIH15, YNT
838	REAL :: WGNUDG, W2NUDG, ALPH1, ALPH2, BET1, BET2, T1P5
839	REAL :: CQ1, CQ2, CQ3, COEFFNP1, COEFFN, TSNEW, TSHLF, T2NEW
840	REAL :: ROFF, WRMAX, PC, DWR, PNET, TENDWR, WRNEW
841	REAL :: COF1, CFNP1WR, CFNWR, PG, FASS
842	REAL :: TENDW2, W2NEW, W2HLF, W2REL, C1, C2, WEQ, CFNP1, CFN, WGNEW
843	REAL :: ALN10, TMP1, TMP2, TMP3, AA, AB, TST, RBH, CTVEG
844	REAL :: QST1,PHIH,PSIOB
845	REAL :: T2NUD, T2NUDF
846	REAL :: VAPPRS, QSBT, RH2MOD
847	REAL :: A,BB,C,D
848
849	!... Parameters
850	REAL :: ZOBS, GAMAH, BETAH, SIGF, BH, CT_SNOW
851	REAL, PARAMETER :: CV = 8.0E-6 ! K-M2/J
852	PARAMETER (ZOBS = 1.5) ! height for observed screen temp., (m)
853	PARAMETER (BH = 15.7)
854	PARAMETER (GAMAH = 16. ) !11.6)
855	PARAMETER (BETAH = 5.0 ) !8.21)
856	PARAMETER (SIGF = 0.5) ! rain interception see LSM (can be 0-1)
857	!PARAMETER (CT_SNOW = 5.54E-5)
858	! New value of CT_SNOW calibrated using multilayer soil model where csnow=6.9E5 J/(m3 K)
859	! from NCAR CSM
860	PARAMETER (CT_SNOW = 2.0E-5)
861	ALN10 = ALOG(10.0)
862	RADNET = SOLDN - (EMISSI (STBOLT TG **4 - LWDN)) ! NET RADIATION
863	!--------------------------------------------------------------------
864	CPOT= (100.0 / PSURF) ** ROVCP ! rcp is global constant(module_model_constants)
865	THETAG = TG * CPOT
866
867	ZOL = Z1/MOL
868	ZOBOL = ZOBS/MOL
869	ZNTOL = ZNT/MOL
870
871	!-----------------------------------------------------------------------------------------
872	IF (MOL .LT. 0.0) THEN
873	Y = ( 1.0 - GAMAH * ZOL ) ** 0.5
874	Y0 = ( 1.0 - GAMAH * ZOBOL ) ** 0.5
875	YNT = ( 1.0 - GAMAH * ZNTOL ) ** 0.5
876	PSIH15 = 2.0 * ALOG((Y + 1.0) / (Y0 + 1.0))
877	PSIH = 2.0 * ALOG((Y + 1.0) / (YNT + 1.0))
878	PSIOB = 2.0 * ALOG((Y0 + 1.0) / (YNT + 1.0))
879	PHIH = 1.0 / Y
880	ELSE
881	IF((ZOL-ZNTOL).LE.1.0) THEN
882	PSIH = -BETAH*(ZOL-ZNTOL)
883	ELSE
884	PSIH = 1.-BETAH-(ZOL-ZNTOL)
885	ENDIF
886	IF ((ZOBOL - ZNTOL) .LE. 1.0) THEN
887	PSIOB = -BETAH * (ZOBOL - ZNTOL)
888	ELSE
889	PSIOB = 1.0 - BETAH - (ZOBOL - ZNTOL)
890	ENDIF
891	PSIH15 = PSIH - PSIOB
892	IF(ZOL.LE.1.0) THEN
893	PHIH = 1.0 + BETAH * ZOL
894	ELSE
895	PHIH = BETAH + ZOL
896	ENDIF
897	ENDIF
898	!-----------------------------------------------------------------------------------------
899	!-- ADD RA AND RB FOR HEAT AND MOISTURE
900	!... RB FOR HEAT = 5 /UST
901	!... RB FOR WATER VAPOR = 5*(0.599/0.709)^2/3 /UST = 4.503/UST
902	RA=PR0* ( ALOG(Z1/ZNT) - PSIH )/(KARMAN*UST)
903	RAH = RA + 5.0 / UST
904	RAW = RA + 4.503 / UST
905	!--------------------------------------------------------------------
906	!-- COMPUTE MOISTURE FLUX
907	CALL QFLUX( DENS1, QV1, TA1, SOLDN, RAW, QSS, &
908	VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, &
909	WG, W2, WR, &
910	RSTMIN, WWLT, WFC, &
911	EG, ER, ETR, CQ4, RS, FASS)
912	!--------------------------------------------------------------------
913
914	!--------------------------------------------------------------------
915	!..........Total evaporation (ET)
916	ET = EG + ER + ETR
917	QST = -ET / (DENS1 * UST)
918
919	LV = 2.83E6 !-- LATENT HEAT OF SUBLIMATION AT 0C FROM STULL(1988)
920	IF (ISNOW .LT. 0.5.AND.TG.GT.273.15) &
921	LV = (2.501 - 0.00237 * (TG - 273.15)) * 1.E6 !-- FROM STULL(1988) in J/KG
922
923	IF (IFLAND .LT. 1.5 ) QFX = ET !-- Recaculate QFX over land to account for P-X LSM EG, ER and ETR
924	LH = LV * QFX
925	!-----------------------------------------------------------------------------------------
926
927	!-----------------------------------------------------------------------------------------
928	! Surface sensible heat flux
929	TST = (THETA1 - THETAG ) / (UST*RAH)
930	HF = UST * TST
931	HFX = AMAX1(-DENS1 * CPAIR * HF, -250.0) ! using -250. from MM5
932	!-----------------------------------------------------------------------------------------
933
934	!-----------------------------------------------------------------------------------------
935	! Compute the diagnosed 2m Q and T consistent with PX LSM
936	QST1 = 0.5*(QST+QST12/PHIH)
937	TA2= (THETAG + TST * (PR0 / KARMAN * (ALOG(ZOBS / ZNT) - PSIOB)+5.))/CPOT
938	QA2= QV1 - QST1 * PR0/ KARMAN * (ALOG(Z1 / ZOBS) - PSIH15)
939	IF (QA2.LE.0.0) QA2 = QV1
940	!-- Relative humidity
941	VAPPRS = SVP1 * EXP(SVP2 * (TA2 - SVPT0) / (TA2 - SVP3))
942	QSBT = EP_2 * VAPPRS / (PSURF - VAPPRS)
943	RH2MOD = QA2 / QSBT
944	!-----------------------------------------------------------------------------------------
945	IF (IFLAND .LT. 1.5 ) THEN
946	W2CG = AMAX1(W2,WWLT)
947	CG = CGSAT * 1.0E-6 * (WSAT/ W2CG) ** &
948	(0.5 * B / ALN10)
949	CT = 1./((1-VEGFRC)/CG + VEGFRC/CV)
950	CT = 1./(SNOW_FRA/CT_SNOW + (1-SNOW_FRA)/CT)
951	CAPG = 1.0/CT
952
953	SOILFLX = 2.0 * PI * TAUINV * (TG - T2)
954	GRDFLX = SOILFLX / CT
955	ENDIF
956	!-----------------------------------------------------------------------------------------
957
958	!--------------------------------------------------------------------
959	!-- ASSIMILATION --- COMPUTE SOIL MOISTURE NUDGING FROM TA2 and RH2
960	!-------COMPUTE ASSIMILATION COEFFICIENTS FOR ALL I
961	IF (IFLAND .LT. 1.5) THEN
962	IF (NUDGEX .EQ. 0) THEN !-- NO NUDGING CASE
963	WGNUDG = 0.0
964	W2NUDG = 0.0
965	T2NUD = 0.0
966	ELSE !-- NUDGING CASE
967
968	CALL SMASS (ISTI, FASS, SOLDN, VEGFRC, RA, WWLT, WFC, &
969	ALPH1, ALPH2, BET1, BET2, T2NUDF)
970
971	!--COMPUTE MODEL RH
972	WGNUDG = ALPH1 * (T2OBS - TA2) + ALPH2 * (RH2OBS - RH2MOD) * 100
973	W2NUDG = BET1 * (T2OBS - TA2) + BET2 * (RH2OBS - RH2MOD) * 100
974	IF (W2 .GE. WFC) W2NUDG = AMIN1(W2NUDG,0.0)
975	IF (W2 .LE. WWLT) W2NUDG = AMAX1(W2NUDG,0.0)
976	T2NUD = T2NUDF * (T2OBS - TA2)
977	!print *, 'T2NUD =',T2NUD,T2NUDF
978	ENDIF
979	ENDIF
980	!-----------------------------------------------------------------------------------------
981
982	!-----------------------------------------------------------------------------------------
983	!-- Compute new values for TS,T2,WG,W2 and WR. No change over ice or water (XLAND > 1)
984	IF (IFLAND .LT. 1.5) THEN
985
986	!-- SOLVE BY CRANK-NIC -- TENDTS=CT*(RADNET-HFX-QFX)-SOILFLX
987	!-- Calculate the coefficients for implicit calculation of TG
988	CQ1 = (1.0 - 0.622 * LV * CRANKP / (r_d * TG)) * QSS
989	CQ2 = 0.622 * LV * QSS * CRANKP / (r_d * TG * TG)
990	CQ3 = DENS1 * BETAP * (1.0 - VEGFRC) / RAW
991	COEFFNP1 = 1.0 + DTPBL * CRANKP * (4.0 * EMISSI * STBOLT * TG ** 3 &
992	* CT + DENS1 * CPAIR / RAH * CPOT * CT + 2.0 * PI &
993	* TAUINV ) + DTPBL * (CT * LV * CQ2 * (CQ3 + CQ4))
994	COEFFN = CT * (GSW + EMISSI * (STBOLT * (4.0 * CRANKP - 1.0) &
995	* TGTGTG*TG + LWDN) & !NET RAD
996	+ DENS1 * CPAIR / RAH * (THETA1 - (1.0 - CRANKP) * THETAG) &
997	- LV * (CQ3 * (CQ1 - QV1) + CQ4 * (CQ1 - QV1))) & !SFC HEAT FLUX
998	- 2.0 * PI * TAUINV * ((1.0 - CRANKP) * TG - T2) !SOIL FLUX
999	TSNEW = (TG + DTPBL * COEFFN) / COEFFNP1
1000	!-- FOR SNOW COVERED SURFACE TEMPERATURE IS NOT MORE THAN ZERO
1001	IF (XICE1 .GT. 0.5) TSNEW = AMIN1(TSNEW,273.15) ! Re-added Jan 2010 to keep ice surface at or below freezing (J. Pleim)
1002	TSHLF = 0.5 * ( TSNEW + TG)
1003	T2NEW = (T2 + DTPBL * TAUINV * T2TFAC * (TSHLF - (1 - CRANKP) * T2) &
1004	+ DTPBL*T2NUD) & ! Added deep temperature nudging
1005	/ (1.0 + DTPBL * TAUINV * T2TFAC * CRANKP)
1006	!-- REPLACE OLD with NEW Value
1007	TG = TSNEW
1008	T2 = T2NEW
1009	ENDIF
1010	!-----------------------------------------------------------------------------------------
1011
1012	!-----------------------------------------------------------------------------------------
1013	! Compute new subsurface soil and canopy moisture values DENS1. No change required over ocean.
1014	IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN
1015	!-- Compute WR
1016	ROFF = 0.0
1017	WRMAX = 0.2E-3 * VEGFRC * LAI ! max. WRMAX IN m
1018
1019	IF(WRMAX.GT.0.0) THEN
1020	!-- PC is precip. intercepted by veg.(M/S)
1021	PC = VEGFRC * SIGF * PRECIP
1022	DWR = (WRMAX - WR) / DTPBL !the tendency to reach max.
1023	PNET = PC - ER/ DENW ! residual of precip. and evap.
1024	IF (PNET .GT. DWR) THEN
1025	ROFF = PNET - DWR
1026	PC = PC - ROFF
1027	ENDIF
1028	IF (QSS .LT. QV1) THEN
1029	TENDWR = PC - ER / DENW
1030	WRNEW = WR + DTPBL * TENDWR
1031	ELSE
1032	COF1 = DENS1 / DENW * VEGFRC * (QSS - QV1) / RAW
1033	!-- using delta=wr/wrmax
1034	CFNP1WR = 1.0 + DTPBL * COF1 * CRANKP / WRMAX
1035	CFNWR = PC - COF1 * (1.0 - CRANKP) * WR / WRMAX
1036	WRNEW = (WR + DTPBL * CFNWR) / CFNP1WR
1037	ENDIF
1038	ELSE
1039	WRNEW=0.0
1040	ENDIF
1041	!---------------------------------------------
1042	!-- Compute W2
1043	PG = DENW * (PRECIP - PC) ! PG is precip. reaching soil (PC already including ROFF)
1044	TENDW2 = 1.0 / (DENW * DS2) * (PG - EG - ETR) &
1045	+ (W2NUDG + WGNUDG) / DS2 ! NUDGING
1046	W2NEW = W2 + DTPBL * TENDW2
1047	W2NEW = AMIN1(W2NEW,WSAT)
1048	W2NEW = AMAX1(W2NEW,0.05)
1049	W2HLF = 0.5 * (W2 + W2NEW)
1050	!.. new values
1051	W2 = W2NEW
1052	WR = AMIN1(WRMAX,WRNEW)
1053	ENDIF
1054	!-----------------------------------------------------------------------------------------
1055
1056	!-----------------------------------------------------------------------------------------
1057	! Compute new surface soil moisture values (WR).
1058	IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN ! over ocean no change to wg w2,wr
1059	!-- FOR SNOW COVERED SURFACE, ASSUME SURFACE IS SATURATED AND
1060	! WG AND W2 ARE NOT CHANGED
1061	IF (ISNOW .GT.0.5) THEN
1062	WG = WSAT
1063	ELSE
1064	W2REL = W2HLF / WSAT
1065	IF (WG .GT. WWLT) THEN
1066	C1 = C1SAT * (WSAT / WG) ** (0.5 * B + 1.0)
1067	ELSE ! elimilate C1 for wg < wilting point
1068	C1 = C1SAT * (WSAT / WWLT) ** (0.5 * B + 1.0)
1069	ENDIF
1070	C2 = C2R * W2HLF / (WSAT - W2HLF + 1.E-11)
1071	IF (W2HLF .EQ. WSAT) THEN
1072	WEQ = WSAT
1073	ELSE
1074	WEQ = W2HLF - AS * WSAT * W2REL ** JP * &
1075	(1.0 - W2REL ** (8.0 * JP))
1076	ENDIF
1077
1078	!.... The beta method, Lee & Pielke (JAM, May 1992)
1079	CFNP1 = 1.0 + DTPBL * C2 * TAUINV * CRANKP
1080	CFN = C1 / (DENW * DS1) * (PG - EG) - C2 * TAUINV * &
1081	((1.0 - CRANKP) * WG - WEQ) + WGNUDG/ DS1
1082
1083	WGNEW = AMAX1((WG + DTPBL * CFN) / CFNP1,0.001)
1084	!-- NEW VALUES
1085	WG = AMIN1(WGNEW,WSAT)
1086
1087	ENDIF !endif for ISNOW
1088	ENDIF !endif for XLAND
1089	!-----------------------------------------------------------------------------------------
1090
1091	END SUBROUTINE surfpx
1092	!***********************************************************************
1093	!***********************************************************************
1094
1095	!***********************************************************************
1096	!***********************************************************************
1097	SUBROUTINE QFLUX (DENS1, QV1, TA1, RG, RAW, QSS, & ! in
1098	VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, & ! in
1099	WG, W2, WR, & ! in
1100	RSTMIN, WWLT, WFC, &
1101	EG, ER, ETR, CQ4, RS, FASS) ! out
1102
1103	!***********************************************************************
1104	!
1105	! FUNCTION:
1106	! THIS SUBROUTINE COMPUTES EVAPORATION FROM BARE SOIL (EG) AND FROM
1107	! THE WET PART OF CANOPY (ER) AND TRANSPIRATION FROM THE DRY PART OF
1108	! CANOPY (ETR).
1109	!
1110	! REVISION HISTORY:
1111	! A. Xiu Oct 2004 - adapted from the PX LSM in MM5 for the WRF system
1112	! R. Gilliam Dec 2006 - Completed WRF V2.1.2 implementation
1113	!
1114	!***********************************************************************
1115	! QFLUX ARGUMENT LIST:
1116	! ----------------------------------------------------------------------
1117	! INPUTS:
1118	!-- DENS1 air density at first layer
1119	!-- QV1 air humidity at first layer
1120	!-- TA1 air temperature at first layer
1121	!-- RG shortwave radition reaching the ground
1122	!-- RAW RA+RB for moisture
1123	!-- QSS saturation mixing ratio at ground
1124	!-- VEGFRC vegetation coverage
1125	!-- ISNOW if snow on the ground
1126	!-- ISTI soil type
1127	!-- IFLAND if land (=1) or water (=2)
1128	!-- LAI leaf area index
1129	!-- BETAP
1130	!-- WG soil moisture at first soil layer
1131	!-- W2 soil moisture at root zone
1132	!-- WR Canopy water
1133	!
1134	! OUTPUTS FROM QFLUX:
1135	!-- EG evaporation from ground (bare soil)
1136	!-- ER evaporation from canopy
1137	!-- ETR transpiration from vegetation
1138	!-- CQ4
1139	!-- RS surface resistence
1140	!-- FASS parameter for soil moisture nudging
1141	! ----------------------------------------------------------------------
1142
1143	IMPLICIT NONE
1144
1145	! ----------------------------------------------------------------------
1146	! DECLARATIONS - INTEGER
1147	INTEGER , INTENT(IN) :: ISTI
1148
1149	! ----------------------------------------------------------------------
1150	! DECLARATIONS - REAL
1151	REAL , INTENT(IN) :: ISNOW, IFLAND
1152	REAL , INTENT(IN) :: DENS1, QV1, TA1, RG, RAW, QSS, &
1153	VEGFRC, LAI, &
1154	WG, W2, WR, RSTMIN
1155	REAL , INTENT(INOUT) :: BETAP
1156	REAL, INTENT(IN) :: WWLT, WFC
1157
1158	REAL , INTENT(OUT) :: EG, ER, ETR, CQ4, RS, FASS
1159
1160	!... Local Variables
1161
1162	!... Real
1163	REAL :: WRMAX, DELTA, SIGG, RADL, RADF, W2AVAIL, W2MXAV
1164	REAL :: FTOT, F1, F2, F3, F4
1165	REAL :: FSHELT, GS, GA, FX
1166
1167
1168
1169	!... Parameters
1170	REAL, PARAMETER :: RSMAX = 5000.0 ! s/m
1171	REAL, PARAMETER :: FTMIN = 0.0000001 ! m/s
1172	REAL, PARAMETER :: F3MIN = 0.25
1173
1174	!
1175	!... for water surface, no canopy evaporation and transpiration
1176	ER = 0.0
1177	ETR = 0.0
1178	CQ4 = 0.0
1179
1180	!... GROUND EVAPORATION (DEPOSITION)
1181	IF (QSS .LT. QV1) BETAP = 1.0
1182	EG = DENS1 * (1.0 - VEGFRC) * BETAP * (QSS - QV1) / RAW
1183
1184	!!---------------------------------------------------------------------
1185	!... CANOPY
1186	IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1187	WRMAX = 0.2E-3 * VEGFRC * LAI ! in unit m
1188	IF (WR .LE. 0.0) THEN
1189	DELTA = 0.0
1190	ELSE
1191	! DELTA = (WR / WRMAX) ** 0.66667
1192	DELTA = WR / WRMAX ! referred to SiB model
1193	ENDIF
1194	IF (QSS .GE. QV1) THEN
1195	SIGG = DELTA
1196	ELSE
1197	SIGG = 1.0
1198	ENDIF
1199
1200	ER = DENS1 * VEGFRC * SIGG * (QSS - QV1) / RAW
1201	! IF(ER.GT.0.5E-3)PRINT *,' WR,SIGG,DENS1,VEGFRC,QSS,QV1,RAW=', &
1202	! WR,SIGG,DENS1,VEGFRC,QSS,QV1,RAW,' ER=',ER
1203	ENDIF
1204	!!---------------------------------------------------------------------
1205
1206	!-- TRANSPIRATION
1207	!!---------------------------------------------------------------------
1208	IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1209
1210	!!!-RADIATION
1211	IF (RSTMIN .GT. 130.0) THEN
1212	RADL = 30.0 ! W/M2
1213	ELSE
1214	RADL = 100.0 ! W/M2
1215	ENDIF
1216	RADF = 1.1 * RG / (RADL * LAI) ! NP89 - EQN34
1217	F1 = (RSTMIN / RSMAX + RADF) / (1.0 + RADF)
1218
1219	!!!-SOIL MOISTURE
1220	W2AVAIL = W2 - WWLT
1221	W2MXAV = WFC - WWLT
1222	F2 = 1.0 / (1.0 + EXP(-5.0 * ( W2AVAIL / W2MXAV - &
1223	(W2MXAV / 3.0 + WWLT)))) ! according JP, 9/94
1224
1225	!-AIR TEMP
1226	!... according to Avissar (1985) and AX 7/95
1227	IF (TA1 .LE. 302.15) THEN
1228	F4 = 1.0 / (1.0 + EXP(-0.41 * (TA1 - 282.05)))
1229	ELSE
1230	F4 = 1.0 / (1.0 + EXP(0.5 * (TA1 - 314.0)))
1231	ENDIF
1232
1233	FTOT = LAI * F1 * F2 * F4
1234	ENDIF
1235
1236	!!---------------------------------------------------------------------
1237	IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN
1238	FSHELT = 1.0 ! go back to NP89
1239	GS = FTOT / (RSTMIN * FSHELT)
1240	GA = 1.0 / RAW
1241	!-- Compute humidity effect according to RH at leaf surf
1242	F3 = 0.5 * (GS - GA + SQRT(GA * GA + GA * GS * &
1243	(4.0 * QV1 / QSS - 2.0) + GS * GS)) / GS
1244	F3 = AMIN1(AMAX1(F3,F3MIN),1.0)
1245
1246	RS = 1.0 / (GS * F3)
1247
1248	!IF(RS .GT. 999999.0) THEN
1249	! print *,'RS,LAI,VEGFRC,ISTI-->',RS,LAI,VEGFRC,ISTI,ISNOW
1250	!ENDIF
1251
1252	!--- Compute Assimilation factor for soil moisture nudging - jp 12/94
1253	!-- Note that the 30 coef is to result in order 1 factor for max
1254	IF (RG .LT. 0.00001) THEN ! do not nudge during night
1255	FX = 0.0
1256	ELSE
1257	FX = 30.0 * F1 * F4 * LAI / (RSTMIN * FSHELT)
1258	ENDIF
1259
1260	FASS = FX
1261	ETR = DENS1 * VEGFRC * (1.0 - SIGG) * (QSS - QV1) / (RAW + RS)
1262	!..... CQ4 is used for the implicit calculation of TG in SURFACE
1263	CQ4 = DENS1 * VEGFRC * ((1.0 - SIGG) / (RAW + RS) + SIGG / RAW)
1264	ENDIF
1265	!!---------------------------------------------------------------------
1266
1267	END SUBROUTINE qflux
1268	!------------------------------------------------------------------------------------------
1269	!------------------------------------------------------------------------------------------
1270
1271	!------------------------------------------------------------------------------------------
1272	SUBROUTINE SMASS (ISTI, FASS, RG, VEGFRC, RA, & !in
1273	WWLT, WFC, &
1274	ALPH1, ALPH2, BET1, BET2, T2NUDF ) !out
1275	!------------------------------------------------------------------------------------------
1276	IMPLICIT NONE
1277	!*******************************************************************
1278	! SMASS COMPUTES SOIL MOISTURE NUDGING COEFFICIENTS
1279	!*******************************************************************
1280	!
1281	!.........Arguments
1282	INTEGER, PARAMETER :: NSCAT = 16 ! max. soil types
1283
1284	INTEGER, INTENT(IN) :: ISTI
1285	REAL, INTENT(IN) :: FASS, RG, VEGFRC, RA
1286	REAL, INTENT(IN) :: WWLT, WFC
1287	REAL, INTENT(OUT) :: ALPH1, ALPH2, BET1, BET2, T2NUDF
1288
1289
1290	!........Local variables
1291
1292	!... Real
1293	REAL :: FBET, FALPH, FRA, FTEXT
1294	REAL, DIMENSION( 1: NSCAT ) :: WFCX, WWLTX
1295
1296	!... Parameters
1297	REAL, PARAMETER :: A1MAX = -10.E-5, A2MAX = 1.E-5 ! m/K, m for 6hr period
1298	REAL, PARAMETER :: B1MAX = -10.E-3, B2MAX = 1.E-3 ! m/K, m (Bouttier et al 1993)
1299	REAL, PARAMETER :: TASSI = 4.6296E-5 ! 1/6hr in 1/sec
1300	REAL, PARAMETER :: RAMIN = 10.0 ! 0.1 s/cm
1301	!
1302	!-- WFC is field capacity (M^3/M^3) (JN90)
1303	DATA WFCX / 0.135, 0.150, 0.195, 0.255, 0.240, 0.255, 0.322, &
1304	0.325, 0.310, 0.370, 0.367, 0.367, 0.367, 0.367, 0.367, 0.367 /
1305	!
1306	!-- WWLT is wilting point (M^3/M^3) (JN90)
1307	DATA WWLTX / 0.068, 0.075, 0.114, 0.179, 0.155, 0.175, 0.218, &
1308	0.250, 0.219, 0.283, 0.286, 0.286, 0.286, 0.286, 0.286, 0.286 /
1309
1310	!
1311	FBET = FASS
1312	FALPH = RG / 1370.0
1313	!--TEXTURE FACTOR NORMALIZED BY LOAM (IST=5)
1314	FRA = RAMIN / RA ! scale by aerodynamic resistance
1315	FTEXT = TASSI * (WWLT + WFC) / (WWLTX(5) + WFCX(5)) * FRA
1316	! write(6,*) ' ftot, fbet=',ftot, fbet,' ftext=',ftext/tassi
1317	!
1318	ALPH1 = A1MAX * FALPH * (1.0 - VEGFRC) * FTEXT
1319	ALPH2 = A2MAX * FALPH * (1.0 - VEGFRC) * FTEXT
1320	BET1 = B1MAX * FBET * VEGFRC * FTEXT
1321	BET2 = B2MAX * FBET * VEGFRC * FTEXT
1322	T2NUDF = 1.0E-5 * MAX((1.0 - 5.0 * FALPH),0.0) ! T2 Nudging at night
1323
1324	END SUBROUTINE smass
1325	!------------------------------------------------------------------------------------------
1326	!------------------------------------------------------------------------------------------
1327	!------------------------------------------------------------------------------------------
1328	SUBROUTINE SOILPROP (SOILCBOT,WEIGHT, ITIMESTEP, MAVAIL, & ! IN
1329	PXLSM_SMOIS_INIT, & ! IN
1330	FWSAT,FWFC,FWWLT,FB,FCGSAT, & ! OUT
1331	FJP,FAS,FC2R,FC1SAT,ISTI, W2 ) ! OUT
1332	IMPLICIT NONE
1333	!*******************************************************************
1334	! SOILPROP COMPUTES SOIL PARAMETERS FOR BOTH BOTTOM AND TOP LAYERS
1335	! USING FRACTIONAL SOIL TYPE. A HARD CODED OPTION IS AVAILIABLE
1336	! TO COMPUTE THE SOIL PARAMETERS USING FRACTIONAL INFORMATION, OR
1337	! TO JUST USE SOIL PARAMETERS OF THE DOMINANT SOIL TYPE
1338	!*******************************************************************
1339	!------------------------------------------------------------------------
1340	!-- SOIL PARAMETERS ARE SPECIFIED BY SOIL TYPE:
1341	! # SOIL TYPE WSAT WFC WWLT B CGSAT JP AS C2R C1SAT
1342	! _ _________ ____ ___ ____ ____ _____ ___ ___ ___ _____
1343	! 1 SAND .395 .135 .068 4.05 3.222 4 .387 3.9 .082
1344	! 2 LOAMY SAND .410 .150 .075 4.38 3.057 4 .404 3.7 .098
1345	! 3 SANDY LOAM .435 .195 .114 4.90 3.560 4 .219 1.8 .132
1346	! 4 SILT LOAM .485 .255 .179 5.30 4.418 6 .105 0.8 .153
1347	! 4 SILT .485 .255 .179 5.30 4.418 6 .105 0.8 .153 NP89 does not have Silt so mapped to Silt Loam
1348	! 5 LOAM .451 .240 .155 5.39 4.111 6 .148 0.8 .191
1349	! 6 SND CLY LM .420 .255 .175 7.12 3.670 6 .135 0.8 .213
1350	! 7 SLT CLY LM .477 .322 .218 7.75 3.593 8 .127 0.4 .385
1351	! 8 CLAY LOAM .476 .325 .250 8.52 3.995 10 .084 0.6 .227
1352	! 9 SANDY CLAY .426 .310 .219 10.40 3.058 8 .139 0.3 .421
1353	! 10 SILTY CLAY .482 .370 .283 10.40 3.729 10 .075 0.3 .375
1354	! 11 CLAY .482 .367 .286 11.40 3.600 12 .083 0.3 .342
1355	!------------------------------------------------------------------------
1356	!
1357	!.........Arguments
1358	INTEGER, PARAMETER :: NSCAT = 16 ! max. soil types
1359	INTEGER, PARAMETER :: NSCATMIN = 16 ! min soil types
1360
1361	INTEGER, INTENT(IN) :: WEIGHT, ITIMESTEP, PXLSM_SMOIS_INIT
1362	REAL, INTENT(IN) :: MAVAIL
1363	REAL, DIMENSION(1:NSCAT), INTENT(IN) :: SOILCBOT
1364	REAL, INTENT(OUT) :: FWSAT,FWFC,FWWLT,FB,FCGSAT, &
1365	FJP,FAS,FC2R,FC1SAT
1366	REAL, INTENT(INOUT) :: W2
1367
1368
1369	INTEGER, INTENT(OUT) :: ISTI
1370	!........Local variables
1371	CHARACTER*4, AVCLASS
1372	CHARACTER*4, DIMENSION( 1: NSCAT ) :: TEXID
1373	!... Integer
1374	INTEGER:: S
1375	!... Real
1376	REAL:: TFRACBOT, CFRAC, SUMSND, SUMCLY, AVS, AVC, AVSLT
1377	REAL, DIMENSION( 1: NSCAT ) :: WSAT, WFC, WWLT, B, CGSAT, AS, &
1378	JP, C2R, C1SAT
1379
1380	REAL, DIMENSION( 1: NSCATMIN ) :: SAND, CLAY
1381	!.......... DATA statement for SOIL PARAMETERS for the 11 soil types
1382	DATA SAND /92.5,80.5,61.1,19.6,4.0,40.0,57.1,11.3,26.8, &
1383	52.0,6.5,10.2,1.0,1.0,1.0,1.0/
1384	DATA CLAY/2.1,4.1,10.9,19.1,7.3,18.8,23.3,32.2,36.6, &
1385	43.0,46.2,58.8,1.0,1.0,1.0,1.0/
1386	DATA TEXID/'Sand','Lsan','Sloa','Sill','Silt','Loam','Sclo', &
1387	'Sicl','Cllo','Sacl','Sicy','Clay','Ormt','Wate', &
1388	'Bedr','Othe'/
1389
1390	!
1391	!-- WSAT is saturated soil moisture (M^3/M^3) (JN90)
1392	DATA WSAT / 0.395, 0.410, 0.435, 0.485, 0.451, 0.420, 0.477, &
1393	0.476, 0.426, 0.482, 0.482, 0.482, 0.482, 0.482, 0.482, 0.482 /
1394	!
1395	!-- WFC is field capacity (M^3/M^3) (JN90)
1396	DATA WFC / 0.135, 0.150, 0.195, 0.255, 0.240, 0.255, 0.322, &
1397	0.325, 0.310, 0.370, 0.367, 0.367, 0.367, 0.367, 0.367, 0.367 /
1398	!
1399	!-- WWLT is wilting point (M^3/M^3) (JN90)
1400	DATA WWLT / 0.068, 0.075, 0.114, 0.179, 0.155, 0.175, 0.218, &
1401	0.250, 0.219, 0.283, 0.286, 0.286, 0.286, 0.286, 0.286, 0.286 /
1402	!
1403	!-- B is slop of the retention curve (NP89)
1404	DATA B / 4.05, 4.38, 4.90, 5.30, 5.39, 7.12, 7.75, &
1405	8.52, 10.40, 10.40, 11.40, 11.40, 11.40, 11.40, 11.40, 11.40 /
1406	!
1407	!-- CGSAT is soil thermal coef. at saturation (10^-6 K M^2 J^-1) (NP89)
1408	DATA CGSAT / 3.222, 3.057, 3.560, 4.418, 4.111, 3.670, 3.593, &
1409	3.995, 3.058, 3.729, 3.600, 3.600, 3.600, 3.600, 3.600, 3.600 /
1410	!
1411	!-- JP is coefficient of WGEQ formulation (NP89)
1412	DATA JP / 4, 4, 4, 6, 6, 6, 8, &
1413	10, 8, 10, 12, 12, 12, 12, 12, 12 /
1414	!
1415	!-- AS is coefficient of WGEQ formulation (NP89)
1416	DATA AS / 0.387, 0.404, 0.219, 0.105, 0.148, 0.135, 0.127, &
1417	0.084, 0.139, 0.075, 0.083, 0.083, 0.083, 0.083, 0.083, 0.083 /
1418	!
1419	!-- C2R is the value of C2 for W2=0.5WSAT (NP89)
1420	DATA C2R / 3.9, 3.7, 1.8, 0.8, 0.8, 0.8, 0.4, &
1421	0.6, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3 /
1422	!
1423	!-- C1SAT is the value of C1 at saturation (NP89)
1424	DATA C1SAT / 0.082, 0.098, 0.132, 0.153, 0.191, 0.213, 0.385, &
1425	0.227, 0.421, 0.375, 0.342, 0.342, 0.342, 0.342, 0.342, 0.342 /
1426	!
1427	!-------------------------------Exicutable starts here--------------------
1428
1429	IF(WEIGHT.GE.1.0) THEN !Compute soil characteristics using weighting determined by fractional soil content
1430	FWSAT =0
1431	FWFC =0
1432	FWWLT =0
1433	FB =0
1434	FCGSAT=0
1435	FJP =0
1436	FAS =0
1437	FC2R =0
1438	FC1SAT=0
1439	TFRACBOT =0
1440	CFRAC=0
1441
1442	DO S=1,NSCAT
1443	IF(SOILCBOT(S).GE.CFRAC) THEN
1444	ISTI=S
1445	CFRAC=SOILCBOT(S)
1446	ENDIF
1447
1448	TFRACBOT=TFRACBOT+SOILCBOT(S)
1449
1450	FWSAT =FWSAT + WSAT(S) *SOILCBOT(S)
1451	FWFC =FWFC + WFC(S) *SOILCBOT(S)
1452	FWWLT =FWWLT + WWLT(S) *SOILCBOT(S)
1453	FB =FB + B(S) *SOILCBOT(S)
1454	FCGSAT=FCGSAT + CGSAT(S) *SOILCBOT(S)
1455	FJP =FJP + JP(S) *SOILCBOT(S)
1456	FAS =FAS + AS(S) *SOILCBOT(S)
1457	FC2R =FC2R + C2R(S) *SOILCBOT(S)
1458	FC1SAT=FC1SAT + C1SAT(S) *SOILCBOT(S)
1459
1460
1461	ENDDO
1462
1463	TFRACBOT = 1/TFRACBOT
1464	FWSAT =FWSAT * TFRACBOT
1465	FWFC =FWFC * TFRACBOT
1466	FWWLT =FWWLT * TFRACBOT
1467	FB =FB * TFRACBOT
1468	FCGSAT=FCGSAT * TFRACBOT
1469	FJP =FJP * TFRACBOT
1470	FAS =FAS * TFRACBOT
1471	FC2R =FC2R * TFRACBOT
1472	FC1SAT=FC1SAT * TFRACBOT
1473
1474	ELSE !Compute soil characteristics by sand and clay fraction
1475
1476	CFRAC = 0.0
1477	SUMSND = 0.0
1478	SUMCLY = 0.0
1479	TFRACBOT = 0.0
1480
1481	DO S = 1,12 !NSCAT
1482	TFRACBOT = TFRACBOT + SOILCBOT(S)
1483	SUMSND = SUMSND + SAND(S) * SOILCBOT(S)
1484	SUMCLY = SUMCLY + CLAY(S) * SOILCBOT(S)
1485
1486	IF(SOILCBOT(S).GE.CFRAC) THEN ! Find Dominant Category and fraction
1487	ISTI=S
1488	CFRAC=SOILCBOT(S)
1489	ENDIF
1490	ENDDO
1491
1492	IF(TFRACBOT.GT.0.001) THEN
1493	AVS = SUMSND / TFRACBOT
1494	AVC = SUMCLY / TFRACBOT
1495	AVSLT = 100 - AVS - AVC
1496
1497	IF(AVS.GT.(85.+ 0.5*AVC)) THEN
1498	AVCLASS= 'Sand'
1499	ISTI = 1
1500	ELSE IF(AVS.GT.(70.+ AVC)) THEN
1501	AVCLASS= 'Lsan'
1502	ISTI = 2
1503	ELSE IF((AVC.LT.20..AND.AVS.GT.52.) &
1504	.OR.(AVC.LE.7.5.AND.AVSLT.LT.50.)) THEN
1505	AVCLASS= 'Sloa'
1506	ISTI = 3
1507	ELSE IF(AVC.LT.35..AND.AVS.GT.45..AND.AVSLT.LT.28.) THEN
1508	AVCLASS= 'Sclo'
1509	ISTI = 6
1510	ELSE IF(AVC.GE.35..AND.AVS.GT.45.) THEN
1511	AVCLASS = 'Sacl'
1512	ISTI = 9
1513	ELSE IF(AVC.LT.27.0.AND.AVSLT.LT.50.) THEN
1514	AVCLASS= 'Loam'
1515	ISTI = 5
1516	ELSE IF(AVC.LT.12..AND.AVSLT.GT.80.) THEN
1517	AVCLASS = 'Silt'
1518	ISTI = 4
1519	ELSE IF(AVC.LT.27.) THEN
1520	AVCLASS = 'Sill'
1521	ISTI = 4
1522	ELSE IF(AVC.LT.40..AND.AVS.GT.20.) THEN
1523	AVCLASS = 'Cllo'
1524	ISTI = 8
1525	ELSE IF(AVC.LT.40.) THEN
1526	AVCLASS = 'Sicl'
1527	ISTI = 7
1528	ELSE IF(AVSLT.GE.40.) THEN
1529	AVCLASS = 'Sicy'
1530	ISTI = 10
1531	ELSE
1532	AVCLASS = 'Clay'
1533	ISTI = 11
1534	ENDIF
1535
1536	ELSE
1537	ISTI=5
1538	AVCLASS = TEXID(ISTI)
1539	ENDIF
1540	FWSAT =WSAT(ISTI)
1541	FWFC =WFC(ISTI)
1542	FWWLT =WWLT(ISTI)
1543	FB =B(ISTI)
1544	FCGSAT=CGSAT(ISTI)
1545	FJP =JP(ISTI)
1546	FAS =AS(ISTI)
1547	FC2R =C2R(ISTI)
1548	FC1SAT=C1SAT(ISTI)
1549
1550
1551	ENDIF
1552
1553	! Compute W2 using soil moisture availiability if pxlsm_smois_init (in namelist) is not zero
1554	IF (ITIMESTEP .EQ. 1 .AND. PXLSM_SMOIS_INIT .GT. 0) THEN
1555	W2 = FWWLT + (0.5(FWSAT+FWFC) - FWWLT) MAVAIL
1556	ENDIF
1557
1558	END SUBROUTINE soilprop
1559	!------------------------------------------------------------------------------------------
1560	!------------------------------------------------------------------------------------------
1561	SUBROUTINE PXSNOW (ITIMESTEP, ASNOW, CSNOW, SNOW, &
1562	SNOWH, SNUP, &
1563	ALB, SNOALB, SHDFAC, SHDMIN, &
1564	HC_SNOW, SNOW_FRA, SNOWC, SNOWALB)
1565	!------------------------------------------------------------------------------------------
1566	IMPLICIT NONE
1567	!*******************************************************************
1568	! Pleim-Xiu LSM Snow model
1569	!*******************************************************************
1570	!---------------------------------------------------------------------------------------------------
1571	! ITIMESTEP -- Model time step index
1572	! ASNOW -- Analyzed snow water equivalent (mm)
1573	! CSNOW -- Current snow water equivalent (mm)
1574	! SNOW -- Snow water equivalent in model (mm)
1575	! SNOWH -- Physical snow depth (m)
1576	! SNUP -- Physical snow depth (landuse dependent) where when below, snow cover is fractional
1577	!
1578	! HC_SNOW -- Heat capacity of snow as a function of depth
1579	! SNOW_FRA-- Factional snow area
1580	! IFSNOW -- Snow flag
1581	! SNOWALB -- Snow albedo
1582	!---------------------------------------------------------------------------------------------------
1583	!--- Arguments
1584	REAL, PARAMETER :: W2SN_CONV = 10.0
1585	REAL, PARAMETER :: CS_SNOWPACK = 2092.0
1586	REAL, PARAMETER :: SALP = 2.6
1587	!-----------------------------------------------
1588	INTEGER, INTENT(IN) :: ITIMESTEP
1589	REAL, INTENT(IN) :: ASNOW, CSNOW, SNUP, ALB, SNOALB, SHDFAC, SHDMIN
1590	REAL, INTENT(INOUT) :: SNOW, SNOWH, SNOWC
1591	REAL, INTENT(OUT) :: HC_SNOW, SNOW_FRA, SNOWALB
1592	!------------------------------------------------------------------------------------
1593
1594
1595	!-----------------------------------------------
1596	! Local variables
1597	!... Real
1598	REAL:: CONV_WAT2SNOW, CSNOWW, RHO_SNOWPACK, &
1599	LIQSN_RATIO, SNEQV, RSNOW
1600	!-----------------------------------------------
1601
1602	SNEQV=ASNOW*0.001 ! Snow water in meters
1603	RHO_SNOWPACK = 450 ! Snowpack density (kg/m3), this should be computed in the future
1604	LIQSN_RATIO = DENW/RHO_SNOWPACK ! Ratio of water density and snowpack density
1605
1606	CONV_WAT2SNOW = LIQSN_RATIO/1000 ! Conversion factor for snow liquid equiv. (mm) to snowpack (m)
1607
1608	SNOW = ASNOW ! Set snow water to analysis value for now, implement a nudging scheme later
1609	SNOWH= SNOW * CONV_WAT2SNOW ! Conversion of snow water (mm) to snow depth (m)
1610
1611
1612	! Is snow cover now present. The limit is 0.45 mm of water eqivalent or about 2 inches snow depth
1613	SNOWC = 0.0
1614	IF (SNOWH .GT. 0.005) SNOWC = 1.0
1615
1616	HC_SNOW = RHO_SNOWPACK * CS_SNOWPACK * SNOWH
1617
1618	IF (SNEQV < SNUP) THEN
1619	RSNOW = SNEQV / SNUP
1620	SNOW_FRA = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
1621	ELSE
1622	SNOW_FRA = 1.0
1623	END IF
1624
1625	SNOWC = SNOW_FRA
1626
1627	SNOWALB = ALB + SNOWC*(SNOALB-ALB)
1628
1629	END SUBROUTINE pxsnow
1630
1631	END MODULE module_sf_pxlsm
1632

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format