Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_pxlsm.F @ 3094

Last change on this file since 3094 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

Rev	Line
[2759]	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