Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_ruclsm.F @ 3567

Last change on this file since 3567 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 162.5 KB

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

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

Download in other formats:

Original Format