Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_ruclsm.F @ 3584

Last change on this file since 3584 was 2759, checked in by aslmd, 3 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 182.3 KB

Rev	Line
[2759]	1	#define LSMRUC_DBG_LVL 3000
	2	!WRF:MODEL_LAYER:PHYSICS
	3	!
	4	MODULE module_sf_ruclsm
	5
	6	USE module_model_constants
	7	USE module_wrf_error
	8
	9	! VEGETATION PARAMETERS
	10	INTEGER :: LUCATS , BARE
	11	integer, PARAMETER :: NLUS=50
	12	CHARACTER*8 LUTYPE
	13	INTEGER, DIMENSION(1:NLUS) :: NROTBL
	14	real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, LAITBL, &
	15	ALBTBL, Z0TBL, LEMITBL, PCTBL, SHDTBL, MAXALB
	16	REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
	17
	18	! SOIL PARAMETERS
	19	INTEGER :: SLCATS
	20	INTEGER, PARAMETER :: NSLTYPE=30
	21	CHARACTER*8 SLTYPE
	22	REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,HC, &
	23	MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
	24
	25	! LSM GENERAL PARAMETERS
	26	INTEGER :: SLPCATS
	27	INTEGER, PARAMETER :: NSLOPE=30
	28	REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
	29	REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
	30	REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
	31	CZIL_DATA
	32
	33	CHARACTER*256 :: err_message
	34
	35
	36	CONTAINS
	37
	38	!-----------------------------------------------------------------
	39	SUBROUTINE LSMRUC( &
	40	DT,KTAU,NSL,ZS, &
	41	RAINBL,SNOW,SNOWH,SNOWC,FRZFRAC,frpcpn, &
	42	Z3D,P8W,T3D,QV3D,QC3D,RHO3D, & !p8W in [PA]
	43	GLW,GSW,EMISS,CHKLOWQ, CHS, &
	44	FLQC,FLHC,MAVAIL,CANWAT,VEGFRA,ALB,ZNT, &
	45	SNOALB,ALBBCK, & !new
	46	QSFC,QSG,QVG,QCG,SOILT1,TSNAV, &
	47	TBOT,IVGTYP,ISLTYP,XLAND,XICE, &
	48	CP,G0,LV,STBOLT, &
	49	SOILMOIS,SH2O,SMAVAIL,SMMAX, &
	50	TSO,SOILT,HFX,QFX,LH, &
	51	SFCRUNOFF,UDRUNOFF,SFCEXC, &
	52	SFCEVP,GRDFLX,ACSNOW, &
	53	SMFR3D,KEEPFR3DFLAG, &
	54	myj, &
	55	ids,ide, jds,jde, kds,kde, &
	56	ims,ime, jms,jme, kms,kme, &
	57	its,ite, jts,jte, kts,kte )
	58	!-----------------------------------------------------------------
	59	IMPLICIT NONE
	60	!-----------------------------------------------------------------
	61	!
	62	! The RUC LSM model is described in:
	63	! Smirnova, T.G., J.M. Brown, and S.G. Benjamin, 1997:
	64	! Performance of different soil model configurations in simulating
	65	! ground surface temperature and surface fluxes.
	66	! Mon. Wea. Rev. 125, 1870-1884.
	67	! Smirnova, T.G., J.M. Brown, and D. Kim, 2000: Parameterization of
	68	! cold-season processes in the MAPS land-surface scheme.
	69	! J. Geophys. Res. 105, 4077-4086.
	70	!-----------------------------------------------------------------
	71	!-- DT time step (second)
	72	! ktau - number of time step
	73	! NSL - number of soil layers
	74	! NZS - number of levels in soil
	75	! ZS - depth of soil levels (m)
	76	!-- RAINBL - accumulated rain in [mm] between the PBL calls
	77	!-- RAINNCV one time step grid scale precipitation (mm/step)
	78	! SNOW - snow water equivalent [mm]
	79	! FRAZFRAC - fraction of frozen precipitation
	80	!-- SNOWC flag indicating snow coverage (1 for snow cover)
	81	!-- Z3D heights (m)
	82	!-- P8W 3D pressure (Pa)
	83	!-- T3D temperature (K)
	84	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	85	! QC3D - 3D cloud water mixing ratio (Kg/Kg)
	86	! RHO3D - 3D air density (kg/m^3)
	87	!-- GLW downward long wave flux at ground surface (W/m^2)
	88	!-- GSW absorbed short wave flux at ground surface (W/m^2)
	89	!-- EMISS surface emissivity (between 0 and 1)
	90	! FLQC - surface exchange coefficient for moisture (kg/m^2/s)
	91	! FLHC - surface exchange coefficient for heat [W/m^2/s/degreeK]
	92	! SFCEXC - surface exchange coefficient for heat [m/s]
	93	! CANWAT - CANOPY MOISTURE CONTENT (mm)
	94	! VEGFRA - vegetation fraction (between 0 and 1)
	95	! ALB - surface albedo (between 0 and 1)
	96	! SNOALB - maximum snow albedo (between 0 and 1)
	97	! ALBBCK - snow-free albedo (between 0 and 1)
	98	! ZNT - roughness length [m]
	99	!-- TBOT soil temperature at lower boundary (K)
	100	! IVGTYP - USGS vegetation type (24 classes)
	101	! ISLTYP - STASGO soil type (16 classes)
	102	!-- XLAND land mask (1 for land, 2 for water)
	103	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
	104	!-- G0 acceleration due to gravity (m/s^2)
	105	!-- LV latent heat of melting (J/kg)
	106	!-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
	107	! SOILMOIS - soil moisture content (volumetric fraction)
	108	! TSO - soil temp (K)
	109	!-- SOILT surface temperature (K)
	110	!-- HFX upward heat flux at the surface (W/m^2)
	111	!-- QFX upward moisture flux at the surface (kg/m^2/s)
	112	!-- LH upward latent heat flux (W/m^2)
	113	! SFCRUNOFF - ground surface runoff [mm]
	114	! UDRUNOFF - underground runoff [mm]
	115	! SFCEVP - total evaporation in [kg/m^2]
	116	! GRDFLX - soil heat flux (W/m^2: negative, if downward from surface)
	117	! ACSNOW - accumulation of snow water [m]
	118	!-- CHKLOWQ - is either 0 or 1 (so far set equal to 1).
	119	!-- used only in MYJPBL.
	120	!-- ims start index for i in memory
	121	!-- ime end index for i in memory
	122	!-- jms start index for j in memory
	123	!-- jme end index for j in memory
	124	!-- kms start index for k in memory
	125	!-- kme end index for k in memory
	126	!-------------------------------------------------------------------------
	127	! INTEGER, PARAMETER :: nzss=5
	128	! INTEGER, PARAMETER :: nddzs=2*(nzss-2)
	129
	130	INTEGER, PARAMETER :: nvegclas=24+3
	131
	132	REAL, INTENT(IN ) :: DT
	133	LOGICAL, INTENT(IN ) :: myj,frpcpn
	134	INTEGER, INTENT(IN ) :: ktau, nsl, &
	135	ims,ime, jms,jme, kms,kme, &
	136	ids,ide, jds,jde, kds,kde, &
	137	its,ite, jts,jte, kts,kte
	138
	139	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	140	INTENT(IN ) :: QV3D, &
	141	QC3D, &
	142	p8w, &
	143	rho3D, &
	144	T3D, &
	145	z3D
	146
	147	REAL, DIMENSION( ims:ime , jms:jme ), &
	148	INTENT(IN ) :: RAINBL, &
	149	GLW, &
	150	GSW, &
	151	SNOALB, &
	152	ALBBCK, &
	153	FLHC, &
	154	FLQC, &
	155	CHS , &
	156	EMISS, &
	157	XICE, &
	158	XLAND, &
	159	VEGFRA, &
	160	TBOT
	161
	162	REAL, DIMENSION( 1:nsl), INTENT(IN ) :: ZS
	163
	164	REAL, DIMENSION( ims:ime , jms:jme ), &
	165	INTENT(INOUT) :: &
	166	SNOW, & !new
	167	SNOWH, &
	168	SNOWC, &
	169	CANWAT, & ! new
	170	ALB, &
	171	MAVAIL, &
	172	SFCEXC, &
	173	ZNT
	174
	175	REAL, DIMENSION( ims:ime , jms:jme ), &
	176	INTENT(IN ) :: &
	177	FRZFRAC
	178
	179	INTEGER, DIMENSION( ims:ime , jms:jme ), &
	180	INTENT(IN ) :: IVGTYP, &
	181	ISLTYP
	182
	183	REAL, INTENT(IN ) :: CP,G0,LV,STBOLT
	184
	185	REAL, DIMENSION( ims:ime , 1:nsl, jms:jme ) , &
	186	INTENT(INOUT) :: SOILMOIS,SH2O,TSO
	187
	188	REAL, DIMENSION( ims:ime, jms:jme ) , &
	189	INTENT(INOUT) :: SOILT, &
	190	HFX, &
	191	QFX, &
	192	LH, &
	193	SFCEVP, &
	194	SFCRUNOFF, &
	195	UDRUNOFF, &
	196	GRDFLX, &
	197	ACSNOW, &
	198	QVG, &
	199	QCG, &
	200	QSFC, &
	201	QSG, &
	202	CHKLOWQ, &
	203	SOILT1, &
	204	TSNAV
	205
	206	REAL, DIMENSION( ims:ime, jms:jme ) , &
	207	INTENT(INOUT) :: SMAVAIL, &
	208	SMMAX
	209
	210	REAL, DIMENSION( its:ite, jts:jte ) :: DEW, &
	211	PC, &
	212	RUNOFF1, &
	213	RUNOFF2, &
	214	EMISSL, &
	215	ZNTL, &
	216	LMAVAIL, &
	217	SMELT, &
	218	SNOH, &
	219	SNFLX, &
	220	SNOM, &
	221	EDIR, &
	222	EC, &
	223	ETT, &
	224	SUBLIM, &
	225	sflx, &
	226	EVAPL, &
	227	PRCPL, &
	228	XICED, &
	229	INFILTR
	230
	231	!--- soil/snow properties
	232	REAL, DIMENSION( ims:ime, 1:nsl, jms:jme) &
	233	:: KEEPFR3DFLAG, &
	234	SMFR3D
	235
	236	REAL &
	237	:: RHOCS, &
	238	RHOSN, &
	239	BCLH, &
	240	DQM, &
	241	KSAT, &
	242	PSIS, &
	243	QMIN, &
	244	QWRTZ, &
	245	REF, &
	246	WILT, &
	247	CANWATR, &
	248	SNOWFRAC, &
	249	SNHEI, &
	250	SNWE
	251
	252	REAL :: CN, &
	253	SAT,CW, &
	254	C1SN, &
	255	C2SN, &
	256	KQWRTZ, &
	257	KICE, &
	258	KWT
	259
	260
	261	REAL, DIMENSION(1:NSL) :: ZSMAIN, &
	262	ZSHALF, &
	263	DTDZS2
	264
	265	REAL, DIMENSION(1:2*(nsl-2)) :: DTDZS
	266
	267	REAL, DIMENSION(1:4001) :: TBQ
	268
	269
	270	REAL, DIMENSION( 1:nsl ) :: SOILM1D, &
	271	TSO1D, &
	272	SOILICE, &
	273	SOILIQW, &
	274	SMFRKEEP
	275
	276	REAL, DIMENSION( 1:nsl ) :: KEEPFR
	277
	278
	279	REAL :: RSM, &
	280	SNWEPRINT, &
	281	SNHEIPRINT
	282
	283	REAL :: PRCPMS, &
	284	NEWSNMS, &
	285	PATM, &
	286	TABS, &
	287	QVATM, &
	288	QCATM, &
	289	Q2SAT, &
	290	SATFLG, &
	291	CONFLX, &
	292	RHO, &
	293	QKMS, &
	294	TKMS, &
	295	INFILTRP
	296	REAL :: cq,r61,r273,arp,brp,x,evs,eis
	297
	298	INTEGER :: NROOT
	299	INTEGER :: ILAND,ISOIL
	300
	301	INTEGER, DIMENSION ( 1:nvegclas ) :: IFOREST
	302
	303	INTEGER :: I,J,K,NZS,NZS1,NDDZS
	304	INTEGER :: k1,l,k2,kp,km
	305
	306
	307	!-----------------------------------------------------------------
	308
	309	NZS=NSL
	310	NDDZS=2*(nzs-2)
	311
	312	!---- table TBQ is for resolution of balance equation in VILKA
	313	CQ=173.15-.05
	314	R273=1./273.15
	315	R61=6.1153*0.62198
	316	ARP=77455.*41.9/461.525
	317	BRP=64.*41.9/461.525
	318
	319	DO K=1,4001
	320	CQ=CQ+.05
	321	! TBQ(K)=R61EXP(ARP(R273-1./CQ)-BRPLOG(CQR273))
	322	EVS=EXP(17.67*(CQ-273.15)/(CQ-29.65))
	323	EIS=EXP(22.514-6.15E3/CQ)
	324	if(CQ.ge.273.15) then
	325	! tbq is in mb
	326	tbq(k) = R61*evs
	327	else
	328	tbq(k) = R61*eis
	329	endif
	330
	331	END DO
	332
	333	!--- Initialize soil/vegetation parameters
	334	!--- This is temporary until SI is added to mass coordinate ---!!!!!
	335
	336	#if ( NMM_CORE == 1 )
	337	if(ktau+1.eq.1) then
	338	#else
	339	if(ktau.eq.1) then
	340	#endif
	341	DO J=jts,jte
	342	DO i=its,ite
	343	do k=1,nsl
	344	! smfr3d (i,k,j)=soilmois(i,k,j)/900.*1.e3
	345	keepfr3dflag(i,k,j)=0.
	346	sh2o (i,k,j)=0.
	347	enddo
	348	!--- initializing to zero snow fraction
	349	snowc(i,j) = min(1.,snowh(i,j)/0.1)
	350	!--- initializing of snow temp
	351	soilt1(i,j)=soilt(i,j)
	352	tsnav(i,j) =0.5*(soilt(i,j)+tso(i,1,j))-273.
	353	qcg (i,j) =0.
	354	patm=P8w(i,kms,j)*1.e-2
	355	QSG (i,j) = QSN(SOILT(i,j),TBQ)/PATM
	356	qvg (i,j) = QSG(i,j)*mavail(i,j)
	357	! qvg (i,j) =qv3d(i,kms,j)
	358	qsfc(i,j) = qsg(i,j)/(1.+qsg(i,j))
	359	SMELT(i,j) = 0.
	360	SNOM (i,j) = 0.
	361	SNFLX(i,j) = 0.
	362	DEW (i,j) = 0.
	363	PC (i,j) = 0.
	364	zntl (i,j) = 0.
	365	RUNOFF1(i,j) = 0.
	366	RUNOFF2(i,j) = 0.
	367	emissl (i,j) = 0.
	368	! Temporarily!!!
	369	canwat(i,j)=0.
	370
	371	! For RUC LSM CHKLOWQ needed for MYJPBL should
	372	! 1 because is actual specific humidity at the surface, and
	373	! not the saturation value
	374	chklowq(i,j) = 1.
	375	infiltr(i,j) = 0.
	376	snoh (i,j) = 0.
	377	edir (i,j) = 0.
	378	ec (i,j) = 0.
	379	ett (i,j) = 0.
	380	sublim(i,j) = 0.
	381	sflx (i,j) = 0.
	382	evapl (i,j) = 0.
	383	prcpl (i,j) = 0.
	384	ENDDO
	385	ENDDO
	386
	387	do k=1,nsl
	388	soilice(k)=0.
	389	soiliqw(k)=0.
	390	enddo
	391	endif
	392
	393	!-----------------------------------------------------------------
	394
	395	PRCPMS = 0.
	396	! NROOT = 4
	397
	398
	399	DO J=jts,jte
	400
	401	DO i=its,ite
	402
	403	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	404	print *,' IN LSMRUC ','ims,ime,jms,jme,its,ite,jts,jte,nzs', &
	405	ims,ime,jms,jme,its,ite,jts,jte,nzs
	406	print *,' IVGTYP, ISLTYP ', ivgtyp(i,j),isltyp(i,j)
	407	print *,' MAVAIL ', mavail(i,j)
	408	print *,' SOILT,QVG,P8w',soilt(i,j),qvg(i,j),p8w(i,1,j)
	409	print *, 'LSMRUC, I,J,xland, QFX,HFX from SFCLAY',i,j,xland(i,j), &
	410	qfx(i,j),hfx(i,j)
	411	print *, ' GSW, GLW =',gsw(i,j),glw(i,j)
	412	print *, 'SOILT, TSO start of time step =',soilt(i,j),(tso(i,k,j),k=1,nsl)
	413	print *, 'SOILMOIS start of time step =',(soilmois(i,k,j),k=1,nsl)
	414	print *, 'SMFROZEN start of time step =',(smfr3d(i,k,j),k=1,nsl)
	415	print *, ' I,J=, after SFCLAY CHS,FLHC ',i,j,chs(i,j),flhc(i,j)
	416	print *, 'LSMRUC, IVGTYP,ISLTYP,ZNT,ALB = ', ivgtyp(i,j),isltyp(i,j),znt(i,j),alb(i,j),i,j
	417	print *, 'LSMRUC I,J,DT,RAINBL =',I,J,dt,RAINBL(i,j)
	418	print *, 'XLAND ---->, ivgtype,isoiltyp,i,j',xland(i,j),ivgtyp(i,j),isltyp(i,j),i,j
	419	ENDIF
	420
	421
	422	ILAND = IVGTYP(i,j)
	423	ISOIL = ISLTYP(I,J)
	424	TABS = T3D(i,kms,j)
	425	QVATM = QV3D(i,kms,j)
	426	QCATM = QC3D(i,kms,j)
	427	PATM = P8w(i,kms,j)*1.e-5
	428	!---- what height is the first level?---- check!!!!!
	429	!-- need to de-stagger from w levels to P levels
	430	CONFLX = Z3D(i,kms,j)
	431	! CONFLX = 0.5*Z3D(i,kms,j)
	432	! CONFLX = 5.
	433	RHO = RHO3D(I,kms,J)
	434	!--- 1*e-3 is to convert from mm/s to m/s
	435	IF(FRPCPN) THEN
	436	PRCPMS = (RAINBL(i,j)/DT1.e-3)(1-FRZFRAC(I,J))
	437	NEWSNMS = (RAINBL(i,j)/DT1.e-3)FRZFRAC(I,J)
	438	ELSE
	439	if (tabs.le.273.15) then
	440	PRCPMS = 0.
	441	NEWSNMS = RAINBL(i,j)/DT*1.e-3
	442	else
	443	PRCPMS = RAINBL(i,j)/DT*1.e-3
	444	NEWSNMS = 0.
	445	endif
	446	ENDIF
	447	!--- rooting depth is 5 levels for forests
	448	! if(iforest(ivgtyp(i,j)).eq.1) nroot=5
	449	!--- convert exchange coeff to [m/s]
	450	! QKMS=FLQC(I,J)/RHO/MAVAIL(I,J)
	451	if (myj) then
	452	QKMS=CHS(i,j)
	453	TKMS=CHS(i,j)
	454	else
	455	QKMS=FLQC(I,J)/RHO/MAVAIL(I,J)
	456	TKMS=FLHC(I,J)/RHO/CP
	457	endif
	458	!--- convert incoming snow and canwat from mm to m
	459	SNWE=SNOW(I,J)*1.E-3
	460	SNHEI=SNOWH(I,J)
	461	CANWATR=CANWAT(I,J)*1.E-3
	462	SNOWFRAC=SNOWC(I,J)
	463
	464	!-----
	465	zsmain(1)=0.
	466	zshalf(1)=0.
	467	do k=2,nzs
	468	zsmain(k)= zs(k)
	469	zshalf(k)=0.5*(zsmain(k-1) + zsmain(k))
	470	enddo
	471
	472
	473	!------------------------------------------------------------
	474	!----- DDZS and DSDZ1 are for implicit soilution of soil eqns.
	475	!-------------------------------------------------------------
	476	NZS1=NZS-1
	477	!-----
	478	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	479	print *,' DT,NZS1, ZSMAIN, ZSHALF --->', dt,nzs1,zsmain,zshalf
	480	ENDIF
	481
	482	DO K=2,NZS1
	483	K1=2*K-3
	484	K2=K1+1
	485	X=DT/2./(ZSHALF(K+1)-ZSHALF(K))
	486	DTDZS(K1)=X/(ZSMAIN(K)-ZSMAIN(K-1))
	487	DTDZS2(K-1)=X
	488	DTDZS(K2)=X/(ZSMAIN(K+1)-ZSMAIN(K))
	489	END DO
	490
	491
	492	CN=0.5 ! exponent
	493	SAT=0.0004 ! canopy water saturated
	494
	495	CW =4.183E6
	496
	497
	498	!--- Constants used in Johansen soil thermal
	499	!--- conductivity method
	500
	501	KQWRTZ=7.7
	502	KICE=2.2
	503	KWT=0.57
	504
	505	!***********************************************************************
	506	!--- Constants for snow density calculations C1SN and C2SN
	507
	508	c1sn=0.026
	509	! c1sn=0.01
	510	c2sn=21.
	511
	512	!***********************************************************************
	513
	514	NROOT= 4
	515	! ! rooting depth
	516
	517	if(SNOWH(i,j).gt.0.) then
	518	RHOSN = SNOW(i,j)/SNOWH(i,j)
	519	else
	520	RHOSN = 300.
	521	endif
	522
	523	!--- initializing soil and surface properties
	524	CALL SOILVEGIN ( ILAND,ISOIL,MYJ,IFOREST, &
	525	EMISSL(I,J),PC(i,j),ZNT(I,J),QWRTZ, &
	526	! EMISSL(I,J),PC(i,j),ZNTL(I,J),QWRTZ, &
	527	RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT )
	528
	529	!-- definition of number of soil levels in the rooting zone
	530	IF(iforest(ivgtyp(i,j)).ne.1) THEN
	531	!---- all vegetation types except evergreen and mixed forests
	532	do k=2,nzs
	533	if(zsmain(k).ge.0.4) then
	534	NROOT=K
	535	goto 111
	536	endif
	537	enddo
	538
	539	ELSE
	540	!---- evergreen and mixed forests
	541	do k=2,nzs
	542	if(zsmain(k).ge.1.1) then
	543	NROOT=K
	544	goto 111
	545	endif
	546	enddo
	547	ENDIF
	548	111 continue
	549
	550	!-----
	551	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	552	print *,' ZS, ZSMAIN, ZSHALF, CONFLX --->', zs,zsmain,zshalf,conflx
	553	print *,'NROOT, iforest, ivgtyp, i,j ', nroot,iforest(ivgtyp(i,j)),ivgtyp(I,J),I,J
	554	ENDIF
	555
	556	!*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
	557
	558
	559	IF((XLAND(I,J)-1.5).GE.0. .or. XICE(I,J).gt.0.5)THEN
	560	!-- Water or ice point
	561	SMAVAIL(I,J)=1.0
	562	SMMAX(I,J)=1.0
	563	! SNOW(I,J)=0.0
	564	LMAVAIL(I,J)=1.0
	565
	566	ILAND=16
	567	ISOIL=14
	568
	569	patm=P8w(i,kms,j)*1.e-2
	570	qvg (i,j) = QSN(SOILT(i,j),TBQ)/PATM
	571	qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
	572	CHKLOWQ(I,J)=1.
	573	Q2SAT=QSN(TABS,TBQ)/PATM
	574
	575	DO K=1,NZS
	576	SOILMOIS(I,K,J)=1.0
	577	TSO(I,K,J)= SOILT(I,J)
	578	ENDDO
	579
	580	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	581	PRINT*,' water point, I=',I, &
	582	'J=',J, 'SOILT=', SOILT(i,j)
	583	ENDIF
	584	!--- decide if this water point is ice:
	585	! if(tabs.le.271.) then
	586	if(xice(i,j).gt.0.5) then
	587	! if(soilt(i,j).le.271.or.xice(i,j).eq.1.) then
	588	! if(tabs.le.271.or.xice(i,j).eq.1.) then
	589	XICED(i,j)=1.
	590	else
	591	XICED(i,j)=0.
	592	endif
	593
	594	IF(XICED(I,J).NE.1.) SNOW(I,J)=0.
	595	IF(XICED(I,J).GT.0.5)THEN
	596	!-- Sea-ice case
	597	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	598	PRINT*,' sea-ice at water point, I=',I, &
	599	'J=',J
	600	ENDIF
	601	ILAND = 24
	602	ISOIL = 16
	603
	604	SMAVAIL(I,J)=1.0
	605	SMMAX(I,J)=1.0
	606	LMAVAIL(I,J)=1.0
	607	! RUC model approach is to keep sea-ice skin temperature equal to temperature
	608	! at first atmospheric level
	609	SOILT(I,J)=MIN(271.,TABS)
	610	ZNT(I,J)=0.011
	611
	612
	613	DO K=1,NZS
	614	SOILMOIS(I,K,J)=1.0
	615	TSO(I,K,J)= MIN(273.15,SOILT(I,J))
	616	ENDDO
	617	ENDIF
	618
	619	! for MYJ surface and PBL scheme
	620	if (myj) then
	621	IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
	622	! IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qsg(I,J))THEN
	623	SATFLG=0.
	624	ELSE
	625	SATFLG=1.0
	626	ENDIF
	627	else
	628	SATFLG=1.0
	629	endif
	630	QFX(I,J)=QFX(I,J)*SATFLG
	631
	632
	633	ELSE
	634
	635	!-- Land point
	636	! Attention!!!! RUC LSM uses soil moisture content minus residual (minimum
	637	! soil moisture content for a given soil type) as a state variable.
	638	! If the WRF model is initialized from the RUC background model, then the
	639	! soil moisture variable is consistent with the RUC LSM.
	640	! If the WRF model is initialized from another background model (ETA, GFS...)
	641	! then the residual value should be subtracted when the 1-d array of soil
	642	! moisture is initialized before the call to SFCTMP, and after SFCTMP qmin
	643	! should be added back in.
	644	!
	645	! soilm1d (k) = min(max(0.,soilmois(i,k,j)-qmin(i,j)),dqm(i,j))
	646
	647
	648	DO k=1,nzs
	649	! soilm1d - soil moisture content minus residual [m3/m3]
	650	soilm1d (k) = min(max(0.,soilmois(i,k,j)),dqm)
	651	tso1d (k) = tso(i,k,j)
	652	ENDDO
	653
	654	do k=1,nzs
	655	smfrkeep(k) = smfr3d(i,k,j)
	656	keepfr (k) = keepfr3dflag(i,k,j)
	657	enddo
	658
	659	! LMAVAIL(I,J)=max(0.00001,min(1.,soilmois(i,1,j)/(REF-QMIN)))
	660	LMAVAIL(I,J)=max(0.00001,min(1.,soilmois(i,1,j)/dqm))
	661
	662	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	663	print *,'LAND, i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO', &
	664	i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO
	665	print *,'CONFLX =',CONFLX
	666	print *,'SMFRKEEP,KEEPFR ',SMFRKEEP,KEEPFR
	667	ENDIF
	668
	669	!-----------------------------------------------------------------
	670	CALL SFCTMP (dt,ktau,conflx,i,j, &
	671	!--- input variables
	672	nzs,nddzs,nroot, &
	673	iland,isoil,xland(i,j),ivgtyp(i,j), &
	674	PRCPMS,NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN, &
	675	PATM,TABS,QVATM,QCATM,RHO, &
	676	GLW(I,J),GSW(I,J),EMISSL(I,J), &
	677	QKMS,TKMS,PC(I,J),LMAVAIL(I,J), &
	678	canwatr,vegfra(I,J),alb(I,J),znt(I,J), &
	679	snoalb(i,j),albbck(i,j), & !new
	680	myj, &
	681	!--- soil fixed fields
	682	QWRTZ, &
	683	rhocs,dqm,qmin,ref, &
	684	wilt,psis,bclh,ksat, &
	685	sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
	686	!--- constants
	687	cp,g0,lv,stbolt,cw,c1sn,c2sn, &
	688	KQWRTZ,KICE,KWT, &
	689	!--- output variables
	690	snweprint,snheiprint,rsm, &
	691	soilm1d,tso1d,smfrkeep,keepfr, &
	692	soilt(I,J),soilt1(i,j),tsnav(i,j),dew(I,J), &
	693	qvg(I,J),qsg(I,J),qcg(I,J),SMELT(I,J), &
	694	SNOH(I,J),SNFLX(I,J),SNOM(I,J),ACSNOW(I,J), &
	695	edir(I,J),ec(I,J),ett(I,J),qfx(I,J), &
	696	lh(I,J),hfx(I,J),sflx(I,J),sublim(I,J), &
	697	evapl(I,J),prcpl(I,J),runoff1(I,J), &
	698	runoff2(I,J),soilice,soiliqw,infiltrp)
	699	!-----------------------------------------------------------------
	700
	701	!*** DIAGNOSTICS
	702	!--- available and maximum soil moisture content in the soil
	703	!--- domain
	704	smavail(i,j) = 0.
	705	smmax (i,j) = 0.
	706
	707	do k=1,nzs-1
	708	smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* &
	709	(zshalf(k+1)-zshalf(k))
	710	smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
	711	(zshalf(k+1)-zshalf(k))
	712	enddo
	713
	714	smavail(i,j)=smavail(i,j)+(qmin+soilm1d(nzs))* &
	715	(zsmain(nzs)-zshalf(nzs))
	716	smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
	717	(zsmain(nzs)-zshalf(nzs))
	718
	719	!--- Convert the water unit into mm
	720	SFCRUNOFF(I,J) = SFCRUNOFF(I,J)+RUNOFF1(I,J)DT1000.0
	721	UDRUNOFF (I,J) = UDRUNOFF(I,J)+RUNOFF2(I,J)*1000.0
	722	SMAVAIL (I,J) = SMAVAIL(I,J) * 1000.
	723	SMMAX (I,J) = SMMAX(I,J) * 1000.
	724	SFCEXC (I,J) = TKMS
	725	! SFCLAY expects QSFC as saturation mixing ration at surface
	726	! QSFC(I,J) = QSG(I,J)
	727	! MYJSFC expects QSFC as actual specific humidity at the surface
	728	QSFC(I,J) = QVG(I,J)/(1.+QVG(I,J))
	729	patm=P8w(i,kms,j)*1.e-2
	730	Q2SAT=QSN(TABS,TBQ)/PATM
	731	! for MYJ surface and PBL scheme
	732	if (myj) then
	733	IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
	734	CHKLOWQ(I,J)=0.
	735	ELSE
	736	CHKLOWQ(I,J)=1.
	737	ENDIF
	738	else
	739	CHKLOWQ(I,J)=1.
	740	endif
	741
	742	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	743	if(CHKLOWQ(I,J).eq.0.) then
	744	print *,'i,j,CHKLOWQ', &
	745	i,j,CHKLOWQ(I,J)
	746	endif
	747	ENDIF
	748
	749	MAVAIL (i,j) = LMAVAIL(I,J)
	750	! SNOW is in [mm], SNWE is in [m]; CANWAT is in mm, CANWATR is in m
	751	SNOW (i,j) = SNWE*1000.
	752	SNOWH (I,J) = SNHEI
	753	CANWAT (I,J) = CANWATR*1000.
	754
	755	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	756	print *,' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
	757	ENDIF
	758	!!! QFX (I,J) = LH(I,J)/LV
	759	SFCEVP (I,J) = SFCEVP (I,J) + QFX (I,J) * DT
	760	GRDFLX (I,J) = -1. * sflx(I,J)
	761
	762	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	763	print *,' QFX after change, LH ', i,j, QFX(i,j),LH(I,J)
	764	ENDIF
	765	!--- SNOWC snow cover flag
	766	SNOWC(I,J)=SNOWFRAC
	767
	768	! IF(SNOWH(I,J).GT.0.02)THEN
	769	! SNOWC(I,J)=1.0
	770	! ELSE
	771	! SNOWC(I,J)=0.0
	772	! ENDIF
	773
	774	INFILTR(I,J) = INFILTRP
	775
	776	!--- get 3d soil fields
	777	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	778	print *,'LAND, i,j,tso1d,soilm1d - end of time step', &
	779	i,j,tso1d,soilm1d
	780	ENDIF
	781
	782	do k=1,nzs
	783
	784	soilmois(i,k,j) = soilm1d(k)
	785	sh2o (i,k,j) = soiliqw(k)
	786	tso(i,k,j) = tso1d(k)
	787	enddo
	788
	789	do k=1,nzs
	790	smfr3d(i,k,j) = smfrkeep(k)
	791	keepfr3dflag(i,k,j) = keepfr (k)
	792	enddo
	793
	794	ENDIF
	795
	796	ENDDO
	797
	798	ENDDO
	799
	800	!-----------------------------------------------------------------
	801	END SUBROUTINE LSMRUC
	802	!-----------------------------------------------------------------
	803
	804
	805
	806	SUBROUTINE SFCTMP (delt,ktau,conflx,i,j, &
	807	!--- input variables
	808	nzs,nddzs,nroot, &
	809	ILAND,ISOIL,XLAND,IVGTYP, &
	810	PRCPMS,NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN, &
	811	PATM,TABS,QVATM,QCATM,rho, &
	812	GLW,GSW,EMISS,QKMS,TKMS,PC, &
	813	MAVAIL,CST,VEGFRA,ALB,ZNT, &
	814	ALB_SNOW,ALB_SNOW_FREE, &
	815	MYJ, &
	816	!--- soil fixed fields
	817	QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
	818	sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
	819	!--- constants
	820	cp,g0,lv,stbolt,cw,c1sn,c2sn, &
	821	KQWRTZ,KICE,KWT, &
	822	!--- output variables
	823	snweprint,snheiprint,rsm, &
	824	soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
	825	tsnav,dew,qvg,qsg,qcg, &
	826	SMELT,SNOH,SNFLX,SNOM,ACSNOW, &
	827	edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
	828	evapl,prcpl,runoff1,runoff2,soilice, &
	829	soiliqw,infiltr)
	830	!-----------------------------------------------------------------
	831	IMPLICIT NONE
	832	!-----------------------------------------------------------------
	833
	834	!--- input variables
	835
	836	INTEGER, INTENT(IN ) :: i,j,nroot,ktau,nzs , &
	837	nddzs !nddzs=2*(nzs-2)
	838
	839	REAL, INTENT(IN ) :: DELT,CONFLX
	840	REAL, INTENT(IN ) :: C1SN,C2SN
	841	LOGICAL, INTENT(IN ) :: myj
	842	!--- 3-D Atmospheric variables
	843	REAL , &
	844	INTENT(IN ) :: PATM, &
	845	TABS, &
	846	QVATM, &
	847	QCATM
	848	REAL , &
	849	INTENT(IN ) :: GLW, &
	850	GSW, &
	851	PC, &
	852	ALB_SNOW, &
	853	ALB_SNOW_FREE, &
	854	VEGFRA, &
	855	XLAND, &
	856	RHO, &
	857	QKMS, &
	858	TKMS
	859
	860	INTEGER, INTENT(IN ) :: IVGTYP
	861	!--- 2-D variables
	862	REAL , &
	863	INTENT(INOUT) :: EMISS, &
	864	MAVAIL, &
	865	SNOWFRAC, &
	866	ALB, &
	867	CST
	868
	869	!--- soil properties
	870	REAL :: &
	871	RHOCS, &
	872	BCLH, &
	873	DQM, &
	874	KSAT, &
	875	PSIS, &
	876	QMIN, &
	877	QWRTZ, &
	878	REF, &
	879	SAT, &
	880	WILT
	881
	882	REAL, INTENT(IN ) :: CN, &
	883	CW, &
	884	CP, &
	885	G0, &
	886	LV, &
	887	STBOLT, &
	888	KQWRTZ, &
	889	KICE, &
	890	KWT
	891
	892	REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
	893	ZSHALF, &
	894	DTDZS2
	895
	896	REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
	897
	898	REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
	899
	900
	901	!--- input/output variables
	902	!-------- 3-d soil moisture and temperature
	903	REAL, DIMENSION( 1:nzs ) , &
	904	INTENT(INOUT) :: TS1D, &
	905	SOILM1D, &
	906	SMFRKEEP
	907	REAL, DIMENSION( 1:nzs ) , &
	908	INTENT(INOUT) :: KEEPFR
	909
	910
	911	INTEGER, INTENT(INOUT) :: ILAND,ISOIL
	912
	913	!-------- 2-d variables
	914	REAL , &
	915	INTENT(INOUT) :: DEW, &
	916	EDIR1, &
	917	EC1, &
	918	ETT1, &
	919	EETA, &
	920	EVAPL, &
	921	INFILTR, &
	922	RHOSN, &
	923	SUBLIM, &
	924	PRCPL, &
	925	QVG, &
	926	QSG, &
	927	QCG, &
	928	QFX, &
	929	HFX, &
	930	S, &
	931	RUNOFF1, &
	932	RUNOFF2, &
	933	ACSNOW, &
	934	SNWE, &
	935	SNHEI, &
	936	SMELT, &
	937	SNOM, &
	938	SNOH, &
	939	SNFLX, &
	940	SOILT, &
	941	SOILT1, &
	942	TSNAV, &
	943	ZNT
	944
	945	!-------- 1-d variables
	946	REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
	947	SOILIQW
	948
	949	REAL, INTENT(OUT) :: RSM, &
	950	SNWEPRINT, &
	951	SNHEIPRINT
	952	!--- Local variables
	953
	954	INTEGER :: K,ILNB
	955
	956	REAL :: BSN, XSN, RHONEWSN , &
	957	RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , &
	958	T3, UPFLUX, XINET
	959	REAL :: snhei_crit, keep_snow_albedo
	960
	961	REAL :: RNET,GSWNEW,EMISSN,ALBSN,ZNTSN
	962	REAL :: VEGFRAC
	963
	964	!-----------------------------------------------------------------
	965	integer, parameter :: ilsnow=99
	966
	967	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	968	print *,' in SFCTMP',i,j,nzs,nddzs,nroot, &
	969	SNWE,RHOSN,SNOM,SMELT,TS1D
	970	ENDIF
	971	! print *,' in SFCTMP',i,j,nzs,nddzs,nroot, &
	972	! IVGTYP,ISOIL,ILAND, &
	973	! PRCPMS,SNWE,RHOSN, &
	974	! PATM,TABS,QVATM,QCATM,rho
	975	! GLW,GSW,EMISS,QKMS,TKMS,PC, &
	976	! cst,vegfrac,alb,znt, &
	977	!--- soil fixed fields
	978	! QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
	979	! sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
	980	!--- constants
	981	! cp,g0,lv,stbolt,cw,c1sn,c2sn, &
	982	! KQWRTZ,KICE,KWT
	983
	984	NEWSN=0.
	985	RAINF = 0.
	986	RSM=0.
	987	INFILTR=0.
	988	VEGFRAC=0.01*VEGFRA
	989
	990	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	991	print *,'I,J,KTAU,QKMS,TKMS', i,j,ktau,qkms,tkms
	992	print *,'GSW, GLW, SOILT, STBOLT, EMISS', &
	993	GSW, GLW, SOILT, STBOLT, EMISS
	994	ENDIF
	995
	996
	997	SNHEI = SNWE * 1000. / RHOSN
	998	!--------------
	999	T3 = STBOLTSOILTSOILT*SOILT
	1000	UPFLUX = T3 *SOILT
	1001	XINET = EMISS*(GLW-UPFLUX)
	1002	RNET = GSW + XINET
	1003
	1004	!Calculate the amount (m) of fresh snow
	1005
	1006	if(snhei.gt.0.0081*1.e3/rhosn) then
	1007	!*** Correct snow density for current temperature (Koren et al. 1999)
	1008	BSN=delt/3600.c1snexp(0.08tsnav-c2snrhosn*1.e-3)
	1009	if(bsnsnwe100..lt.1.e-4) goto 777
	1010	XSN=rhosn(exp(bsnsnwe100.)-1.)/(bsnsnwe*100.)
	1011	rhosn=MIN(MAX(100.,XSN),400.)
	1012	! rhosn=MIN(MAX(50.,XSN),400.)
	1013	777 continue
	1014
	1015	else
	1016	rhosn =200.
	1017	rhonewsn =100.
	1018	endif
	1019
	1020	! IF(TABS.LE.273.15)THEN
	1021
	1022	newsn=newsnms*delt
	1023	!--- consider for now that all PRCPMS went into snow
	1024	! prcpms = 0.
	1025	!---- ACSNOW - accumulation of snow water [m]
	1026	acsnow=acsnow+newsn
	1027
	1028	IF(NEWSN.GE.1.E-8) THEN
	1029	!*** Calculate fresh snow density (t > -15C, else MIN value)
	1030	!*** Eq. 10 from Koren et al. (1999)
	1031
	1032	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	1033	print *, 'THERE IS NEW SNOW, newsn', newsn
	1034	ENDIF
	1035	if(tabs.lt.258.15) then
	1036	! rhonewsn=50.
	1037	rhonewsn=100.
	1038
	1039	else
	1040	rhonewsn=1.e3max((0.10+0.0017(Tabs-273.15+15.)**1.5) &
	1041	, 0.10)
	1042	! rhonewsn=1.e3max((0.05+0.0017(Tabs-273.15+15.)**1.5) &
	1043	! , 0.05)
	1044	rhonewsn=MIN(rhonewsn,400.)
	1045	! rhonewsn=100.
	1046	endif
	1047
	1048
	1049	!*** Define average snow density of the snow pack considering
	1050	!*** the amount of fresh snow (eq. 9 in Koren et al.(1999)
	1051	!*** without snow melt )
	1052	xsn=(rhosnsnwe+rhonewsnnewsn)/ &
	1053	(snwe+newsn)
	1054	rhosn=MIN(MAX(100.,XSN),400.)
	1055	! rhosn=MIN(MAX(50.,XSN),400.)
	1056
	1057	snwe=snwe+newsn
	1058	snhei=snwe*1.E3/rhosn
	1059	NEWSN=NEWSN*1.E3/rhosn
	1060	endif
	1061
	1062	! ELSE
	1063	!--- TABS is above freezing. Needed precip rates from microphysics
	1064	!--- to do a better job with mixed phase precip.
	1065
	1066	! NEWSN = 0.
	1067	!
	1068	! ENDIF
	1069
	1070	IF(PRCPMS.NE.0.) THEN
	1071
	1072	! PRCPMS is liquid precipitation rate
	1073	! RAINF is a flag used for calculation of rain water
	1074	! heat content contribution into heat budget equation. Rain's temperature
	1075	! is set equal to air temperature at the first atmospheric
	1076	! level.
	1077
	1078	RAINF=1.
	1079	ENDIF
	1080
	1081	! IF((XLAND-1.5).GE.0.)THEN
	1082	! IF(ILAND.EQ.16) THEN
	1083	! SNHEI=0.
	1084	! SNWE=0.
	1085	! ELSE
	1086
	1087	IF(SNHEI.GT.0.0) THEN
	1088	!--- Set of surface parameters should be changed to snow values for grid
	1089	!--- points where the snow cover exceeds snow threshold of 2 cm
	1090	EMISS = 0.91
	1091
	1092	! GSWNEW = GSW
	1093	! The following lines compute albedo depending on snow
	1094	! depth. For now commented out.
	1095	! alb_snow_free=0.2
	1096	! alb_snow=0.70
	1097	! SNHEI_CRIT=0.05
	1098
	1099	SNHEI_CRIT=0.01601*1.e3/rhosn
	1100	SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
	1101
	1102	KEEP_SNOW_ALBEDO = 0.
	1103	IF (NEWSN.GT.0.) KEEP_SNOW_ALBEDO = 1.
	1104
	1105	!--- GSW in-coming solar
	1106	GSWNEW=GSW/(1.-ALB)
	1107
	1108	ALB = MAX(keep_snow_albedo*alb_snow, &
	1109	MIN((alb_snow_free + &
	1110	(alb_snow - alb_snow_free) * &
	1111	(snhei/(2.*SNHEI_CRIT))), alb_snow))
	1112	!--- recompute absorbed solar radiation and net radiation
	1113	!--- for new value of albedo
	1114	gswnew=gswnew*(1.-alb)
	1115
	1116	XINET = EMISS*(GLW-UPFLUX)
	1117	RNET = GSWnew + XINET
	1118
	1119	CALL SNOWSOIL ( & !--- input variables
	1120	i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
	1121	ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
	1122	RHOSN,PATM,QVATM,QCATM, &
	1123	GLW,GSWnew,EMISS,RNET,IVGTYP, &
	1124	QKMS,TKMS,PC,CST, &
	1125	RHO,VEGFRAC,ALB,ZNT, &
	1126	MYJ, &
	1127	!--- soil fixed fields
	1128	QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
	1129	sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
	1130	!--- constants
	1131	lv,CP,G0,cw,stbolt,tabs, &
	1132	KQWRTZ,KICE,KWT, &
	1133	!--- output variables
	1134	ilnb,snweprint,snheiprint,rsm, &
	1135	soilm1d,ts1d,smfrkeep,keepfr, &
	1136	dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
	1137	SMELT,SNOH,SNFLX,SNOM,edir1,ec1,ett1,eeta, &
	1138	qfx,hfx,s,sublim,prcpl,runoff1,runoff2, &
	1139	mavail,soilice,soiliqw,infiltr )
	1140
	1141	if(snhei.eq.0.) then
	1142
	1143	! if(snhei.le.2.e-2) then
	1144	!--- all snow is melted
	1145	! gswnew=gswnew/(1.-alb)
	1146	alb=alb_snow_free
	1147	! gswnew=gswnew*(1.-alb)
	1148	endif
	1149
	1150	ELSE
	1151
	1152	snheiprint=0.
	1153	snweprint=0.
	1154
	1155	CALL SOIL( &
	1156	!--- input variables
	1157	i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
	1158	PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSW, &
	1159	EMISS,RNET,QKMS,TKMS,PC,cst,rho,vegfrac, &
	1160	!--- soil fixed fields
	1161	QWRTZ,rhocs,dqm,qmin,ref,wilt, &
	1162	psis,bclh,ksat,sat,cn, &
	1163	zsmain,zshalf,DTDZS,DTDZS2,tbq, &
	1164	!--- constants
	1165	lv,CP,G0,cw,stbolt,tabs, &
	1166	KQWRTZ,KICE,KWT, &
	1167	!--- output variables
	1168	soilm1d,ts1d,smfrkeep,keepfr, &
	1169	dew,soilt,qvg,qsg,qcg,edir1,ec1, &
	1170	ett1,eeta,qfx,hfx,s,evapl,prcpl,runoff1, &
	1171	runoff2,mavail,soilice,soiliqw, &
	1172	infiltr)
	1173
	1174	ENDIF
	1175	! ENDIF
	1176
	1177	!
	1178
	1179	! RETURN
	1180	! END
	1181	!---------------------------------------------------------------
	1182	END SUBROUTINE SFCTMP
	1183	!---------------------------------------------------------------
	1184
	1185
	1186	FUNCTION QSN(TN,T)
	1187	!****************************************************************
	1188	REAL, DIMENSION(1:4001), INTENT(IN ) :: T
	1189	REAL, INTENT(IN ) :: TN
	1190
	1191	REAL QSN, R,R1,R2
	1192	INTEGER I
	1193
	1194	R=(TN-173.15)/.05+1.
	1195	I=INT(R)
	1196	IF(I.GE.1) goto 10
	1197	I=1
	1198	R=1.
	1199	10 IF(I.LE.4000) GOTO 20
	1200	I=4000
	1201	R=4001.
	1202	20 R1=T(I)
	1203	R2=R-I
	1204	QSN=(T(I+1)-R1)*R2 + R1
	1205	! print *,' in QSN, I,R,R1,R2,T(I+1),TN, QSN', I,R,r1,r2,t(i+1),tn,QSN
	1206	! RETURN
	1207	! END
	1208	!-----------------------------------------------------------------------
	1209	END FUNCTION QSN
	1210	!------------------------------------------------------------------------
	1211
	1212
	1213	SUBROUTINE SOIL ( &
	1214	!--- input variables
	1215	i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,&
	1216	PRCPMS,RAINF,PATM,QVATM,QCATM, &
	1217	GLW,GSW,EMISS,RNET, &
	1218	QKMS,TKMS,PC,cst,rho,vegfrac, &
	1219	!--- soil fixed fields
	1220	QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
	1221	sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
	1222	!--- constants
	1223	xlv,CP,G0_P,cw,stbolt,TABS, &
	1224	KQWRTZ,KICE,KWT, &
	1225	!--- output variables
	1226	soilmois,tso,smfrkeep,keepfr, &
	1227	dew,soilt,qvg,qsg,qcg, &
	1228	edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
	1229	prcpl,runoff1,runoff2,mavail,soilice, &
	1230	soiliqw,infiltrp)
	1231
	1232	!*************************************************************
	1233	! Energy and moisture budget for vegetated surfaces
	1234	! without snow, heat diffusion amf Richards eqns. in
	1235	! soil
	1236	!
	1237	! DELT - time step (s)
	1238	! ktau - numver of time step
	1239	! CONFLX - depth of constant flux layer (m)
	1240	! J,I - the location of grid point
	1241	! IME, JME, KME, NZS - dimensions of the domain
	1242	! NROOT - number of levels within the root zone
	1243	! PRCPMS - precipitation rate in m/s
	1244	! PATM - pressure [bar]
	1245	! QVATM,QCATM - cloud and water vapor mixing ratio (kg/kg)
	1246	! at the first atm. level
	1247	! GLW, GSW - incoming longwave and absorbed shortwave
	1248	! radiation at the surface (W/m^2)
	1249	! EMISS,RNET - emissivity of the ground surface (0-1) and net
	1250	! radiation at the surface (W/m^2)
	1251	! QKMS - exchange coefficient for water vapor in the
	1252	! surface layer (m/s)
	1253	! TKMS - exchange coefficient for heat in the surface
	1254	! layer (m/s)
	1255	! PC - plant coefficient (resistance) (0-1)
	1256	! RHO - density of atmosphere near sueface (kg/m^3)
	1257	! VEGFRAC - greeness fraction
	1258	! RHOCS - volumetric heat capacity of dry soil
	1259	! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
	1260	! REF, WILT - field capacity soil moisture and the
	1261	! wilting point (m^3/m^3)
	1262	! PSIS - matrix potential at saturation (m)
	1263	! BCLH - exponent for Clapp-Hornberger parameterization
	1264	! KSAT - saturated hydraulic conductivity (m/s)
	1265	! SAT - maximum value of water intercepted by canopy (m)
	1266	! CN - exponent for calculation of canopy water
	1267	! ZSMAIN - main levels in soil (m)
	1268	! ZSHALF - middle of the soil layers (m)
	1269	! DTDZS,DTDZS2 - dt/(2.dzshalfdzmain) and dt/dzshalf in soil
	1270	! TBQ - table to define saturated mixing ration
	1271	! of water vapor for given temperature and pressure
	1272	! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
	1273	! DEW - dew in kg/m^2s
	1274	! SOILT - skin temperature (K)
	1275	! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
	1276	! water vapor and cloud at the ground
	1277	! surface, respectively (kg/kg)
	1278	! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
	1279	! canopy water, transpiration in kg m-2 s-1 and total
	1280	! evaporation in m s-1.
	1281	! QFX, HFX - latent and sensible heat fluxes (W/m^2)
	1282	! S - soil heat flux in the top layer (W/m^2)
	1283	! RUNOFF - surface runoff (m/s)
	1284	! RUNOFF2 - underground runoff (m)
	1285	! MAVAIL - moisture availability in the top soil layer (0-1)
	1286	! INFILTRP - infiltration flux from the top of soil domain (m/s)
	1287	!
	1288	!*****************************************************************
	1289	IMPLICIT NONE
	1290	!-----------------------------------------------------------------
	1291
	1292	!--- input variables
	1293
	1294	INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
	1295	nddzs !nddzs=2*(nzs-2)
	1296	INTEGER, INTENT(IN ) :: i,j,iland,isoil
	1297	REAL, INTENT(IN ) :: DELT,CONFLX
	1298	!--- 3-D Atmospheric variables
	1299	REAL, &
	1300	INTENT(IN ) :: PATM, &
	1301	QVATM, &
	1302	QCATM
	1303	!--- 2-D variables
	1304	REAL, &
	1305	INTENT(IN ) :: GLW, &
	1306	GSW, &
	1307	EMISS, &
	1308	RHO, &
	1309	PC, &
	1310	VEGFRAC, &
	1311	QKMS, &
	1312	TKMS
	1313
	1314	!--- soil properties
	1315	REAL, &
	1316	INTENT(IN ) :: RHOCS, &
	1317	BCLH, &
	1318	DQM, &
	1319	KSAT, &
	1320	PSIS, &
	1321	QMIN, &
	1322	QWRTZ, &
	1323	REF, &
	1324	WILT
	1325
	1326	REAL, INTENT(IN ) :: CN, &
	1327	CW, &
	1328	KQWRTZ, &
	1329	KICE, &
	1330	KWT, &
	1331	XLV, &
	1332	g0_p
	1333
	1334
	1335	REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
	1336	ZSHALF, &
	1337	DTDZS2
	1338
	1339	REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
	1340
	1341	REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
	1342
	1343
	1344	!--- input/output variables
	1345	!-------- 3-d soil moisture and temperature
	1346	REAL, DIMENSION( 1:nzs ) , &
	1347	INTENT(INOUT) :: TSO, &
	1348	SOILMOIS, &
	1349	SMFRKEEP
	1350
	1351	REAL, DIMENSION( 1:nzs ) , &
	1352	INTENT(INOUT) :: KEEPFR
	1353
	1354	!-------- 2-d variables
	1355	REAL, &
	1356	INTENT(INOUT) :: DEW, &
	1357	CST, &
	1358	EDIR1, &
	1359	EC1, &
	1360	ETT1, &
	1361	EETA, &
	1362	EVAPL, &
	1363	PRCPL, &
	1364	MAVAIL, &
	1365	QVG, &
	1366	QSG, &
	1367	QCG, &
	1368	RNET, &
	1369	QFX, &
	1370	HFX, &
	1371	S, &
	1372	SAT, &
	1373	RUNOFF1, &
	1374	RUNOFF2, &
	1375	SOILT
	1376
	1377	!-------- 1-d variables
	1378	REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
	1379	SOILIQW
	1380
	1381	!--- Local variables
	1382
	1383	REAL :: INFILTRP, transum , &
	1384	RAINF, PRCPMS , &
	1385	TABS, T3, UPFLUX, XINET
	1386	REAL :: CP,G0,LV,STBOLT,xlmelt,dzstop , &
	1387	can,epot,fac,fltot,ft,fq,hft , &
	1388	q1,ras,rhoice,sph , &
	1389	trans,zn,ci,cvw,tln,tavln,pi , &
	1390	DD1,CMC2MS,DRYCAN,WETCAN , &
	1391	INFMAX,RIW
	1392	REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
	1393	thdif,tranf,tav,soilmoism , &
	1394	soilicem,soiliqwm,detal , &
	1395	fwsat,lwsat,told,smold
	1396
	1397	REAL :: drip
	1398
	1399	INTEGER :: nzs1,nzs2,k
	1400
	1401	!-----------------------------------------------------------------
	1402
	1403	!-- define constants
	1404	! STBOLT=5.670151E-8
	1405	RHOICE=900.
	1406	CI=RHOICE*2100.
	1407	XLMELT=3.335E+5
	1408	cvw=cw
	1409
	1410	SAT=0.0004
	1411	prcpl=prcpms
	1412
	1413	!--- Initializing local arrays
	1414	DO K=1,NZS
	1415	TRANSP (K)=0.
	1416	soilmoism(k)=0.
	1417	soilice (k)=0.
	1418	soiliqw (k)=0.
	1419	soilicem (k)=0.
	1420	soiliqwm (k)=0.
	1421	lwsat (k)=0.
	1422	fwsat (k)=0.
	1423	tav (k)=0.
	1424	cap (k)=0.
	1425	thdif (k)=0.
	1426	diffu (k)=0.
	1427	hydro (k)=0.
	1428	tranf (k)=0.
	1429	detal (k)=0.
	1430	told (k)=0.
	1431	smold (k)=0.
	1432	ENDDO
	1433
	1434	NZS1=NZS-1
	1435	NZS2=NZS-2
	1436	dzstop=1./(zsmain(2)-zsmain(1))
	1437	RAS=RHO*1.E-3
	1438	RIW=rhoice*1.e-3
	1439
	1440	!--- Computation of volumetric content of ice in soil
	1441
	1442	DO K=1,NZS
	1443	!- main levels
	1444	tln=log(tso(k)/273.15)
	1445	if(tln.lt.0.) then
	1446	soiliqw(k)=(dqm+qmin)(XLMELT &
	1447	(tso(k)-273.15)/tso(k)/9.81/psis) &
	1448	**(-1./bclh)-qmin
	1449	soiliqw(k)=max(0.,soiliqw(k))
	1450	soiliqw(k)=min(soiliqw(k),soilmois(k))
	1451	soilice(k)=(soilmois(k)-soiliqw(k))/RIW
	1452
	1453	!---- melting and freezing is balanced, soil ice cannot increase
	1454	if(keepfr(k).eq.1.) then
	1455	soilice(k)=min(soilice(k),smfrkeep(k))
	1456	soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
	1457	endif
	1458
	1459	else
	1460	soilice(k)=0.
	1461	soiliqw(k)=soilmois(k)
	1462	endif
	1463
	1464	ENDDO
	1465
	1466	DO K=1,NZS1
	1467	!- middle of soil layers
	1468	tav(k)=0.5*(tso(k)+tso(k+1))
	1469	soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
	1470	tavln=log(tav(k)/273.15)
	1471
	1472	if(tavln.lt.0.) then
	1473	soiliqwm(k)=(dqm+qmin)(XLMELT &
	1474	(tav(k)-273.15)/tav(k)/9.81/psis) &
	1475	**(-1./bclh)-qmin
	1476	fwsat(k)=dqm-soiliqwm(k)
	1477	lwsat(k)=soiliqwm(k)+qmin
	1478	soiliqwm(k)=max(0.,soiliqwm(k))
	1479	soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
	1480	soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
	1481	!---- melting and freezing is balanced, soil ice cannot increase
	1482	if(keepfr(k).eq.1.) then
	1483	soilicem(k)=min(soilicem(k), &
	1484	0.5*(smfrkeep(k)+smfrkeep(k+1)))
	1485	soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
	1486	fwsat(k)=dqm-soiliqwm(k)
	1487	lwsat(k)=soiliqwm(k)+qmin
	1488	endif
	1489
	1490	else
	1491	soilicem(k)=0.
	1492	soiliqwm(k)=soilmoism(k)
	1493	lwsat(k)=dqm+qmin
	1494	fwsat(k)=0.
	1495	endif
	1496
	1497	ENDDO
	1498
	1499	do k=1,nzs
	1500	if(soilice(k).gt.0.) then
	1501	smfrkeep(k)=soilice(k)
	1502	else
	1503	smfrkeep(k)=soilmois(k)/riw
	1504	endif
	1505	enddo
	1506
	1507	!******************************************************************
	1508	! SOILPROP computes thermal diffusivity, and diffusional and
	1509	! hydraulic condeuctivities
	1510	!******************************************************************
	1511	CALL SOILPROP( &
	1512	!--- input variables
	1513	nzs,fwsat,lwsat,tav,keepfr, &
	1514	soilmois,soiliqw,soilice, &
	1515	soilmoism,soiliqwm,soilicem, &
	1516	!--- soil fixed fields
	1517	QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
	1518	!--- constants
	1519	riw,xlmelt,CP,G0_P,cvw,ci, &
	1520	kqwrtz,kice,kwt, &
	1521	!--- output variables
	1522	thdif,diffu,hydro,cap)
	1523
	1524	!********************************************************************
	1525	!--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
	1526
	1527	DRIP=0.
	1528	DD1=0.
	1529
	1530	FQ=QKMS
	1531
	1532	Q1=-QKMSRAS(QVATM - QSG)
	1533
	1534	DEW=0.
	1535	IF(QVATM.GE.QSG)THEN
	1536	DEW=FQ*(QVATM-QSG)
	1537	ENDIF
	1538	IF(DEW.NE.0.)THEN
	1539	DD1=CST+DELT(PRCPMS +DEWRAS)*vegfrac
	1540	ELSE
	1541	DD1=CST+ &
	1542	DELT(PRCPMS+RASFQ*(QVATM-QSG) &
	1543	(CST/SAT)CN)vegfrac
	1544	ENDIF
	1545
	1546	IF(DD1.LT.0.) DD1=0.
	1547	if(vegfrac.eq.0.)then
	1548	cst=0.
	1549	drip=0.
	1550	endif
	1551	IF (vegfrac.GT.0.) THEN
	1552	CST=DD1
	1553	IF(CST.GT.SAT) THEN
	1554	CST=SAT
	1555	DRIP=DD1-SAT
	1556	ENDIF
	1557	ENDIF
	1558
	1559	!--- WETCAN is the fraction of vegetated area covered by canopy
	1560	!--- water, and DRYCAN is the fraction of vegetated area where
	1561	!--- transpiration may take place.
	1562
	1563	WETCAN=(CST/SAT)**CN
	1564	DRYCAN=1.-WETCAN
	1565
	1566	! print *,'CST,DRIP',cst,drip
	1567
	1568	!**************************************************************
	1569	! TRANSF computes transpiration function
	1570	!**************************************************************
	1571	CALL TRANSF( &
	1572	!--- input variables
	1573	nzs,nroot,soiliqw,tabs, &
	1574	!--- soil fixed fields
	1575	dqm,qmin,ref,wilt,zshalf, &
	1576	!--- output variables
	1577	tranf,transum)
	1578
	1579
	1580	!--- Save soil temp and moisture from the beginning of time step
	1581	do k=1,nzs
	1582	told(k)=tso(k)
	1583	smold(k)=soilmois(k)
	1584	enddo
	1585
	1586	!**************************************************************
	1587	! SOILTEMP soilves heat budget and diffusion eqn. in soil
	1588	!**************************************************************
	1589
	1590	CALL SOILTEMP( &
	1591	!--- input variables
	1592	i,j,iland,isoil, &
	1593	delt,ktau,conflx,nzs,nddzs,nroot, &
	1594	PRCPMS,RAINF, &
	1595	PATM,TABS,QVATM,QCATM,EMISS,RNET, &
	1596	QKMS,TKMS,PC,rho,vegfrac, &
	1597	thdif,cap,drycan,wetcan, &
	1598	transum,dew,mavail, &
	1599	!--- soil fixed fields
	1600	dqm,qmin,bclh,zsmain,zshalf,DTDZS,tbq, &
	1601	!--- constants
	1602	xlv,CP,G0_P,cvw,stbolt, &
	1603	!--- output variables
	1604	tso,soilt,qvg,qsg,qcg)
	1605
	1606	!************************************************************************
	1607
	1608	!--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
	1609	ETT1=0.
	1610	DEW=0.
	1611
	1612	IF(QVATM.GE.QSG)THEN
	1613	DEW=QKMS*(QVATM-QSG)
	1614	DO K=1,NZS
	1615	TRANSP(K)=0.
	1616	ENDDO
	1617	ELSE
	1618	DO K=1,NROOT
	1619	TRANSP(K)=VEGFRACRASQKMS* &
	1620	(QVATM-QSG)* &
	1621	PCTRANF(K)DRYCAN/ZSHALF(NROOT+1)
	1622	IF(TRANSP(K).GT.0.) TRANSP(K)=0.
	1623	ETT1=ETT1-TRANSP(K)
	1624	ENDDO
	1625	DO k=nroot+1,nzs
	1626	transp(k)=0.
	1627	enddo
	1628	ENDIF
	1629
	1630	!-- Recalculating of volumetric content of frozen water in soil
	1631	DO K=1,NZS
	1632	!- main levels
	1633	tln=log(tso(k)/273.15)
	1634	if(tln.lt.0.) then
	1635	soiliqw(k)=(dqm+qmin)(XLMELT &
	1636	(tso(k)-273.15)/tso(k)/9.81/psis) &
	1637	**(-1./bclh)-qmin
	1638	soiliqw(k)=max(0.,soiliqw(k))
	1639	soiliqw(k)=min(soiliqw(k),soilmois(k))
	1640	soilice(k)=(soilmois(k)-soiliqw(k))/riw
	1641	!---- melting and freezing is balanced, soil ice cannot increase
	1642	if(keepfr(k).eq.1.) then
	1643	soilice(k)=min(soilice(k),smfrkeep(k))
	1644	soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
	1645	endif
	1646
	1647	else
	1648	soilice(k)=0.
	1649	soiliqw(k)=soilmois(k)
	1650	endif
	1651	ENDDO
	1652
	1653	INFMAX=999.
	1654	!--- The threshold when the infiltration stops is:
	1655	!--- volumetric content of unfrozen pores < 0.12
	1656	if((dqm+qmin-riw*soilicem(1)).lt.0.12) &
	1657	INFMAX=0.
	1658
	1659	!*************************************************************************
	1660	! SOILMOIST solves moisture budget (Smirnova et al., 1996, EQ.22,28)
	1661	! and Richards eqn.
	1662	!*************************************************************************
	1663	CALL SOILMOIST ( &
	1664	!-- input
	1665	delt,nzs,nddzs,DTDZS,DTDZS2, &
	1666	zsmain,zshalf,diffu,hydro, &
	1667	QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
	1668	QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
	1669	!-- soil properties
	1670	DQM,QMIN,REF,KSAT,RAS,INFMAX, &
	1671	!-- output
	1672	SOILMOIS,MAVAIL,RUNOFF1, &
	1673	RUNOFF2,INFILTRP)
	1674
	1675	!--- KEEPFR is 1 when the temperature and moisture in soil
	1676	!--- are both increasing. In this case soil ice should not
	1677	!--- be increasing according to the freezing curve.
	1678	!--- Some part of ice is melted, but additional water is
	1679	!--- getting frozen. Thus, only structure of frozen soil is
	1680	!--- changed, and phase changes are not affecting the heat
	1681	!--- transfer. This situation may happen when it rains on the
	1682	!--- frozen soil.
	1683
	1684	do k=1,nzs
	1685	if (soilice(k).gt.0.) then
	1686	if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
	1687	keepfr(k)=1.
	1688	else
	1689	keepfr(k)=0.
	1690	endif
	1691	endif
	1692	enddo
	1693	!--- THE DIAGNOSTICS OF SURFACE FLUXES
	1694
	1695	T3 = STBOLTSOILTSOILT*SOILT
	1696	UPFLUX = T3 *SOILT
	1697	XINET = EMISS*(GLW-UPFLUX)
	1698	RNET = GSW + XINET
	1699	HFT=-TKMSCPRHO*(TABS-SOILT)
	1700	Q1=-QKMSRAS(QVATM - QSG)
	1701	EDIR1 =-(1.-vegfrac)QKMSRAS* &
	1702	(QVATM-QVG)
	1703	IF (Q1.LE.0.) THEN
	1704	! --- condensation
	1705	EC1=0.
	1706	EDIR1=0.
	1707	ETT1=0.
	1708	! EETA=0.
	1709	QFX=- XLVRHODEW
	1710	EETA= QFX/XLV
	1711	ELSE
	1712	! --- evaporation
	1713	EC1 = Q1 * WETCAN
	1714	CMC2MS=CST/DELT
	1715	if(EC1.gt.CMC2MS) cst=0.
	1716	EC1=MIN(CMC2MS,EC1)*vegfrac
	1717	EETA = (EDIR1 + EC1 + ETT1)*1.E3
	1718	! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
	1719	QFX= XLV * EETA
	1720	ENDIF
	1721	EVAPL=QFX/XLV
	1722	S=THDIF(1)CAP(1)DZSTOP*(TSO(1)-TSO(2))
	1723	HFX=HFT
	1724	FLTOT=RNET-HFT-QFX-S
	1725
	1726	222 CONTINUE
	1727
	1728	1123 FORMAT(I5,8F12.3)
	1729	1133 FORMAT(I7,8E12.4)
	1730	123 format(i6,f6.2,7f8.1)
	1731	122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
	1732
	1733
	1734	! RETURN
	1735	! END
	1736	!-------------------------------------------------------------------
	1737	END SUBROUTINE SOIL
	1738	!-------------------------------------------------------------------
	1739
	1740
	1741	SUBROUTINE SNOWSOIL ( &
	1742	!--- input variables
	1743	i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
	1744	ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
	1745	RHOSN, &
	1746	PATM,QVATM,QCATM, &
	1747	GLW,GSW,EMISS,RNET,IVGTYP, &
	1748	QKMS,TKMS,PC,cst,rho,vegfrac,alb,znt, &
	1749	MYJ, &
	1750	!--- soil fixed fields
	1751	QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
	1752	sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
	1753	!--- constants
	1754	xlv,CP,G0_P,cw,stbolt,TABS, &
	1755	KQWRTZ,KICE,KWT, &
	1756	!--- output variables
	1757	ilnb,snweprint,snheiprint,rsm, &
	1758	soilmois,tso,smfrkeep,keepfr, &
	1759	dew,soilt,soilt1,tsnav, &
	1760	qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, &
	1761	edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
	1762	prcpl,runoff1,runoff2,mavail,soilice, &
	1763	soiliqw,infiltrp )
	1764
	1765	!***************************************************************
	1766	! Energy and moisture budget for snow, heat diffusion eqns.
	1767	! in snow and soil, Richards eqn. for soil covered with snow
	1768	!
	1769	! DELT - time step (s)
	1770	! ktau - numver of time step
	1771	! CONFLX - depth of constant flux layer (m)
	1772	! J,I - the location of grid point
	1773	! IME, JME, NZS - dimensions of the domain
	1774	! NROOT - number of levels within the root zone
	1775	! PRCPMS - precipitation rate in m/s
	1776	! NEWSNOW - pcpn in soilid form (m)
	1777	! SNHEI, SNWE - snow height and snow water equivalent (m)
	1778	! RHOSN - snow density (kg/m-3)
	1779	! PATM - pressure (bar)
	1780	! QVATM,QCATM - cloud and water vapor mixing ratio
	1781	! at the first atm. level (kg/kg)
	1782	! GLW, GSW - incoming longwave and absorbed shortwave
	1783	! radiation at the surface (W/m^2)
	1784	! EMISS,RNET - emissivity (0-1) of the ground surface and net
	1785	! radiation at the surface (W/m^2)
	1786	! QKMS - exchange coefficient for water vapor in the
	1787	! surface layer (m/s)
	1788	! TKMS - exchange coefficient for heat in the surface
	1789	! layer (m/s)
	1790	! PC - plant coefficient (resistance) (0-1)
	1791	! RHO - density of atmosphere near surface (kg/m^3)
	1792	! VEGFRAC - greeness fraction (0-1)
	1793	! RHOCS - volumetric heat capacity of dry soil (J/m^3/K)
	1794	! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
	1795	! REF, WILT - field capacity soil moisture and the
	1796	! wilting point (m^3/m^3)
	1797	! PSIS - matrix potential at saturation (m)
	1798	! BCLH - exponent for Clapp-Hornberger parameterization
	1799	! KSAT - saturated hydraulic conductivity (m/s)
	1800	! SAT - maximum value of water intercepted by canopy (m)
	1801	! CN - exponent for calculation of canopy water
	1802	! ZSMAIN - main levels in soil (m)
	1803	! ZSHALF - middle of the soil layers (m)
	1804	! DTDZS,DTDZS2 - dt/(2.dzshalfdzmain) and dt/dzshalf in soil
	1805	! TBQ - table to define saturated mixing ration
	1806	! of water vapor for given temperature and pressure
	1807	! ilnb - number of layers in snow
	1808	! rsm - liquid water inside snow pack (m)
	1809	! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
	1810	! DEW - dew in (kg/m^2 s)
	1811	! SOILT - skin temperature (K)
	1812	! SOILT1 - snow temperature at 7.5 cm depth (K)
	1813	! TSNAV - average temperature of snow pack (C)
	1814	! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
	1815	! water vapor and cloud at the ground
	1816	! surface, respectively (kg/kg)
	1817	! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
	1818	! canopy water, transpiration (kg m-2 s-1) and total
	1819	! evaporation in (m s-1).
	1820	! QFX, HFX - latent and sensible heat fluxes (W/m^2)
	1821	! S - soil heat flux in the top layer (W/m^2)
	1822	! SUBLIM - snow sublimation (kg/m^2/s)
	1823	! RUNOFF1 - surface runoff (m/s)
	1824	! RUNOFF2 - underground runoff (m)
	1825	! MAVAIL - moisture availability in the top soil layer (0-1)
	1826	! SOILICE - content of soil ice in soil layers (m^3/m^3)
	1827	! SOILIQW - lliquid water in soil layers (m^3/m^3)
	1828	! INFILTRP - infiltration flux from the top of soil domain (m/s)
	1829	! XINET - net long-wave radiation (W/m^2)
	1830	!
	1831	!*******************************************************************
	1832
	1833	IMPLICIT NONE
	1834	!-------------------------------------------------------------------
	1835	!--- input variables
	1836
	1837	INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
	1838	nddzs !nddzs=2*(nzs-2)
	1839	INTEGER, INTENT(IN ) :: i,j,isoil
	1840
	1841	REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
	1842	RAINF,NEWSNOW
	1843
	1844	LOGICAL, INTENT(IN ) :: myj
	1845
	1846	!--- 3-D Atmospheric variables
	1847	REAL, &
	1848	INTENT(IN ) :: PATM, &
	1849	QVATM, &
	1850	QCATM
	1851	!--- 2-D variables
	1852	REAL , &
	1853	INTENT(IN ) :: GLW, &
	1854	GSW, &
	1855	RHO, &
	1856	PC, &
	1857	VEGFRAC, &
	1858	QKMS, &
	1859	TKMS
	1860
	1861	INTEGER, INTENT(IN ) :: IVGTYP
	1862	!--- soil properties
	1863	REAL , &
	1864	INTENT(IN ) :: RHOCS, &
	1865	BCLH, &
	1866	DQM, &
	1867	KSAT, &
	1868	PSIS, &
	1869	QMIN, &
	1870	QWRTZ, &
	1871	REF, &
	1872	SAT, &
	1873	WILT
	1874
	1875	REAL, INTENT(IN ) :: CN, &
	1876	CW, &
	1877	XLV, &
	1878	G0_P, &
	1879	KQWRTZ, &
	1880	KICE, &
	1881	KWT
	1882
	1883
	1884	REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
	1885	ZSHALF, &
	1886	DTDZS2
	1887
	1888	REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
	1889
	1890	REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
	1891
	1892
	1893	!--- input/output variables
	1894	!-------- 3-d soil moisture and temperature
	1895	REAL, DIMENSION( 1:nzs ) , &
	1896	INTENT(INOUT) :: TSO, &
	1897	SOILMOIS, &
	1898	SMFRKEEP
	1899
	1900	REAL, DIMENSION( 1:nzs ) , &
	1901	INTENT(INOUT) :: KEEPFR
	1902
	1903
	1904	INTEGER, INTENT(INOUT) :: ILAND
	1905
	1906
	1907	!-------- 2-d variables
	1908	REAL , &
	1909	INTENT(INOUT) :: DEW, &
	1910	CST, &
	1911	EDIR1, &
	1912	EC1, &
	1913	ETT1, &
	1914	EETA, &
	1915	RHOSN, &
	1916	SUBLIM, &
	1917	PRCPL, &
	1918	ALB, &
	1919	EMISS, &
	1920	ZNT, &
	1921	MAVAIL, &
	1922	QVG, &
	1923	QSG, &
	1924	QCG, &
	1925	QFX, &
	1926	HFX, &
	1927	S, &
	1928	RUNOFF1, &
	1929	RUNOFF2, &
	1930	SNWE, &
	1931	SNHEI, &
	1932	SMELT, &
	1933	SNOM, &
	1934	SNOH, &
	1935	SNFLX, &
	1936	SOILT, &
	1937	SOILT1, &
	1938	SNOWFRAC, &
	1939	TSNAV
	1940
	1941	INTEGER, INTENT(INOUT) :: ILNB
	1942
	1943	!-------- 1-d variables
	1944	REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
	1945	SOILIQW
	1946
	1947	REAL, INTENT(OUT) :: RSM, &
	1948	SNWEPRINT, &
	1949	SNHEIPRINT
	1950	!--- Local variables
	1951
	1952
	1953	INTEGER :: nzs1,nzs2,k
	1954
	1955	REAL :: INFILTRP, RHONEWSN,TRANSUM , &
	1956	SNTH, NEWSN , &
	1957	TABS, T3, UPFLUX, XINET , &
	1958	BETA, SNWEPR,EPDT,PP
	1959	REAL :: CP,G0,LV,xlvm,STBOLT,xlmelt,dzstop , &
	1960	can,epot,fac,fltot,ft,fq,hft , &
	1961	q1,ras,rhoice,sph , &
	1962	trans,zn,ci,cvw,tln,tavln,pi , &
	1963	DD1,CMC2MS,DRYCAN,WETCAN , &
	1964	INFMAX,RIW,DELTSN,H,UMVEG
	1965
	1966	REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
	1967	thdif,tranf,tav,soilmoism , &
	1968	soilicem,soiliqwm,detal , &
	1969	fwsat,lwsat,told,smold
	1970	REAL :: drip
	1971
	1972	REAL :: RNET
	1973
	1974	!-----------------------------------------------------------------
	1975
	1976	cvw=cw
	1977	XLMELT=3.335E+5
	1978	!-- the next line calculates heat of sublimation of water vapor
	1979	XLVm=XLV+XLMELT
	1980	! STBOLT=5.670151E-8
	1981
	1982	!--- SNOW flag -- 99
	1983	ILAND=99
	1984
	1985	!--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
	1986	!--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
	1987	!--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
	1988	!--- computed using SNWE=0.03 m and current snow density.
	1989	!--- SNTH - the threshold below which the snow layer is combined with
	1990	!--- the top soil layer. SNTH is computed using snwe=0.016 m, and
	1991	!--- equals 4 cm for snow density 400 kg/m^3.
	1992
	1993	DELTSN=0.0301*1.e3/rhosn
	1994	snth=0.01601*1.e3/rhosn
	1995
	1996	RHOICE=900.
	1997	CI=RHOICE*2100.
	1998	RAS=RHO*1.E-3
	1999	RIW=rhoice*1.e-3
	2000	MAVAIL=1.
	2001	RSM=0.
	2002
	2003	DO K=1,NZS
	2004	TRANSP (K)=0.
	2005	soilmoism (k)=0.
	2006	soiliqwm (k)=0.
	2007	soilice (k)=0.
	2008	soilicem (k)=0.
	2009	lwsat (k)=0.
	2010	fwsat (k)=0.
	2011	tav (k)=0.
	2012	cap (k)=0.
	2013	diffu (k)=0.
	2014	hydro (k)=0.
	2015	thdif (k)=0.
	2016	tranf (k)=0.
	2017	detal (k)=0.
	2018	told (k)=0.
	2019	smold (k)=0.
	2020	ENDDO
	2021
	2022	snweprint=0.
	2023	snheiprint=0.
	2024	prcpl=prcpms
	2025
	2026	!*** DELTSN is the depth of the top layer of snow where
	2027	!*** there is a temperature gradient, the rest of the snow layer
	2028	!*** is considered to have constant temperature
	2029
	2030
	2031	NZS1=NZS-1
	2032	NZS2=NZS-2
	2033	DZSTOP=1./(zsmain(2)-zsmain(1))
	2034
	2035	!----- THE CALCULATION OF THERMAL DIFFUSIVITY, DIFFUSIONAL AND ---
	2036	!----- HYDRAULIC CONDUCTIVITY (SMIRNOVA ET AL. 1996, EQ.2,5,6) ---
	2037	!tgs - the following loop is added to define the amount of frozen
	2038	!tgs - water in soil if there is any
	2039	DO K=1,NZS
	2040
	2041	tln=log(tso(k)/273.15)
	2042	if(tln.lt.0.) then
	2043	soiliqw(k)=(dqm+qmin)(XLMELT &
	2044	(tso(k)-273.15)/tso(k)/9.81/psis) &
	2045	**(-1./bclh)-qmin
	2046	soiliqw(k)=max(0.,soiliqw(k))
	2047	soiliqw(k)=min(soiliqw(k),soilmois(k))
	2048	soilice(k)=(soilmois(k)-soiliqw(k))/riw
	2049
	2050	!---- melting and freezing is balanced, soil ice cannot increase
	2051	if(keepfr(k).eq.1.) then
	2052	soilice(k)=min(soilice(k),smfrkeep(k))
	2053	soiliqw(k)=max(0.,soilmois(k)-soilice(k)rhoice1.e-3)
	2054	endif
	2055
	2056	else
	2057	soilice(k)=0.
	2058	soiliqw(k)=soilmois(k)
	2059	endif
	2060
	2061	ENDDO
	2062
	2063	DO K=1,NZS1
	2064
	2065	tav(k)=0.5*(tso(k)+tso(k+1))
	2066	soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
	2067	tavln=log(tav(k)/273.15)
	2068
	2069	if(tavln.lt.0.) then
	2070	soiliqwm(k)=(dqm+qmin)(XLMELT &
	2071	(tav(k)-273.15)/tav(k)/9.81/psis) &
	2072	**(-1./bclh)-qmin
	2073	fwsat(k)=dqm-soiliqwm(k)
	2074	lwsat(k)=soiliqwm(k)+qmin
	2075	soiliqwm(k)=max(0.,soiliqwm(k))
	2076	soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
	2077	soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
	2078	!---- melting and freezing is balanced, soil ice cannot increase
	2079	if(keepfr(k).eq.1.) then
	2080	soilicem(k)=min(soilicem(k), &
	2081	0.5*(smfrkeep(k)+smfrkeep(k+1)))
	2082	soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
	2083	fwsat(k)=dqm-soiliqwm(k)
	2084	lwsat(k)=soiliqwm(k)+qmin
	2085	endif
	2086
	2087	else
	2088	soilicem(k)=0.
	2089	soiliqwm(k)=soilmoism(k)
	2090	lwsat(k)=dqm+qmin
	2091	fwsat(k)=0.
	2092
	2093	endif
	2094	ENDDO
	2095
	2096	do k=1,nzs
	2097	if(soilice(k).gt.0.) then
	2098	smfrkeep(k)=soilice(k)
	2099	else
	2100	smfrkeep(k)=soilmois(k)/riw
	2101	endif
	2102	enddo
	2103
	2104
	2105	! print *,'etaf,etal,etamf,etaml,lwsat,fwsat',
	2106	! 1 soilice,soiliqw,soilicem,soiliqwm,lwsat,fwsat
	2107
	2108	!******************************************************************
	2109	! SOILPROP computes thermal diffusivity, and diffusional and
	2110	! hydraulic condeuctivities
	2111	!******************************************************************
	2112	CALL SOILPROP( &
	2113	!--- input variables
	2114	nzs,fwsat,lwsat,tav,keepfr, &
	2115	soilmois,soiliqw,soilice, &
	2116	soilmoism,soiliqwm,soilicem, &
	2117	!--- soil fixed fields
	2118	QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
	2119	!--- constants
	2120	riw,xlmelt,CP,G0_P,cvw,ci, &
	2121	kqwrtz,kice,kwt, &
	2122	!--- output variables
	2123	thdif,diffu,hydro,cap)
	2124
	2125	!********************************************************************
	2126	!--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
	2127
	2128	DRIP=0.
	2129	SMELT=0.
	2130	DD1=0.
	2131	H=1.
	2132
	2133	FQ=QKMS
	2134
	2135
	2136	!--- If vegfrac.ne.0. then part of falling snow can be
	2137	!--- intercepted by the canopy.
	2138
	2139	DEW=0.
	2140	UMVEG=1.-vegfrac
	2141	EPOT = -FQ*(QVATM-QSG)
	2142
	2143	IF(vegfrac.EQ.0.) then
	2144	cst=0.
	2145	drip=0.
	2146	ELSE
	2147	IF(EPOT.GE.0.) THEN
	2148	! Evaporation
	2149	DD1=CST+(NEWSNOWRHOSN1.E-3 &
	2150	!-- need to think more if we want this change....
	2151	-DELT(RASEPOT &
	2152	! -DELT(-PRCPMS+RASEPOT &
	2153	(CST/SAT)CN)) vegfrac
	2154	ELSE
	2155	! Sublimation
	2156	DEW = - EPOT
	2157	DD1=CST+(NEWSNOWRHOSN1.E-3+delt*( &
	2158	! DD1=CST+(NEWSNOWRHOSN1.E-3+delt*(PRCPMS &
	2159	+DEWRAS)) vegfrac
	2160	ENDIF
	2161
	2162	IF(DD1.LT.0.) DD1=0.
	2163	IF (vegfrac.GT.0.) THEN
	2164	CST=DD1
	2165	IF(CST.GT.SAT*vegfrac) THEN
	2166	CST=SAT*vegfrac
	2167	DRIP=DD1-SAT*vegfrac
	2168	ENDIF
	2169	ENDIF
	2170
	2171
	2172	!--- In SFCTMP NEWSNOW is added to SNHEI as if there is no vegetation
	2173	!--- With vegetation part of NEWSNOW can be intercepted by canopy until
	2174	!--- the saturation is reached. After the canopy saturation is reached
	2175	!--- DRIP in the solid form will be added to SNOW cover.
	2176
	2177	!!! 4 Nov 07 SNWE=(SNHEI-vegfracNEWSNOW)RHOSN*1.E-3 &
	2178	!!! 4 Nov 07 + DRIP &
	2179	! - 10% of liquid precip could be added to snow water
	2180	! - this is based on SnowMIP2.
	2181	! - something more intelligent should be done to liquid water
	2182	SNWE = SNWE &
	2183	+0.10prcpmsdelt
	2184
	2185
	2186	ENDIF
	2187
	2188	DRIP=0.
	2189	SNHEI=SNWE*1.e3/RHOSN
	2190	SNWEPR=SNWE
	2191
	2192	! check if all snow can evaporate during DT
	2193	BETA=1.
	2194	EPDT = EPOT * RAS DELTUMVEG
	2195	IF(SNWEPR.LE.EPDT) THEN
	2196	BETA=SNWEPR/max(1.e-8,EPDT)
	2197	SNWE=0.
	2198	SNHEI=0.
	2199	ENDIF
	2200
	2201	WETCAN=(CST/SAT)**CN
	2202	DRYCAN=1.-WETCAN
	2203
	2204	!**************************************************************
	2205	! TRANSF computes transpiration function
	2206	!**************************************************************
	2207	CALL TRANSF( &
	2208	!--- input variables
	2209	nzs,nroot,soiliqw,tabs, &
	2210	!--- soil fixed fields
	2211	dqm,qmin,ref,wilt,zshalf, &
	2212	!--- output variables
	2213	tranf,transum)
	2214
	2215	!--- Save soil temp and moisture from the beginning of time step
	2216	do k=1,nzs
	2217	told(k)=tso(k)
	2218	smold(k)=soilmois(k)
	2219	enddo
	2220
	2221	!**************************************************************
	2222	! SOILTEMP soilves heat budget and diffusion eqn. in soil
	2223	!**************************************************************
	2224
	2225	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	2226	print *, 'TSO before calling SNOWTEMP: ', tso
	2227	ENDIF
	2228	CALL SNOWTEMP( &
	2229	!--- input variables
	2230	i,j,iland,isoil, &
	2231	delt,ktau,conflx,nzs,nddzs,nroot, &
	2232	snwe,snwepr,snhei,newsnow,snowfrac, &
	2233	beta,deltsn,snth,rhosn, &
	2234	PRCPMS,RAINF, &
	2235	PATM,TABS,QVATM,QCATM, &
	2236	GLW,GSW,EMISS,RNET, &
	2237	QKMS,TKMS,PC,rho,vegfrac, &
	2238	thdif,cap,drycan,wetcan,cst, &
	2239	tranf,transum,dew,mavail, &
	2240	!--- soil fixed fields
	2241	dqm,qmin,psis,bclh, &
	2242	zsmain,zshalf,DTDZS,tbq, &
	2243	!--- constants
	2244	xlvm,CP,G0_P,cvw,stbolt, &
	2245	!--- output variables
	2246	snweprint,snheiprint,rsm, &
	2247	tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
	2248	smelt,snoh,snflx,ilnb)
	2249
	2250	!************************************************************************
	2251	!--- RECALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
	2252	DEW=0.
	2253	ETT1=0.
	2254	PP=PATM*1.E3
	2255	QSG= QSN(SOILT,TBQ)/PP
	2256	EPOT = -FQ*(QVATM-QSG)
	2257	IF(EPOT.GE.0.) THEN
	2258	! Evaporation
	2259	DO K=1,NROOT
	2260	TRANSP(K)=vegfracRASFQ*(QVATM-QSG) &
	2261	PCtranf(K)*DRYCAN/zshalf(NROOT+1)
	2262	IF(TRANSP(K).GT.0.) TRANSP(K)=0.
	2263	ETT1=ETT1-TRANSP(K)
	2264	ENDDO
	2265	DO k=nroot+1,nzs
	2266	transp(k)=0.
	2267	enddo
	2268
	2269	ELSE
	2270	! Sublimation
	2271	DEW=-EPOT
	2272	DO K=1,NZS
	2273	TRANSP(K)=0.
	2274	ENDDO
	2275	ETT1=0.
	2276	ENDIF
	2277
	2278	!-- recalculating of frozen water in soil
	2279	DO K=1,NZS
	2280	tln=log(tso(k)/273.15)
	2281	if(tln.lt.0.) then
	2282	soiliqw(k)=(dqm+qmin)(XLMELT &
	2283	(tso(k)-273.15)/tso(k)/9.81/psis) &
	2284	**(-1./bclh)-qmin
	2285	soiliqw(k)=max(0.,soiliqw(k))
	2286	soiliqw(k)=min(soiliqw(k),soilmois(k))
	2287	soilice(k)=(soilmois(k)-soiliqw(k))/riw
	2288	!---- melting and freezing is balanced, soil ice cannot increase
	2289	if(keepfr(k).eq.1.) then
	2290	soilice(k)=min(soilice(k),smfrkeep(k))
	2291	soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
	2292	endif
	2293
	2294	else
	2295	soilice(k)=0.
	2296	soiliqw(k)=soilmois(k)
	2297	endif
	2298	ENDDO
	2299
	2300	INFMAX=999.
	2301	!--- The threshold when the infiltration stops is:
	2302	!--- volumetric content of unfrozen pores < 0.12
	2303	soilicem(1)=0.5*(soilice(1)+soilice(2))
	2304	if((dqm+qmin-riw*soilicem(1)).lt.0.12) &
	2305	INFMAX=0.
	2306
	2307	!*************************************************************************
	2308	!--- TQCAN FOR SOLUTION OF MOISTURE BALANCE (Smirnova et al. 1996, EQ.22,28)
	2309	! AND TSO,ETA PROFILES
	2310	!*************************************************************************
	2311	CALL SOILMOIST ( &
	2312	!-- input
	2313	delt,nzs,nddzs,DTDZS,DTDZS2, &
	2314	zsmain,zshalf,diffu,hydro, &
	2315	! 4 Nov 07 QSG,QVG,QCG,QCATM,QVATM,-0.9*PRCPMS/(1.-vegfrac), &
	2316	QSG,QVG,QCG,QCATM,QVATM,-0.9*PRCPMS, &
	2317	! QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
	2318	0.,TRANSP,0., &
	2319	0.,SMELT,soilice,vegfrac, &
	2320	!-- soil properties
	2321	DQM,QMIN,REF,KSAT,RAS,INFMAX, &
	2322	!-- output
	2323	soilmois,MAVAIL,RUNOFF1, &
	2324	RUNOFF2,infiltrp)
	2325
	2326	! 4 Nov 07 - update CST for snow melt
	2327	if(snwe.ne.0.) then
	2328	CST=(1.-min(1.,smelt/snwe))*CST
	2329	else
	2330	CST=0.
	2331	endif
	2332
	2333	!-- Restore land-use parameters if snow is less than threshold
	2334	IF(SNHEI.EQ.0.) then
	2335	tsnav=soilt-273.15
	2336	CALL SNOWFREE(ivgtyp,myj,emiss, &
	2337	znt,iland)
	2338	smelt=smelt+snwe/delt
	2339	rsm=0.
	2340	! snwe=0.
	2341	ENDIF
	2342
	2343	SNOM=SNOM+SMELT*DELT
	2344
	2345	!--- KEEPFR is 1 when the temperature and moisture in soil
	2346	!--- are both increasing. In this case soil ice should not
	2347	!--- be increasing according to the freezing curve.
	2348	!--- Some part of ice is melted, but additional water is
	2349	!--- getting frozen. Thus, only structure of frozen soil is
	2350	!--- changed, and phase changes are not affecting the heat
	2351	!--- transfer. This situation may happen when it rains on the
	2352	!--- frozen soil.
	2353
	2354	do k=1,nzs
	2355	if (soilice(k).gt.0.) then
	2356	if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
	2357	keepfr(k)=1.
	2358	else
	2359	keepfr(k)=0.
	2360	endif
	2361	endif
	2362	enddo
	2363	!--- THE DIAGNOSTICS OF SURFACE FLUXES
	2364
	2365	T3 = STBOLTSOILTSOILT*SOILT
	2366	UPFLUX = T3 *SOILT
	2367	XINET = EMISS*(GLW-UPFLUX)
	2368	RNET = GSW + XINET
	2369	HFT=- TKMSCPRHO*(TABS-SOILT)
	2370	Q1 = - FQRAS (QVATM - QSG)
	2371	EDIR1 = Q1UMVEG BETA
	2372
	2373	IF (Q1.LT.0.) THEN
	2374	! --- condensation
	2375	EC1=0.
	2376	EDIR1=0.
	2377	ETT1=0.
	2378	! EETA=0.
	2379	DEW=FQ*(QVATM-QSG)
	2380	QFX= -XLVmRHODEW
	2381	sublim=QFX/XLVm
	2382	eeta=QFX/XLVm
	2383	ELSE
	2384	! --- evaporation
	2385	EC1 = Q1 * WETCAN
	2386	CMC2MS=CST/DELT
	2387	if(EC1.gt.CMC2MS) cst=0.
	2388	EC1=MIN(CMC2MS,EC1)*vegfrac
	2389	EETA = (EDIR1 + EC1 + ETT1)*1.E3
	2390	! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
	2391	QFX= XLVm * EETA
	2392	sublim=(EDIR1 + EC1)*1.E3
	2393	ENDIF
	2394	s=THDIF(1)CAP(1)dzstop*(tso(1)-tso(2))
	2395	HFX=HFT
	2396	FLTOT=RNET-HFT-QFX-S
	2397
	2398	222 CONTINUE
	2399
	2400	1123 FORMAT(I5,8F12.3)
	2401	1133 FORMAT(I7,8E12.4)
	2402	123 format(i6,f6.2,7f8.1)
	2403	122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
	2404
	2405
	2406	! RETURN
	2407	! END
	2408	!-------------------------------------------------------------------
	2409	END SUBROUTINE SNOWSOIL
	2410	!-------------------------------------------------------------------
	2411
	2412
	2413	SUBROUTINE SOILTEMP( &
	2414	!--- input variables
	2415	i,j,iland,isoil, &
	2416	delt,ktau,conflx,nzs,nddzs,nroot, &
	2417	PRCPMS,RAINF,PATM,TABS,QVATM,QCATM, &
	2418	EMISS,RNET, &
	2419	QKMS,TKMS,PC,RHO,VEGFRAC, &
	2420	THDIF,CAP,DRYCAN,WETCAN, &
	2421	TRANSUM,DEW,MAVAIL, &
	2422	!--- soil fixed fields
	2423	DQM,QMIN,BCLH, &
	2424	ZSMAIN,ZSHALF,DTDZS,TBQ, &
	2425	!--- constants
	2426	XLV,CP,G0_P,CVW,STBOLT, &
	2427	!--- output variables
	2428	TSO,SOILT,QVG,QSG,QCG)
	2429
	2430	!*************************************************************
	2431	! Energy budget equation and heat diffusion eqn are
	2432	! solved here and
	2433	!
	2434	! DELT - time step (s)
	2435	! ktau - numver of time step
	2436	! CONFLX - depth of constant flux layer (m)
	2437	! IME, JME, KME, NZS - dimensions of the domain
	2438	! NROOT - number of levels within the root zone
	2439	! PRCPMS - precipitation rate in m/s
	2440	! COTSO, RHTSO - coefficients for implicit solution of
	2441	! heat diffusion equation
	2442	! THDIF - thermal diffusivity (m^2/s)
	2443	! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
	2444	! water vapor and cloud at the ground
	2445	! surface, respectively (kg/kg)
	2446	! PATM - pressure [bar]
	2447	! QC3D,QV3D - cloud and water vapor mixing ratio
	2448	! at the first atm. level (kg/kg)
	2449	! EMISS,RNET - emissivity (0-1) of the ground surface and net
	2450	! radiation at the surface (W/m^2)
	2451	! QKMS - exchange coefficient for water vapor in the
	2452	! surface layer (m/s)
	2453	! TKMS - exchange coefficient for heat in the surface
	2454	! layer (m/s)
	2455	! PC - plant coefficient (resistance)
	2456	! RHO - density of atmosphere near surface (kg/m^3)
	2457	! VEGFRAC - greeness fraction (0-1)
	2458	! CAP - volumetric heat capacity (J/m^3/K)
	2459	! DRYCAN - dry fraction of vegetated area where
	2460	! transpiration may take place (0-1)
	2461	! WETCAN - fraction of vegetated area covered by canopy
	2462	! water (0-1)
	2463	! TRANSUM - transpiration function integrated over the
	2464	! rooting zone (m)
	2465	! DEW - dew in kg/m^2s
	2466	! MAVAIL - fraction of maximum soil moisture in the top
	2467	! layer (0-1)
	2468	! ZSMAIN - main levels in soil (m)
	2469	! ZSHALF - middle of the soil layers (m)
	2470	! DTDZS - dt/(2.dzshalfdzmain)
	2471	! TBQ - table to define saturated mixing ration
	2472	! of water vapor for given temperature and pressure
	2473	! TSO - soil temperature (K)
	2474	! SOILT - skin temperature (K)
	2475	!
	2476	!****************************************************************
	2477
	2478	IMPLICIT NONE
	2479	!-----------------------------------------------------------------
	2480
	2481	!--- input variables
	2482
	2483	INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
	2484	nddzs !nddzs=2*(nzs-2)
	2485	INTEGER, INTENT(IN ) :: i,j,iland,isoil
	2486	REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS, RAINF
	2487	REAL, INTENT(INOUT) :: DRYCAN,WETCAN,TRANSUM
	2488	!--- 3-D Atmospheric variables
	2489	REAL, &
	2490	INTENT(IN ) :: PATM, &
	2491	QVATM, &
	2492	QCATM
	2493	!--- 2-D variables
	2494	REAL , &
	2495	INTENT(IN ) :: &
	2496	EMISS, &
	2497	RHO, &
	2498	RNET, &
	2499	PC, &
	2500	VEGFRAC, &
	2501	DEW, &
	2502	QKMS, &
	2503	TKMS
	2504
	2505	!--- soil properties
	2506	REAL , &
	2507	INTENT(IN ) :: &
	2508	BCLH, &
	2509	DQM, &
	2510	QMIN
	2511
	2512	REAL, INTENT(IN ) :: CP, &
	2513	CVW, &
	2514	XLV, &
	2515	STBOLT, &
	2516	TABS, &
	2517	G0_P
	2518
	2519
	2520	REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
	2521	ZSHALF, &
	2522	THDIF, &
	2523	CAP
	2524
	2525	REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
	2526
	2527	REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
	2528
	2529
	2530	!--- input/output variables
	2531	!-------- 3-d soil moisture and temperature
	2532	REAL, DIMENSION( 1:nzs ) , &
	2533	INTENT(INOUT) :: TSO
	2534
	2535	!-------- 2-d variables
	2536	REAL , &
	2537	INTENT(INOUT) :: &
	2538	MAVAIL, &
	2539	QVG, &
	2540	QSG, &
	2541	QCG, &
	2542	SOILT
	2543
	2544
	2545	!--- Local variables
	2546
	2547	REAL :: x,x1,x2,x4,dzstop,can,ft,sph , &
	2548	tn,trans,umveg,denom
	2549
	2550	REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
	2551	PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
	2552	TDENOM
	2553
	2554	REAL :: C,CC,AA1,RHCS,H1
	2555
	2556	REAL, DIMENSION(1:NZS) :: cotso,rhtso
	2557
	2558	INTEGER :: nzs1,nzs2,k,k1,kn,kk
	2559
	2560	!-----------------------------------------------------------------
	2561
	2562
	2563	NZS1=NZS-1
	2564	NZS2=NZS-2
	2565	dzstop=1./(ZSMAIN(2)-ZSMAIN(1))
	2566
	2567	do k=1,nzs
	2568	cotso(k)=0.
	2569	rhtso(k)=0.
	2570	enddo
	2571	!******************************************************************************
	2572	! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
	2573	!******************************************************************************
	2574	! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
	2575	! h1=DTDZS(8)THDIF(nzs-1)(ZSHALF(nzs)-ZSHALF(nzs-1))/did
	2576	! cotso(1)=h1/(1.+h1)
	2577	! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
	2578	! 1 (1.+h1)
	2579	cotso(1)=0.
	2580	rhtso(1)=TSO(NZS)
	2581	DO 33 K=1,NZS2
	2582	KN=NZS-K
	2583	K1=2*KN-3
	2584	X1=DTDZS(K1)*THDIF(KN-1)
	2585	X2=DTDZS(K1+1)*THDIF(KN)
	2586	FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
	2587	-X2*(TSO(KN)-TSO(KN+1))
	2588	DENOM=1.+X1+X2-X2*cotso(K)
	2589	cotso(K+1)=X1/DENOM
	2590	rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
	2591	33 CONTINUE
	2592
	2593	!************************************************************************
	2594	!--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
	2595
	2596	RHCS=CAP(1)
	2597	H=MAVAIL
	2598	IF(DEW.NE.0.)THEN
	2599	DRYCAN=0.
	2600	WETCAN=1.
	2601	ENDIf
	2602	TRANS=PCTRANSUMDRYCAN/ZSHALF(NROOT+1)
	2603	CAN=WETCAN+TRANS
	2604	UMVEG=1.-VEGFRAC
	2605	FKT=TKMS
	2606	D1=cotso(NZS1)
	2607	D2=rhtso(NZS1)
	2608	TN=SOILT
	2609	D9=THDIF(1)RHCSdzstop
	2610	D10=TKMSCPRHO
	2611	R211=.5*CONFLX/DELT
	2612	R21=R211CPRHO
	2613	R22=.5/(THDIF(1)DELTdzstop**2)
	2614	R6=EMISS STBOLT.5TN*4
	2615	R7=R6/TN
	2616	D11=RNET+R6
	2617	TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
	2618	+RAINFCVWPRCPMS
	2619	FKQ=QKMS*RHO
	2620	R210=R211*RHO
	2621	C=VEGFRACFKQCAN
	2622	CC=C*XLV/TDENOM
	2623	AA=XLV(FKQUMVEG+R210)/TDENOM
	2624	BB=(D10TABS+R21TN+XLV(QVATM &
	2625	(FKQ*UMVEG+C) &
	2626	+R210QVG)+D11+D9(D2+R22*TN) &
	2627	+RAINFCVWPRCPMS*max(273.15,TABS) &
	2628	)/TDENOM
	2629	AA1=AA+CC
	2630	PP=PATM*1.E3
	2631	AA1=AA1/PP
	2632	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	2633	PRINT *,' VILKA-1'
	2634	print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
	2635	D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
	2636	print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT
	2637	print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
	2638	R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
	2639	print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
	2640	tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
	2641	ENDIF
	2642	CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
	2643	TQ2=QVATM+QCATM
	2644	TX2=TQ2*(1.-H)
	2645	Q1=TX2+H*QS1
	2646	IF(Q1.LT.QS1) GOTO 100
	2647	!--- if no saturation - goto 100
	2648	!--- if saturation - goto 90
	2649	90 QVG=QS1
	2650	QSG=QS1
	2651	TSO(1)=TS1
	2652	QCG=Q1-QS1
	2653	GOTO 200
	2654	100 BB=BB-AA*TX2
	2655	AA=(AA*H+CC)/PP
	2656	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	2657	PRINT *,' VILKA-2'
	2658	print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
	2659	D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
	2660	print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
	2661	R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
	2662
	2663	print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
	2664	tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
	2665	ENDIF
	2666
	2667	CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
	2668	Q1=TX2+H*QS1
	2669	IF(Q1.GT.QS1) GOTO 90
	2670	QSG=QS1
	2671	QVG=Q1
	2672	TSO(1)=TS1
	2673	QCG=0.
	2674	200 CONTINUE
	2675
	2676	!--- SOILT - skin temperature
	2677	SOILT=TS1
	2678
	2679	!---- Final solution for soil temperature - TSO
	2680	DO K=2,NZS
	2681	KK=NZS-K+1
	2682	TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
	2683	END DO
	2684
	2685	! return
	2686	! end
	2687	!--------------------------------------------------------------------
	2688	END SUBROUTINE SOILTEMP
	2689	!--------------------------------------------------------------------
	2690
	2691
	2692	SUBROUTINE SNOWTEMP( &
	2693	!--- input variables
	2694	i,j,iland,isoil, &
	2695	delt,ktau,conflx,nzs,nddzs,nroot, &
	2696	snwe,snwepr,snhei,newsnow,snowfrac, &
	2697	beta,deltsn,snth,rhosn, &
	2698	PRCPMS,RAINF, &
	2699	PATM,TABS,QVATM,QCATM, &
	2700	GLW,GSW,EMISS,RNET, &
	2701	QKMS,TKMS,PC,RHO,VEGFRAC, &
	2702	THDIF,CAP,DRYCAN,WETCAN,CST, &
	2703	TRANF,TRANSUM,DEW,MAVAIL, &
	2704	!--- soil fixed fields
	2705	DQM,QMIN,PSIS,BCLH, &
	2706	ZSMAIN,ZSHALF,DTDZS,TBQ, &
	2707	!--- constants
	2708	XLVM,CP,G0_P,CVW,STBOLT, &
	2709	!--- output variables
	2710	SNWEPRINT,SNHEIPRINT,RSM, &
	2711	TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, &
	2712	SMELT,SNOH,SNFLX,ILNB)
	2713
	2714	!********************************************************************
	2715	! Energy budget equation and heat diffusion eqn are
	2716	! solved here to obtain snow and soil temperatures
	2717	!
	2718	! DELT - time step (s)
	2719	! ktau - numver of time step
	2720	! CONFLX - depth of constant flux layer (m)
	2721	! IME, JME, KME, NZS - dimensions of the domain
	2722	! NROOT - number of levels within the root zone
	2723	! PRCPMS - precipitation rate in m/s
	2724	! COTSO, RHTSO - coefficients for implicit solution of
	2725	! heat diffusion equation
	2726	! THDIF - thermal diffusivity (W/m/K)
	2727	! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
	2728	! water vapor and cloud at the ground
	2729	! surface, respectively (kg/kg)
	2730	! PATM - pressure [bar]
	2731	! QCATM,QVATM - cloud and water vapor mixing ratio
	2732	! at the first atm. level (kg/kg)
	2733	! EMISS,RNET - emissivity (0-1) of the ground surface and net
	2734	! radiation at the surface (W/m^2)
	2735	! QKMS - exchange coefficient for water vapor in the
	2736	! surface layer (m/s)
	2737	! TKMS - exchange coefficient for heat in the surface
	2738	! layer (m/s)
	2739	! PC - plant coefficient (resistance)
	2740	! RHO - density of atmosphere near surface (kg/m^3)
	2741	! VEGFRAC - greeness fraction (0-1)
	2742	! CAP - volumetric heat capacity (J/m^3/K)
	2743	! DRYCAN - dry fraction of vegetated area where
	2744	! transpiration may take place (0-1)
	2745	! WETCAN - fraction of vegetated area covered by canopy
	2746	! water (0-1)
	2747	! TRANSUM - transpiration function integrated over the
	2748	! rooting zone (m)
	2749	! DEW - dew in kg/m^2/s
	2750	! MAVAIL - fraction of maximum soil moisture in the top
	2751	! layer (0-1)
	2752	! ZSMAIN - main levels in soil (m)
	2753	! ZSHALF - middle of the soil layers (m)
	2754	! DTDZS - dt/(2.dzshalfdzmain)
	2755	! TBQ - table to define saturated mixing ration
	2756	! of water vapor for given temperature and pressure
	2757	! TSO - soil temperature (K)
	2758	! SOILT - skin temperature (K)
	2759	!
	2760	!*********************************************************************
	2761
	2762	IMPLICIT NONE
	2763	!---------------------------------------------------------------------
	2764	!--- input variables
	2765
	2766	INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
	2767	nddzs !nddzs=2*(nzs-2)
	2768
	2769	INTEGER, INTENT(IN ) :: i,j,iland,isoil
	2770	REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
	2771	RAINF,NEWSNOW,DELTSN,SNTH , &
	2772	TABS,TRANSUM,SNWEPR
	2773
	2774	!--- 3-D Atmospheric variables
	2775	REAL, &
	2776	INTENT(IN ) :: PATM, &
	2777	QVATM, &
	2778	QCATM
	2779	!--- 2-D variables
	2780	REAL , &
	2781	INTENT(IN ) :: GLW, &
	2782	GSW, &
	2783	RHO, &
	2784	PC, &
	2785	VEGFRAC, &
	2786	QKMS, &
	2787	TKMS
	2788
	2789	!--- soil properties
	2790	REAL , &
	2791	INTENT(IN ) :: &
	2792	BCLH, &
	2793	DQM, &
	2794	PSIS, &
	2795	QMIN
	2796
	2797	REAL, INTENT(IN ) :: CP, &
	2798	CVW, &
	2799	STBOLT, &
	2800	XLVM, &
	2801	G0_P
	2802
	2803
	2804	REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
	2805	ZSHALF, &
	2806	THDIF, &
	2807	CAP, &
	2808	TRANF
	2809
	2810	REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
	2811
	2812	REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
	2813
	2814
	2815	!--- input/output variables
	2816	!-------- 3-d soil moisture and temperature
	2817	REAL, DIMENSION( 1:nzs ) , &
	2818	INTENT(INOUT) :: TSO
	2819
	2820
	2821	!-------- 2-d variables
	2822	REAL , &
	2823	INTENT(INOUT) :: DEW, &
	2824	CST, &
	2825	RHOSN, &
	2826	EMISS, &
	2827	MAVAIL, &
	2828	QVG, &
	2829	QSG, &
	2830	QCG, &
	2831	SNWE, &
	2832	SNHEI, &
	2833	SNOWFRAC, &
	2834	SMELT, &
	2835	SNOH, &
	2836	SNFLX, &
	2837	SOILT, &
	2838	SOILT1, &
	2839	TSNAV
	2840
	2841	REAL, INTENT(INOUT) :: DRYCAN, WETCAN
	2842
	2843	REAL, INTENT(OUT) :: RSM, &
	2844	SNWEPRINT, &
	2845	SNHEIPRINT
	2846	INTEGER, INTENT(OUT) :: ilnb
	2847	!--- Local variables
	2848
	2849
	2850	INTEGER :: nzs1,nzs2,k,k1,kn,kk
	2851
	2852	REAL :: x,x1,x2,x4,dzstop,can,ft,sph, &
	2853	tn,trans,umveg,denom
	2854
	2855	REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
	2856
	2857	REAL :: t3,upflux,xinet,ras, &
	2858	xlmelt,rhocsn,thdifsn, &
	2859	beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
	2860
	2861	REAL :: fso,fsn, &
	2862	FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
	2863	PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2, &
	2864	TDENOM,C,CC,AA1,RHCS,H1, &
	2865	tsob, snprim, sh1, sh2, &
	2866	smeltg,snohg,snodif,soh, &
	2867	CMC2MS,TNOLD,QGOLD,SNOHGNEW
	2868
	2869	REAL, DIMENSION(1:NZS) :: transp,cotso,rhtso
	2870	REAL :: edir1, &
	2871	ec1, &
	2872	ett1, &
	2873	eeta, &
	2874	s, &
	2875	qfx, &
	2876	hfx
	2877
	2878	REAL :: RNET,rsmfrac,soiltfrac,hsn
	2879
	2880	!-----------------------------------------------------------------
	2881
	2882	do k=1,nzs
	2883	transp (k)=0.
	2884	cotso (k)=0.
	2885	rhtso (k)=0.
	2886	enddo
	2887
	2888	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	2889	print *, 'SNOWTEMP: SNHEI,SNTH,SOILT1: ',SNHEI,SNTH,SOILT1,soilt
	2890	ENDIF
	2891	XLMELT=3.335E+5
	2892	RHOCSN=2090.* RHOSN
	2893	THDIFSN = 0.265/RHOCSN
	2894	RAS=RHO*1.E-3
	2895
	2896	SOILTFRAC=SOILT
	2897
	2898	SMELT=0.
	2899	SOH=0.
	2900	SMELTG=0.
	2901	SNOHG=0.
	2902	SNODIF=0.
	2903	RSM = 0.
	2904	RSMFRAC = 0.
	2905	fsn=1.
	2906	fso=0.
	2907	hsn=snhei
	2908
	2909	NZS1=NZS-1
	2910	NZS2=NZS-2
	2911
	2912	QGOLD=QVG
	2913	TNOLD=SOILT
	2914	DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
	2915
	2916	!******************************************************************************
	2917	! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
	2918	!******************************************************************************
	2919	! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
	2920	! h1=DTDZS(8)THDIF(nzs-1)(ZSHALF(nzs)-ZSHALF(nzs-1))/did
	2921	! cotso(1)=h1/(1.+h1)
	2922	! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
	2923	! 1 (1.+h1)
	2924
	2925	cotso(1)=0.
	2926	rhtso(1)=TSO(NZS)
	2927	DO 33 K=1,NZS2
	2928	KN=NZS-K
	2929	K1=2*KN-3
	2930	X1=DTDZS(K1)*THDIF(KN-1)
	2931	X2=DTDZS(K1+1)*THDIF(KN)
	2932	FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
	2933	-X2*(TSO(KN)-TSO(KN+1))
	2934	DENOM=1.+X1+X2-X2*cotso(K)
	2935	cotso(K+1)=X1/DENOM
	2936	rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
	2937	33 CONTINUE
	2938	!--- THE NZS element in COTSO and RHTSO will be for snow
	2939	!--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
	2940	IF(SNHEI.GE.SNTH) then
	2941	! if(snhei.le.DELTSN+DELTSN) then
	2942	if(snhei.le.DELTSN+SNTH) then
	2943	!-- 1-layer snow model
	2944	ilnb=1
	2945	snprim=snhei
	2946	soilt1=tso(1)
	2947	tsob=tso(1)
	2948	XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
	2949	DDZSN = XSN / SNPRIM
	2950	X1SN = DDZSN * thdifsn
	2951	X2 = DTDZS(1)*THDIF(1)
	2952	FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
	2953	-X2*(TSO(1)-TSO(2))
	2954	DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
	2955	cotso(NZS)=X1SN/DENOM
	2956	rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
	2957	cotsn=cotso(NZS)
	2958	rhtsn=rhtso(NZS)
	2959	!*** Average temperature of snow pack (C)
	2960	tsnav=0.5*(soilt+tso(1)) &
	2961	-273.15
	2962
	2963	else
	2964	!-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
	2965	ilnb=2
	2966	snprim=deltsn
	2967	tsob=soilt1
	2968	XSN = DELT/2./(0.5*SNHEI)
	2969	XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
	2970	DDZSN = XSN / DELTSN
	2971	DDZSN1 = XSN1 / (SNHEI-DELTSN)
	2972	X1SN = DDZSN * thdifsn
	2973	X1SN1 = DDZSN1 * thdifsn
	2974	X2 = DTDZS(1)*THDIF(1)
	2975	FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
	2976	-X2*(TSO(1)-TSO(2))
	2977	DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
	2978	cotso(nzs)=x1sn1/denom
	2979	rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
	2980	ftsnow = soilt1+x1sn*(soilt-soilt1) &
	2981	-x1sn1*(soilt1-tso(1))
	2982	denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
	2983	cotsn=x1sn/denomsn
	2984	rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
	2985	!*** Average temperature of snow pack (C)
	2986	tsnav=0.5/snhei((soilt+soilt1)deltsn &
	2987	+(soilt1+tso(1))*(SNHEI-DELTSN)) &
	2988	-273.15
	2989	endif
	2990	ENDIF
	2991
	2992	IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
	2993	!--- snow is too thin to be treated separately, therefore it
	2994	!--- is combined with the first soil layer.
	2995	fsn=SNHEI/(SNHEI+zsmain(2))
	2996	fso=1.-fsn
	2997	soilt1=tso(1)
	2998	tsob=tso(2)
	2999	snprim=SNHEI+zsmain(2)
	3000	XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
	3001	DDZSN = XSN /snprim
	3002	X1SN = DDZSN * (fsnthdifsn+fsothdif(1))
	3003	X2=DTDZS(2)*THDIF(2)
	3004	FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
	3005	X2*(TSO(2)-TSO(3))
	3006	denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
	3007	cotso(nzs1) = x1sn/denom
	3008	rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
	3009	tsnav=0.5*(soilt+tso(1)) &
	3010	-273.15
	3011	ENDIF
	3012
	3013	!************************************************************************
	3014	!--- THE HEAT BALANCE EQUATION (Smirnova et al. 1996, EQ. 21,26)
	3015
	3016	ETT1=0.
	3017	EPOT=-QKMS*(QVATM-QSG)
	3018	RHCS=CAP(1)
	3019	H=MAVAIL
	3020	IF(DEW.NE.0.)THEN
	3021	DRYCAN=0.
	3022	WETCAN=1.
	3023	ENDIF
	3024	TRANS=PCTRANSUMDRYCAN/ZSHALF(NROOT+1)
	3025	CAN=WETCAN+TRANS
	3026	UMVEG=1.-VEGFRAC
	3027	FKT=TKMS
	3028	D1=cotso(NZS1)
	3029	D2=rhtso(NZS1)
	3030	TN=SOILT
	3031	D9=THDIF(1)RHCSdzstop
	3032	D10=TKMSCPRHO
	3033	R211=.5*CONFLX/DELT
	3034	R21=R211CPRHO
	3035	R22=.5/(THDIF(1)DELTdzstop**2)
	3036	R6=EMISS STBOLT.5TN*4
	3037	R7=R6/TN
	3038	D11=RNET+R6
	3039
	3040	IF(SNHEI.GE.SNTH) THEN
	3041	! if(snhei.le.DELTSN+DELTSN) then
	3042	if(snhei.le.DELTSN+SNTH) then
	3043	!--- 1-layer snow
	3044	D1SN = cotso(NZS)
	3045	D2SN = rhtso(NZS)
	3046	else
	3047	!--- 2-layer snow
	3048	D1SN = cotsn
	3049	D2SN = rhtsn
	3050	endif
	3051	D9SN= THDIFSN*RHOCSN / SNPRIM
	3052	R22SN = SNPRIMSNPRIM0.5/(THDIFSN*DELT)
	3053	ENDIF
	3054
	3055	IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
	3056	!--- thin snow is combined with soil
	3057	D1SN = D1
	3058	D2SN = D2
	3059	D9SN = (fsnTHDIFSNRHOCSN+fsoTHDIF(1)RHCS)/ &
	3060	snprim
	3061	R22SN = snprimsnprim0.5 &
	3062	/((fsnTHDIFSN+fsoTHDIF(1))*delt)
	3063	ENDIF
	3064
	3065	IF(SNHEI.eq.0.)then
	3066	!--- all snow is sublimated
	3067	D9SN = D9
	3068	R22SN = R22
	3069	D1SN = D1
	3070	D2SN = D2
	3071	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	3072	print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
	3073	ENDIF
	3074	ENDIF
	3075
	3076	!---- TDENOM for snow
	3077
	3078	TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
	3079	+RAINFCVWPRCPMS &
	3080	+RHOCSN*NEWSNOW/DELT
	3081
	3082	FKQ=QKMS*RHO
	3083	R210=R211*RHO
	3084	C=VEGFRACFKQCAN
	3085	CC=C*XLVM/TDENOM
	3086	AA=XLVM(BETAFKQ*UMVEG+R210)/TDENOM
	3087	BB=(D10TABS+R21TN+XLVM(QVATM &
	3088	(BETAFKQUMVEG+C) &
	3089	+R210QVG)+D11+D9SN(D2SN+R22SN*TN) &
	3090	+RAINFCVWPRCPMS*max(273.15,TABS) &
	3091	+ RHOCSNNEWSNOW/DELTmin(273.15,TABS) &
	3092	)/TDENOM
	3093	AA1=AA+CC
	3094	PP=PATM*1.E3
	3095	AA1=AA1/PP
	3096	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	3097	print *,'VILKA-SNOW'
	3098	print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
	3099	tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
	3100	ENDIF
	3101
	3102	CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
	3103	TQ2=QVATM+QCATM
	3104	TX2=TQ2*(1.-H)
	3105	Q1=TX2+H*QS1
	3106	!--- it is saturation over snow
	3107	90 QVG=QS1
	3108	QSG=QS1
	3109	QCG=Q1-QS1
	3110
	3111	!--- SOILT - skin temperature
	3112	SOILT=TS1
	3113
	3114	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	3115	print *,' AFTER VILKA-SNOW'
	3116	print *,' TS1,QS1: ', ts1,qs1
	3117	ENDIF
	3118
	3119	! Solution for temperature at 7.5 cm depth and snow-soil interface
	3120	IF(SNHEI.GE.SNTH) THEN
	3121	! if(snhei.gt.DELTSN+DELTSN) then
	3122	if(snhei.gt.DELTSN+SNTH) then
	3123	!-- 2-layer snow model
	3124	SOILT1=rhtsn+cotsn*SOILT
	3125	TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT1
	3126	tsob=soilt1
	3127	else
	3128	!-- 1 layer in snow
	3129	TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT
	3130	SOILT1=TSO(1)
	3131	tsob=tso(1)
	3132	endif
	3133	ELSE
	3134	!-- all snow is evaporated
	3135	TSO(1)=SOILT
	3136	SOILT1=SOILT
	3137	tsob=SOILT
	3138	ENDIF
	3139
	3140	!---- Final solution for TSO
	3141	DO K=2,NZS
	3142	KK=NZS-K+1
	3143	TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
	3144	END DO
	3145	!--- For thin snow layer combined with the top soil layer
	3146	!--- TSO is computed by linear inmterpolation between SOILT
	3147	!--- and TSO(2)
	3148
	3149	if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then
	3150	tso(1)=tso(2)+(soilt-tso(2))*fso
	3151	SOILT1=TSO(1)
	3152	tsob=tso(2)
	3153	!!! tsob=tso(1)
	3154	endif
	3155
	3156	!--- IF SOILT > 273.15 F then melting of snow can happen
	3157	IF(SOILT.GE.273.15.AND.SNHEI.GT.0.) THEN
	3158	!!! SOILT=273.15
	3159	soiltfrac=snowfrac273.15+(1.-snowfrac)SOILT
	3160	soilt=soiltfrac
	3161	QSG= QSN(soilt,TBQ)/PP
	3162	!!! QSG= QSN(273.15,TBQ)/PP
	3163	QVG=QSG
	3164	T3 = STBOLTSOILTSOILT*SOILT
	3165	UPFLUX = T3 * SOILT
	3166	XINET = EMISS*(GLW-UPFLUX)
	3167	RNET = GSW + XINET
	3168	EPOT = -QKMS*(QVATM-QSG)
	3169	Q1=EPOT*RAS
	3170
	3171	IF (Q1.LE.0.) THEN
	3172	! --- condensation
	3173	DEW=-EPOT
	3174	DO K=1,NZS
	3175	TRANSP(K)=0.
	3176	ENDDO
	3177
	3178	QFX= XLVMRHODEW
	3179	EETA=QFX/XLVM
	3180	ELSE
	3181	! --- evaporation
	3182	DO K=1,NROOT
	3183	TRANSP(K)=-VEGFRAC*q1 &
	3184	PCTRANF(K)*DRYCAN/zshalf(NROOT+1)
	3185	IF(TRANSP(K).GT.0.) TRANSP(K)=0.
	3186	ETT1=ETT1-TRANSP(K)
	3187	ENDDO
	3188	DO k=nroot+1,nzs
	3189	transp(k)=0.
	3190	enddo
	3191
	3192	EDIR1 = Q1UMVEG BETA
	3193	EC1 = Q1 * WETCAN *VEGFRAC
	3194	CMC2MS=CST/DELT
	3195	EC1=MIN(CMC2MS,EC1)
	3196	EETA = (EDIR1 + EC1 + ETT1)*1.E3
	3197	! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
	3198	QFX= - XLVM * EETA
	3199	ENDIF
	3200
	3201	HFX=D10*(TABS-soilt)
	3202	!!! HFX=D10*(TABS-273.15)
	3203
	3204	IF(SNHEI.GE.SNTH)then
	3205	SOH=thdifsnRHOCSN(soilt-TSOB)/SNPRIM
	3206	! SOH=thdifsnRHOCSN(273.15-TSOB)/SNPRIM
	3207	SNFLX=SOH
	3208	ELSE
	3209	SOH=(fsnthdifsnrhocsn+fsothdif(1)rhcs)* &
	3210	(soilt-TSOB)/snprim
	3211	!!! (273.15-TSOB)/snprim
	3212	SNFLX=SOH
	3213	ENDIF
	3214
	3215	X= (R21+D9SNR22SN)(soilt-TNOLD) + &
	3216	!!! X= (R21+D9SNR22SN)(273.15-TNOLD) + &
	3217	XLVMR210(QSG-QGOLD)
	3218	!-- SNOH is energy flux of snow phase change
	3219	SNOH=RNET+QFX +HFX &
	3220	+RHOCSNNEWSNOW/DELT(min(273.15,TABS)-TN) &
	3221	-SOH-X+RAINFCVWPRCPMS* &
	3222	(max(273.15,TABS)-TN)
	3223	SNOH=AMAX1(0.,SNOH)
	3224	!-- SMELT is speed of melting in M/S
	3225	SMELT= SNOH /XLMELT*1.E-3
	3226	SMELT=AMIN1(SMELT,SNWEPR/DELT-BETAEPOTRAS)
	3227	! SMELT=AMIN1(SMELT,SNWEPR/DELT-BETAEPOTRAS*UMVEG)
	3228
	3229	SNOHGNEW=SMELTXLMELT1.E3
	3230	SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
	3231
	3232	SNOH=SNOHGNEW
	3233	! SNOHSMELTXLMELT1.E3
	3234
	3235	!*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
	3236	!!! rsm=0.13smeltdelt
	3237	if(snwe.gt.0.) then
	3238	rsmfrac=min(0.18,(max(0.08,0.10/snwe*0.13)))
	3239	else
	3240	rsmfrac=0.13
	3241	endif
	3242
	3243	rsm=rsmfracsmeltdelt
	3244	SMELT=SMELT-rsm/delt
	3245
	3246	!-- correction of liquid equivalent of snow depth
	3247	!-- due to evaporation and snow melt
	3248	SNWE = AMAX1(0.,(SNWEPR- &
	3249	(SMELT+BETAEPOTRAS)*DELT &
	3250	! (SMELT+BETAEPOTRASUMVEG)DELT &
	3251	) )
	3252
	3253	!--- If all snow melts, then 13% of snow melt we kept in the
	3254	!--- snow pack should be added back to snow melt and infiltrate
	3255	!--- into soil.
	3256	if(snwe.le.rsm) then
	3257	smelt=smelt+rsm/delt
	3258	snwe=0.
	3259	rsm=0.
	3260	SOILT=SNODIFDELT/RHCSZSHALF(2) &
	3261	+soiltfrac
	3262	!!! +273.15
	3263	else
	3264	!*** Correct snow density on effect of snow melt, melted
	3265	!*** from the top of the snow. 13% of melted water
	3266	!*** remains in the pack and changes its density.
	3267	!*** Eq. 9 (with my correction) in Koren et al. (1999)
	3268
	3269	if(snwe.gt.snthrhosn1.e-3) then
	3270	xsn=(rhosn(snwe-rsm)+1.e3rsm)/ &
	3271	snwe
	3272	rhosn=MIN(XSN,400.)
	3273
	3274	RHOCSN=2090.* RHOSN
	3275	thdifsn = 0.265/RHOCSN
	3276	endif
	3277
	3278	SOILT=SOILTFRAC
	3279
	3280	endif
	3281
	3282	!--- If there is no snow melting then just evaporation
	3283	!--- or condensation cxhanges SNWE
	3284	ELSE
	3285	EPOT=-QKMS*(QVATM-QSG)
	3286	SNWE = AMAX1(0.,(SNWEPR- &
	3287	BETAEPOTRAS*DELT))
	3288	! BETAEPOTRASUMVEGDELT))
	3289
	3290	ENDIF
	3291	!*** Correct snow density on effect of snow melt, melted
	3292	!*** from the top of the snow. 13% of melted water
	3293	!*** remains in the pack and changes its density.
	3294	!*** Eq. 9 (with my correction) in Koren et al. (1999)
	3295
	3296	SNHEI=SNWE *1.E3 / RHOSN
	3297
	3298	!-- Snow melt from the top is done. But if ground surface temperature
	3299	!-- is above freezing snow can melt from the bottom. The following
	3300	!-- piece of code will check if bottom melting is possible.
	3301
	3302	IF(TSO(1).GE.273.15.AND.SNHEI.GT.0.) THEN
	3303	soiltfrac=snowfrac273.15+(1.-snowfrac)TSO(1)
	3304
	3305	SNOHG=(TSO(1)-soiltfrac)(RHCSzshalf(2)+ &
	3306	RHOCSN0.5SNHEI) / DELT
	3307	SNOHG=AMAX1(0.,SNOHG)
	3308	SNODIF=0.
	3309	! TSO(1)=273.15
	3310	SMELTG=SNOHG/XLMELT*1.E-3
	3311	! SMELTG=AMIN1(SMELTG,SNWE/DELT)
	3312	if(SNWE-SMELTG*DELT.ge.rsm) then
	3313	! SNWE = SNWE-SMELTG*DELT
	3314	SNWE = AMAX1(0.,SNWE-SMELTG*DELT)
	3315	else
	3316	smeltg=snwe/delt
	3317	snwe=0.
	3318	rsm=0.
	3319	endif
	3320
	3321	SNOHGNEW=SMELTGXLMELT1.e3
	3322	SNODIF=AMAX1(0.,(SNOHG-SNOHGNEW))
	3323	TSO(1)=soiltfrac
	3324	if(snwe.eq.0.)then
	3325	TSO(1)=SNODIFDELT/RHCSzshalf(2) + soiltfrac
	3326	!!! TSO(1)=SNODIFDELT/RHCSzshalf(2) + 273.15
	3327	endif
	3328
	3329	SMELT=SMELT+SMELTG
	3330	SNOH=SNOH+SNOHGNEW
	3331
	3332	ENDIF
	3333
	3334	SNHEI=SNWE *1.E3 / RHOSN
	3335
	3336	snweprint=snwe
	3337	! &
	3338	!--- if VEGFRAC.ne.0. then some snow stays on the canopy
	3339	!--- and should be added to SNWE for water conservation
	3340	! 4 Nov 07 +VEGFRAC*cst
	3341	snheiprint=snweprint*1.E3 / RHOSN
	3342
	3343	IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
	3344	print *, 'snweprint : ',snweprint
	3345	print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
	3346	ENDIF
	3347	!--- Compute flux in the top snow layer
	3348	SNFLX=D9SN*(SOILT-TSOB)
	3349	IF(SNHEI.GT.0.) THEN
	3350	if(ilnb.gt.1) then
	3351	tsnav=0.5/snhei((soilt+soilt1)deltsn &
	3352	+(soilt1+tso(1))*(SNHEI-DELTSN)) &
	3353	-273.15
	3354	else
	3355	tsnav=0.5*(soilt+tso(1)) - 273.15
	3356	endif
	3357	ENDIF
	3358
	3359	! return
	3360	! end
	3361	!------------------------------------------------------------------------
	3362	END SUBROUTINE SNOWTEMP
	3363	!------------------------------------------------------------------------
	3364
	3365
	3366	SUBROUTINE SOILMOIST ( &
	3367	!--input parameters
	3368	DELT,NZS,NDDZS,DTDZS,DTDZS2, &
	3369	ZSMAIN,ZSHALF,DIFFU,HYDRO, &
	3370	QSG,QVG,QCG,QCATM,QVATM,PRCP, &
	3371	QKMS,TRANSP,DRIP, &
	3372	DEW,SMELT,SOILICE,VEGFRAC, &
	3373	!--soil properties
	3374	DQM,QMIN,REF,KSAT,RAS,INFMAX, &
	3375	!--output
	3376	SOILMOIS,MAVAIL,RUNOFF,RUNOFF2,INFILTRP)
	3377	!*************************************************************************
	3378	! moisture balance equation and Richards eqn.
	3379	! are solved here
	3380	!
	3381	! DELT - time step (s)
	3382	! IME,JME,NZS - dimensions of soil domain
	3383	! ZSMAIN - main levels in soil (m)
	3384	! ZSHALF - middle of the soil layers (m)
	3385	! DTDZS - dt/(2.dzshalfdzmain)
	3386	! DTDZS2 - dt/(2.*dzshalf)
	3387	! DIFFU - diffusional conductivity (m^2/s)
	3388	! HYDRO - hydraulic conductivity (m/s)
	3389	! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
	3390	! water vapor and cloud at the ground
	3391	! surface, respectively (kg/kg)
	3392	! QCATM,QVATM - cloud and water vapor mixing ratio
	3393	! at the first atm. level (kg/kg)
	3394	! PRCP - precipitation rate in m/s
	3395	! QKMS - exchange coefficient for water vapor in the
	3396	! surface layer (m/s)
	3397	! TRANSP - transpiration from the soil layers (m/s)
	3398	! DRIP - liquid water dripping from the canopy to soil (m)
	3399	! DEW - dew in kg/m^2s
	3400	! SMELT - melting rate in m/s
	3401	! SOILICE - volumetric content of ice in soil (m^3/m^3)
	3402	! VEGFRAC - greeness fraction (0-1)
	3403	! RAS - ration of air density to soil density
	3404	! INFMAX - maximum infiltration rate (kg/m^2/s)
	3405	!
	3406	! SOILMOIS - volumetric soil moisture, 6 levels (m^3/m^3)
	3407	! MAVAIL - fraction of maximum soil moisture in the top
	3408	! layer (0-1)
	3409	! RUNOFF - surface runoff (m/s)
	3410	! RUNOFF2 - underground runoff (m)
	3411	! INFILTRP - point infiltration flux into soil (m/s)
	3412	! /(snow bottom runoff) (mm/s)
	3413	!
	3414	! COSMC, RHSMC - coefficients for implicit solution of
	3415	! Richards equation
	3416	!******************************************************************
	3417	IMPLICIT NONE
	3418	!------------------------------------------------------------------
	3419	!--- input variables
	3420	REAL, INTENT(IN ) :: DELT
	3421	INTEGER, INTENT(IN ) :: NZS,NDDZS
	3422
	3423	! input variables
	3424
	3425	REAL, DIMENSION(1:NZS), INTENT(IN ) :: ZSMAIN, &
	3426	ZSHALF, &
	3427	DIFFU, &
	3428	HYDRO, &
	3429	TRANSP, &
	3430	SOILICE, &
	3431	DTDZS2
	3432
	3433	REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
	3434
	3435	REAL, INTENT(IN ) :: QSG,QVG,QCG,QCATM,QVATM , &
	3436	QKMS,VEGFRAC,DRIP,PRCP , &
	3437	DEW,SMELT , &
	3438	DQM,QMIN,REF,KSAT,RAS
	3439
	3440	! output
	3441
	3442	REAL, DIMENSION( 1:nzs ) , &
	3443
	3444	INTENT(INOUT) :: SOILMOIS
	3445
	3446	REAL, INTENT(INOUT) :: MAVAIL,RUNOFF,RUNOFF2,INFILTRP, &
	3447	INFMAX
	3448
	3449	! local variables
	3450
	3451	REAL, DIMENSION( 1:nzs ) :: COSMC,RHSMC
	3452
	3453	REAL :: DZS,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10
	3454	REAL :: REFKDT,REFDK,DELT1,F1MAX,F2MAX
	3455	REAL :: F1,F2,FD,KDT,VAL,DDT,PX
	3456	REAL :: QQ,UMVEG,INFMAX1,TRANS
	3457	REAL :: TOTLIQ,FLX,FLXSAT,QTOT
	3458	REAL :: DID,X1,X2,X4,DENOM,Q2,Q4
	3459	REAL :: dice,fcr,acrt,frzx,sum,cvfrz
	3460
	3461	INTEGER :: NZS1,NZS2,K,KK,K1,KN,ialp1,jj,jk
	3462
	3463	!******************************************************************************
	3464	! COEFFICIENTS FOR THOMAS ALGORITHM FOR SOILMOIS
	3465	!******************************************************************************
	3466	NZS1=NZS-1
	3467	NZS2=NZS-2
	3468
	3469	118 format(6(10Pf23.19))
	3470
	3471	do k=1,nzs
	3472	cosmc(k)=0.
	3473	rhsmc(k)=0.
	3474	enddo
	3475
	3476	DID=(ZSMAIN(NZS)-ZSHALF(NZS))*2.
	3477	X1=ZSMAIN(NZS)-ZSMAIN(NZS1)
	3478	! DENOM=DID/DELT+DIFFU(NZS1)/X1
	3479	! COSMC(1)=DIFFU(NZS1)/X1/DENOM
	3480	! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT
	3481	! 1 +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS)
	3482	! 1 -HYDRO(NZS1)SOILMOIS(NZS1))DID
	3483	! 1 /X1) /DENOM
	3484
	3485	DENOM=(1.+DIFFU(nzs1)/X1/DIDDELT+HYDRO(NZS)/(2.DID)*DELT)
	3486	COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 &
	3487	+HYDRO(NZS1)/2./DID)/DENOM
	3488	RHSMC(1)=(SOILMOIS(NZS)+TRANSP(NZS)*DELT/ &
	3489	DID)/DENOM
	3490
	3491	DO 330 K=1,NZS2
	3492	KN=NZS-K
	3493	K1=2*KN-3
	3494	X4=2.DTDZS(K1)DIFFU(KN-1)
	3495	X2=2.DTDZS(K1+1)DIFFU(KN)
	3496	Q4=X4+HYDRO(KN-1)*DTDZS2(KN-1)
	3497	Q2=X2-HYDRO(KN+1)*DTDZS2(KN-1)
	3498	DENOM=1.+X2+X4-Q2*COSMC(K)
	3499	COSMC(K+1)=Q4/DENOM
	3500	330 RHSMC(K+1)=(SOILMOIS(KN)+Q2*RHSMC(K) &
	3501	+TRANSP(KN) &
	3502	/(ZSHALF(KN+1)-ZSHALF(KN)) &
	3503	*DELT)/DENOM
	3504
	3505	! --- MOISTURE BALANCE BEGINS HERE
	3506
	3507	TRANS=TRANSP(1)
	3508	UMVEG=1.-VEGFRAC
	3509
	3510	RUNOFF=0.
	3511	RUNOFF2=0.
	3512	DZS=ZSMAIN(2)
	3513	R1=COSMC(NZS1)
	3514	R2= RHSMC(NZS1)
	3515	R3=DIFFU(1)/DZS
	3516	R4=R3+HYDRO(1)*.5
	3517	R5=R3-HYDRO(2)*.5
	3518	R6=QKMS*RAS
	3519	!-- Total liquid water available on the top of soil domain
	3520	!-- Without snow - 3 sources of water: precipitation,
	3521	!-- water dripping from the canopy and dew
	3522	!-- With snow - only one source of water - snow melt
	3523
	3524	! print *,'PRCP,DRIP,DEW,umveg,ras,smelt',
	3525	! 1 PRCP,DRIP,DEW,umveg,ras,smelt
	3526	! if (drip.ne.0.) then
	3527	! print *,'DRIP non-zero'
	3528	! write(6,191) drip
	3529	! write (6,191)soilmois(1)
	3530	! write (6,191)soilmois(2)
	3531	! endif
	3532	191 format (f23.19)
	3533
	3534	TOTLIQ=UMVEGPRCP-DRIP/DELT-UMVEGDEW*RAS-SMELT
	3535
	3536
	3537	FLX=TOTLIQ
	3538	INFILTRP=TOTLIQ
	3539
	3540	! ----------- FROZEN GROUND VERSION -------------------------
	3541	! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
	3542	! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
	3543	! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT.
	3544	! BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT
	3545	! CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM.
	3546	! THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6})
	3547	!
	3548	! Current logic doesn't allow CVFRZ be bigger than 3
	3549	CVFRZ = 3.
	3550
	3551	!-- SCHAAKE/KOREN EXPRESSION for calculation of max infiltration
	3552	REFKDT=3.
	3553	REFDK=3.4341E-6
	3554	DELT1=DELT/86400.
	3555	F1MAX=DQM*ZSHALF(2)
	3556	F2MAX=DQM*(ZSHALF(3)-ZSHALF(2))
	3557	F1=F1MAX*(1.-SOILMOIS(1)/DQM)
	3558	F2=F2MAX*(1.-SOILMOIS(2)/DQM)
	3559	FD=F1+F2
	3560	KDT=REFKDT*KSAT/REFDK
	3561	VAL=(1.-EXP(-KDT*DELT1))
	3562	DDT = FD*VAL
	3563	PX= - TOTLIQ * DELT
	3564	IF(PX.LT.0.0) PX = 0.0
	3565	if(ddt.eq.0.) then
	3566	infmax1=ksat
	3567	else
	3568	INFMAX1 = (PX*(DDT/(PX+DDT)))/DELT
	3569	INFMAX1 = MIN(INFMAX1, KSAT)
	3570	endif
	3571	! print *,'INFMAX1=,ksat',infmax1,ksat,f1,f2,kdt,val,ddt,px
	3572	! ----------- FROZEN GROUND VERSION --------------------------
	3573	! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
	3574	!
	3575	! ------------------------------------------------------------------
	3576
	3577	DICE = soilice(1)*zshalf(2)
	3578	DO K=2,NZS1
	3579	DICE = DICE + ( ZSHALF(K+1) - ZSHALF(K) ) * soilice(k)
	3580	ENDDO
	3581	FRZX= 0.28((dqm+qmin)/ref) (0.400 / 0.482)
	3582	FCR = 1.
	3583	IF ( DICE .GT. 1.E-2) THEN
	3584	ACRT = CVFRZ * FRZX / DICE
	3585	SUM = 1.
	3586	IALP1 = CVFRZ - 1
	3587	DO JK = 1,IALP1
	3588	K = 1
	3589	DO JJ = JK+1, IALP1
	3590	K = K * JJ
	3591	END DO
	3592	SUM = SUM + (ACRT ** ( CVFRZ-JK)) / FLOAT (K)
	3593	END DO
	3594	FCR = 1. - EXP(-ACRT) * SUM
	3595	END IF
	3596	! print *,'FCR--------',fcr
	3597	INFMAX1 = INFMAX1* FCR
	3598	INFMAX1 = MIN(INFMAX1, KSAT)
	3599	! -------------------------------------------------------------------
	3600
	3601	INFMAX = MIN(INFMAX,INFMAX1)
	3602	!----
	3603	IF (-TOTLIQ.GE.INFMAX)THEN
	3604	RUNOFF=-TOTLIQ-INFMAX
	3605	FLX=-INFMAX
	3606	ENDIF
	3607	! INFILTRP is total infiltration flux in M/S
	3608	INFILTRP=FLX
	3609	! print *,'PRCIP',infiltrp,flx,infmax
	3610	! Solution of moisture budget
	3611	R7=.5*DZS/DELT
	3612	R4=R4+R7
	3613	FLX=FLX-SOILMOIS(1)*R7
	3614	R8=UMVEG*R6
	3615	QTOT=QVATM+QCATM
	3616	R9=TRANS
	3617	R10=QTOT-QSG
	3618	!-- evaporation regime
	3619	IF(R10.LE.0.) THEN
	3620	QQ=(R5R2-FLX+R9)/(R4-R5R1-R10*R8/(REF-QMIN))
	3621	FLXSAT=-DQM(R4-R5R1-R10*R8/(REF-QMIN)) &
	3622	+R5*R2+R9
	3623	ELSE
	3624	!-- dew formation regime
	3625	QQ=(R2R5-FLX+R8(QTOT-QCG-QVG)+R9)/(R4-R1*R5)
	3626	FLXSAT=-DQM(R4-R1R5)+R2R5+R8(QTOT-QVG-QCG)+R9
	3627	END IF
	3628
	3629	IF(QQ.LT.0.) THEN
	3630	SOILMOIS(1)=0.
	3631
	3632	ELSE IF(QQ.GT.DQM) THEN
	3633	!-- saturation
	3634	SOILMOIS(1)=DQM
	3635	RUNOFF2=runoff2+(FLXSAT-FLX)*DELT
	3636	RUNOFF=RUNOFF+(FLXSAT-FLX)
	3637	ELSE
	3638	SOILMOIS(1)=max(1.e-8,QQ)
	3639	END IF
	3640
	3641	!--- FINAL SOLUTION FOR SOILMOIS
	3642	! DO K=2,NZS
	3643	DO K=2,NZS-1
	3644	KK=NZS-K+1
	3645	QQ=COSMC(KK)*SOILMOIS(K-1)+RHSMC(KK)
	3646
	3647	IF (QQ.LT.0.) THEN
	3648	SOILMOIS(K)=0.
	3649
	3650	ELSE IF(QQ.GT.DQM) THEN
	3651	!-- saturation
	3652	SOILMOIS(K)=DQM
	3653	IF(K.EQ.NZS)THEN
	3654	RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSMAIN(K)-ZSHALF(K))
	3655	ELSE
	3656	RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSMAIN(K+1)-ZSHALF(K))
	3657	ENDIF
	3658	ELSE
	3659	SOILMOIS(K)=max(1.e-8,QQ)
	3660	END IF
	3661	END DO
	3662
	3663	! MAVAIL=min(1.,SOILMOIS(1)/(REF-QMIN))
	3664	MAVAIL=min(1.,SOILMOIS(1)/DQM)
	3665	if (MAVAIL.EQ.0.) MAVAIL=.00001
	3666
	3667	! RETURN
	3668	! END
	3669	!-------------------------------------------------------------------
	3670	END SUBROUTINE SOILMOIST
	3671	!-------------------------------------------------------------------
	3672
	3673
	3674	SUBROUTINE SOILPROP( &
	3675	!--- input variables
	3676	nzs,fwsat,lwsat,tav,keepfr, &
	3677	soilmois,soiliqw,soilice, &
	3678	soilmoism,soiliqwm,soilicem, &
	3679	!--- soil fixed fields
	3680	QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
	3681	!--- constants
	3682	riw,xlmelt,CP,G0_P,cvw,ci, &
	3683	kqwrtz,kice,kwt, &
	3684	!--- output variables
	3685	thdif,diffu,hydro,cap)
	3686
	3687	!******************************************************************
	3688	! SOILPROP computes thermal diffusivity, and diffusional and
	3689	! hydraulic condeuctivities
	3690	!******************************************************************
	3691	! NX,NY,NZS - dimensions of soil domain
	3692	! FWSAT, LWSAT - volumetric content of frozen and liquid water
	3693	! for saturated condition at given temperatures (m^3/m^3)
	3694	! TAV - temperature averaged for soil layers (K)
	3695	! SOILMOIS - volumetric soil moisture at the main soil levels (m^3/m^3)
	3696	! SOILMOISM - volumetric soil moisture averaged for layers (m^3/m^3)
	3697	! SOILIQWM - volumetric liquid soil moisture averaged for layers (m^3/m^3)
	3698	! SOILICEM - volumetric content of soil ice averaged for layers (m^3/m^3)
	3699	! THDIF - thermal diffusivity for soil layers (W/m/K)
	3700	! DIFFU - diffusional conductivity (m^2/s)
	3701	! HYDRO - hydraulic conductivity (m/s)
	3702	! CAP - volumetric heat capacity (J/m^3/K)
	3703	!
	3704	!******************************************************************
	3705
	3706	IMPLICIT NONE
	3707	!-----------------------------------------------------------------
	3708
	3709	!--- soil properties
	3710	INTEGER, INTENT(IN ) :: NZS
	3711	REAL , &
	3712	INTENT(IN ) :: RHOCS, &
	3713	BCLH, &
	3714	DQM, &
	3715	KSAT, &
	3716	PSIS, &
	3717	QWRTZ, &
	3718	QMIN
	3719
	3720	REAL, DIMENSION( 1:nzs ) , &
	3721	INTENT(IN ) :: SOILMOIS, &
	3722	keepfr
	3723
	3724
	3725	REAL, INTENT(IN ) :: CP, &
	3726	CVW, &
	3727	RIW, &
	3728	kqwrtz, &
	3729	kice, &
	3730	kwt, &
	3731	XLMELT, &
	3732	G0_P
	3733
	3734
	3735
	3736	!--- output variables
	3737	REAL, DIMENSION(1:NZS) , &
	3738	INTENT(INOUT) :: cap,diffu,hydro , &
	3739	thdif,tav , &
	3740	soilmoism , &
	3741	soiliqw,soilice , &
	3742	soilicem,soiliqwm , &
	3743	fwsat,lwsat
	3744
	3745	!--- local variables
	3746	REAL, DIMENSION(1:NZS) :: hk,detal,kasat,kjpl
	3747
	3748	REAL :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci
	3749	REAL :: tln,tavln,tn,pf,a,am,ame,h
	3750	INTEGER :: nzs1,k
	3751
	3752	!-- for Johansen thermal conductivity
	3753	REAL :: kzero,gamd,kdry,kas,x5,sr,ke
	3754
	3755
	3756	nzs1=nzs-1
	3757
	3758	!-- Constants for Johansen (1975) thermal conductivity
	3759	kzero =2. ! if qwrtz > 0.2
	3760
	3761
	3762	do k=1,nzs
	3763	detal (k)=0.
	3764	kasat (k)=0.
	3765	kjpl (k)=0.
	3766	hk (k)=0.
	3767	enddo
	3768
	3769	ws=dqm+qmin
	3770	x1=xlmelt/(g0_p*psis)
	3771	x2=x1/bclh*ws
	3772	x4=(bclh+1.)/bclh
	3773	!--- Next 3 lines are for Johansen thermal conduct.
	3774	gamd=(1.-ws)*2700.
	3775	kdry=(0.135gamd+64.7)/(2700.-0.947gamd)
	3776	kas=kqwrtz*qwrtzkzero**(1.-qwrtz)
	3777
	3778	DO K=1,NZS1
	3779	tn=tav(k) - 273.15
	3780	wd=ws - riw*soilicem(k)
	3781	psif=psis100.(wd/(soiliqwm(k)+qmin))**bclh &
	3782	* (ws/wd)**3.
	3783	!--- PSIF should be in [CM] to compute PF
	3784	pf=log10(abs(psif))
	3785	fact=1.+riw*soilicem(k)
	3786	!--- HK is for McCumber thermal conductivity
	3787	IF(PF.LE.5.2) THEN
	3788	HK(K)=420.EXP(-(PF+2.7))fact
	3789	ELSE
	3790	HK(K)=.1744*fact
	3791	END IF
	3792
	3793	IF(soilicem(k).NE.0.AND.TN.LT.0.) then
	3794	!--- DETAL is taking care of energy spent on freezing or released from
	3795	! melting of soil water
	3796
	3797	DETAL(K)=273.15X2/(TAV(K)TAV(K))* &
	3798	(TAV(K)/(X1TN))*X4
	3799
	3800	if(keepfr(k).eq.1.) then
	3801	detal(k)=0.
	3802	endif
	3803
	3804	ENDIF
	3805
	3806	!--- Next 10 lines calculate Johansen thermal conductivity KJPL
	3807	kasat(k)=kas*(1.-ws)kice**fwsat(k) &
	3808	kwt*lwsat(k)
	3809
	3810	X5=(soilmoism(k)+qmin)/ws
	3811	if(soilicem(k).eq.0.) then
	3812	sr=max(0.101,x5)
	3813	ke=log10(sr)+1.
	3814	!--- next 2 lines - for coarse soils
	3815	! sr=max(0.0501,x5)
	3816	! ke=0.7*log10(sr)+1.
	3817	else
	3818	ke=x5
	3819	endif
	3820
	3821	kjpl(k)=ke*(kasat(k)-kdry)+kdry
	3822
	3823	!--- CAP -volumetric heat capacity
	3824	CAP(K)=(1.-WS)*RHOCS &
	3825	+ (soiliqwm(K)+qmin)*CVW &
	3826	+ soilicem(K)*CI &
	3827	+ (dqm-soilmoism(k))CP1.2 &
	3828	- DETAL(K)1.e3xlmelt
	3829
	3830	a=RIW*soilicem(K)
	3831
	3832	if((ws-a).lt.0.12)then
	3833	diffu(K)=0.
	3834	else
	3835	H=max(0.,(soilmoism(K)-a)/(max(1.e-8,(dqm-a))))
	3836	facd=1.
	3837	if(a.ne.0.)facd=1.-a/max(1.e-8,soilmoism(K))
	3838	ame=max(1.e-8,dqm-riw*soilicem(K))
	3839	!--- DIFFU is diffusional conductivity of soil water
	3840	diffu(K)=-BCLHKSATPSIS/ame* &
	3841	(dqm/ame)**3. &
	3842	H(BCLH+2.)facd
	3843	endif
	3844
	3845	! diffu(K)=-BCLHKSATPSIS/dqm &
	3846	! H*(BCLH+2.)
	3847
	3848
	3849	!--- thdif - thermal diffusivity
	3850	! thdif(K)=HK(K)/CAP(K)
	3851	!--- Use thermal conductivity from Johansen (1975)
	3852	thdif(K)=KJPL(K)/CAP(K)
	3853
	3854	END DO
	3855
	3856	DO K=1,NZS
	3857
	3858	if((ws-riw*soilice(k)).lt.0.12)then
	3859	hydro(k)=0.
	3860	else
	3861	fach=1.
	3862	if(soilice(k).ne.0.) &
	3863	fach=1.-riw*soilice(k)/max(1.e-8,soilmois(k))
	3864	am=max(1.e-8,dqm-riw*soilice(k))
	3865	!--- HYDRO is hydraulic conductivity of soil water
	3866	hydro(K)=KSAT/am* &
	3867	(soiliqw(K)/am) &
	3868	*(2.BCLH+2.) &
	3869	* fach
	3870	endif
	3871
	3872	ENDDO
	3873
	3874	! RETURN
	3875	! END
	3876
	3877	!-----------------------------------------------------------------------
	3878	END SUBROUTINE SOILPROP
	3879	!-----------------------------------------------------------------------
	3880
	3881
	3882	SUBROUTINE TRANSF( &
	3883	!--- input variables
	3884	nzs,nroot,soiliqw,tabs, &
	3885	!--- soil fixed fields
	3886	dqm,qmin,ref,wilt,zshalf, &
	3887	!--- output variables
	3888	tranf,transum)
	3889
	3890	!-------------------------------------------------------------------
	3891	!--- TRANF(K) - THE TRANSPIRATION FUNCTION (Smirnova et al. 1996, EQ. 18,19)
	3892	!*******************************************************************
	3893	! NX,NY,NZS - dimensions of soil domain
	3894	! SOILIQW - volumetric liquid soil moisture at the main levels (m^3/m^3)
	3895	! TRANF - the transpiration function at levels (m)
	3896	! TRANSUM - transpiration function integrated over the rooting zone (m)
	3897	!
	3898	!*******************************************************************
	3899	IMPLICIT NONE
	3900	!-------------------------------------------------------------------
	3901
	3902	!--- input variables
	3903
	3904	INTEGER, INTENT(IN ) :: nroot,nzs
	3905
	3906	REAL , &
	3907	INTENT(IN ) :: TABS
	3908	!--- soil properties
	3909	REAL , &
	3910	INTENT(IN ) :: DQM, &
	3911	QMIN, &
	3912	REF, &
	3913	WILT
	3914
	3915	REAL, DIMENSION(1:NZS), INTENT(IN) :: soiliqw, &
	3916	ZSHALF
	3917
	3918	!-- output
	3919	REAL, DIMENSION(1:NZS), INTENT(OUT) :: TRANF
	3920	REAL, INTENT(OUT) :: TRANSUM
	3921
	3922	!-- local variables
	3923	REAL :: totliq, did
	3924	INTEGER :: k
	3925
	3926	!-- for non-linear root distribution
	3927	REAL :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4
	3928	REAL :: FTEM
	3929	REAL, DIMENSION(1:NZS) :: PART
	3930	!--------------------------------------------------------------------
	3931
	3932	do k=1,nzs
	3933	part(k)=0.
	3934	enddo
	3935
	3936	transum=0.
	3937	totliq=soiliqw(1)+qmin
	3938	sm1=totliq
	3939	sm2=sm1*sm1
	3940	sm3=sm2*sm1
	3941	sm4=sm3*sm1
	3942	ap0=0.299
	3943	ap1=-8.152
	3944	ap2=61.653
	3945	ap3=-115.876
	3946	ap4=59.656
	3947	gx=ap0+ap1sm1+ap2sm2+ap3sm3+ap4sm4
	3948	if(totliq.ge.ref) gx=1.
	3949	if(totliq.le.0.) gx=0.
	3950	if(gx.gt.1.) gx=1.
	3951	if(gx.lt.0.) gx=0.
	3952	DID=zshalf(2)
	3953	part(1)=DID*gx
	3954	!--- air temperature function
	3955	! Avissar (1985) and AX 7/95
	3956	IF (TABS .LE. 302.15) THEN
	3957	FTEM = 1.0 / (1.0 + EXP(-0.41 * (TABS - 282.05)))
	3958	ELSE
	3959	FTEM = 1.0 / (1.0 + EXP(0.5 * (TABS - 314.0)))
	3960	ENDIF
	3961	!---
	3962	IF(TOTLIQ.GT.REF) THEN
	3963	TRANF(1)=DID
	3964	ELSE IF(TOTLIQ.LE.WILT) THEN
	3965	TRANF(1)=0.
	3966	ELSE
	3967	TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID
	3968	ENDIF
	3969	!-- uncomment next line for non-linear root distribution
	3970	!cc TRANF(1)=part(1)
	3971	TRANF(1)=TRANF(1)*FTEM
	3972
	3973	DO K=2,NROOT
	3974	totliq=soiliqw(k)+qmin
	3975	sm1=totliq
	3976	sm2=sm1*sm1
	3977	sm3=sm2*sm1
	3978	sm4=sm3*sm1
	3979	gx=ap0+ap1sm1+ap2sm2+ap3sm3+ap4sm4
	3980	if(totliq.ge.ref) gx=1.
	3981	if(totliq.le.0.) gx=0.
	3982	if(gx.gt.1.) gx=1.
	3983	if(gx.lt.0.) gx=0.
	3984	DID=zshalf(K+1)-zshalf(K)
	3985	part(k)=did*gx
	3986	IF(totliq.GE.REF) THEN
	3987	TRANF(K)=DID
	3988	ELSE IF(totliq.LE.WILT) THEN
	3989	TRANF(K)=0.
	3990	ELSE
	3991	TRANF(K)=(totliq-WILT) &
	3992	/(REF-WILT)*DID
	3993	ENDIF
	3994	!-- uncomment next line for non-linear root distribution
	3995	!cc TRANF(k)=part(k)
	3996	TRANF(k)=TRANF(k)*FTEM
	3997	END DO
	3998
	3999	!-- TRANSUM - total for the rooting zone
	4000	transum=0.
	4001	DO K=1,NROOT
	4002	transum=transum+tranf(k)
	4003	END DO
	4004
	4005	! RETURN
	4006	! END
	4007	!-----------------------------------------------------------------
	4008	END SUBROUTINE TRANSF
	4009	!-----------------------------------------------------------------
	4010
	4011
	4012	SUBROUTINE VILKA(TN,D1,D2,PP,QS,TS,TT,NSTEP,ii,j,iland,isoil)
	4013	!--------------------------------------------------------------
	4014	!--- VILKA finds the solution of energy budget at the surface
	4015	!--- using table T,QS computed from Clausius-Klapeiron
	4016	!--------------------------------------------------------------
	4017	REAL, DIMENSION(1:4001), INTENT(IN ) :: TT
	4018	REAL, INTENT(IN ) :: TN,D1,D2,PP
	4019	INTEGER, INTENT(IN ) :: NSTEP,ii,j,iland,isoil
	4020
	4021	REAL, INTENT(OUT ) :: QS, TS
	4022
	4023	REAL :: F1,T1,T2,RN
	4024	INTEGER :: I,I1
	4025
	4026	I=(TN-1.7315E2)/.05+1
	4027	T1=173.1+FLOAT(I)*.05
	4028	F1=T1+D1*TT(I)-D2
	4029	I1=I-F1/(.05+D1*(TT(I+1)-TT(I)))
	4030	I=I1
	4031	IF(I.GT.4000.OR.I.LT.1) GOTO 1
	4032	10 I1=I
	4033	T1=173.1+FLOAT(I)*.05
	4034	F1=T1+D1*TT(I)-D2
	4035	RN=F1/(.05+D1*(TT(I+1)-TT(I)))
	4036	I=I-INT(RN)
	4037	IF(I.GT.4000.OR.I.LT.1) GOTO 1
	4038	IF(I1.NE.I) GOTO 10
	4039	TS=T1-.05*RN
	4040	QS=(TT(I)+(TT(I)-TT(I+1))*RN)/PP
	4041	GOTO 20
	4042	1 PRINT *,' AVOST IN VILKA '
	4043	! WRITE(12,*)'AVOST',TN,D1,D2,PP,NSTEP
	4044	PRINT *,TN,D1,D2,PP,NSTEP,I,TT(i),ii,j,iland,isoil
	4045	CALL wrf_error_fatal (' AVOST IN VILKA ' )
	4046	20 CONTINUE
	4047	! RETURN
	4048	! END
	4049	!-----------------------------------------------------------------------
	4050	END SUBROUTINE VILKA
	4051	!-----------------------------------------------------------------------
	4052
	4053
	4054	SUBROUTINE SOILVEGIN ( IVGTYP,ISLTYP,MYJ, &
	4055	IFOREST,EMISS,PC,ZNT,QWRTZ, &
	4056	RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT )
	4057
	4058	!************************************************************************
	4059	! Set-up soil and vegetation Parameters in the case when
	4060	! snow disappears during the forecast and snow parameters
	4061	! shold be replaced by surface parameters according to
	4062	! soil and vegetation types in this point.
	4063	!
	4064	! Output:
	4065	!
	4066	!
	4067	! Soil parameters:
	4068	! DQM: MAX soil moisture content - MIN (m^3/m^3)
	4069	! REF: Reference soil moisture (m^3/m^3)
	4070	! WILT: Wilting PT soil moisture contents (m^3/m^3)
	4071	! QMIN: Air dry soil moist content limits (m^3/m^3)
	4072	! PSIS: SAT soil potential coefs. (m)
	4073	! KSAT: SAT soil diffusivity/conductivity coefs. (m/s)
	4074	! BCLH: Soil diffusivity/conductivity exponent.
	4075	!
	4076	! ************************************************************************
	4077	IMPLICIT NONE
	4078	!---------------------------------------------------------------------------
	4079	integer, parameter :: nsoilclas=19
	4080	integer, parameter :: nvegclas=24
	4081	integer, parameter :: iwater=16
	4082	integer, parameter :: ilsnow=99
	4083
	4084
	4085	!--- soiltyp classification according to STATSGO(nclasses=16)
	4086	!
	4087	! 1 SAND SAND
	4088	! 2 LOAMY SAND LOAMY SAND
	4089	! 3 SANDY LOAM SANDY LOAM
	4090	! 4 SILT LOAM SILTY LOAM
	4091	! 5 SILT SILTY LOAM
	4092	! 6 LOAM LOAM
	4093	! 7 SANDY CLAY LOAM SANDY CLAY LOAM
	4094	! 8 SILTY CLAY LOAM SILTY CLAY LOAM
	4095	! 9 CLAY LOAM CLAY LOAM
	4096	! 10 SANDY CLAY SANDY CLAY
	4097	! 11 SILTY CLAY SILTY CLAY
	4098	! 12 CLAY LIGHT CLAY
	4099	! 13 ORGANIC MATERIALS LOAM
	4100	! 14 WATER
	4101	! 15 BEDROCK
	4102	! Bedrock is reclassified as class 14
	4103	! 16 OTHER (land-ice)
	4104	! 17 Playa
	4105	! 18 Lava
	4106	! 19 White Sand
	4107	!
	4108	!----------------------------------------------------------------------
	4109	REAL LQMA(nsoilclas),LRHC(nsoilclas), &
	4110	LPSI(nsoilclas),LQMI(nsoilclas), &
	4111	LBCL(nsoilclas),LKAS(nsoilclas), &
	4112	LWIL(nsoilclas),LREF(nsoilclas), &
	4113	DATQTZ(nsoilclas)
	4114	!-- LQMA Rawls et al.[1982]
	4115	! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
	4116	! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
	4117	!---
	4118	!-- Clapp, R. and G. Hornberger, 1978: Empirical equations for some soil
	4119	! hydraulic properties, Water Resour. Res., 14, 601-604.
	4120
	4121	!-- Clapp et al. [1978]
	4122	DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
	4123	0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
	4124	0.20, 0.435, 0.468, 0.200, 0.339/
	4125
	4126	!-- LREF Rawls et al.[1982]
	4127	! DATA LREF /0.091, 0.125, 0.207, 0.330, 0.360, 0.270, 0.255,
	4128	! & 0.366, 0.318, 0.339, 0.387, 0.396, 0.329, 1.0, 0.108, 0.283/
	4129
	4130	!-- Clapp et al. [1978]
	4131	DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
	4132	0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
	4133	0.1, 0.249, 0.454, 0.17, 0.236/
	4134
	4135	!-- LWIL Rawls et al.[1982]
	4136	! DATA LWIL/0.033, 0.055, 0.095, 0.133, 0.133, 0.117, 0.148,
	4137	! & 0.208, 0.197, 0.239, 0.250, 0.272, 0.066, 0.0, 0.006, 0.029/
	4138
	4139	!-- Clapp et al. [1978]
	4140	DATA LWIL/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
	4141	0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
	4142	0.006, 0.114, 0.030, 0.006, 0.01/
	4143
	4144	! DATA LQMI/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, 0.067,
	4145	! & 0.120, 0.103, 0.100, 0.126, 0.138, 0.066, 0.0, 0.006, 0.028/
	4146
	4147	!-- Carsel and Parrish [1988]
	4148	DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
	4149	0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
	4150	0.004, 0.065, 0.020, 0.004, 0.008/
	4151
	4152	!-- LPSI Cosby et al[1984]
	4153	! DATA LPSI/0.060, 0.036, 0.141, 0.759, 0.759, 0.355, 0.135,
	4154	! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
	4155	! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
	4156
	4157	!-- Clapp et al. [1978]
	4158	DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
	4159	0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
	4160	0.121, 0.218, 0.468, 0.069, 0.069/
	4161
	4162	!-- LKAS Rawls et al.[1982]
	4163	! DATA LKAS/5.83E-5, 1.70E-5, 7.19E-6, 1.89E-6, 1.89E-6,
	4164	! & 3.67E-6, 1.19E-6, 4.17E-7, 6.39E-7, 3.33E-7, 2.50E-7,
	4165	! & 1.67E-7, 3.38E-6, 0.0, 1.41E-4, 1.41E-5/
	4166
	4167	!-- Clapp et al. [1978]
	4168	DATA LKAS/1.76E-4, 1.56E-4, 3.47E-5, 7.20E-6, 7.20E-6, &
	4169	6.95E-6, 6.30E-6, 1.70E-6, 2.45E-6, 2.17E-6, &
	4170	1.03E-6, 1.28E-6, 6.95E-6, 0.0, 1.41E-4, &
	4171	3.47E-5, 1.28E-6, 1.41E-4, 1.76E-4/
	4172
	4173	!-- LBCL Cosby et al [1984]
	4174	! DATA LBCL/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, 6.66,
	4175	! & 8.72, 8.17, 10.73, 10.39, 11.55, 5.25, 0.0, 2.79, 4.26/
	4176
	4177	!-- Clapp et al. [1978]
	4178	DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
	4179	7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
	4180	4.05, 4.90, 11.55, 2.79, 2.79/
	4181
	4182	DATA LRHC /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
	4183	1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
	4184
	4185	DATA DATQTZ/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
	4186	0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
	4187
	4188	!--------------------------------------------------------------------------
	4189	!
	4190	! USGS Vegetation Types
	4191	!
	4192	! 1: Urban and Built-Up Land
	4193	! 2: Dryland Cropland and Pasture
	4194	! 3: Irrigated Cropland and Pasture
	4195	! 4: Mixed Dryland/Irrigated Cropland and Pasture
	4196	! 5: Cropland/Grassland Mosaic
	4197	! 6: Cropland/Woodland Mosaic
	4198	! 7: Grassland
	4199	! 8: Shrubland
	4200	! 9: Mixed Shrubland/Grassland
	4201	! 10: Savanna
	4202	! 11: Deciduous Broadleaf Forest
	4203	! 12: Deciduous Needleleaf Forest
	4204	! 13: Evergreen Broadleaf Forest
	4205	! 14: Evergreen Needleleaf Fores
	4206	! 15: Mixed Forest
	4207	! 16: Water Bodies
	4208	! 17: Herbaceous Wetland
	4209	! 18: Wooded Wetland
	4210	! 19: Barren or Sparsely Vegetated
	4211	! 20: Herbaceous Tundra
	4212	! 21: Wooded Tundra
	4213	! 22: Mixed Tundra
	4214	! 23: Bare Ground Tundra
	4215	! 24: Snow or Ice
	4216	!
	4217	! 25: Playa
	4218	! 26: Lava
	4219	! 27: White Sand
	4220
	4221
	4222	!---- Below are the arrays for the vegetation parameters
	4223	REAL LALB(nvegclas),LMOI(nvegclas),LEMI(nvegclas), &
	4224	LROU(nvegclas),LTHI(nvegclas),LSIG(nvegclas), &
	4225	LPC(nvegclas), NROTBL(nvegclas)
	4226
	4227	!************************************************************************
	4228	!---- vegetation parameters
	4229	!
	4230	!-- USGS model
	4231	!
	4232	DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
	4233	.12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55/
	4234	DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
	4235	.95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95/
	4236	!-- Roughness length is changed for forests and some others
	4237	! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
	4238	! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
	4239	DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
	4240	.5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
	4241
	4242	DATA LMOI/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
	4243	.5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95/
	4244	!
	4245	!---- still needs to be corrected
	4246	!
	4247	! DATA LPC/ 15*.8,0.,.8,.8,.5,.5,.5,.5,.5,.0/
	4248	DATA LPC /0.4,0.3,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,5*0.55,0.,0.55,0.55, &
	4249	0.3,0.3,0.4,0.4,0.3,0./
	4250
	4251	! DATA LPC /0.6,60.8,0.7,0.75,60.8,0.,0.8,0.8, &
	4252	! 0.5,0.7,0.6,0.7,0.5,0./
	4253
	4254
	4255	!***************************************************************************
	4256
	4257
	4258	INTEGER :: &
	4259	IVGTYP, &
	4260	ISLTYP
	4261
	4262	LOGICAL, INTENT(IN ) :: myj
	4263
	4264	REAL , &
	4265	INTENT ( OUT) :: pc
	4266
	4267	REAL , &
	4268	INTENT (INOUT ) :: emiss, &
	4269	znt
	4270	!--- soil properties
	4271	REAL , &
	4272	INTENT( OUT) :: RHOCS, &
	4273	BCLH, &
	4274	DQM, &
	4275	KSAT, &
	4276	PSIS, &
	4277	QMIN, &
	4278	QWRTZ, &
	4279	REF, &
	4280	WILT
	4281
	4282	INTEGER, DIMENSION( 1:(nvegclas+3) ) , &
	4283	INTENT ( OUT) :: iforest
	4284
	4285
	4286
	4287	INTEGER, DIMENSION( 1:(nvegclas+3) ) :: if1
	4288	INTEGER :: kstart, kfin, lstart, lfin
	4289	INTEGER :: i,j,k
	4290
	4291	!***********************************************************************
	4292	! DATA ZS1/0.0,0.05,0.20,0.40,1.6,3.0/ ! o - levels in soil
	4293	! DATA ZS2/0.0,0.025,0.125,0.30,1.,2.3/ ! x - levels in soil
	4294	DATA IF1/120,1,1,1,120/
	4295
	4296	do k=1,nvegclas+3
	4297	iforest(k)=if1(k)
	4298	enddo
	4299
	4300
	4301	! EMISS = LEMI(IVGTYP)
	4302	EMISS = LEMITBL(IVGTYP)
	4303	! When MYJ sfc scheme is used - better use recommended in MYJSFCINIT
	4304	! values of roughness length, and not redefine it here.
	4305	! The table in this routine is the one we use in RUC with RUC LSM.
	4306
	4307	!tgs 3 Oct 2007 - MYJSFCINIT roughness produced unrealistic fluxes,
	4308	!tgs - will not use it any more!
	4309	!tgs if (.not. myj) then
	4310	! ZNT = LROU(IVGTYP)
	4311	ZNT = Z0TBL(IVGTYP)
	4312	!tgs endif
	4313
	4314	! PC = LPC (IVGTYP)
	4315	PC = PCTBL(IVGTYP)
	4316	! RHOCS = LRHC(ISLTYP)*1.E6
	4317	RHOCS = HC(ISLTYP)*1.E6
	4318
	4319	!tgs - parameters from SOILPARM.TBL
	4320	BCLH = BB(ISLTYP)
	4321	DQM = MAXSMC(ISLTYP)- &
	4322	DRYSMC(ISLTYP)
	4323	KSAT = SATDK(ISLTYP)
	4324	PSIS = - SATPSI(ISLTYP)
	4325	QMIN = DRYSMC(ISLTYP)
	4326	REF = REFSMC(ISLTYP)
	4327	WILT = WLTSMC(ISLTYP)
	4328	QWRTZ = QTZ(ISLTYP)
	4329
	4330	! BCLH = LBCL(ISLTYP)
	4331	! DQM = LQMA(ISLTYP)- &
	4332	! LQMI(ISLTYP)
	4333	! KSAT = LKAS(ISLTYP)
	4334	! PSIS = - LPSI(ISLTYP)
	4335	! QMIN = LQMI(ISLTYP)
	4336	! REF = LREF(ISLTYP)
	4337	! WILT = LWIL(ISLTYP)
	4338	! QWRTZ = DATQTZ(ISLTYP)
	4339
	4340	!--------------------------------------------------------------------------
	4341	END SUBROUTINE SOILVEGIN
	4342	!--------------------------------------------------------------------------
	4343
	4344
	4345	SUBROUTINE SNOWFREE (ivgtyp,myj,emiss,znt,iland)
	4346	!************************************************************************
	4347	! Set-up soil and vegetation Parameters in the case when
	4348	! snow disappears during the forecast and snow parameters
	4349	! shold be replaced by surface parameters according to
	4350	! soil and vegetation types in this point.
	4351	!
	4352	!***************************************************************************
	4353	IMPLICIT NONE
	4354	!---------------------------------------------------------------------------
	4355	integer, parameter :: nvegclas=24
	4356
	4357
	4358	INTEGER :: IVGTYP
	4359
	4360	LOGICAL, INTENT(IN ) :: myj
	4361
	4362	REAL, INTENT(INOUT) :: &
	4363	emiss, &
	4364	znt
	4365	INTEGER, INTENT(INOUT) :: ILAND
	4366
	4367	!---- Below are the arrays for the vegetation parameters
	4368	REAL, DIMENSION( 1:nvegclas ) :: LALB, &
	4369	LEMI, &
	4370	LROU_MYJ,&
	4371	LROU
	4372
	4373	!************************************************************************
	4374	!-- USGS model
	4375	!
	4376	DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
	4377	.12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55/
	4378	DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
	4379	.95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95/
	4380	!-- Roughness length is changed for forests and some others
	4381	! next 2 lines - table from RUC
	4382	! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
	4383	! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
	4384
	4385	DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
	4386	.5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
	4387
	4388	! With MYJSFC better use the table from MYJSFCINIT
	4389	DATA LROU_MYJ/1.0,.07,.07,.07,.07,.15,.08,.03,.05,.86,.8,.85, &
	4390	2.65,1.09,.8,.001,.04,.05,.01,.04,.06,.05,.03,.001/
	4391
	4392
	4393
	4394	!--------------------------------------------------------------------------
	4395
	4396	EMISS = LEMITBL(IVGTYP)
	4397	!tgs 3 Oct 2007 - LROU_MYJ gives unrealistic surface fluxes with RUC LSM with RUC LSM
	4398	! if(myj) then
	4399	! ZNT = LROU_MYJ(IVGTYP)
	4400	! else
	4401	ZNT = Z0TBL(IVGTYP)
	4402	!!! ZNT = LROU(IVGTYP)
	4403	! endif
	4404	ILAND = IVGTYP
	4405	! ---
	4406
	4407	! RETURN
	4408	! END
	4409	!--------------------------------------------------------------------------
	4410	END SUBROUTINE SNOWFREE
	4411
	4412	SUBROUTINE LSMRUCINIT( SMFR3D,TSLB,SMOIS,ISLTYP,mavail, &
	4413	nzs, restart, &
	4414	allowed_to_read , &
	4415	ids,ide, jds,jde, kds,kde, &
	4416	ims,ime, jms,jme, kms,kme, &
	4417	its,ite, jts,jte, kts,kte )
	4418	#ifdef WRF_CHEM
	4419	USE module_data_gocart_dust
	4420	#endif
	4421	IMPLICIT NONE
	4422
	4423
	4424	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	4425	ims,ime, jms,jme, kms,kme, &
	4426	its,ite, jts,jte, kts,kte, &
	4427	nzs
	4428
	4429	REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
	4430	INTENT(IN) :: TSLB, &
	4431	SMOIS
	4432
	4433	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
	4434	INTENT(INOUT) :: ISLTYP
	4435
	4436	REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
	4437	INTENT(INOUT) :: SMFR3D
	4438
	4439	REAL, DIMENSION( ims:ime, jms:jme ) , &
	4440	INTENT(INOUT) :: MAVAIL
	4441
	4442	REAL, DIMENSION ( 1:nzs ) :: SOILIQW
	4443
	4444	LOGICAL , INTENT(IN) :: restart, allowed_to_read
	4445
	4446	!
	4447	INTEGER :: I,J,L,itf,jtf
	4448	REAL :: RIW,XLMELT,TLN,DQM,REF,PSIS,QMIN,BCLH
	4449
	4450	INTEGER :: errflag
	4451
	4452	! itf=min0(ite,ide-1)
	4453	! jtf=min0(jte,jde-1)
	4454
	4455
	4456	RIW=900.*1.e-3
	4457	XLMELT=3.335E+5
	4458
	4459	! initialize three LSM related tables
	4460	IF ( allowed_to_read ) THEN
	4461	CALL wrf_message( 'INITIALIZE THREE LSM RELATED TABLES' )
	4462	CALL RUCLSM_PARM_INIT
	4463	ENDIF
	4464
	4465	IF(.not.restart)THEN
	4466
	4467	itf=min0(ite,ide-1)
	4468	jtf=min0(jte,jde-1)
	4469
	4470	errflag = 0
	4471	DO j = jts,jtf
	4472	DO i = its,itf
	4473	IF ( ISLTYP( i,j ) .LT. 1 ) THEN
	4474	errflag = 1
	4475	WRITE(err_message,*)"module_sf_ruclsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
	4476	CALL wrf_message(err_message)
	4477	ENDIF
	4478	ENDDO
	4479	ENDDO
	4480	IF ( errflag .EQ. 1 ) THEN
	4481	CALL wrf_error_fatal( "module_sf_ruclsm.F: lsminit: out of range value "// &
	4482	"of ISLTYP. Is this field in the input?" )
	4483	ENDIF
	4484
	4485	#ifdef WRF_CHEM
	4486	!
	4487	! need this parameter for dust parameterization in wrf/chem
	4488	!
	4489	do I=1,NSLTYPE
	4490	porosity(i)=maxsmc(i)
	4491	enddo
	4492	#endif
	4493	DO J=jts,jtf
	4494	DO I=its,itf
	4495
	4496	! CALL SOILIN ( ISLTYP(I,J), DQM, REF, PSIS, QMIN, BCLH )
	4497
	4498
	4499	!--- Computation of volumetric content of ice in soil
	4500	!--- and initialize MAVAIL
	4501	DQM = MAXSMC (ISLTYP(I,J)) - &
	4502	DRYSMC (ISLTYP(I,J))
	4503	REF = REFSMC (ISLTYP(I,J))
	4504	PSIS = - SATPSI (ISLTYP(I,J))
	4505	QMIN = DRYSMC (ISLTYP(I,J))
	4506	BCLH = BB (ISLTYP(I,J))
	4507
	4508
	4509	!!! IF (.not.restart) THEN
	4510
	4511	if(isltyp(i,j).ne.14) then
	4512	mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/dqm))
	4513	! mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/(ref-qmin)))
	4514	else
	4515	mavail(i,j) = 1.
	4516	endif
	4517	DO L=1,NZS
	4518	if(isltyp(i,j).ne.14) then
	4519	!-- for land points initialize soil ice
	4520	tln=log(TSLB(i,l,j)/273.15)
	4521
	4522	if(tln.lt.0.) then
	4523	soiliqw(l)=(dqm+qmin)(XLMELT &
	4524	(tslb(i,l,j)-273.15)/tslb(i,l,j)/9.81/psis) &
	4525	**(-1./bclh)-qmin
	4526	soiliqw(l)=max(0.,soiliqw(l))
	4527	soiliqw(l)=min(soiliqw(l),smois(i,l,j))
	4528	smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/RIW
	4529
	4530	else
	4531	smfr3d(i,l,j)=0.
	4532	endif
	4533	else
	4534	!-- for water points
	4535	smfr3d(i,l,j)=0.
	4536	endif
	4537
	4538	ENDDO
	4539	! ENDIF
	4540
	4541	ENDDO
	4542	ENDDO
	4543
	4544	ENDIF
	4545
	4546	END SUBROUTINE lsmrucinit
	4547	!
	4548	!-----------------------------------------------------------------
	4549	SUBROUTINE RUCLSM_PARM_INIT
	4550	!-----------------------------------------------------------------
	4551
	4552	character*8 :: MMINLU, MMINSL
	4553
	4554	MMINLU='USGS-RUC'
	4555	MMINSL='STAS-RUC'
	4556	call RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
	4557
	4558	!-----------------------------------------------------------------
	4559	END SUBROUTINE RUCLSM_PARM_INIT
	4560	!-----------------------------------------------------------------
	4561
	4562	!-----------------------------------------------------------------
	4563	SUBROUTINE RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
	4564	!-----------------------------------------------------------------
	4565
	4566	IMPLICIT NONE
	4567
	4568	integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
	4569	integer :: ierr
	4570	INTEGER , PARAMETER :: OPEN_OK = 0
	4571
	4572	character*8 :: MMINLU, MMINSL
	4573	character*128 :: mess , message
	4574	logical, external :: wrf_dm_on_monitor
	4575
	4576
	4577	!-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
	4578	! ALBBCK: SFC albedo (in percentage)
	4579	! Z0: Roughness length (m)
	4580	! LEMI: Emissivity
	4581	! PC: Plant coefficient for transpiration function
	4582	! -- the rest of the parameters are read in but not used currently
	4583	! SHDFAC: Green vegetation fraction (in percentage)
	4584	! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
	4585	! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
	4586	! the monthly green vegetation data
	4587	! CMXTBL: MAX CNPY Capacity (m)
	4588	! NROTBL: Rooting depth (layer)
	4589	! RSMIN: Mimimum stomatal resistance (s m-1)
	4590	! RSMAX: Max. stomatal resistance (s m-1)
	4591	! RGL: Parameters used in radiation stress function
	4592	! HS: Parameter used in vapor pressure deficit functio
	4593	! TOPT: Optimum transpiration air temperature. (K)
	4594	! CMCMAX: Maximum canopy water capacity
	4595	! CFACTR: Parameter used in the canopy inteception calculati
	4596	! SNUP: Threshold snow depth (in water equivalent m) that
	4597	! implies 100% snow cover
	4598	! LAI: Leaf area index (dimensionless)
	4599	! MAXALB: Upper bound on maximum albedo over deep snow
	4600	!
	4601	!-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
	4602	!
	4603
	4604	IF ( wrf_dm_on_monitor() ) THEN
	4605
	4606	OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
	4607	IF(ierr .NE. OPEN_OK ) THEN
	4608	WRITE(message,FMT='(A)') &
	4609	'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
	4610	CALL wrf_error_fatal ( message )
	4611	END IF
	4612
	4613	WRITE ( mess, * ) 'INPUT VEGPARM FOR ',MMINLU
	4614	CALL wrf_message( mess )
	4615
	4616	LUMATCH=0
	4617
	4618	READ (19,*)
	4619	READ (19,2000,END=2002)LUTYPE
	4620	READ (19,*)LUCATS,IINDEX
	4621	2000 FORMAT (A8)
	4622
	4623	IF(LUTYPE.NE.MMINLU)THEN
	4624	DO LC=1,LUCATS
	4625	READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC), &
	4626	SHDTBL(LC),NROTBL(LC),RSTBL(LC),RGLTBL(LC), &
	4627	HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
	4628	ENDDO
	4629	READ (19,*)
	4630	READ (19,*)TOPT_DATA
	4631	READ (19,*)
	4632	READ (19,*)CMCMAX_DATA
	4633	READ (19,*)
	4634	READ (19,*)CFACTR_DATA
	4635	READ (19,*)
	4636	READ (19,*)RSMAX_DATA
	4637	READ (19,*)
	4638	READ (19,*)BARE
	4639	ENDIF
	4640
	4641	READ (19,*)
	4642	READ (19,2000,END=2002)LUTYPE
	4643	READ (19,*)LUCATS,IINDEX
	4644
	4645
	4646	IF(LUTYPE.EQ.MMINLU)THEN
	4647	WRITE( mess , * ) 'VEGPARM FOR ',LUTYPE,' FOUND', &
	4648	LUCATS,' CATEGORIES'
	4649	CALL wrf_message( mess )
	4650	LUMATCH=1
	4651	ENDIF
	4652
	4653	IF(LUTYPE.EQ.MMINLU)THEN
	4654	DO LC=1,LUCATS
	4655	READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),LEMITBL(LC),PCTBL(LC), &
	4656	SHDTBL(LC),NROTBL(LC),RSTBL(LC),RGLTBL(LC), &
	4657	HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
	4658	ENDDO
	4659	!
	4660	READ (19,*)
	4661	READ (19,*)TOPT_DATA
	4662	READ (19,*)
	4663	READ (19,*)CMCMAX_DATA
	4664	READ (19,*)
	4665	READ (19,*)CFACTR_DATA
	4666	READ (19,*)
	4667	READ (19,*)RSMAX_DATA
	4668	READ (19,*)
	4669	READ (19,*)BARE
	4670	ENDIF
	4671	!
	4672	2002 CONTINUE
	4673
	4674	CLOSE (19)
	4675	ENDIF
	4676
	4677	CALL wrf_dm_bcast_string ( LUTYPE , 8 )
	4678	CALL wrf_dm_bcast_integer ( LUCATS , 1 )
	4679	CALL wrf_dm_bcast_integer ( IINDEX , 1 )
	4680	CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
	4681	CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
	4682	CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
	4683	CALL wrf_dm_bcast_real ( LEMITBL , NLUS )
	4684	CALL wrf_dm_bcast_real ( PCTBL , NLUS )
	4685	CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
	4686	CALL wrf_dm_bcast_real ( NROTBL , NLUS )
	4687	CALL wrf_dm_bcast_real ( RSTBL , NLUS )
	4688	CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
	4689	CALL wrf_dm_bcast_real ( HSTBL , NLUS )
	4690	CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
	4691	CALL wrf_dm_bcast_real ( LAITBL , NLUS )
	4692	CALL wrf_dm_bcast_real ( MAXALB , NLUS )
	4693	CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
	4694	CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
	4695	CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
	4696	CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
	4697	CALL wrf_dm_bcast_integer ( BARE , 1 )
	4698
	4699	!
	4700	!-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
	4701	!
	4702	IF ( wrf_dm_on_monitor() ) THEN
	4703	OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
	4704	IF(ierr .NE. OPEN_OK ) THEN
	4705	WRITE(message,FMT='(A)') &
	4706	'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
	4707	CALL wrf_error_fatal ( message )
	4708	END IF
	4709
	4710	WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICAION = ',MMINSL
	4711	CALL wrf_message( mess )
	4712
	4713	LUMATCH=0
	4714
	4715	READ (19,*)
	4716	READ (19,2000,END=2003)SLTYPE
	4717	READ (19,*)SLCATS,IINDEX
	4718	IF(SLTYPE.NE.MMINSL)THEN
	4719	DO LC=1,SLCATS
	4720	READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
	4721	REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
	4722	WLTSMC(LC), QTZ(LC)
	4723	ENDDO
	4724	ENDIF
	4725	READ (19,*)
	4726	READ (19,2000,END=2003)SLTYPE
	4727	READ (19,*)SLCATS,IINDEX
	4728
	4729	IF(SLTYPE.EQ.MMINSL)THEN
	4730	WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
	4731	SLCATS,' CATEGORIES'
	4732	CALL wrf_message ( mess )
	4733	LUMATCH=1
	4734	ENDIF
	4735	IF(SLTYPE.EQ.MMINSL)THEN
	4736	DO LC=1,SLCATS
	4737	READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
	4738	REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
	4739	WLTSMC(LC), QTZ(LC)
	4740	ENDDO
	4741	ENDIF
	4742
	4743	2003 CONTINUE
	4744
	4745	CLOSE (19)
	4746	ENDIF
	4747
	4748	CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
	4749	CALL wrf_dm_bcast_string ( SLTYPE , 8 )
	4750	CALL wrf_dm_bcast_string ( MMINSL , 8 ) ! since this is reset above, see oct2 ^
	4751	CALL wrf_dm_bcast_integer ( SLCATS , 1 )
	4752	CALL wrf_dm_bcast_integer ( IINDEX , 1 )
	4753	CALL wrf_dm_bcast_real ( BB , NSLTYPE )
	4754	CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
	4755	CALL wrf_dm_bcast_real ( HC , NSLTYPE )
	4756	CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
	4757	CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
	4758	CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
	4759	CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
	4760	CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
	4761	CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
	4762	CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
	4763
	4764	IF(LUMATCH.EQ.0)THEN
	4765	CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
	4766	CALL wrf_message( 'MATCH SOILPARM TABLE' )
	4767	CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
	4768	ENDIF
	4769
	4770	!
	4771	!-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
	4772	!
	4773	IF ( wrf_dm_on_monitor() ) THEN
	4774	OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
	4775	IF(ierr .NE. OPEN_OK ) THEN
	4776	WRITE(message,FMT='(A)') &
	4777	'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
	4778	CALL wrf_error_fatal ( message )
	4779	END IF
	4780
	4781	READ (19,*)
	4782	READ (19,*)
	4783	READ (19,*) NUM_SLOPE
	4784
	4785	SLPCATS=NUM_SLOPE
	4786
	4787	DO LC=1,SLPCATS
	4788	READ (19,*)SLOPE_DATA(LC)
	4789	ENDDO
	4790
	4791	READ (19,*)
	4792	READ (19,*)SBETA_DATA
	4793	READ (19,*)
	4794	READ (19,*)FXEXP_DATA
	4795	READ (19,*)
	4796	READ (19,*)CSOIL_DATA
	4797	READ (19,*)
	4798	READ (19,*)SALP_DATA
	4799	READ (19,*)
	4800	READ (19,*)REFDK_DATA
	4801	READ (19,*)
	4802	READ (19,*)REFKDT_DATA
	4803	READ (19,*)
	4804	READ (19,*)FRZK_DATA
	4805	READ (19,*)
	4806	READ (19,*)ZBOT_DATA
	4807	READ (19,*)
	4808	READ (19,*)CZIL_DATA
	4809	READ (19,*)
	4810	READ (19,*)SMLOW_DATA
	4811	READ (19,*)
	4812	READ (19,*)SMHIGH_DATA
	4813	CLOSE (19)
	4814	ENDIF
	4815
	4816	CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
	4817	CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
	4818	CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
	4819	CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
	4820	CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
	4821	CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
	4822	CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
	4823	CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
	4824	CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
	4825	CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
	4826	CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
	4827	CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
	4828	CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
	4829	CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
	4830
	4831
	4832	!-----------------------------------------------------------------
	4833	END SUBROUTINE RUCLSM_SOILVEGPARM
	4834	!-----------------------------------------------------------------
	4835
	4836
	4837	SUBROUTINE SOILIN (ISLTYP, DQM, REF, PSIS, QMIN, BCLH )
	4838
	4839	!--- soiltyp classification according to STATSGO(nclasses=16)
	4840	!
	4841	! 1 SAND SAND
	4842	! 2 LOAMY SAND LOAMY SAND
	4843	! 3 SANDY LOAM SANDY LOAM
	4844	! 4 SILT LOAM SILTY LOAM
	4845	! 5 SILT SILTY LOAM
	4846	! 6 LOAM LOAM
	4847	! 7 SANDY CLAY LOAM SANDY CLAY LOAM
	4848	! 8 SILTY CLAY LOAM SILTY CLAY LOAM
	4849	! 9 CLAY LOAM CLAY LOAM
	4850	! 10 SANDY CLAY SANDY CLAY
	4851	! 11 SILTY CLAY SILTY CLAY
	4852	! 12 CLAY LIGHT CLAY
	4853	! 13 ORGANIC MATERIALS LOAM
	4854	! 14 WATER
	4855	! 15 BEDROCK
	4856	! Bedrock is reclassified as class 14
	4857	! 16 OTHER (land-ice)
	4858	! extra classes from Fei Chen
	4859	! 17 Playa
	4860	! 18 Lava
	4861	! 19 White Sand
	4862	!
	4863	!----------------------------------------------------------------------
	4864	integer, parameter :: nsoilclas=19
	4865
	4866	integer, intent ( in) :: isltyp
	4867	real, intent ( out) :: dqm,ref,qmin,psis
	4868
	4869	REAL LQMA(nsoilclas),LREF(nsoilclas),LBCL(nsoilclas), &
	4870	LPSI(nsoilclas),LQMI(nsoilclas)
	4871
	4872	!-- LQMA Rawls et al.[1982]
	4873	! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
	4874	! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
	4875	!---
	4876	!-- Clapp, R. and G. Hornberger, Empirical equations for some soil
	4877	! hydraulic properties, Water Resour. Res., 14,601-604,1978.
	4878	!-- Clapp et al. [1978]
	4879	DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
	4880	0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
	4881	0.20, 0.435, 0.468, 0.200, 0.339/
	4882
	4883	!-- Clapp et al. [1978]
	4884	DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
	4885	0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
	4886	0.1, 0.249, 0.454, 0.17, 0.236/
	4887
	4888	!-- Carsel and Parrish [1988]
	4889	DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
	4890	0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
	4891	0.004, 0.065, 0.020, 0.004, 0.008/
	4892
	4893	!-- Clapp et al. [1978]
	4894	DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
	4895	0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
	4896	0.121, 0.218, 0.468, 0.069, 0.069/
	4897
	4898	!-- Clapp et al. [1978]
	4899	DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
	4900	7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
	4901	4.05, 4.90, 11.55, 2.79, 2.79/
	4902
	4903
	4904	DQM = LQMA(ISLTYP)- &
	4905	LQMI(ISLTYP)
	4906	REF = LREF(ISLTYP)
	4907	PSIS = - LPSI(ISLTYP)
	4908	QMIN = LQMI(ISLTYP)
	4909	BCLH = LBCL(ISLTYP)
	4910
	4911	END SUBROUTINE SOILIN
	4912
	4913	END MODULE module_sf_ruclsm

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

Download in other formats:

Original Format