Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_pxlsm.F @ 3567

Last change on this file since 3567 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
Property svn:executable set to ``*
File size: 75.5 KB

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

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

Download in other formats:

Original Format