Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_ruclsm.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

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

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

Download in other formats:

Original Format