Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_ruclsm.F @ 3431

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

Line
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