Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_lsm_nmm.F @ 3094

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

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3	MODULE MODULE_SF_LSM_NMM
4
5	USE MODULE_MPP
6	USE MODULE_MODEL_CONSTANTS
7
8	REAL, SAVE :: SCFX(30)
9
10	INTEGER, SAVE :: ISEASON
11	CHARACTER*256 :: errmess
12
13	CONTAINS
14
15	!-----------------------------------------------------------------------
16	SUBROUTINE NMMLSM(DZ8W,QV3D,P8W3D,RHO3D, &
17	& T3D,TH3D,TSK,CHS, &
18	& HFX,QFX,QGH,GSW,GLW,ELFLX,RMOL, & ! added for WRF CHEM
19	& SMSTAV,SMSTOT,SFCRUNOFF, &
20	& UDRUNOFF,IVGTYP,ISLTYP,VEGFRA,SFCEVP,POTEVP, &
21	& GRDFLX,SFCEXC,ACSNOW,ACSNOM,SNOPCX, &
22	& ALBSF,TMN,XLAND,XICE,QZ0, &
23	& TH2,Q2,SNOWC,CHS2,QSFC,TBOT,CHKLOWQ,RAINBL, &
24	& NUM_SOIL_LAYERS,DT,DZS,ITIMESTEP, &
25	& SMOIS,TSLB,SNOW,CANWAT,CPM,ROVCP,SR, &
26	& ALB,SNOALB,SMLIQ,SNOWH, &
27	& IDS,IDE, JDS,JDE, KDS,KDE, &
28	& IMS,IME, JMS,JME, KMS,KME, &
29	& ITS,ITE, JTS,JTE, KTS,KTE )
30	!-----------------------------------------------------------------------
31	!-----------------------------------------------------------------------
32	IMPLICIT NONE
33	!-----------------------------------------------------------------------
34	!-----------------------------------------------------------------------
35	!-- DZ8W thickness of layers (m)
36	!-- T3D temperature (K)
37	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
38	!-- P8W3D 3D pressure on layer interfaces (Pa)
39	!-- FLHC exchange coefficient for heat (m/s)
40	!-- FLQC exchange coefficient for moisture (m/s)
41	!-- PSFC surface pressure (Pa)
42	!-- XLAND land mask (1 for land, 2 for water)
43	!-- TMN soil temperature at lower boundary (K)
44	!-- HFX upward heat flux at the surface (W/m^2)
45	!-- QFX upward moisture flux at the surface (kg/m^2/s)
46	!-- TSK surface temperature (K)
47	!-- GSW NET downward short wave flux at ground surface (W/m^2)
48	!-- GLW downward long wave flux at ground surface (W/m^2)
49	!-- ELFLX actual latent heat flux (w m-2: positive, if up from surface)
50	!-- SFCEVP accumulated surface evaporation (W/m^2)
51	!-- POTEVP accumulated potential evaporation (W/m^2)
52	!-- CAPG heat capacity for soil (J/K/m^3)
53	!-- THC thermal inertia (Cal/cm/K/s^0.5)
54	!-- TBOT bottom soil temperature (local yearly-mean sfc air temperature)
55	!-- SNOWC flag indicating snow coverage (1 for snow cover)
56	!-- EMISS surface emissivity (between 0 and 1)
57	!-- DELTSM time step (second)
58	!-- ROVCP R/CP
59	!-- SR fraction of frozen precip (0.0 to 1.0)
60	!-- XLV latent heat of melting (J/kg)
61	!-- DTMIN time step (minute)
62	!-- IFSNOW ifsnow=1 for snow-cover effects
63	!-- SVP1 constant for saturation vapor pressure (kPa)
64	!-- SVP2 constant for saturation vapor pressure (dimensionless)
65	!-- SVP3 constant for saturation vapor pressure (K)
66	!-- SVPT0 constant for saturation vapor pressure (K)
67	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
68	!-- EP2 constant for specific humidity calculation
69	! (R_d/R_v) (dimensionless)
70	!-- KARMAN Von Karman constant
71	!-- EOMEG angular velocity of earth's rotation (rad/s)
72	!-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
73	!-- STEM soil temperature in 5-layer model
74	!-- ZS depths of centers of soil layers
75	!-- DZS thicknesses of soil layers
76	!-- num_soil_layers the number of soil layers
77	!-- ACSNOW accumulated snowfall (water equivalent) (mm)
78	!-- ACSNOM accumulated snowmelt (water equivalent) (mm)
79	!-- SNOPCX snow phase change heat flux (W/m^2)
80	!-- ids start index for i in domain
81	!-- ide end index for i in domain
82	!-- jds start index for j in domain
83	!-- jde end index for j in domain
84	!-- kds start index for k in domain
85	!-- kde end index for k in domain
86	!-- ims start index for i in memory
87	!-- ime end index for i in memory
88	!-- jms start index for j in memory
89	!-- jme end index for j in memory
90	!-- kms start index for k in memory
91	!-- kme end index for k in memory
92	!-- its start index for i in tile
93	!-- ite end index for i in tile
94	!-- jts start index for j in tile
95	!-- jte end index for j in tile
96	!-- kts start index for k in tile
97	!-- kte end index for k in tile
98	!-----------------------------------------------------------------------
99	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
100	& IMS,IME,JMS,JME,KMS,KME, &
101	& ITS,ITE,JTS,JTE,KTS,KTE
102	!
103	INTEGER,INTENT(IN) :: NUM_SOIL_LAYERS,ITIMESTEP
104	!
105	REAL,INTENT(IN) :: DT,ROVCP
106	!
107	REAL,DIMENSION(IMS:IME,1:NUM_SOIL_LAYERS,JMS:JME), &
108	& INTENT(INOUT) :: SMOIS, & ! new
109	SMLIQ, & ! new
110	TSLB !
111
112	REAL,DIMENSION(1:NUM_SOIL_LAYERS),INTENT(IN) :: DZS
113	!
114	REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
115	& TSK, & !was TGB (temperature)
116	& HFX, &
117	& QFX, &
118	& QSFC,&
119	& SNOW, & !new
120	& SNOWH, & !new
121	& ALB, &
122	& SNOALB, &
123	& ALBSF, &
124	& SNOWC, &
125	& CANWAT, & ! new
126	& SMSTAV, &
127	& SMSTOT, &
128	& SFCRUNOFF, &
129	& UDRUNOFF, &
130	& SFCEVP, &
131	& POTEVP, &
132	& GRDFLX, &
133	& ACSNOW, &
134	& ACSNOM, &
135	& SNOPCX, &
136	& Q2, &
137	& TH2, &
138	& SFCEXC
139
140	INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IVGTYP, &
141	ISLTYP
142
143	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: TMN, &
144	XLAND, &
145	XICE, &
146	VEGFRA, &
147	GSW, &
148	GLW, &
149	QZ0, &
150	SR
151
152	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: QV3D, &
153	P8W3D, &
154	RHO3D, &
155	TH3D, &
156	T3D, &
157	DZ8W
158
159	!
160	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: RAINBL
161	!
162	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CHS2, &
163	CHS, &
164	QGH, &
165	CPM
166	!
167	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: TBOT
168	!
169	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CHKLOWQ, &
170	ELFLX
171	! added for WRF-CHEM, 20041205, JM -- not used in this routine as yet
172	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: RMOL
173
174	! LOCAL VARS
175
176	REAL,DIMENSION(ITS:ITE) :: QV1D, &
177	& T1D, &
178	& TH1D, &
179	& ZA1D, &
180	& P8W1D, &
181	& PSFC1D, &
182	& RHO1D, &
183	& PREC1D
184
185	INTEGER :: I,J
186	REAL :: RATIOMX
187	!-----------------------------------------------------------------------
188	!-----------------------------------------------------------------------
189
190	DO J=JTS,JTE
191
192	DO I=ITS,ITE
193	T1D(I) = T3D(I,1,J)
194	TH1D(I) = TH3D(I,1,J)
195	!!! QV1D(I) = QV3D(I,1,J)
196	RATIOMX = QV3D(I,1,J)
197	QV1D(I) = RATIOMX/(1.+RATIOMX)
198	P8W1D(I) = (P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
199	PSFC1D(I) = P8W3D(I,1,J)
200	ZA1D(I) = 0.5*DZ8W(I,1,J)
201	RHO1D(I) = RHO3D(I,1,J)
202	PREC1D(I) = RAINBL(I,J)/DT
203	ENDDO
204
205	!FLHC = SFCEXC
206
207	!-----------------------------------------------------------------------
208	CALL SURFCE(J,ZA1D,QV1D,P8W1D,PSFC1D,RHO1D,T1D,TH1D,TSK, &
209	CHS(IMS,J),PREC1D,HFX,QFX,QGH(IMS,J),GSW,GLW, &
210	SMSTAV,SMSTOT,SFCRUNOFF, &
211	UDRUNOFF,IVGTYP,ISLTYP,VEGFRA,SFCEVP,POTEVP,GRDFLX, &
212	ELFLX,SFCEXC,ACSNOW,ACSNOM,SNOPCX, &
213	ALBSF,TMN,XLAND,XICE,QZ0, &
214	TH2,Q2,SNOWC,CHS2(IMS,J),QSFC,TBOT,CHKLOWQ, &
215	NUM_SOIL_LAYERS,DT,DZS,ITIMESTEP, &
216	SMOIS,TSLB,SNOW,CANWAT,CPM(IMS,J),ROVCP,SR, &
217	ALB,SNOALB,SMLIQ,SNOWH, &
218	IMS,IME,JMS,JME,KMS,KME, &
219	ITS,ITE,JTS,JTE,KTS,KTE )
220	!
221	ENDDO
222
223	END SUBROUTINE NMMLSM
224
225	!-----------------------------------------------------------------------
226	!-----------------------------------------------------------------------
227	SUBROUTINE SURFCE(J,ZA,QV,P8W,PSFC,RHO,T,TH,TSK,CHS,PREC,HFX,QFX, &
228	QGH,GSW,GLW,SMSTAV,SMSTOT,SFCRUNOFF,UDRUNOFF, &
229	IVGTYP,ISLTYP,VEGFRA,SFCEVP,POTEVP,GRDFLX, &
230	ELFLX,SFCEXC,ACSNOW,ACSNOM,SNOPCX, &
231	ALBSF,TMN,XLAND,XICE,QZ0, &
232	TH2,Q2,SNOWC,CHS2,QSFC,TBOT,CHKLOWQ, &
233	NUM_SOIL_LAYERS,DT,DZS,ITIMESTEP, &
234	SMOIS,TSLB,SNOW,CANWAT,CPM,ROVCP,SR, &
235	ALB,SNOALB,SMLIQ,SNOWH, &
236	IMS,IME,JMS,JME,KMS,KME, &
237	ITS,ITE,JTS,JTE,KTS,KTE )
238	!------------------------------------------------------------------------
239	IMPLICIT NONE
240	!------------------------------------------------------------------------
241	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
242	! . . .
243	! SUBPROGRAM: SURFCE CALCULATE SURFACE CONDITIONS
244	! PRGRMMR: F. CHEN DATE: 97-12-06
245	!
246	! ABSTRACT:
247	! THIS ROUTINE IS THE DRIVER FOR COMPUTATION OF GROUND CONDITIONS
248	! BY USING A LAND SURFACE MODEL (LSM).
249	!
250	! PROGRAM HISTORY LOG:
251	! 97-12-06 CHEN - ORIGINATOR
252	!
253	! REFERENCES:
254	! PAN AND MAHRT (1987) BOUN. LAYER METEOR.
255	! CHEN ET AL. (1996) J. GEOPHYS. RES.
256	! CHEN ET AL. (1997) BOUN. LAYER METEOR.
257	! CHEN and Dudhia (2000) Mon. Wea. Rev.
258	!
259	! SUBPROGRAMS CALLED:
260	! SFLX
261	!
262	! SET LOCAL PARAMETERS.
263	!----------------------------------------------------------------------
264	INTEGER, INTENT(IN ) :: IMS,IME, JMS,JME, KMS,KME, &
265	ITS,ITE, JTS,JTE, KTS,KTE, &
266	J,ITIMESTEP
267
268	INTEGER , INTENT(IN) :: NUM_SOIL_LAYERS
269
270	REAL, INTENT(IN ) :: DT,ROVCP
271
272	REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
273
274
275	REAL, PARAMETER :: PQ0=379.90516
276	REAL, PARAMETER :: TRESH=.95E0,A2=17.2693882,A3=273.16,A4=35.86, &
277	T0=273.16E0,T1=274.16E0,ROW=1.E3, &
278	ELWV=2.50E6,ELIV=XLS,ELIW=XLF, &
279	A23M4=A2*(A3-A4), RLIVWV=ELIV/ELWV, &
280	ROWLIW=ROWELIW,ROWLIV=ROWELIV,CAPA=R_D/CP
281
282	INTEGER, PARAMETER :: NROOT=3
283	!
284	REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
285	INTENT(INOUT) :: SMOIS, & ! new
286	SMLIQ, & ! new
287	TSLB ! new !STEMP
288
289
290	REAL, DIMENSION( ims:ime, jms:jme ) , &
291	INTENT(INOUT) :: TSK, & !was TGB (temperature)
292	HFX, & !new
293	QFX, & !new
294	QSFC,& !new
295	SNOW, & !new
296	SNOWH, & !new
297	ALB, &
298	SNOALB, &
299	ALBSF, &
300	SNOWC, &
301	CANWAT, & ! new
302	SMSTAV, &
303	SMSTOT, &
304	SFCRUNOFF, &
305	UDRUNOFF, &
306	SFCEVP, &
307	POTEVP, &
308	GRDFLX, &
309	ACSNOW, &
310	ACSNOM, &
311	SNOPCX
312
313	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
314	INTENT(IN ) :: IVGTYP, &
315	ISLTYP
316
317	REAL, DIMENSION( ims:ime, jms:jme ) , &
318	INTENT(IN ) :: TMN, &
319	XLAND, &
320	XICE, &
321	VEGFRA, &
322	GSW, &
323	GLW, &
324	QZ0, &
325	SR
326
327	REAL, DIMENSION( ims:ime, jms:jme ) , &
328	INTENT(INOUT) :: Q2, &
329	TH2, &
330	SFCEXC
331
332	REAL, DIMENSION( ims:ime, jms:jme ) , &
333	INTENT(OUT) :: TBOT
334
335
336	REAL, DIMENSION( ims:ime, jms:jme ) , &
337	INTENT(OUT) :: CHKLOWQ, &
338	ELFLX
339
340	REAL, DIMENSION( ims:ime ) , &
341	INTENT(IN ) :: QGH, &
342	CHS, &
343	CPM, &
344	CHS2
345
346	! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE LSM
347	REAL, DIMENSION( its:ite ) , &
348	INTENT(IN ) :: ZA, &
349	TH, &
350	QV, &
351	T, &
352	p8w, &
353	PSFC, &
354	rho, &
355	PREC ! one time step in mm
356
357	REAL, DIMENSION( its:ite ) :: TGDSA
358
359	! LOCAL VARS
360
361	REAL, DIMENSION(1:num_soil_layers) :: SMLIQ1D,SMOIS1D,STEMP1D
362
363	!----------------------------------------------------------------------
364	!*** DECLARATIONS FOR IMPLICIT NONE
365
366	REAL :: APELM,APES,FDTLIW,FDTW,Q2SAT,Z,FK,SOLDN,SFCTMP,SFCTH2, &
367	SFCPRS,PRCP,Q2K,DQSDTK,SATFLG,TBOTK,CHK,VGFRCK,T1K,LWDN, &
368	CMCK,Q2M,SNODPK,PLFLX,HFLX,GFLX,RNOF1K, &
369	RNOF2K,Q1K,SMELTK,SOILQW,SOILQM,T2K,PRESK,CHFF,STIMESTEP, &
370	ALB1D,SNOALB1D,SNOWH1D,ALBSF1D,SOLNET,FFROZP, &
371	DUM1,DUM2,DUM3,DUM4,DUM5,DUM6,DUM7
372
373	INTEGER :: I,K,NS,ICE,IVGTPK,ISLTPK,ISPTPK,NOOUT,NSOIL,LZ
374
375	!----------------------------------------------------------------------
376	!***********************************************************************
377	! START SURFCE HERE
378	!***
379	!*** SET CONSTANTS CALCULATED HERE FOR CLARITY.
380	!***
381	FDTLIW=DT/ROWLIW
382	! FDTLIV=DT/ROWLIV
383	FDTW=DT/(XLV*RHOWATER)
384	!***
385	!*** SET LSM CONSTANTS AND TIME INDEPENDENT VARIABLES
386	!*** INITIALIZE LSM HISTORICAL VARIABLES
387	!***
388	!-----------------------------------------------------------------------
389
390	NSOIL=num_soil_layers
391
392	IF(ITIMESTEP.EQ.1)THEN
393	DO 50 I=its,ite
394	!*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
395	IF((XLAND(I,J)-1.5).GE.0.)THEN
396	! check sea-ice point
397	IF(XICE(I,J).EQ.1.)PRINT*,' sea-ice at water point, I=',I, &
398	'J=',J
399	!*** Open Water Case
400	SMSTAV(I,J)=1.0
401	SMSTOT(I,J)=1.0
402	DO NS=1,NSOIL
403	SMOIS(I,NS,J)=1.0
404	TSLB(I,NS,J)=273.16 !STEMP
405	ENDDO
406	ELSE
407	IF(XICE(I,J).EQ.1.)THEN
408	!*** SEA-ICE CASE
409	SMSTAV(I,J)=1.0
410	SMSTOT(I,J)=1.0
411	DO NS=1,NSOIL
412	SMOIS(I,NS,J)=1.0
413	ENDDO
414	ENDIF
415	ENDIF
416	!
417	50 CONTINUE
418	ENDIF
419	!-----------------------------------------------------------------------
420	DO 100 I=its,ite
421	! SFCPRS=(A(KL)PSB(I,J)+PTOP+PP3D(I,J,KL)0.001)*1.E3
422	SFCPRS=p8w(I) !Pressure in middle of lowest layer
423	Q2SAT=QGH(I)
424	! CHKLOWQ(I,J)=1.
425	CHFF=CHS(I)RHO(I)CPM(I)
426	!CHKRHOCP
427	! TGDSA: potential T
428	TGDSA(I)=TSK(I,J)(1.E5/SFCPRS)*ROVCP
429	!
430	!*** CHECK FOR SATURATION AT THE LOWEST MODEL LEVEL
431	!
432	Q2K=QV(I)
433	APES=(1.E5/PSFC(I))**CAPA
434	!
435	IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
436	SATFLG=0.
437	CHKLOWQ(I,J)=0.
438	ELSE
439	SATFLG=1.0
440	CHKLOWQ(I,J)=1.
441	ENDIF
442	!
443	TBOT(I,J)=273.16
444	!***
445	!*** LOADING AND UNLOADING MM5/LSM LAND SOIL VARIABLES
446	!***
447	IF((XLAND(I,J)-1.5).GE.0.)THEN
448	!*** Water
449	HFX(I,J)=HFX(I,J)/APES
450	QFX(I,J)=QFX(I,J)*SATFLG
451	SFCEVP(I,J)=SFCEVP(I,J)+QFX(I,J)*DT
452	ELSE
453	!*** LAND OR SEA-ICE
454	!ATEC ICE=INT(XICE(I,J)+0.3)
455	IF (XICE(I,J) .GT. 0.5) THEN
456	ICE=1
457	ELSE
458	ICE=0
459	ENDIF
460	!
461	Q2K=MIN(QV(I),Q2SAT)
462	Z=ZA(I)
463	! FK=GSW(I,J)+GLW(I,J)
464	LWDN=GLW(I,J)
465	!
466	!*** GSW is net downward shortwave
467	!
468	! SOLNET=GSW(I,J)
469	!
470	!*** GSW is total downward shortwave
471	!
472	SOLDN=GSW(I,J)
473	!
474	!*** Simple use of albedo to determine total incoming solar shortwave SOLDN
475	!*** (no solar zenith angle correction)
476	!
477	! SOLDN=SOLNET/(1.-ALB(I,J))
478	SOLNET=SOLDN*(1.-ALB(I,J))
479	!
480	ALBSF1D=ALBSF(I,J)
481	SNOALB1D=SNOALB(I,J)
482	SFCTMP=T(I)
483	!!! SFCTH2=SFCTMP+(0.0097545*Z)
484	APELM=(1.E5/SFCPRS)**CAPA
485	SFCTH2=SFCTMP*APELM
486	SFCTH2=SFCTH2/APES
487	PRCP=PREC(I)
488	!!! Q2K=QV(I)
489	!!! Q2SAT=PQ0/SFCPRSEXP(A2(SFCTMP-A3)/(SFCTMP-A4))
490	DQSDTK=Q2SATA23M4/(SFCTMP-A4)*2
491	IF(ICE.EQ.0)THEN
492	TBOTK=TMN(I,J)
493	ELSE
494	TBOTK=271.16
495	ENDIF
496	CHK=CHS(I)
497	IVGTPK=IVGTYP(I,J)
498	IF(IVGTPK.EQ.0)IVGTPK=13
499	ISLTPK=ISLTYP(I,J)
500	IF(ISLTPK.EQ.0)ISLTPK=9
501	! hardwire slope type (ISPTPK)=1
502	ISPTPK=1
503	VGFRCK=VEGFRA(I,J)/100.
504	IF(IVGTPK.EQ.25) VGFRCK=0.0001
505	IF(ISLTPK.EQ.14.AND.XICE(I,J).EQ.0.)THEN
506	PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
507	PRINT*,i,j,'RESET SOIL in surfce.F'
508	! ISLTYP(I,J)=7
509	ISLTPK=7
510	ENDIF
511	T1K=TSK(I,J)
512	CMCK=CANWAT(I,J)
513	!*** convert snow depth from mm to meter
514	SNODPK=SNOW(I,J)*0.001
515	SNOWH1D=SNOWH(I,J)*0.001
516	!
517	!*** fraction of frozen precip
518	!
519	FFROZP=SR(I,J)
520	!
521	DO 70 NS=1,NSOIL
522	SMOIS1D(NS)=SMOIS(I,NS,J)
523	SMLIQ1D(NS)=SMLIQ(I,NS,J)
524	STEMP1D(NS)=TSLB(I,NS,J) !STEMP
525	70 CONTINUE
526
527	!
528	! print*,'BF SFLX','ISLTPK',ISLTPK,'IVGTPK=',IVGTPK,'SMOIS1D',&
529	! SMOIS1D,'STEMP1',STEMP1D,'VGFRCK',VGFRCK
530	!-----------------------------------------------------------------------
531	! old WRF call to SFLX
532	! CALL SFLX(ICE,SATFLG,DT,Z,NSOIL,NROOT,DZS,FK,SOLDN,SFCPRS, &
533	! PRCP,SFCTMP,SFCTH2,Q2K,Q2SAT,DQSDTK,TBOTK,CHK,CHFF, &
534	! IVGTPK,ISLTPK,VGFRCK,PLFLX,ELFLX,HFLX,GFLX,RNOF1K,RNOF2K,&
535	! Q1K,SMELTK,T1K,CMCK,SMOIS1D,STEMP1D,SNODPK,SOILQW,SOILQM)
536	!-----------------------------------------------------------------------
537	! ----------------------------------------------------------------------
538	! Ek 12 June 2002 - NEW CALL SFLX
539	! ops Eta call to SFLX ...'tailor' this to WRF
540	! CALL SFLX
541	! I (ICE,DTK,Z,NSOIL,SLDPTH,
542	! I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,SFCTH2,Q2K,SFCSPD,Q2SAT,DQSDTK,
543	! I IVGTPK,ISLTPK,ISPTPK,
544	! I VGFRCK,PTU,TBOT,ALB,SNOALB,
545	! 2 CMCK,T1K,STCK,SMCK,SH2OK,SNOWH,SNODPK,ALB2D,CHK,CMK,
546	! O PLFLX,ELFLX,HFLX,GFLX,RNOF1K,RNOF2K,Q1K,SMELTK,
547	! O SOILQW,SOILQM,DUM1,DUM2,DUM3,DUM4)
548	!-----------------------------------------------------------------------
549	CALL SFLX &
550	(FFROZP,ICE,DT,Z,NSOIL,DZS, &
551	LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,SFCTH2,Q2K,DUM5,Q2SAT, &
552	DQSDTK,IVGTPK,ISLTPK,ISPTPK, &
553	VGFRCK,DUM6,TBOTK,ALBSF1D,SNOALB1D, &
554	CMCK,T1K,STEMP1D,SMOIS1D,SMLIQ1D,SNOWH1D,SNODPK,ALB1D,CHK,DUM7, &
555	PLFLX,ELFLX(I,J),HFLX,GFLX,RNOF1K,RNOF2K,Q1K,SMELTK, &
556	SOILQW,SOILQM,DUM1,DUM2,DUM3,DUM4,I,J)
557	!-----------------------------------------------------------------------
558	!*** DIAGNOSTICS
559	! Convert the water unit into mm
560	SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RNOF1KDT1000.0
561	UDRUNOFF(I,J)=UDRUNOFF(I,J)+RNOF2KDT1000.0
562	SMSTAV(I,J)=SOILQW
563
564	!mp
565	if (abs(SMSTAV(I,J)) .lt. 3.5) then
566	else
567	write(errmess,*) 'bad SMSTAV: ', I,J,SMSTAV(I,J)
568	CALL wrf_message( errmess )
569	endif
570	!mp
571
572	SMSTOT(I,J)=SOILQM*1000.
573	SFCEXC(I,J)=CHK
574	! IF(SNOB(I,J).GT.0..OR.SICE(I,J).GT.0.)THEN
575	! QFC1(I,J)=QFC1(I,J)*RLIVWV
576	! ENDIF
577	IF(FFROZP.GT.0.5)THEN
578	ACSNOW(I,J)=ACSNOW(I,J)+PREC(I)*DT
579	ENDIF
580	IF(SNOW(I,J).GT.0.)THEN
581	ACSNOM(I,J)=ACSNOM(I,J)+SMELTK*1000.
582	SNOPCX(I,J)=SNOPCX(I,J)-SMELTK/FDTLIW
583	ENDIF
584	POTEVP(I,J)=POTEVP(I,J)+PLFLX*FDTW
585	! POTFLX(I,J)=POTFLX(I,J)-PLFLX
586	!*** WRF LOWER BOUNDARY CONDITIONS
587	GRDFLX(I,J)=GFLX
588	HFX(I,J)=HFLX
589	QFX(I,J)=ELFLX(I,J)/ELWV
590	SFCEVP(I,J)=SFCEVP(I,J)+QFX(I,J)*DT
591	TSK(I,J)=T1K
592	T2K=T1K-HFX(I,J)/(RHO(I)CPM(I)CHS2(I))
593	TH2(I,J)=T2K(1.E5/SFCPRS)*ROVCP
594	Q2M=Q1K-QFX(I,J)/(RHO(I)*CHS2(I))
595	!!!!!! Q2(I,J)=Q2M
596	!!!!!! Q2(I,J)=Q2K
597	! t2k=th2k/(1.E5/SFCPRS)**ROVCP
598	! QS(I,J)=Q1K
599	!!! QSFC(I,J)=Q1K
600	!*** UPDATE STATE VARIABLES
601	SNOW(I,J)=SNODPK*1000.0
602	SNOWH(I,J)=SNOWH1D*1000.0
603	CANWAT(I,J)=CMCK
604	IF(SNOW(I,J).GT.1.0)THEN
605	! ALB(I,J)=0.01ALBD(IVGTPK,ISEASON)(1.+SCFX(IVGTPK))
606	SNOWC(I,J)=1.0
607	ELSE
608	! ALB(I,J)=0.01*ALBD(IVGTPK,ISEASON)
609	SNOWC(I,J)=0.0
610	ENDIF
611	! update albedo
612	ALB(I,J)=ALB1D
613	! update bottom soil temperature
614	TBOT(I,J)=TBOTK
615
616	DO 80 NS=1,NSOIL
617	SMOIS(I,NS,J)=SMOIS1D(NS)
618	SMLIQ(I,NS,J)=SMLIQ1D(NS)
619	TSLB(I,NS,J)=STEMP1D(NS) ! STEMP
620	80 CONTINUE
621	ENDIF
622	#if 0
623	NOOUT=0
624
625	IF((ITIMESTEP.EQ.1.OR.MOD(ITIMESTEP,1).EQ.0) &
626	.AND. I .EQ.29.AND.J.EQ.23) THEN
627	! print*, 'GLW',GLW(I,J),'GSW',GSW(I,J)
628	print*, 'T2K',T2K,'T1K',T1K,'HFX',HFX(I,J),'RHO',RHO(I),'CPM',CPM(I), &
629	'CHS2',CHS2(I),'soil T',STEMP1D,'soil m', SMOIS1D
630	! print*,'Q2M',Q2M,'Q1K',Q1K,'QFX',QFX(I,J),'RHO',RHO(I),'CHS2',CHS2(I),'latent',ELFLX
631	ENDIF
632
633	IF(NOOUT.EQ.1)GOTO 100
634	! write output to 29
635	IF(ITIMESTEP.EQ.1.AND.I.EQ.1.AND.J.EQ.1) &
636	WRITE (29,*)&
637	'itimestep ',' FK ',' SOLDN ',' SFCPR ', &
638	' SFCTMP ',' Q2K ',' TBOTK ', &
639	' CHK ',' ELFLX ',' HFLX ',' GFLX ', &
640	' RNOF1K ',' RNOF2K ',' T1K ',' CMCK ', &
641	' SMCK1 ',' SMCK2 ',' SMCK3 ',' SMCK4 ', &
642	' STCK1 ',' STCK2 ',' STCK3 ',' STCK4 ', &
643	' SNODPK ',' T2 ', &
644	' Q2 ',' SMSTOT ',' SFCEVP ', ' RAIN'
645	IF((ITIMESTEP.EQ.1.OR.MOD(ITIMESTEP,1).EQ.0) &
646	.AND. I .EQ.29.AND.J.EQ.23) THEN
647	print *,'outputting at itimestep =', itimestep
648	STIMESTEP=FLOAT(itimestep)
649	WRITE (29,1029)STIMESTEP,FK,SOLDN,SFCPRS/100.,SFCTMP,1000.* &
650	Q2K,TBOTK,1000.*CHK,ELFLX(i,j),HFLX,GFLX,SFCRUNOFF(I,J)&
651	,UDRUNOFF(I,J),T1K,CMCK,SMOIS1D,STEMP1D,SNODPK, &
652	! ,UDRUNOFF(I,J),T1K,CMCK,SMOIS1D(3),SMOIS1D(7),SMOIS1D(11),&
653	! SMOIS1D(14),STEMP1D(3), STEMP1D(7),STEMP1D(11), &
654	! STEMP1D(14), SNODPK, &
655	T2K,1000.*Q2M,SMSTOT(I,J),SFCEVP(I,J),PRCP
656	1029 FORMAT (29F10.4)
657	! IF(ITIMESTEP.EQ.0)WRITE (39,*)' P ',' T ', &
658	! ' TH ',' Q ',' U ',' V ', &
659	! ' QC '
660	! WRITE (39,1039)itimestep
661	! DO K=kts,kte
662	! WRITE (39,1039)PRESK,TX(I,K),THX(I,K),1000.*QX(I,K),UX(I,K),&
663	! VX(I,K),1000.*QCX(I,K)
664	1039 FORMAT (7F10.5)
665	! ENDDO
666	ENDIF
667	!
668	#endif
669	100 CONTINUE
670	!
671	!-----------------------------------------------------------------------
672	END SUBROUTINE SURFCE
673	!-----------------------------------------------------------------------
674
675	SUBROUTINE SFLX ( &
676	FFROZP,ICE,DT,ZLVL,NSOIL,SLDPTH, &
677	LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,TH2,Q2,SFCSPD,Q2SAT, &
678	DQSDT2,VEGTYP,SOILTYP,SLOPETYP, &
679	SHDFAC,PTU,TBOT,ALB,SNOALB, &
680	CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM, &
681	ETP,ETA,SHEAT,SSOIL,RUNOFF1,RUNOFF2,Q1,SNOMLT, &
682	SOILW,SOILM,SMCWLT,SMCDRY,SMCREF,SMCMAX,I,J)
683	! ----------------------------------------------------------------------
684	! & ETA,SHEAT, &
685	! ----------------------------------------------------------------------
686	! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
687	! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
688	! MODEL). OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
689	! ----------------------------------------------------------------------
690	! & EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, &
691	! & BETA,ETP,SSOIL, &
692	! & FLX1,FLX2,FLX3, &
693	! & SNOMLT,SNCOVR, &
694	! & RUNOFF1,RUNOFF2,RUNOFF3, &
695	! & RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, &
696	! & SOILW,SOILM, &
697	! & SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT,I,J)
698
699	IMPLICIT NONE
700
701	! ----------------------------------------------------------------------
702	! SUBROUTINE SFLX - VERSION 2.7 - June 2nd 2003
703	! ----------------------------------------------------------------------
704	! SUB-DRIVER FOR "NOAH/OSU LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
705	! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
706	! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
707	! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
708	! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
709	! RADIATION AND PRECIP)
710	! ----------------------------------------------------------------------
711	! SFLX ARGUMENT LIST KEY:
712	! ----------------------------------------------------------------------
713	! C CONFIGURATION INFORMATION
714	! F FORCING DATA
715	! I OTHER (INPUT) FORCING DATA
716	! S SURFACE CHARACTERISTICS
717	! H HISTORY (STATE) VARIABLES
718	! O OUTPUT VARIABLES
719	! D DIAGNOSTIC OUTPUT
720	! ----------------------------------------------------------------------
721	! 1. CONFIGURATION INFORMATION (C):
722	! ----------------------------------------------------------------------
723	! ICE SEA-ICE FLAG (=1: SEA-ICE, =0: LAND)
724	! DT TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
725	! 1800 SECS OR LESS)
726	! ZLVL HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
727	! NSOIL NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
728	! PARAMETER NSOLD SET BELOW)
729	! SLDPTH THE THICKNESS OF EACH SOIL LAYER (M)
730	! ----------------------------------------------------------------------
731	! 2. FORCING DATA (F):
732	! ----------------------------------------------------------------------
733	! LWDN LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
734	! SOLDN SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
735	! SFCPRS PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
736	! PRCP PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
737	! SFCTMP AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
738	! TH2 AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
739	! Q2 MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
740	! ----------------------------------------------------------------------
741	! 3. OTHER FORCING (INPUT) DATA (I):
742	! ----------------------------------------------------------------------
743	! SFCSPD WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
744	! Q2SAT SAT MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
745	! DQSDT2 SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
746	! (KG KG-1 K-1)
747	! ----------------------------------------------------------------------
748	! 4. CANOPY/SOIL CHARACTERISTICS (S):
749	! ----------------------------------------------------------------------
750	! VEGTYP VEGETATION TYPE (INTEGER INDEX)
751	! SOILTYP SOIL TYPE (INTEGER INDEX)
752	! SLOPETYP CLASS OF SFC SLOPE (INTEGER INDEX)
753	! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
754	! (FRACTION= 0.0-1.0)
755	! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
756	! (FRACTION= 0.0-1.0) <= SHDFAC
757	! PTU PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
758	! (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
759	! VEG PARMS)
760	! ALB BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
761	! DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
762	! MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
763	! INCLUDE DIURNAL SUN ANGLE EFFECT)
764	! SNOALB UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
765	! ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
766	! TBOT BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
767	! TEMPERATURE)
768	! ----------------------------------------------------------------------
769	! 5. HISTORY (STATE) VARIABLES (H):
770	! ----------------------------------------------------------------------
771	! CMC CANOPY MOISTURE CONTENT (M)
772	! T1 GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
773	! STC(NSOIL) SOIL TEMP (K)
774	! SMC(NSOIL) TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
775	! SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
776	! NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
777	! SNOWH ACTUAL SNOW DEPTH (M)
778	! SNEQV LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
779	! NOTE: SNOW DENSITY = SNEQV/SNOWH
780	! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
781	! =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
782	! =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
783	! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
784	! (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
785	! IT HAS BEEN MULTIPLIED BY WIND SPEED.
786	! CM SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
787	! CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
788	! MULTIPLIED BY WIND SPEED. CM IS NOT NEEDED IN SFLX
789	! ----------------------------------------------------------------------
790	! 6. OUTPUT (O):
791	! ----------------------------------------------------------------------
792	! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
793	! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL. FOR THIS APPLICATION,
794	! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
795	! NECESSARY. OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
796	! ETA ACTUAL LATENT HEAT FLUX (W M-2: NEGATIVE, IF UP FROM
797	! SURFACE)
798	! SHEAT SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
799	! SURFACE)
800	! ----------------------------------------------------------------------
801	! EC CANOPY WATER EVAPORATION (W M-2)
802	! EDIR DIRECT SOIL EVAPORATION (W M-2)
803	! ET(NSOIL) PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
804	! (W M-2)
805	! ETT TOTAL PLANT TRANSPIRATION (W M-2)
806	! ESNOW SUBLIMATION FROM SNOWPACK (W M-2)
807	! DRIP THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
808	! WATER-HOLDING CAPACITY (M)
809	! DEW DEWFALL (OR FROSTFALL FOR T<273.15) (M)
810	! ----------------------------------------------------------------------
811	! BETA RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
812	! ETP POTENTIAL EVAPORATION (W M-2)
813	! SSOIL SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
814	! ----------------------------------------------------------------------
815	! FLX1 PRECIP-SNOW SFC (W M-2)
816	! FLX2 FREEZING RAIN LATENT HEAT FLUX (W M-2)
817	! FLX3 PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
818	! ----------------------------------------------------------------------
819	! SNOMLT SNOW MELT (M) (WATER EQUIVALENT)
820	! SNCOVR FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
821	! ----------------------------------------------------------------------
822	! RUNOFF1 SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
823	! RUNOFF2 SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
824	! SOIL LAYER
825	! RUNOFF3 NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
826	! FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP
827	! ----------------------------------------------------------------------
828	! RC CANOPY RESISTANCE (S M-1)
829	! PC PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
830	! = ACTUAL TRANSPIRATION
831	! XLAI LEAF AREA INDEX (DIMENSIONLESS)
832	! RSMIN MINIMUM CANOPY RESISTANCE (S M-1)
833	! RCS INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
834	! RCT AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
835	! RCQ ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
836	! RCSOIL SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
837	! ----------------------------------------------------------------------
838	! 7. DIAGNOSTIC OUTPUT (D):
839	! ----------------------------------------------------------------------
840	! SOILW AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
841	! BETWEEN SMCWLT AND SMCMAX)
842	! SOILM TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
843	! ----------------------------------------------------------------------
844	! 8. PARAMETERS (P):
845	! ----------------------------------------------------------------------
846	! SMCWLT WILTING POINT (VOLUMETRIC)
847	! SMCDRY DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
848	! LAYER ENDS (VOLUMETRIC)
849	! SMCREF SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
850	! STRESS (VOLUMETRIC)
851	! SMCMAX POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
852	! (VOLUMETRIC)
853	! NROOT NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
854	! IN SUBROUTINE REDPRM.
855	! ----------------------------------------------------------------------
856	INTEGER NSOLD
857	PARAMETER(NSOLD = 20)
858
859	! ----------------------------------------------------------------------
860	! DECLARATIONS - LOGICAL
861	! ----------------------------------------------------------------------
862	LOGICAL FRZGRA
863	LOGICAL SATURATED
864	LOGICAL SNOWNG
865
866	! ----------------------------------------------------------------------
867	! DECLARATIONS - INTEGER
868	! ----------------------------------------------------------------------
869	INTEGER ICE
870	INTEGER K
871	INTEGER KZ
872	INTEGER NSOIL
873	INTEGER NROOT
874	INTEGER SLOPETYP
875	INTEGER SOILTYP
876	INTEGER VEGTYP
877	INTEGER I
878	INTEGER J
879
880	! ----------------------------------------------------------------------
881	! DECLARATIONS - REAL
882	! ----------------------------------------------------------------------
883	REAL ALBEDO
884	REAL ALB
885	REAL BEXP
886	REAL BETA
887	REAL CFACTR
888	REAL CH
889	REAL CM
890	REAL CMC
891	REAL CMCMAX
892	REAL CP
893	! REAL CSNOW
894	REAL CSOIL
895	REAL CZIL
896	REAL DEW
897	REAL DF1
898	REAL DF1H
899	REAL DF1A
900	REAL DKSAT
901	REAL DT
902	REAL DWSAT
903	REAL DQSDT2
904	REAL DSOIL
905	REAL DTOT
906	REAL DRIP
907	REAL EC
908	REAL EDIR
909	REAL ESNOW
910	REAL ET(NSOIL)
911	REAL ETT
912	REAL FRCSNO
913	REAL FRCSOI
914	REAL EPSCA
915	REAL ETA
916	REAL ETP
917	REAL FDOWN
918	REAL F1
919	REAL FLX1
920	REAL FLX2
921	REAL FLX3
922	REAL FXEXP
923	REAL FRZX
924	REAL SHEAT
925	REAL HS
926	REAL KDT
927	REAL LWDN
928	REAL LVH2O
929	REAL PC
930	REAL PRCP
931	REAL PTU
932	REAL PRCP1
933	REAL PSISAT
934	REAL Q2
935	REAL Q2SAT
936	REAL QUARTZ
937	REAL R
938	REAL RCH
939	REAL REFKDT
940	REAL RR
941	REAL RTDIS(NSOLD)
942	REAL RUNOFF1
943	REAL RUNOFF2
944	REAL RGL
945	REAL RUNOFF3
946	REAL RSMAX
947	REAL RC
948	REAL RSMIN
949	REAL RCQ
950	REAL RCS
951	REAL RCSOIL
952	REAL RCT
953	REAL RSNOW
954	REAL SNDENS
955	REAL SNCOND
956	REAL SSOIL
957	REAL SBETA
958	REAL SFCPRS
959	REAL SFCSPD
960	REAL SFCTMP
961	REAL SHDFAC
962	REAL SHDMIN
963	REAL SH2O(NSOIL)
964	REAL SLDPTH(NSOIL)
965	REAL SMCDRY
966	REAL SMCMAX
967	REAL SMCREF
968	REAL SMCWLT
969	REAL SMC(NSOIL)
970	REAL SNEQV
971	REAL SNCOVR
972	REAL SNOWH
973	REAL SN_NEW
974	REAL SLOPE
975	REAL SNUP
976	REAL SALP
977	REAL SNOALB
978	REAL STC(NSOIL)
979	REAL SNOMLT
980	REAL SOLDN
981	REAL SOILM
982	REAL SOILW
983	REAL SOILWM
984	REAL SOILWW
985	REAL T1
986	REAL T1V
987	REAL T24
988	REAL T2V
989	REAL TBOT
990	REAL TH2
991	REAL TH2V
992	REAL TOPT
993	REAL TFREEZ
994	REAL TSNOW
995	REAL XLAI
996	REAL ZLVL
997	REAL ZBOT
998	REAL Z0
999	REAL ZSOIL(NSOLD)
1000
1001	REAL FFROZP
1002	REAL SOLNET
1003	REAL LSUBS
1004
1005	REAL Q1
1006
1007	! ----------------------------------------------------------------------
1008	! DECLARATIONS - PARAMETERS
1009	! ----------------------------------------------------------------------
1010	PARAMETER(TFREEZ = 273.15)
1011	PARAMETER(LVH2O = 2.501E+6)
1012	PARAMETER(LSUBS = 2.83E+6)
1013	PARAMETER(R = 287.04)
1014	PARAMETER(CP = 1004.5)
1015
1016	! ----------------------------------------------------------------------
1017	! INITIALIZATION
1018	! ----------------------------------------------------------------------
1019	RUNOFF1 = 0.0
1020	RUNOFF2 = 0.0
1021	RUNOFF3 = 0.0
1022	SNOMLT = 0.0
1023
1024	! ----------------------------------------------------------------------
1025	! THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE CASE
1026	! ----------------------------------------------------------------------
1027	IF (ICE .EQ. 1) THEN
1028
1029	! ----------------------------------------------------------------------
1030	! SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
1031	! ----------------------------------------------------------------------
1032	DO KZ = 1,NSOIL
1033	ZSOIL(KZ) = -3.*FLOAT(KZ)/FLOAT(NSOIL)
1034	END DO
1035
1036	ELSE
1037
1038	! ----------------------------------------------------------------------
1039	! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
1040	! EACH SOIL LAYER. NOTE: SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
1041	! GROUND)
1042	! ----------------------------------------------------------------------
1043	ZSOIL(1) = -SLDPTH(1)
1044	DO KZ = 2,NSOIL
1045	ZSOIL(KZ) = -SLDPTH(KZ)+ZSOIL(KZ-1)
1046	END DO
1047
1048	ENDIF
1049
1050	! ----------------------------------------------------------------------
1051	! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
1052	! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
1053	! ----------------------------------------------------------------------
1054	CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP, &
1055	& CFACTR,CMCMAX,RSMAX,TOPT,REFKDT,KDT,SBETA, &
1056	& SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,PSISAT,SLOPE, &
1057	& SNUP,SALP,BEXP,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF, &
1058	& SMCDRY,F1,QUARTZ,FXEXP,RTDIS,SLDPTH,ZSOIL, &
1059	& NROOT,NSOIL,Z0,CZIL,XLAI,CSOIL,PTU)
1060
1061	! ----------------------------------------------------------------------
1062	! INITIALIZE PRECIPITATION LOGICALS.
1063	! ----------------------------------------------------------------------
1064	SNOWNG = .FALSE.
1065	FRZGRA = .FALSE.
1066
1067	! ----------------------------------------------------------------------
1068	! IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP
1069	! ----------------------------------------------------------------------
1070	IF (ICE .EQ. 1) THEN
1071	SNEQV = 0.01
1072	SNOWH = 0.05
1073	SNDENS = SNEQV/SNOWH
1074	ENDIF
1075
1076	! ----------------------------------------------------------------------
1077	! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
1078	! SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
1079	! SUBROUTINE)
1080	! ----------------------------------------------------------------------
1081	IF (SNEQV .EQ. 0.0) THEN
1082	SNDENS = 0.0
1083	SNOWH = 0.0
1084	SNCOND = 1.0
1085	ELSE
1086	SNDENS = SNEQV/SNOWH
1087	SNCOND = CSNOW(SNDENS)
1088	ENDIF
1089
1090	! ----------------------------------------------------------------------
1091	! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
1092	! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
1093	! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
1094	! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
1095	! ----------------------------------------------------------------------
1096	IF (PRCP .GT. 0.0) THEN
1097	! IF (SFCTMP .LE. TFREEZ) THEN
1098	IF (FFROZP .GT. 0.5) THEN
1099	SNOWNG = .TRUE.
1100	ELSE
1101	IF (T1 .LE. TFREEZ) FRZGRA = .TRUE.
1102	ENDIF
1103	ENDIF
1104
1105	! ----------------------------------------------------------------------
1106	! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
1107	! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
1108	! IT TO THE EXISTING SNOWPACK.
1109	! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
1110	! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
1111	! ----------------------------------------------------------------------
1112	IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
1113	SN_NEW = PRCP * DT * 0.001
1114	SNEQV = SNEQV + SN_NEW
1115	PRCP1 = 0.0
1116
1117	! ----------------------------------------------------------------------
1118	! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
1119	! UPDATE SNOW THERMAL CONDUCTIVITY
1120	! ----------------------------------------------------------------------
1121	CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
1122	SNCOND = CSNOW (SNDENS)
1123	ELSE
1124
1125	! ----------------------------------------------------------------------
1126	! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
1127	! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
1128	! ANY CANOPY "DRIP" ADDED TO THIS LATER)
1129	! ----------------------------------------------------------------------
1130	PRCP1 = PRCP
1131
1132	ENDIF
1133
1134	! ----------------------------------------------------------------------
1135	! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
1136	! ----------------------------------------------------------------------
1137	IF (ICE .EQ. 0) THEN
1138
1139	! ----------------------------------------------------------------------
1140	! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
1141	! ----------------------------------------------------------------------
1142	IF (SNEQV .EQ. 0.0) THEN
1143	SNCOVR = 0.0
1144	ALBEDO = ALB
1145
1146	ELSE
1147	! ----------------------------------------------------------------------
1148	! DETERMINE SNOW FRACTIONAL COVERAGE.
1149	! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
1150	! ----------------------------------------------------------------------
1151	CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
1152	! MEK JAN 2006, LIMIT SNOW COVER TO A MAXIMUM FRACTION OF 0.98
1153	SNCOVR = MIN(SNCOVR,0.98)
1154	CALL ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1155	ENDIF
1156
1157	ELSE
1158	! ----------------------------------------------------------------------
1159	! SNOW COVER, ALBEDO OVER SEA-ICE
1160	! ----------------------------------------------------------------------
1161	SNCOVR = 1.0
1162	! changed in version 2.6 on June 2nd 2003
1163	! ALBEDO = 0.60
1164	ALBEDO = 0.65
1165	ENDIF
1166
1167	! ----------------------------------------------------------------------
1168	! THERMAL CONDUCTIVITY FOR SEA-ICE CASE
1169	! ----------------------------------------------------------------------
1170	IF (ICE .EQ. 1) THEN
1171	DF1 = 2.2
1172
1173	ELSE
1174
1175	! ----------------------------------------------------------------------
1176	! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
1177	! CALCULATION OF THE THERMAL DIFFUSIVITY. TREATMENT OF THE
1178	! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
1179	! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
1180	! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
1181	! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
1182	! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
1183	! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
1184	! LIMIT OF VERY THIN SNOWPACK. THIS TREATMENT ALSO ELIMINATES
1185	! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
1186	! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
1187	! ----------------------------------------------------------------------
1188	! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
1189	! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
1190	! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
1191	! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
1192	! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
1193	! ----------------------------------------------------------------------
1194	CALL TDFCND (DF1,SMC(1),QUARTZ,SMCMAX,SH2O(1))
1195
1196	! ----------------------------------------------------------------------
1197	! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
1198	! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
1199	! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
1200	! ----------------------------------------------------------------------
1201	DF1 = DF1 * EXP(SBETA*SHDFAC)
1202	ENDIF
1203
1204	! ----------------------------------------------------------------------
1205	! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
1206	! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
1207	! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
1208	! ----------------------------------------------------------------------
1209	DSOIL = -(0.5 * ZSOIL(1))
1210
1211	IF (SNEQV .EQ. 0.) THEN
1212	SSOIL = DF1 * (T1 - STC(1) ) / DSOIL
1213	ELSE
1214	DTOT = SNOWH + DSOIL
1215	FRCSNO = SNOWH/DTOT
1216	FRCSOI = DSOIL/DTOT
1217	!
1218	! 1. HARMONIC MEAN (SERIES FLOW)
1219	! DF1 = (SNCONDDF1)/(FRCSOISNCOND+FRCSNO*DF1)
1220	DF1H = (SNCONDDF1)/(FRCSOISNCOND+FRCSNO*DF1)
1221	! 2. ARITHMETIC MEAN (PARALLEL FLOW)
1222	! DF1 = FRCSNOSNCOND + FRCSOIDF1
1223	DF1A = FRCSNOSNCOND + FRCSOIDF1
1224	!
1225	! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
1226	! DF1 = (SNCOND*FRCSNO)(DF1**FRCSOI)
1227	! TEST - MBEK, 10 Jan 2002
1228	! weigh DF by snow fraction
1229	! DF1 = DF1HSNCOVR + DF1A(1.0-SNCOVR)
1230	! DF1 = DF1HSNCOVR + DF1(1.0-SNCOVR)
1231	DF1 = DF1ASNCOVR + DF1(1.0-SNCOVR)
1232
1233	! ----------------------------------------------------------------------
1234	! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
1235	! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
1236	! MID-LAYER SOIL TEMPERATURE
1237	! ----------------------------------------------------------------------
1238	SSOIL = DF1 * (T1 - STC(1) ) / DTOT
1239	ENDIF
1240
1241	! MEK -- DEBUG -- AUG 2005
1242	! WRITE(,) 'T1,STC(1),DSOIL=',T1,STC(1),DSOIL
1243	! WRITE(,) 'DF1,SBETA,SHDFAC=',DF1,SBETA,SHDFAC
1244	! WRITE(,) 'SSOIL=',SSOIL
1245
1246	! ----------------------------------------------------------------------
1247	! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
1248	! THE PREVIOUS TIMESTEP.
1249	! ----------------------------------------------------------------------
1250	IF (SNCOVR .GT. 0.) THEN
1251	CALL SNOWZ0 (SNCOVR,Z0)
1252	ENDIF
1253
1254	! ----------------------------------------------------------------------
1255	! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
1256	! HEAT AND MOISTURE.
1257	!
1258	! NOTE !!!
1259	! COMMENT OUT CALL SFCDIF, IF SFCDIF ALREADY CALLED IN CALLING PROGRAM
1260	! (SUCH AS IN COUPLED ATMOSPHERIC MODEL).
1261	!
1262	! NOTE !!!
1263	! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
1264	! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
1265	! (CZIL) ARE SET THERE VIA NAMELIST I/O.
1266	!
1267	! NOTE !!!
1268	! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
1269	! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. HENCE THE CH
1270	! RETURNED FROM SFCDIF HAS UNITS OF M/S. THE IMPORTANT COMPANION
1271	! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
1272	! AIR DENSITY AND PARAMETER "CP". "RCH" IS COMPUTED IN "CALL PENMAN".
1273	! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
1274	!
1275	! NOTE !!!
1276	! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
1277	! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT, BUT CM IS NOT USED HERE.
1278	! ----------------------------------------------------------------------
1279	! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
1280	! SFCDIF AND PENMAN.
1281	! ----------------------------------------------------------------------
1282	T2V = SFCTMP * (1.0 + 0.61 * Q2 )
1283	! ----------------------------------------------------------------------
1284	! COMMENT OUT BELOW 2 LINES IF CALL SFCDIF IS COMMENTED OUT, I.E. IN THE
1285	! COUPLED MODEL.
1286	! ----------------------------------------------------------------------
1287	! T1V = T1 * (1.0 + 0.61 * Q2)
1288	! TH2V = TH2 * (1.0 + 0.61 * Q2)
1289	!
1290	! CALL SFCDIF (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
1291
1292	! ----------------------------------------------------------------------
1293	! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
1294	! PENMAN EP SUBROUTINE THAT FOLLOWS
1295	! ----------------------------------------------------------------------
1296	! FDOWN = SOLDN*(1.0-ALBEDO) + LWDN
1297	FDOWN = SOLNET + LWDN
1298
1299	! ----------------------------------------------------------------------
1300	! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
1301	! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
1302	! CALCULATIONS.
1303	! ----------------------------------------------------------------------
1304	CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
1305	& Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
1306	& DQSDT2,FLX2)
1307
1308	! ----------------------------------------------------------------------
1309	! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
1310	! IF NONZERO GREENNESS FRACTION
1311	! ----------------------------------------------------------------------
1312	IF (SHDFAC .GT. 0.) THEN
1313
1314	! ----------------------------------------------------------------------
1315	! FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
1316	! BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
1317	! ----------------------------------------------------------------------
1318	CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, &
1319	& SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1320	& TOPT,RSMAX,RGL,HS,XLAI, &
1321	& RCS,RCT,RCQ,RCSOIL)
1322
1323	ENDIF
1324
1325	! ----------------------------------------------------------------------
1326	! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
1327	! EXISTS OR NOT:
1328	! ----------------------------------------------------------------------
1329	ESNOW = 0.0
1330	IF (SNEQV .EQ. 0.0) THEN
1331	CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
1332	& SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
1333	& SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL, &
1334	& STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
1335	& SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL, &
1336	& DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
1337	& RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
1338	& QUARTZ,FXEXP,CSOIL, &
1339	& BETA,DRIP,DEW,FLX1,FLX2,FLX3)
1340	ELSE
1341	CALL SNOPAC (ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, &
1342	& SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1343	& SBETA,DF1, &
1344	& Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
1345	& SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS, &
1346	& SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT,SNUP, &
1347	& ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
1348	& RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
1349	& ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
1350	& BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW)
1351	! ESNOW = ETA
1352	ENDIF
1353
1354	! ----------------------------------------------------------------------
1355	! PREPARE SENSIBLE HEAT (H) FOR RETURN TO PARENT MODEL
1356	! ----------------------------------------------------------------------
1357	SHEAT = -(CH * CP * SFCPRS)/(R * T2V) * ( TH2 - T1 )
1358
1359	! ----------------------------------------------------------------------
1360	! CONVERT UNITS AND/OR SIGN OF TOTAL EVAP (ETA), POTENTIAL EVAP (ETP),
1361	! SUBSURFACE HEAT FLUX (S), AND RUNOFFS FOR WHAT PARENT MODEL EXPECTS
1362	! CONVERT ETA FROM KG M-2 S-1 TO W M-2
1363	! ----------------------------------------------------------------------
1364	! ETA = ETA*LVH2O
1365	! ETP = ETP*LVH2O
1366
1367	! ----------------------------------------------------------------------
1368	EDIR = EDIR * LVH2O
1369	EC = EC * LVH2O
1370	DO K=1,4
1371	ET(K) = ET(K) * LVH2O
1372	ENDDO
1373	ETT = ETT * LVH2O
1374	ESNOW = ESNOW * LSUBS
1375	ETP = ETP((1.-SNCOVR)LVH2O + SNCOVR*LSUBS)
1376	IF (ETP .GT. 0.) THEN
1377	ETA = EDIR + EC + ETT + ESNOW
1378	ELSE
1379	ETA = ETP
1380	ENDIF
1381	! ----------------------------------------------------------------------
1382	! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
1383	! ----------------------------------------------------------------------
1384	IF (ETP == 0.0) THEN
1385	BETA = 0.0
1386	ELSE
1387	BETA = ETA/ETP
1388	ENDIF
1389
1390	! ----------------------------------------------------------------------
1391
1392	! ----------------------------------------------------------------------
1393	! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
1394	! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
1395	! SSOIL<0: COOL THE SURFACE (DAY TIME)
1396	! ----------------------------------------------------------------------
1397	SSOIL = -1.0*SSOIL
1398
1399	! ----------------------------------------------------------------------
1400	! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
1401	! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW
1402	! ----------------------------------------------------------------------
1403	RUNOFF3 = RUNOFF3/DT
1404	RUNOFF2 = RUNOFF2+RUNOFF3
1405
1406	! ----------------------------------------------------------------------
1407	! TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE
1408	! SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION
1409	! ----------------------------------------------------------------------
1410	SOILM = -1.0SMC(1)ZSOIL(1)
1411	DO K = 2,NSOIL
1412	SOILM = SOILM+SMC(K)*(ZSOIL(K-1)-ZSOIL(K))
1413	END DO
1414	SOILWM = -1.0(SMCMAX-SMCWLT)ZSOIL(1)
1415	SOILWW = -1.0(SMC(1)-SMCWLT)ZSOIL(1)
1416	DO K = 2,NROOT
1417	SOILWM = SOILWM+(SMCMAX-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1418	SOILWW = SOILWW+(SMC(K)-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1419	END DO
1420	SOILW = SOILWW/SOILWM
1421
1422	! ----------------------------------------------------------------------
1423	! END SUBROUTINE SFLX
1424	! ----------------------------------------------------------------------
1425	END SUBROUTINE SFLX
1426
1427	SUBROUTINE ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1428
1429	IMPLICIT NONE
1430
1431	! ----------------------------------------------------------------------
1432	! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
1433	! ALB SNOWFREE ALBEDO
1434	! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
1435	! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1436	! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1437	! SNCOVR FRACTIONAL SNOW COVER
1438	! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
1439	! TSNOW SNOW SURFACE TEMPERATURE (K)
1440	! ----------------------------------------------------------------------
1441	REAL ALB, SNOALB, SHDFAC, SHDMIN, SNCOVR, ALBEDO, TSNOW
1442
1443	! ----------------------------------------------------------------------
1444	! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
1445	! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
1446	! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
1447	! (1985, JCAM, VOL 24, 402-411)
1448	! ----------------------------------------------------------------------
1449	! changed in version 2.6 on June 2nd 2003
1450	! ALBEDO = ALB + (1.0-(SHDFAC-SHDMIN))SNCOVR(SNOALB-ALB)
1451	ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
1452	IF (ALBEDO .GT. SNOALB) ALBEDO=SNOALB
1453
1454	! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
1455	! IF (TSNOW.LE.263.16) THEN
1456	! ALBEDO=SNOALB
1457	! ELSE
1458	! IF (TSNOW.LT.273.16) THEN
1459	! TM=0.1*(TSNOW-263.16)
1460	! ALBEDO=0.5((0.9-0.2(TM*3))+(0.8-0.16(TM**3)))
1461	! ELSE
1462	! ALBEDO=0.67
1463	! ENDIF
1464	! ENDIF
1465
1466	! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
1467	! IF (TSNOW.LT.273.16) THEN
1468	! ALBEDO=SNOALB-0.008*DT/86400
1469	! ELSE
1470	! ALBEDO=(SNOALB-0.5)EXP(-0.24DT/86400)+0.5
1471	! ENDIF
1472
1473	! ----------------------------------------------------------------------
1474	! END SUBROUTINE ALCALC
1475	! ----------------------------------------------------------------------
1476	END SUBROUTINE ALCALC
1477
1478	SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
1479	& SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1480	& TOPT,RSMAX,RGL,HS,XLAI, &
1481	& RCS,RCT,RCQ,RCSOIL)
1482
1483	IMPLICIT NONE
1484
1485	! ----------------------------------------------------------------------
1486	! SUBROUTINE CANRES
1487	! ----------------------------------------------------------------------
1488	! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
1489	! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
1490	! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
1491	! MOISTURE RATHER THAN TOTAL)
1492	! ----------------------------------------------------------------------
1493	! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
1494	! NOILHAN (1990, BLM)
1495	! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
1496	! AND TABLE 2 OF SEC. 3.1.2
1497	! ----------------------------------------------------------------------
1498	! INPUT:
1499	! SOLAR INCOMING SOLAR RADIATION
1500	! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1501	! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
1502	! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1503	! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1504	! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
1505	! SFCPRS SURFACE PRESSURE
1506	! SMC VOLUMETRIC SOIL MOISTURE
1507	! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
1508	! NSOIL NO. OF SOIL LAYERS
1509	! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
1510	! XLAI LEAF AREA INDEX
1511	! SMCWLT WILTING POINT
1512	! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
1513	! SETS IN)
1514	! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
1515	! SURBOUTINE REDPRM
1516	! OUTPUT:
1517	! PC PLANT COEFFICIENT
1518	! RC CANOPY RESISTANCE
1519	! ----------------------------------------------------------------------
1520	INTEGER NSOLD
1521	PARAMETER(NSOLD = 20)
1522
1523	INTEGER K
1524	INTEGER NROOT
1525	INTEGER NSOIL
1526
1527	REAL CH
1528	REAL CP
1529	REAL DELTA
1530	REAL DQSDT2
1531	REAL FF
1532	REAL GX
1533	REAL HS
1534	REAL P
1535	REAL PART(NSOLD)
1536	REAL PC
1537	REAL Q2
1538	REAL Q2SAT
1539	REAL RC
1540	REAL RSMIN
1541	REAL RCQ
1542	REAL RCS
1543	REAL RCSOIL
1544	REAL RCT
1545	REAL RD
1546	REAL RGL
1547	REAL RR
1548	REAL RSMAX
1549	REAL SFCPRS
1550	REAL SFCTMP
1551	REAL SIGMA
1552	REAL SLV
1553	REAL SMC(NSOIL)
1554	REAL SMCREF
1555	REAL SMCWLT
1556	REAL SOLAR
1557	REAL TOPT
1558	REAL SLVCP
1559	REAL ST1
1560	REAL TAIR4
1561	REAL XLAI
1562	REAL ZSOIL(NSOIL)
1563
1564	PARAMETER(CP = 1004.5)
1565	PARAMETER(RD = 287.04)
1566	PARAMETER(SIGMA = 5.67E-8)
1567	PARAMETER(SLV = 2.501000E6)
1568
1569	! ----------------------------------------------------------------------
1570	! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
1571	! ----------------------------------------------------------------------
1572	RCS = 0.0
1573	RCT = 0.0
1574	RCQ = 0.0
1575	RCSOIL = 0.0
1576	RC = 0.0
1577
1578	! ----------------------------------------------------------------------
1579	! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
1580	! ----------------------------------------------------------------------
1581	FF = 0.552.0SOLAR/(RGL*XLAI)
1582	RCS = (FF + RSMIN/RSMAX) / (1.0 + FF)
1583	RCS = MAX(RCS,0.0001)
1584
1585	! ----------------------------------------------------------------------
1586	! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
1587	! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
1588	! ----------------------------------------------------------------------
1589	RCT = 1.0 - 0.0016((TOPT-SFCTMP)*2.0)
1590	RCT = MAX(RCT,0.0001)
1591
1592	! ----------------------------------------------------------------------
1593	! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
1594	! RCQ EXPRESSION FROM SSIB
1595	! ----------------------------------------------------------------------
1596	RCQ = 1.0/(1.0+HS*(Q2SAT-Q2))
1597	RCQ = MAX(RCQ,0.01)
1598
1599	! ----------------------------------------------------------------------
1600	! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
1601	! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
1602	! ----------------------------------------------------------------------
1603	GX = (SMC(1) - SMCWLT) / (SMCREF - SMCWLT)
1604	IF (GX .GT. 1.) GX = 1.
1605	IF (GX .LT. 0.) GX = 0.
1606
1607	! ----------------------------------------------------------------------
1608	! USE SOIL DEPTH AS WEIGHTING FACTOR
1609	! ----------------------------------------------------------------------
1610	PART(1) = (ZSOIL(1)/ZSOIL(NROOT)) * GX
1611	! ----------------------------------------------------------------------
1612	! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1613	! PART(1) = RTDIS(1) * GX
1614	! ----------------------------------------------------------------------
1615	IF (NROOT .GT. 1) THEN
1616	DO K = 2,NROOT
1617	GX = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT)
1618	IF (GX .GT. 1.) GX = 1.
1619	IF (GX .LT. 0.) GX = 0.
1620	! ----------------------------------------------------------------------
1621	! USE SOIL DEPTH AS WEIGHTING FACTOR
1622	! ----------------------------------------------------------------------
1623	PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT)) * GX
1624	! ----------------------------------------------------------------------
1625	! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1626	! PART(K) = RTDIS(K) * GX
1627	! ----------------------------------------------------------------------
1628	END DO
1629	ENDIF
1630
1631	DO K = 1,NROOT
1632	RCSOIL = RCSOIL+PART(K)
1633	END DO
1634	RCSOIL = MAX(RCSOIL,0.0001)
1635
1636	! ----------------------------------------------------------------------
1637	! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
1638	! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
1639	! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
1640	! PC * LINERIZED PENMAN POTENTIAL EVAP =
1641	! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
1642	! ----------------------------------------------------------------------
1643	RC = RSMIN/(XLAIRCSRCTRCQRCSOIL)
1644
1645	! TAIR4 = SFCTMP**4.
1646	! ST1 = (4.SIGMARD)/CP
1647	! SLVCP = SLV/CP
1648	! RR = ST1TAIR4/(SFCPRSCH) + 1.0
1649	RR = (4.SIGMARD/CP)(SFCTMP4.)/(SFCPRSCH) + 1.0
1650	DELTA = (SLV/CP)*DQSDT2
1651
1652	PC = (RR+DELTA)/(RR(1.+RCCH)+DELTA)
1653
1654	! ----------------------------------------------------------------------
1655	! END SUBROUTINE CANRES
1656	! ----------------------------------------------------------------------
1657	END SUBROUTINE CANRES
1658
1659	SUBROUTINE DEVAP (EDIR1,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
1660	& DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1661
1662	IMPLICIT NONE
1663
1664	! ----------------------------------------------------------------------
1665	! SUBROUTINE DEVAP
1666	! ----------------------------------------------------------------------
1667	! CALCULATE DIRECT SOIL EVAPORATION
1668	! ----------------------------------------------------------------------
1669	REAL BEXP
1670	! REAL DEVAP
1671	REAL EDIR1
1672	REAL DKSAT
1673	REAL DWSAT
1674	REAL ETP1
1675	REAL FX
1676	REAL FXEXP
1677	REAL SHDFAC
1678	REAL SMC
1679	REAL SMCDRY
1680	REAL SMCMAX
1681	REAL ZSOIL
1682	REAL SMCREF
1683	REAL SMCWLT
1684	REAL SRATIO
1685
1686	! ----------------------------------------------------------------------
1687	! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
1688	! WHEN FXEXP=1.
1689	! FX > 1 REPRESENTS DEMAND CONTROL
1690	! FX < 1 REPRESENTS FLUX CONTROL
1691	! ----------------------------------------------------------------------
1692	SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
1693	IF (SRATIO .GT. 0.) THEN
1694	FX = SRATIO**FXEXP
1695	FX = MAX ( MIN ( FX, 1. ) ,0. )
1696	ELSE
1697	FX = 0.
1698	ENDIF
1699
1700	! ----------------------------------------------------------------------
1701	! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
1702	! ----------------------------------------------------------------------
1703	! DEVAP = FX * ( 1.0 - SHDFAC ) * ETP1
1704	EDIR1 = FX * ( 1.0 - SHDFAC ) * ETP1
1705
1706	! ----------------------------------------------------------------------
1707	! END SUBROUTINE DEVAP
1708	! ----------------------------------------------------------------------
1709	END SUBROUTINE DEVAP
1710
1711	SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
1712	& SH2O, &
1713	& SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
1714	& SMCREF,SHDFAC,CMCMAX, &
1715	& SMCDRY,CFACTR, &
1716	& EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1717
1718	IMPLICIT NONE
1719
1720	! ----------------------------------------------------------------------
1721	! SUBROUTINE EVAPO
1722	! ----------------------------------------------------------------------
1723	! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
1724	! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
1725	! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
1726	! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
1727	! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
1728	! ----------------------------------------------------------------------
1729	INTEGER NSOLD
1730	PARAMETER(NSOLD = 20)
1731
1732	INTEGER I
1733	INTEGER K
1734	INTEGER NSOIL
1735	INTEGER NROOT
1736
1737	REAL BEXP
1738	REAL CFACTR
1739	REAL CMC
1740	REAL CMC2MS
1741	REAL CMCMAX
1742	! REAL DEVAP
1743	REAL DKSAT
1744	REAL DT
1745	REAL DWSAT
1746	REAL EC1
1747	REAL EDIR1
1748	REAL ET1(NSOIL)
1749	REAL ETA1
1750	REAL ETP1
1751	REAL ETT1
1752	REAL FXEXP
1753	REAL PC
1754	REAL Q2
1755	REAL RTDIS(NSOIL)
1756	REAL SFCTMP
1757	REAL SHDFAC
1758	REAL SMC(NSOIL)
1759	REAL SH2O(NSOIL)
1760	REAL SMCDRY
1761	REAL SMCMAX
1762	REAL SMCREF
1763	REAL SMCWLT
1764	REAL ZSOIL(NSOIL)
1765
1766	! ----------------------------------------------------------------------
1767	! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
1768	! GREATER THAN ZERO.
1769	! ----------------------------------------------------------------------
1770	EDIR1 = 0.
1771	EC1 = 0.
1772	DO K = 1,NSOIL
1773	ET1(K) = 0.
1774	END DO
1775	ETT1 = 0.
1776
1777	IF (ETP1 .GT. 0.0) THEN
1778
1779	! ----------------------------------------------------------------------
1780	! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
1781	! ONLY IF VEG COVER NOT COMPLETE.
1782	! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
1783	! ----------------------------------------------------------------------
1784	IF (SHDFAC .LT. 1.) THEN
1785	CALL DEVAP (EDIR1,ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1786	! EDIR = DEVAP(ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1787	& BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1788	ENDIF
1789
1790	! ----------------------------------------------------------------------
1791	! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
1792	! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
1793	! ----------------------------------------------------------------------
1794	IF (SHDFAC.GT.0.0) THEN
1795
1796	CALL TRANSP (ET1,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
1797	& CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
1798
1799	DO K = 1,NSOIL
1800	ETT1 = ETT1 + ET1(K)
1801	END DO
1802
1803	! ----------------------------------------------------------------------
1804	! CALCULATE CANOPY EVAPORATION.
1805	! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
1806	! ----------------------------------------------------------------------
1807	IF (CMC .GT. 0.0) THEN
1808	EC1 = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
1809	ELSE
1810	EC1 = 0.0
1811	ENDIF
1812
1813	! ----------------------------------------------------------------------
1814	! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
1815	! CANOPY. -F.CHEN, 18-OCT-1994
1816	! ----------------------------------------------------------------------
1817	CMC2MS = CMC / DT
1818	EC1 = MIN ( CMC2MS, EC1 )
1819	ENDIF
1820	ENDIF
1821
1822	! ----------------------------------------------------------------------
1823	! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
1824	! ----------------------------------------------------------------------
1825	ETA1 = EDIR1 + ETT1 + EC1
1826
1827	! ----------------------------------------------------------------------
1828	! END SUBROUTINE EVAPO
1829	! ----------------------------------------------------------------------
1830	END SUBROUTINE EVAPO
1831
1832	SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
1833	& TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
1834	& F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
1835
1836	IMPLICIT NONE
1837
1838	! ----------------------------------------------------------------------
1839	! SUBROUTINE HRT
1840	! ----------------------------------------------------------------------
1841	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1842	! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
1843	! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1844	! ----------------------------------------------------------------------
1845	INTEGER NSOLD
1846	PARAMETER(NSOLD = 20)
1847
1848	LOGICAL ITAVG
1849
1850	INTEGER I
1851	INTEGER K
1852	INTEGER NSOIL
1853
1854	! ----------------------------------------------------------------------
1855	! DECLARE WORK ARRAYS NEEDED IN TRI-DIAGONAL IMPLICIT SOLVER
1856	! ----------------------------------------------------------------------
1857	REAL AI(NSOLD)
1858	REAL BI(NSOLD)
1859	REAL CI(NSOLD)
1860
1861	! ----------------------------------------------------------------------
1862	! DECLARATIONS
1863	! ----------------------------------------------------------------------
1864	REAL BEXP
1865	REAL CAIR
1866	REAL CH2O
1867	REAL CICE
1868	REAL CSOIL
1869	REAL DDZ
1870	REAL DDZ2
1871	REAL DENOM
1872	REAL DF1
1873	REAL DF1N
1874	REAL DF1K
1875	REAL DT
1876	REAL DTSDZ
1877	REAL DTSDZ2
1878	REAL F1
1879	REAL HCPCT
1880	REAL PSISAT
1881	REAL QUARTZ
1882	REAL QTOT
1883	REAL RHSTS(NSOIL)
1884	REAL SSOIL
1885	REAL SICE
1886	REAL SMC(NSOIL)
1887	REAL SH2O(NSOIL)
1888	REAL SMCMAX
1889	! REAL SNKSRC
1890	REAL STC(NSOIL)
1891	REAL T0
1892	REAL TAVG
1893	REAL TBK
1894	REAL TBK1
1895	REAL TBOT
1896	REAL ZBOT
1897	REAL TSNSR
1898	REAL TSURF
1899	REAL YY
1900	REAL ZSOIL(NSOIL)
1901	REAL ZZ1
1902
1903	PARAMETER(T0 = 273.15)
1904
1905	! ----------------------------------------------------------------------
1906	! SET SPECIFIC HEAT CAPACITIES OF AIR, WATER, ICE, SOIL MINERAL
1907	! ----------------------------------------------------------------------
1908	PARAMETER(CAIR = 1004.0)
1909	PARAMETER(CH2O = 4.2E6)
1910	PARAMETER(CICE = 2.106E6)
1911	! NOTE: CSOIL NOW SET IN ROUTINE REDPRM AND PASSED IN
1912	! PARAMETER(CSOIL = 1.26E6)
1913
1914	! ----------------------------------------------------------------------
1915	! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
1916	! ----------------------------------------------------------------------
1917	ITAVG = .TRUE.
1918	! ITAVG = .FALSE.
1919
1920	! ----------------------------------------------------------------------
1921	! BEGIN SECTION FOR TOP SOIL LAYER
1922	! ----------------------------------------------------------------------
1923	! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
1924	! ----------------------------------------------------------------------
1925	HCPCT = SH2O(1)CH2O + (1.0-SMCMAX)CSOIL + (SMCMAX-SMC(1))*CAIR &
1926	& + ( SMC(1) - SH2O(1) )*CICE
1927
1928	! ----------------------------------------------------------------------
1929	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1930	! ----------------------------------------------------------------------
1931	DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
1932	AI(1) = 0.0
1933	CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
1934	BI(1) = -CI(1) + DF1 / (0.5 * ZSOIL(1) * ZSOIL(1)HCPCTZZ1)
1935
1936	! ----------------------------------------------------------------------
1937	! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
1938	! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
1939	! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
1940	! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
1941	! ----------------------------------------------------------------------
1942	DTSDZ = (STC(1) - STC(2)) / (-0.5 * ZSOIL(2))
1943	SSOIL = DF1 * (STC(1) - YY) / (0.5 * ZSOIL(1) * ZZ1)
1944	RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
1945
1946	! ----------------------------------------------------------------------
1947	! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
1948	! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
1949	! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
1950	! ----------------------------------------------------------------------
1951	QTOT = SSOIL - DF1*DTSDZ
1952
1953	! ----------------------------------------------------------------------
1954	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
1955	! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
1956	! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
1957	! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
1958	! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
1959	! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
1960	! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
1961	! LATER IN FUNCTION SUBROUTINE SNKSRC
1962	! ----------------------------------------------------------------------
1963	IF (ITAVG) THEN
1964	TSURF = (YY + (ZZ1-1) * STC(1)) / ZZ1
1965	CALL TBND (STC(1),STC(2),ZSOIL,ZBOT,1,NSOIL,TBK)
1966	ENDIF
1967
1968	! ----------------------------------------------------------------------
1969	! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
1970	! ----------------------------------------------------------------------
1971	SICE = SMC(1) - SH2O(1)
1972
1973	! ----------------------------------------------------------------------
1974	! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
1975	! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
1976	! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
1977	! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
1978	! ----------------------------------------------------------------------
1979	IF ( (SICE .GT. 0.) .OR. (TSURF .LT. T0) .OR. &
1980	& (STC(1) .LT. T0) .OR. (TBK .LT. T0) ) THEN
1981
1982	IF (ITAVG) THEN
1983	CALL TMPAVG(TAVG,TSURF,STC(1),TBK,ZSOIL,NSOIL,1)
1984	ELSE
1985	TAVG = STC(1)
1986	ENDIF
1987	TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1), &
1988	& ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1989
1990	RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1991	ENDIF
1992
1993	! ----------------------------------------------------------------------
1994	! THIS ENDS SECTION FOR TOP SOIL LAYER.
1995	! ----------------------------------------------------------------------
1996	! INITIALIZE DDZ2
1997	! ----------------------------------------------------------------------
1998	DDZ2 = 0.0
1999
2000	! ----------------------------------------------------------------------
2001	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2002	! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
2003	! ----------------------------------------------------------------------
2004	DF1K = DF1
2005	DO K = 2,NSOIL
2006
2007	! ----------------------------------------------------------------------
2008	! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
2009	! ----------------------------------------------------------------------
2010	HCPCT = SH2O(K)CH2O +(1.0-SMCMAX)CSOIL +(SMCMAX-SMC(K))*CAIR &
2011	& + ( SMC(K) - SH2O(K) )*CICE
2012
2013	IF (K .NE. NSOIL) THEN
2014	! ----------------------------------------------------------------------
2015	! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
2016	! ----------------------------------------------------------------------
2017	! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
2018	! ----------------------------------------------------------------------
2019	CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2020
2021	! ----------------------------------------------------------------------
2022	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
2023	! ----------------------------------------------------------------------
2024	DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2025	DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2026
2027	! ----------------------------------------------------------------------
2028	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
2029	! ----------------------------------------------------------------------
2030	DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2031	CI(K) = -DF1N * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2032
2033	! ----------------------------------------------------------------------
2034	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2035	! TEMP AT BOTTOM OF LAYER.
2036	! ----------------------------------------------------------------------
2037	IF (ITAVG) THEN
2038	CALL TBND (STC(K),STC(K+1),ZSOIL,ZBOT,K,NSOIL,TBK1)
2039	ENDIF
2040	ELSE
2041
2042	! ----------------------------------------------------------------------
2043	! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
2044	! BOTTOM LAYER.
2045	! ----------------------------------------------------------------------
2046	CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2047
2048	! ----------------------------------------------------------------------
2049	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
2050	! ----------------------------------------------------------------------
2051	DENOM = .5 * (ZSOIL(K-1) + ZSOIL(K)) - ZBOT
2052	DTSDZ2 = (STC(K)-TBOT) / DENOM
2053
2054	! ----------------------------------------------------------------------
2055	! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
2056	! ----------------------------------------------------------------------
2057	CI(K) = 0.
2058
2059	! ----------------------------------------------------------------------
2060	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2061	! TEMP AT BOTTOM OF LAST LAYER.
2062	! ----------------------------------------------------------------------
2063	IF (ITAVG) THEN
2064	CALL TBND (STC(K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
2065	ENDIF
2066
2067	ENDIF
2068	! ----------------------------------------------------------------------
2069	! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
2070	! ----------------------------------------------------------------------
2071	! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2072	! ----------------------------------------------------------------------
2073	DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2074	RHSTS(K) = ( DF1N * DTSDZ2 - DF1K * DTSDZ ) / DENOM
2075	QTOT = -1.0DENOMRHSTS(K)
2076	SICE = SMC(K) - SH2O(K)
2077
2078	IF ( (SICE .GT. 0.) .OR. (TBK .LT. T0) .OR. &
2079	& (STC(K) .LT. T0) .OR. (TBK1 .LT. T0) ) THEN
2080
2081	IF (ITAVG) THEN
2082	CALL TMPAVG(TAVG,TBK,STC(K),TBK1,ZSOIL,NSOIL,K)
2083	ELSE
2084	TAVG = STC(K)
2085	ENDIF
2086	TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL, &
2087	& SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2088	RHSTS(K) = RHSTS(K) - TSNSR / DENOM
2089	ENDIF
2090
2091	! ----------------------------------------------------------------------
2092	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2093	! ----------------------------------------------------------------------
2094	AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2095	BI(K) = -(AI(K) + CI(K))
2096
2097	! ----------------------------------------------------------------------
2098	! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
2099	! ----------------------------------------------------------------------
2100	TBK = TBK1
2101	DF1K = DF1N
2102	DTSDZ = DTSDZ2
2103	DDZ = DDZ2
2104	END DO
2105
2106	! ----------------------------------------------------------------------
2107	! END SUBROUTINE HRT
2108	! ----------------------------------------------------------------------
2109	END SUBROUTINE HRT
2110
2111	SUBROUTINE HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2112
2113	IMPLICIT NONE
2114
2115	! ----------------------------------------------------------------------
2116	! SUBROUTINE HRTICE
2117	! ----------------------------------------------------------------------
2118	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2119	! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE PACK. ALSO TO
2120	! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
2121	! OF THE IMPLICIT TIME SCHEME.
2122	! ----------------------------------------------------------------------
2123	INTEGER NSOLD
2124	PARAMETER(NSOLD = 20)
2125
2126	INTEGER K
2127	INTEGER NSOIL
2128
2129	REAL AI(NSOLD)
2130	REAL BI(NSOLD)
2131	REAL CI(NSOLD)
2132
2133	REAL DDZ
2134	REAL DDZ2
2135	REAL DENOM
2136	REAL DF1
2137	REAL DTSDZ
2138	REAL DTSDZ2
2139	REAL HCPCT
2140	REAL RHSTS(NSOIL)
2141	REAL SSOIL
2142	REAL STC(NSOIL)
2143	REAL TBOT
2144	REAL YY
2145	REAL ZBOT
2146	REAL ZSOIL(NSOIL)
2147	REAL ZZ1
2148
2149	DATA TBOT /271.16/
2150
2151	! ----------------------------------------------------------------------
2152	! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
2153	! HCPCT = 1880.0*917.0.
2154	! ----------------------------------------------------------------------
2155	PARAMETER(HCPCT = 1.72396E+6)
2156
2157	! ----------------------------------------------------------------------
2158	! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
2159	! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
2160	! ----------------------------------------------------------------------
2161	! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
2162	! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
2163	! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
2164	! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
2165	! ----------------------------------------------------------------------
2166	ZBOT = ZSOIL(NSOIL)
2167
2168	! ----------------------------------------------------------------------
2169	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2170	! ----------------------------------------------------------------------
2171	DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
2172	AI(1) = 0.0
2173	CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
2174	BI(1) = -CI(1) + DF1/(0.5 * ZSOIL(1) * ZSOIL(1) * HCPCT * ZZ1)
2175
2176	! ----------------------------------------------------------------------
2177	! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
2178	! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
2179	! RHSTS FOR THE TOP SOIL LAYER.
2180	! ----------------------------------------------------------------------
2181	DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) )
2182	SSOIL = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 )
2183	RHSTS(1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL(1) * HCPCT )
2184
2185	! ----------------------------------------------------------------------
2186	! INITIALIZE DDZ2
2187	! ----------------------------------------------------------------------
2188	DDZ2 = 0.0
2189
2190	! ----------------------------------------------------------------------
2191	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2192	! ----------------------------------------------------------------------
2193	DO K = 2,NSOIL
2194	IF (K .NE. NSOIL) THEN
2195
2196	! ----------------------------------------------------------------------
2197	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
2198	! ----------------------------------------------------------------------
2199	DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2200	DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2201
2202	! ----------------------------------------------------------------------
2203	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
2204	! ----------------------------------------------------------------------
2205	DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2206	CI(K) = -DF1 * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2207	ELSE
2208
2209	! ----------------------------------------------------------------------
2210	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
2211	! ----------------------------------------------------------------------
2212	DTSDZ2 = (STC(K)-TBOT)/(.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT)
2213
2214	! ----------------------------------------------------------------------
2215	! SET MATRIX COEF, CI TO ZERO.
2216	! ----------------------------------------------------------------------
2217	CI(K) = 0.
2218	ENDIF
2219
2220	! ----------------------------------------------------------------------
2221	! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2222	! ----------------------------------------------------------------------
2223	DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2224	RHSTS(K) = ( DF1 * DTSDZ2 - DF1 * DTSDZ ) / DENOM
2225
2226	! ----------------------------------------------------------------------
2227	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2228	! ----------------------------------------------------------------------
2229	AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2230	BI(K) = -(AI(K) + CI(K))
2231
2232	! ----------------------------------------------------------------------
2233	! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
2234	! ----------------------------------------------------------------------
2235	DTSDZ = DTSDZ2
2236	DDZ = DDZ2
2237
2238	END DO
2239	! ----------------------------------------------------------------------
2240	! END SUBROUTINE HRTICE
2241	! ----------------------------------------------------------------------
2242	END SUBROUTINE HRTICE
2243
2244	SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
2245
2246	IMPLICIT NONE
2247
2248	! ----------------------------------------------------------------------
2249	! SUBROUTINE HSTEP
2250	! ----------------------------------------------------------------------
2251	! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
2252	! ----------------------------------------------------------------------
2253	INTEGER NSOLD
2254	PARAMETER(NSOLD = 20)
2255
2256	INTEGER K
2257	INTEGER NSOIL
2258
2259	REAL AI(NSOLD)
2260	REAL BI(NSOLD)
2261	REAL CI(NSOLD)
2262	REAL CIin(NSOLD)
2263	REAL DT
2264	REAL RHSTS(NSOIL)
2265	REAL RHSTSin(NSOIL)
2266	REAL STCIN(NSOIL)
2267	REAL STCOUT(NSOIL)
2268
2269	! ----------------------------------------------------------------------
2270	! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
2271	! ----------------------------------------------------------------------
2272	DO K = 1,NSOIL
2273	RHSTS(K) = RHSTS(K) * DT
2274	AI(K) = AI(K) * DT
2275	BI(K) = 1. + BI(K) * DT
2276	CI(K) = CI(K) * DT
2277	END DO
2278
2279	! ----------------------------------------------------------------------
2280	! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
2281	! ----------------------------------------------------------------------
2282	DO K = 1,NSOIL
2283	RHSTSin(K) = RHSTS(K)
2284	END DO
2285	DO K = 1,NSOIL
2286	CIin(K) = CI(K)
2287	END DO
2288
2289	! ----------------------------------------------------------------------
2290	! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
2291	! ----------------------------------------------------------------------
2292	CALL ROSR12(CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
2293
2294	! ----------------------------------------------------------------------
2295	! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
2296	! ----------------------------------------------------------------------
2297	DO K = 1,NSOIL
2298	STCOUT(K) = STCIN(K) + CI(K)
2299	END DO
2300
2301	! ----------------------------------------------------------------------
2302	! END SUBROUTINE HSTEP
2303	! ----------------------------------------------------------------------
2304	END SUBROUTINE HSTEP
2305
2306	SUBROUTINE NOPAC(ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
2307	& SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
2308	& SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL, &
2309	& STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
2310	& SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
2311	& DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
2312	& RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
2313	& QUARTZ,FXEXP,CSOIL, &
2314	& BETA,DRIP,DEW,FLX1,FLX2,FLX3)
2315
2316	IMPLICIT NONE
2317
2318	! ----------------------------------------------------------------------
2319	! SUBROUTINE NOPAC
2320	! ----------------------------------------------------------------------
2321	! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
2322	! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
2323	! PRESENT.
2324	! ----------------------------------------------------------------------
2325	INTEGER ICE
2326	INTEGER NROOT
2327	INTEGER NSOIL
2328
2329	REAL BEXP
2330	REAL BETA
2331	REAL CFACTR
2332	REAL CMC
2333	REAL CMCMAX
2334	REAL CP
2335	REAL CSOIL
2336	REAL DEW
2337	REAL DF1
2338	REAL DKSAT
2339	REAL DRIP
2340	REAL DT
2341	REAL DWSAT
2342	REAL EC
2343	REAL EDIR
2344	REAL EPSCA
2345	REAL ETA
2346	REAL ETA1
2347	REAL ETP
2348	REAL ETP1
2349	REAL ET(NSOIL)
2350	REAL ETT
2351	REAL FDOWN
2352	REAL F1
2353	REAL FXEXP
2354	REAL FLX1
2355	REAL FLX2
2356	REAL FLX3
2357	REAL FRZFACT
2358	REAL KDT
2359	REAL PC
2360	REAL PRCP
2361	REAL PRCP1
2362	REAL PSISAT
2363	REAL Q2
2364	REAL QUARTZ
2365	REAL RCH
2366	REAL RR
2367	REAL RTDIS(NSOIL)
2368	REAL RUNOFF1
2369	REAL RUNOFF2
2370	REAL RUNOFF3
2371	REAL SSOIL
2372	REAL SBETA
2373	REAL SFCTMP
2374	REAL SHDFAC
2375	REAL SH2O(NSOIL)
2376	REAL SIGMA
2377	REAL SLOPE
2378	REAL SMC(NSOIL)
2379	REAL SMCDRY
2380	REAL SMCMAX
2381	REAL SMCREF
2382	REAL SMCWLT
2383	REAL STC(NSOIL)
2384	REAL T1
2385	REAL T24
2386	REAL TBOT
2387	REAL TH2
2388	REAL YY
2389	REAL YYNUM
2390	REAL ZBOT
2391	REAL ZSOIL(NSOIL)
2392	REAL ZZ1
2393
2394	REAL EC1
2395	REAL EDIR1
2396	REAL ET1(NSOIL)
2397	REAL ETT1
2398
2399	INTEGER K
2400
2401	PARAMETER(CP = 1004.5)
2402	PARAMETER(SIGMA = 5.67E-8)
2403
2404	! ----------------------------------------------------------------------
2405	! EXECUTABLE CODE BEGINS HERE:
2406	! CONVERT ETP FROM KG M-2 S-1 TO MS-1 AND INITIALIZE DEW.
2407	! ----------------------------------------------------------------------
2408	PRCP1 = PRCP * 0.001
2409	ETP1 = ETP * 0.001
2410	DEW = 0.0
2411
2412	EDIR = 0.
2413	EDIR1 = 0.
2414	EC = 0.
2415	EC1 = 0.
2416	DO K = 1,NSOIL
2417	ET(K) = 0.
2418	ET1(K) = 0.
2419	END DO
2420	ETT = 0.
2421	ETT1 = 0.
2422
2423	IF (ETP .GT. 0.0) THEN
2424
2425	! ----------------------------------------------------------------------
2426	! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1'.
2427	! ----------------------------------------------------------------------
2428	CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2429	& SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2430	& SMCREF,SHDFAC,CMCMAX, &
2431	& SMCDRY,CFACTR, &
2432	& EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2433	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2434	& SH2O,SLOPE,KDT,FRZFACT, &
2435	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2436	& SHDFAC,CMCMAX, &
2437	& RUNOFF1,RUNOFF2,RUNOFF3, &
2438	& EDIR1,EC1,ET1, &
2439	& DRIP)
2440
2441	! ----------------------------------------------------------------------
2442	! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2443	! ----------------------------------------------------------------------
2444	ETA = ETA1 * 1000.0
2445
2446	! ----------------------------------------------------------------------
2447	! EDIR = EDIR1 * 1000.0
2448	! EC = EC1 * 1000.0
2449	! ETT = ETT1 * 1000.0
2450	! ET(1) = ET1(1) * 1000.0
2451	! ET(2) = ET1(2) * 1000.0
2452	! ET(3) = ET1(3) * 1000.0
2453	! ET(4) = ET1(4) * 1000.0
2454	! ----------------------------------------------------------------------
2455
2456	ELSE
2457
2458	! ----------------------------------------------------------------------
2459	! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
2460	! ETP1 TO ZERO).
2461	! ----------------------------------------------------------------------
2462	DEW = -ETP1
2463	! ETP1 = 0.0
2464
2465	! ----------------------------------------------------------------------
2466	! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
2467	! ----------------------------------------------------------------------
2468	PRCP1 = PRCP1 + DEW
2469	!
2470	! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
2471	! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
2472	! & SMCREF,SHDFAC,CMCMAX,
2473	! & SMCDRY,CFACTR,
2474	! & EDIR1,EC1,ET1,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2475	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2476	& SH2O,SLOPE,KDT,FRZFACT, &
2477	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2478	& SHDFAC,CMCMAX, &
2479	& RUNOFF1,RUNOFF2,RUNOFF3, &
2480	& EDIR1,EC1,ET1, &
2481	& DRIP)
2482
2483	! ----------------------------------------------------------------------
2484	! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
2485	! ----------------------------------------------------------------------
2486	! ETA = ETA1 * 1000.0
2487
2488	! ----------------------------------------------------------------------
2489	! EDIR = EDIR1 * 1000.0
2490	! EC = EC1 * 1000.0
2491	! ETT = ETT1 * 1000.0
2492	! ET(1) = ET1(1) * 1000.0
2493	! ET(2) = ET1(2) * 1000.0
2494	! ET(3) = ET1(3) * 1000.0
2495	! ET(4) = ET1(4) * 1000.0
2496	! ----------------------------------------------------------------------
2497
2498	ENDIF
2499
2500	! ----------------------------------------------------------------------
2501	! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2502	! ----------------------------------------------------------------------
2503	! ETA = ETA1 * 1000.0
2504
2505	! ----------------------------------------------------------------------
2506	EDIR = EDIR1 * 1000.0
2507	EC = EC1 * 1000.0
2508	DO K = 1,NSOIL
2509	ET(K) = ET1(K) * 1000.0
2510	! ET(1) = ET1(1) * 1000.0
2511	! ET(2) = ET1(2) * 1000.0
2512	! ET(3) = ET1(3) * 1000.0
2513	! ET(4) = ET1(4) * 1000.0
2514	ENDDO
2515	ETT = ETT1 * 1000.0
2516	! ----------------------------------------------------------------------
2517
2518	! ----------------------------------------------------------------------
2519	! BASED ON ETP AND E VALUES, DETERMINE BETA
2520	! ----------------------------------------------------------------------
2521	IF ( ETP .LE. 0.0 ) THEN
2522	BETA = 0.0
2523	IF ( ETP .LT. 0.0 ) THEN
2524	BETA = 1.0
2525	ETA = ETP
2526	ENDIF
2527	ELSE
2528	BETA = ETA / ETP
2529	ENDIF
2530
2531	! ----------------------------------------------------------------------
2532	! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
2533	! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
2534	! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
2535	! ----------------------------------------------------------------------
2536	CALL TDFCND (DF1,SMC(1),QUARTZ,SMCMAX,SH2O(1))
2537
2538	! ----------------------------------------------------------------------
2539	! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
2540	! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
2541	! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
2542	! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
2543	! ROUTINE SFLX)
2544	! ----------------------------------------------------------------------
2545	DF1 = DF1 * EXP(SBETA*SHDFAC)
2546
2547	! ----------------------------------------------------------------------
2548	! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
2549	! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
2550	! ----------------------------------------------------------------------
2551	YYNUM = FDOWN - SIGMA * T24
2552	YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-BETA*EPSCA) / RR
2553	ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0
2554
2555	CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2556	& TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2557	& QUARTZ,CSOIL)
2558
2559	! ----------------------------------------------------------------------
2560	! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
2561	! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
2562	! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
2563	! ----------------------------------------------------------------------
2564	FLX1 = 0.0
2565	FLX3 = 0.0
2566
2567	! ----------------------------------------------------------------------
2568	! END SUBROUTINE NOPAC
2569	! ----------------------------------------------------------------------
2570	END SUBROUTINE NOPAC
2571
2572	SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
2573	& Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
2574	& DQSDT2,FLX2)
2575
2576	IMPLICIT NONE
2577
2578	! ----------------------------------------------------------------------
2579	! SUBROUTINE PENMAN
2580	! ----------------------------------------------------------------------
2581	! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
2582	! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
2583	! CALLING ROUTINE FOR LATER USE.
2584	! ----------------------------------------------------------------------
2585	LOGICAL SNOWNG
2586	LOGICAL FRZGRA
2587
2588	REAL A
2589	REAL BETA
2590	REAL CH
2591	REAL CP
2592	REAL CPH2O
2593	REAL CPICE
2594	REAL DELTA
2595	REAL DQSDT2
2596	REAL ELCP
2597	REAL EPSCA
2598	REAL ETP
2599	REAL FDOWN
2600	REAL FLX2
2601	REAL FNET
2602	REAL LSUBC
2603	REAL LSUBF
2604	REAL PRCP
2605	REAL Q2
2606	REAL Q2SAT
2607	REAL R
2608	REAL RAD
2609	REAL RCH
2610	REAL RHO
2611	REAL RR
2612	REAL SSOIL
2613	REAL SFCPRS
2614	REAL SFCTMP
2615	REAL SIGMA
2616	REAL T24
2617	REAL T2V
2618	REAL TH2
2619
2620	PARAMETER(CP = 1004.6)
2621	PARAMETER(CPH2O = 4.218E+3)
2622	PARAMETER(CPICE = 2.106E+3)
2623	PARAMETER(R = 287.04)
2624	PARAMETER(ELCP = 2.4888E+3)
2625	PARAMETER(LSUBF = 3.335E+5)
2626	PARAMETER(LSUBC = 2.501000E+6)
2627	PARAMETER(SIGMA = 5.67E-8)
2628
2629	! ----------------------------------------------------------------------
2630	! EXECUTABLE CODE BEGINS HERE:
2631	! ----------------------------------------------------------------------
2632	FLX2 = 0.0
2633
2634	! ----------------------------------------------------------------------
2635	! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
2636	! ----------------------------------------------------------------------
2637	DELTA = ELCP * DQSDT2
2638	T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
2639	RR = T24 * 6.48E-8 /(SFCPRS * CH) + 1.0
2640	RHO = SFCPRS / (R * T2V)
2641	RCH = RHO * CP * CH
2642
2643	! ----------------------------------------------------------------------
2644	! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
2645	! EFFECTS CAUSED BY FALLING PRECIPITATION.
2646	! ----------------------------------------------------------------------
2647	IF (.NOT. SNOWNG) THEN
2648	IF (PRCP .GT. 0.0) RR = RR + CPH2O*PRCP/RCH
2649	ELSE
2650	RR = RR + CPICE*PRCP/RCH
2651	ENDIF
2652
2653	FNET = FDOWN - SIGMA*T24 - SSOIL
2654
2655	! ----------------------------------------------------------------------
2656	! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
2657	! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
2658	! ----------------------------------------------------------------------
2659	IF (FRZGRA) THEN
2660	FLX2 = -LSUBF * PRCP
2661	FNET = FNET - FLX2
2662	ENDIF
2663
2664	! ----------------------------------------------------------------------
2665	! FINISH PENMAN EQUATION CALCULATIONS.
2666	! ----------------------------------------------------------------------
2667	RAD = FNET/RCH + TH2 - SFCTMP
2668	A = ELCP * (Q2SAT - Q2)
2669	EPSCA = (ARR + RADDELTA) / (DELTA + RR)
2670	ETP = EPSCA * RCH / LSUBC
2671
2672	! ----------------------------------------------------------------------
2673	! END SUBROUTINE PENMAN
2674	! ----------------------------------------------------------------------
2675	END SUBROUTINE PENMAN
2676
2677	SUBROUTINE REDPRM ( &
2678	& VEGTYP,SOILTYP,SLOPETYP, &
2679	& CFACTR,CMCMAX,RSMAX,TOPT,REFKDT,KDT,SBETA, &
2680	& SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,PSISAT,SLOPE, &
2681	& SNUP,SALP,BEXP,DKSAT,DWSAT, &
2682	& SMCMAX,SMCWLT,SMCREF, &
2683	& SMCDRY,F1,QUARTZ,FXEXP,RTDIS,SLDPTH,ZSOIL, &
2684	& NROOT,NSOIL,Z0,CZIL,LAI,CSOIL,PTU)
2685
2686	USE module_wrf_error
2687	IMPLICIT NONE
2688	! ----------------------------------------------------------------------
2689	! SUBROUTINE REDPRM
2690	! ----------------------------------------------------------------------
2691	! INTERNALLY SET (DEFAULT VALUESS), OR OPTIONALLY READ-IN VIA NAMELIST
2692	! I/O, ALL SOIL AND VEGETATION PARAMETERS REQUIRED FOR THE EXECUSION OF
2693	! THE NOAH LSM.
2694	!
2695	! OPTIONAL NON-DEFAULT PARAMETERS CAN BE READ IN, ACCOMMODATING UP TO 30
2696	! SOIL, VEG, OR SLOPE CLASSES, IF THE DEFAULT MAX NUMBER OF SOIL, VEG,
2697	! AND/OR SLOPE TYPES IS RESET.
2698	!
2699	! FUTURE UPGRADES OF ROUTINE REDPRM MUST EXPAND TO INCORPORATE SOME OF
2700	! THE EMPIRICAL PARAMETERS OF THE FROZEN SOIL AND SNOWPACK PHYSICS (SUCH
2701	! AS IN ROUTINES FRH2O, SNOWPACK, AND SNOW_NEW) NOT YET SET IN THIS
2702	! REDPRM ROUTINE, BUT RATHER SET IN LOWER LEVEL SUBROUTINES.
2703	!
2704	! SET MAXIMUM NUMBER OF SOIL-, VEG-, AND SLOPETYP IN DATA STATEMENT.
2705	! ----------------------------------------------------------------------
2706	INTEGER MAX_SLOPETYP
2707	INTEGER MAX_SOILTYP
2708	INTEGER MAX_VEGTYP
2709
2710	PARAMETER(MAX_SLOPETYP = 30)
2711	PARAMETER(MAX_SOILTYP = 30)
2712	PARAMETER(MAX_VEGTYP = 30)
2713
2714	! ----------------------------------------------------------------------
2715	! NUMBER OF DEFINED SOIL-, VEG-, AND SLOPETYPS USED.
2716	! ----------------------------------------------------------------------
2717	INTEGER DEFINED_VEG
2718	INTEGER DEFINED_SOIL
2719	INTEGER DEFINED_SLOPE
2720
2721	DATA DEFINED_VEG/27/
2722	DATA DEFINED_SOIL/19/
2723	DATA DEFINED_SLOPE/9/
2724
2725	! ----------------------------------------------------------------------
2726	! SET-UP SOIL PARAMETERS FOR GIVEN SOIL TYPE
2727	! INPUT: SOLTYP: SOIL TYPE (INTEGER INDEX)
2728	! OUTPUT: SOIL PARAMETERS:
2729	! MAXSMC: MAX SOIL MOISTURE CONTENT (POROSITY)
2730	! REFSMC: REFERENCE SOIL MOISTURE (ONSET OF SOIL MOISTURE
2731	! STRESS IN TRANSPIRATION)
2732	! WLTSMC: WILTING PT SOIL MOISTURE CONTENTS
2733	! DRYSMC: AIR DRY SOIL MOIST CONTENT LIMITS
2734	! SATPSI: SATURATED SOIL POTENTIAL
2735	! SATDK: SATURATED SOIL HYDRAULIC CONDUCTIVITY
2736	! BB: THE 'B' PARAMETER
2737	! SATDW: SATURATED SOIL DIFFUSIVITY
2738	! F11: USED TO COMPUTE SOIL DIFFUSIVITY/CONDUCTIVITY
2739	! QUARTZ: SOIL QUARTZ CONTENT
2740	! ----------------------------------------------------------------------
2741	! SOIL STATSGO
2742	! TYPE CLASS
2743	! ---- -------
2744	! 1 SAND
2745	! 2 LOAMY SAND
2746	! 3 SANDY LOAM
2747	! 4 SILT LOAM
2748	! 5 SILT
2749	! 6 LOAM
2750	! 7 SANDY CLAY LOAM
2751	! 8 SILTY CLAY LOAM
2752	! 9 CLAY LOAM
2753	! 10 SANDY CLAY
2754	! 11 SILTY CLAY
2755	! 12 CLAY
2756	! 13 ORGANIC MATERIAL
2757	! 14 WATER
2758	! 15 BEDROCK
2759	! 16 OTHER(land-ice)
2760	! 17 PLAYA
2761	! 18 LAVA
2762	! 19 WHITE SAND
2763	! ----------------------------------------------------------------------
2764
2765	REAL BB(MAX_SOILTYP)
2766	REAL DRYSMC(MAX_SOILTYP)
2767	REAL F11(MAX_SOILTYP)
2768	REAL MAXSMC(MAX_SOILTYP)
2769	REAL REFSMC(MAX_SOILTYP)
2770	REAL SATPSI(MAX_SOILTYP)
2771	REAL SATDK(MAX_SOILTYP)
2772	REAL SATDW(MAX_SOILTYP)
2773	REAL WLTSMC(MAX_SOILTYP)
2774	REAL QTZ(MAX_SOILTYP)
2775
2776	REAL BEXP
2777	REAL DKSAT
2778	REAL DWSAT
2779	REAL F1
2780	REAL PTU
2781	REAL QUARTZ
2782	REAL REFSMC1
2783	REAL SMCDRY
2784	REAL SMCMAX
2785	REAL SMCREF
2786	REAL SMCWLT
2787	REAL WLTSMC1
2788
2789	! ----------------------------------------------------------------------
2790	! SOIL TEXTURE-RELATED ARRAYS.
2791	! ----------------------------------------------------------------------
2792	DATA MAXSMC/0.395, 0.421, 0.434, 0.476, 0.476, 0.439, &
2793	& 0.404, 0.464, 0.465, 0.406, 0.468, 0.457, &
2794	& 0.464, 0.464, 0.200, 0.421, 0.457, 0.200, &
2795	& 0.395, 0.000, 0.000, 0.000, 0.000, 0.000, &
2796	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2797	! ----------------------------------------------------------------------
2798	DATA SATPSI/0.0350, 0.0363, 0.1413, 0.7586, 0.7586, 0.3548, &
2799	& 0.1349, 0.6166, 0.2630, 0.0977, 0.3236, 0.4677, &
2800	& 0.3548, 0.3548, 0.0350, 0.0363, 0.4677, 0.0350, &
2801	& 0.0350, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, &
2802	& 0.000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000/
2803	! ----------------------------------------------------------------------
2804	DATA SATDK /1.7600E-4, 1.4078E-5, 5.2304E-6, 2.8089E-6, 2.8089E-6,&
2805	& 3.3770E-6, 4.4518E-6, 2.0348E-6, 2.4464E-6, 7.2199E-6,&
2806	& 1.3444E-6, 9.7394E-7, 3.3770E-6, 3.3770E-6, 1.4078E-5,&
2807	& 1.4078E-5, 9.7394E-7, 1.4078E-5, 1.7600E-4, 0.0,&
2808	& 0.0, 0.0, 0.0, 0.0, 0.0,&
2809	& 0.0, 0.0, 0.0, 0.0, 0.0/
2810	! ----------------------------------------------------------------------
2811	DATA BB /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
2812	& 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
2813	& 5.25, 5.25, 4.05, 4.26, 11.55, 4.05, &
2814	& 4.05, 0.00, 0.00, 0.00, 0.00, 0.00, &
2815	& 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2816	! ----------------------------------------------------------------------
2817	DATA QTZ /0.92, 0.82, 0.60, 0.25, 0.10, 0.40, &
2818	& 0.60, 0.10, 0.35, 0.52, 0.10, 0.25, &
2819	& 0.05, 0.05, 0.07, 0.25, 0.60, 0.52, &
2820	& 0.92, 0.00, 0.00, 0.00, 0.00, 0.00, &
2821	& 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2822
2823	! ----------------------------------------------------------------------
2824	! THE FOLLOWING 5 PARAMETERS ARE DERIVED LATER IN REDPRM.F FROM THE SOIL
2825	! DATA, AND ARE JUST GIVEN HERE FOR REFERENCE AND TO FORCE STATIC
2826	! STORAGE ALLOCATION. -DAG LOHMANN, FEB. 2001
2827	! ----------------------------------------------------------------------
2828	! !!!!!!!!!!!!!! The following values in the table are NOT used
2829	! !!!!!!!!!!!!!! and are just given for reference
2830	DATA REFSMC/0.196, 0.248, 0.282, 0.332, 0.332, 0.301, &
2831	& 0.293, 0.368, 0.361, 0.320, 0.388, 0.389, &
2832	& 0.319, 0.000, 0.116, 0.248, 0.389, 0.116, &
2833	& 0.196, 0.000, 0.000, 0.000, 0.000, 0.000, &
2834	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2835	! !!!!!!!!!!!!!! The following values in the table are NOT used
2836	! !!!!!!!!!!!!!! and are just given for reference
2837	DATA WLTSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2838	& 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2839	& 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2840	& 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2841	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2842	! !!!!!!!!!!!!!! The following values in the table are NOT used
2843	! !!!!!!!!!!!!!! and are just given for reference
2844	DATA DRYSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2845	& 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2846	& 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2847	& 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2848	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2849
2850	! !!!!!!!!!!!!!! The following values in the table are NOT used
2851	! !!!!!!!!!!!!!! and are just given for reference
2852	DATA SATDW /0.632E-4, 0.517E-5, 0.807E-5, 0.239E-4, 0.239E-4, &
2853	& 0.143E-4, 0.101E-4, 0.236E-4, 0.113E-4, 0.186E-4, &
2854	& 0.966E-5, 0.115E-4, 0.136E-4, 0.0, 0.998E-5, &
2855	& 0.517E-5, 0.115E-4, 0.998E-5, 0.632E-4, 0.0, &
2856	& 0.0, 0.0, 0.0, 0.0, 0.0, &
2857	& 0.0, 0.0, 0.0, 0.0, 0.0/
2858	! !!!!!!!!!!!!!! The following values in the table are NOT used
2859	! !!!!!!!!!!!!!! and are just given for reference
2860	DATA F11 /-1.090, -1.041, -0.568, 0.162, 0.162, -0.327, &
2861	& -1.535, -1.118, -1.297, -3.211, -1.916, -2.258, &
2862	& -0.201, 0.000, -2.287, -1.041, -2.258, -2.287, &
2863	& -1.090, 0.000, 0.000, 0.000, 0.000, 0.000, &
2864	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2865
2866	! ----------------------------------------------------------------------
2867	! SET-UP VEGETATION PARAMETERS FOR A GIVEN VEGETAION TYPE:
2868	! INPUT: VEGTYP = VEGETATION TYPE (INTEGER INDEX)
2869	! OUPUT: VEGETATION PARAMETERS
2870	! SHDFAC: VEGETATION GREENNESS FRACTION
2871	! RSMIN: MIMIMUM STOMATAL RESISTANCE
2872	! RGL: PARAMETER USED IN SOLAR RAD TERM OF
2873	! CANOPY RESISTANCE FUNCTION
2874	! HS: PARAMETER USED IN VAPOR PRESSURE DEFICIT TERM OF
2875	! CANOPY RESISTANCE FUNCTION
2876	! SNUP: THRESHOLD SNOW DEPTH (IN WATER EQUIVALENT M) THAT
2877	! IMPLIES 100% SNOW COVER
2878	! ----------------------------------------------------------------------
2879	! CLASS USGS-WRF VEGETATION/SURFACE TYPE
2880	! 1 Urban and Built-Up Land
2881	! 2 Dryland Cropland and Pasture
2882	! 3 Irrigated Cropland and Pasture
2883	! 4 Mixed Dryland/Irrigated Cropland and Pasture
2884	! 5 Cropland/Grassland Mosaic
2885	! 6 Cropland/Woodland Mosaic
2886	! 7 Grassland
2887	! 8 Shrubland
2888	! 9 Mixed Shrubland/Grassland
2889	! 10 Savanna
2890	! 11 Deciduous Broadleaf Forest
2891	! 12 Deciduous Needleleaf Forest
2892	! 13 Evergreen Broadleaf Forest
2893	! 14 Evergreen Needleleaf Forest
2894	! 15 Mixed Forest
2895	! 16 Water Bodies
2896	! 17 Herbaceous Wetland
2897	! 18 Wooded Wetland
2898	! 19 Barren or Sparsely Vegetated
2899	! 20 Herbaceous Tundra
2900	! 21 Wooded Tundra
2901	! 22 Mixed Tundra
2902	! 23 Bare Ground Tundra
2903	! 24 Snow or Ice
2904	! 25 Playa
2905	! 26 Lava
2906	! 27 White Sand
2907	! ----------------------------------------------------------------------
2908
2909	INTEGER NROOT
2910	INTEGER NROOT_DATA(MAX_VEGTYP)
2911
2912	REAL FRZFACT
2913	REAL HS
2914	REAL HSTBL(MAX_VEGTYP)
2915	REAL LAI
2916	REAL LAI_DATA(MAX_VEGTYP)
2917	REAL PSISAT
2918	REAL RSMIN
2919	REAL RGL
2920	REAL RGLTBL(MAX_VEGTYP)
2921	REAL RSMTBL(MAX_VEGTYP)
2922	REAL SHDFAC
2923	REAL SNUP
2924	REAL SNUPX(MAX_VEGTYP)
2925	REAL Z0
2926	REAL Z0_DATA(MAX_VEGTYP)
2927
2928	! ----------------------------------------------------------------------
2929	! VEGETATION CLASS-RELATED ARRAYS
2930	! ----------------------------------------------------------------------
2931	! DATA NROOT_DATA /2,3,3,3,3,3,3,3,3,3,
2932	! & 4,4,4,4,4,2,2,2,2,3,
2933	! & 3,3,2,2,2,2,2,0,0,0/
2934	DATA NROOT_DATA /1,3,3,3,3,3,3,3,3,3, &
2935	& 4,4,4,4,4,0,2,2,1,3, &
2936	& 3,3,2,1,1,1,1,0,0,0/
2937	DATA RSMTBL /200.0, 70.0, 70.0, 70.0, 70.0, 70.0, &
2938	& 70.0, 300.0, 170.0, 70.0, 100.0, 150.0, &
2939	& 150.0, 125.0, 125.0, 100.0, 40.0, 100.0, &
2940	& 300.0, 150.0, 150.0, 150.0, 200.0, 200.0, &
2941	& 40.0, 100.0, 300.0, 0.0, 0.0, 0.0/
2942	DATA RGLTBL /100.0, 100.0, 100.0, 100.0, 100.0, 65.0, &
2943	& 100.0, 100.0, 100.0, 65.0, 30.0, 30.0, &
2944	& 30.0, 30.0, 30.0, 30.0, 100.0, 30.0, &
2945	& 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, &
2946	& 100.0, 100.0, 100.0, 0.0, 0.0, 0.0/
2947	DATA HSTBL /42.00, 36.25, 36.25, 36.25, 36.25, 44.14, &
2948	& 36.35, 42.00, 39.18, 54.53, 54.53, 47.35, &
2949	& 41.69, 47.35, 51.93, 51.75, 60.00, 51.93, &
2950	& 42.00, 42.00, 42.00, 42.00, 42.00, 42.00, &
2951	& 36.25, 42.00, 42.00, 0.00, 0.00, 0.00/
2952	DATA SNUPX /0.020, 0.020, 0.020, 0.020, 0.020, 0.020, &
2953	& 0.020, 0.020, 0.020, 0.040, 0.040, 0.040, &
2954	& 0.040, 0.040, 0.040, 0.010, 0.013, 0.020, &
2955	& 0.013, 0.020, 0.020, 0.020, 0.020, 0.013, &
2956	& 0.013, 0.013, 0.013, 0.000, 0.000, 0.000/
2957	DATA Z0_DATA / 1.00, 0.07, 0.07, 0.07, 0.07, 0.15, &
2958	& 0.08, 0.03, 0.05, 0.86, 0.80, 0.85, &
2959	& 2.65, 1.09, 0.80, 0.001, 0.04, 0.05, &
2960	& 0.01, 0.04, 0.06, 0.05, 0.03, 0.001, &
2961	& 0.01, 0.15, 0.01, 0.00, 0.00, 0.00/
2962	DATA LAI_DATA /4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2963	& 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2964	& 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2965	& 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2966	& 4.0, 4.0, 4.0, 0.0, 0.0, 0.0/
2967
2968	! ----------------------------------------------------------------------
2969	! CLASS PARAMETER 'SLOPETYP' WAS INCLUDED TO ESTIMATE LINEAR RESERVOIR
2970	! COEFFICIENT 'SLOPE' TO THE BASEFLOW RUNOFF OUT OF THE BOTTOM LAYER.
2971	! LOWEST CLASS (SLOPETYP=0) MEANS HIGHEST SLOPE PARAMETER = 1.
2972	! DEFINITION OF SLOPETYP FROM 'ZOBLER' SLOPE TYPE:
2973	! SLOPE CLASS PERCENT SLOPE
2974	! 1 0-8
2975	! 2 8-30
2976	! 3 > 30
2977	! 4 0-30
2978	! 5 0-8 & > 30
2979	! 6 8-30 & > 30
2980	! 7 0-8, 8-30, > 30
2981	! 9 GLACIAL ICE
2982	! BLANK OCEAN/SEA
2983	! ----------------------------------------------------------------------
2984	! NOTE:
2985	! CLASS 9 FROM 'ZOBLER' FILE SHOULD BE REPLACED BY 8 AND 'BLANK' 9
2986	! ----------------------------------------------------------------------
2987	REAL SLOPE
2988	REAL SLOPE_DATA(MAX_SLOPETYP)
2989
2990	DATA SLOPE_DATA /0.1, 0.6, 1.0, 0.35, 0.55, 0.8, &
2991	& 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, &
2992	& 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2993	& 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2994	& 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/
2995
2996	! ----------------------------------------------------------------------
2997	! SET NAMELIST FILE NAME
2998	! ----------------------------------------------------------------------
2999	CHARACTER*50 NAMELIST_NAME
3000
3001	! ----------------------------------------------------------------------
3002	! SET UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOIL, VEG, SLOPE TYPE)
3003	! ----------------------------------------------------------------------
3004	INTEGER I
3005	INTEGER NSOIL
3006	INTEGER SLOPETYP
3007	INTEGER SOILTYP
3008	INTEGER VEGTYP
3009
3010	INTEGER BARE
3011	! DATA BARE /11/
3012	DATA BARE /19/
3013
3014	LOGICAL LPARAM
3015	DATA LPARAM /.TRUE./
3016
3017	LOGICAL LFIRST
3018	DATA LFIRST /.TRUE./
3019
3020	! ----------------------------------------------------------------------
3021	! PARAMETER USED TO CALCULATE ROUGHNESS LENGTH OF HEAT.
3022	! ----------------------------------------------------------------------
3023	REAL CZIL
3024	REAL CZIL_DATA
3025	! changed in version 2.6 June 2nd 2003
3026	! DATA CZIL_DATA /0.2/
3027	DATA CZIL_DATA /0.1/
3028
3029	! ----------------------------------------------------------------------
3030	! PARAMETER USED TO CALUCULATE VEGETATION EFFECT ON SOIL HEAT FLUX.
3031	! ----------------------------------------------------------------------
3032	REAL SBETA
3033	REAL SBETA_DATA
3034	DATA SBETA_DATA /-2.0/
3035
3036	! ----------------------------------------------------------------------
3037	! BARE SOIL EVAPORATION EXPONENT USED IN DEVAP.
3038	! ----------------------------------------------------------------------
3039	REAL FXEXP
3040	REAL FXEXP_DATA
3041	DATA FXEXP_DATA /2.0/
3042
3043	! ----------------------------------------------------------------------
3044	! SOIL HEAT CAPACITY [J M-3 K-1]
3045	! ----------------------------------------------------------------------
3046	REAL CSOIL
3047	REAL CSOIL_DATA
3048	! DATA CSOIL_DATA /1.26E+6/
3049	DATA CSOIL_DATA /2.00E+6/
3050
3051	! ----------------------------------------------------------------------
3052	! SPECIFY SNOW DISTRIBUTION SHAPE PARAMETER SALP - SHAPE PARAMETER OF
3053	! DISTRIBUTION FUNCTION OF SNOW COVER. FROM ANDERSON'S DATA (HYDRO-17)
3054	! BEST FIT IS WHEN SALP = 2.6
3055	! ----------------------------------------------------------------------
3056	REAL SALP
3057	REAL SALP_DATA
3058	! changed for version 2.6 June 2nd 2003
3059	! DATA SALP_DATA /2.6/
3060	DATA SALP_DATA /4.0/
3061
3062	! ----------------------------------------------------------------------
3063	! KDT IS DEFINED BY REFERENCE REFKDT AND DKSAT; REFDK=2.E-6 IS THE SAT.
3064	! DK. VALUE FOR THE SOIL TYPE 2
3065	! ----------------------------------------------------------------------
3066	REAL REFDK
3067	REAL REFDK_DATA
3068	DATA REFDK_DATA /2.0E-6/
3069
3070	REAL REFKDT
3071	REAL REFKDT_DATA
3072	DATA REFKDT_DATA /3.0/
3073
3074	REAL FRZX
3075	REAL KDT
3076
3077	! ----------------------------------------------------------------------
3078	! FROZEN GROUND PARAMETER, FRZK, DEFINITION: ICE CONTENT THRESHOLD ABOVE
3079	! WHICH FROZEN SOIL IS IMPERMEABLE REFERENCE VALUE OF THIS PARAMETER FOR
3080	! THE LIGHT CLAY SOIL (TYPE=3) FRZK = 0.15 M.
3081	! ----------------------------------------------------------------------
3082	REAL FRZK
3083	REAL FRZK_DATA
3084	DATA FRZK_DATA /0.15/
3085
3086	REAL RTDIS(NSOIL)
3087	REAL SLDPTH(NSOIL)
3088	REAL ZSOIL(NSOIL)
3089
3090	! ----------------------------------------------------------------------
3091	! SET TWO CANOPY WATER PARAMETERS.
3092	! ----------------------------------------------------------------------
3093	REAL CFACTR
3094	REAL CFACTR_DATA
3095	DATA CFACTR_DATA /0.5/
3096
3097	REAL CMCMAX
3098	REAL CMCMAX_DATA
3099	DATA CMCMAX_DATA /0.5E-3/
3100
3101	! ----------------------------------------------------------------------
3102	! SET MAX. STOMATAL RESISTANCE.
3103	! ----------------------------------------------------------------------
3104	REAL RSMAX
3105	REAL RSMAX_DATA
3106	DATA RSMAX_DATA /5000.0/
3107
3108	! ----------------------------------------------------------------------
3109	! SET OPTIMUM TRANSPIRATION AIR TEMPERATURE.
3110	! ----------------------------------------------------------------------
3111	REAL TOPT
3112	REAL TOPT_DATA
3113	DATA TOPT_DATA /298.0/
3114
3115	! ----------------------------------------------------------------------
3116	! SPECIFY DEPTH[M] OF LOWER BOUNDARY SOIL TEMPERATURE.
3117	! ----------------------------------------------------------------------
3118	REAL ZBOT
3119	REAL ZBOT_DATA
3120	! changed for version 2.5.2
3121	! DATA ZBOT_DATA /-3.0/
3122	DATA ZBOT_DATA /-8.0/
3123
3124	! ----------------------------------------------------------------------
3125	! SET TWO SOIL MOISTURE WILT, SOIL MOISTURE REFERENCE PARAMETERS
3126	! ----------------------------------------------------------------------
3127	REAL SMLOW
3128	REAL SMLOW_DATA
3129	DATA SMLOW_DATA /0.5/
3130
3131	REAL SMHIGH
3132	REAL SMHIGH_DATA
3133	! changed in 2.6 from 3 to 6 on June 2nd 2003
3134	DATA SMHIGH_DATA /3.0/
3135	! DATA SMHIGH_DATA /6.0/
3136
3137	! ----------------------------------------------------------------------
3138	! NAMELIST DEFINITION:
3139	! ----------------------------------------------------------------------
3140	NAMELIST /SOIL_VEG/ SLOPE_DATA, RSMTBL, RGLTBL, HSTBL, SNUPX, &
3141	& BB, DRYSMC, F11, MAXSMC, REFSMC, SATPSI, SATDK, SATDW, &
3142	& WLTSMC, QTZ, LPARAM, ZBOT_DATA, SALP_DATA, CFACTR_DATA, &
3143	& CMCMAX_DATA, SBETA_DATA, RSMAX_DATA, TOPT_DATA, &
3144	& REFDK_DATA, FRZK_DATA, BARE, DEFINED_VEG, DEFINED_SOIL, &
3145	& DEFINED_SLOPE, FXEXP_DATA, NROOT_DATA, REFKDT_DATA, Z0_DATA, &
3146	& CZIL_DATA, LAI_DATA, CSOIL_DATA
3147
3148	! ----------------------------------------------------------------------
3149	! READ NAMELIST FILE TO OVERRIDE DEFAULT PARAMETERS ONLY ONCE.
3150	! NAMELIST_NAME must be 50 characters or less.
3151	! ----------------------------------------------------------------------
3152	IF (LFIRST) THEN
3153	! WRITE(,) 'READ NAMELIST'
3154	! OPEN(58, FILE = 'namelist_filename.txt')
3155	! READ(58,'(A)') NAMELIST_NAME
3156	! CLOSE(58)
3157	! WRITE(,) 'Namelist Filename is ', NAMELIST_NAME
3158	! OPEN(59, FILE = NAMELIST_NAME)
3159	! 50 CONTINUE
3160	! READ(59, SOIL_VEG, END=100)
3161	! IF (LPARAM) GOTO 50
3162	! 100 CONTINUE
3163	! CLOSE(59)
3164	! WRITE(*,NML=SOIL_VEG)
3165	LFIRST = .FALSE.
3166	IF (DEFINED_SOIL .GT. MAX_SOILTYP) THEN
3167	WRITE(wrf_err_message,*) 'Warning: DEFINED_SOIL too large in namelist'
3168	CALL wrf_error_fatal( wrf_err_message )
3169	! STOP 222
3170	ENDIF
3171	IF (DEFINED_VEG .GT. MAX_VEGTYP) THEN
3172	WRITE(wrf_err_message,*) 'Warning: DEFINED_VEG too large in namelist'
3173	CALL wrf_error_fatal( wrf_err_message )
3174	! STOP 222
3175	ENDIF
3176	IF (DEFINED_SLOPE .GT. MAX_SLOPETYP) THEN
3177	WRITE(wrf_err_message,*) 'Warning: DEFINED_SLOPE too large in namelist'
3178	CALL wrf_error_fatal( wrf_err_message )
3179	! STOP 222
3180	ENDIF
3181
3182	SMLOW = SMLOW_DATA
3183	SMHIGH = SMHIGH_DATA
3184
3185	DO I = 1,DEFINED_SOIL
3186	SATDW(I) = BB(I)SATDK(I)(SATPSI(I)/MAXSMC(I))
3187	F11(I) = ALOG10(SATPSI(I)) + BB(I)*ALOG10(MAXSMC(I)) + 2.0
3188	REFSMC1 = MAXSMC(I)*(5.79E-9/SATDK(I)) &
3189	& *(1.0/(2.0BB(I)+3.0))
3190	REFSMC(I) = REFSMC1 + (MAXSMC(I)-REFSMC1) / SMHIGH
3191	WLTSMC1 = MAXSMC(I) * (200.0/SATPSI(I))**(-1.0/BB(I))
3192	WLTSMC(I) = WLTSMC1 - SMLOW * WLTSMC1
3193
3194	! ----------------------------------------------------------------------
3195	! CURRENT VERSION DRYSMC VALUES THAT EQUATE TO WLTSMC.
3196	! FUTURE VERSION COULD LET DRYSMC BE INDEPENDENTLY SET VIA NAMELIST.
3197	! ----------------------------------------------------------------------
3198	DRYSMC(I) = WLTSMC(I)
3199	END DO
3200
3201	! ----------------------------------------------------------------------
3202	! END LFIRST BLOCK
3203	! ----------------------------------------------------------------------
3204	ENDIF
3205
3206	IF (SOILTYP .GT. DEFINED_SOIL) THEN
3207	WRITE(wrf_err_message,*) 'Warning: too many soil types'
3208	CALL wrf_error_fatal( wrf_err_message )
3209	! STOP 333
3210	ENDIF
3211	IF (VEGTYP .GT. DEFINED_VEG) THEN
3212	WRITE(wrf_err_message,*) 'Warning: too many veg types'
3213	! STOP 333
3214	ENDIF
3215	IF (SLOPETYP .GT. DEFINED_SLOPE) THEN
3216	WRITE(wrf_err_message,*) 'Warning: too many slope types'
3217	! STOP 333
3218	ENDIF
3219
3220	! ----------------------------------------------------------------------
3221	! SET-UP UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOILTYP, VEGTYP OR
3222	! SLOPETYP)
3223	! ----------------------------------------------------------------------
3224	ZBOT = ZBOT_DATA
3225	SALP = SALP_DATA
3226	CFACTR = CFACTR_DATA
3227	CMCMAX = CMCMAX_DATA
3228	SBETA = SBETA_DATA
3229	RSMAX = RSMAX_DATA
3230	TOPT = TOPT_DATA
3231	REFDK = REFDK_DATA
3232	FRZK = FRZK_DATA
3233	FXEXP = FXEXP_DATA
3234	REFKDT = REFKDT_DATA
3235	CZIL = CZIL_DATA
3236	CSOIL = CSOIL_DATA
3237
3238	! ----------------------------------------------------------------------
3239	! SET-UP SOIL PARAMETERS
3240	! ----------------------------------------------------------------------
3241	BEXP = BB(SOILTYP)
3242	DKSAT = SATDK(SOILTYP)
3243	DWSAT = SATDW(SOILTYP)
3244	F1 = F11(SOILTYP)
3245	KDT = REFKDT * DKSAT/REFDK
3246	PSISAT = SATPSI(SOILTYP)
3247	QUARTZ = QTZ(SOILTYP)
3248	SMCDRY = DRYSMC(SOILTYP)
3249	SMCMAX = MAXSMC(SOILTYP)
3250	SMCREF = REFSMC(SOILTYP)
3251	SMCWLT = WLTSMC(SOILTYP)
3252	FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
3253
3254	! ----------------------------------------------------------------------
3255	! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
3256	! ----------------------------------------------------------------------
3257	FRZX = FRZK * FRZFACT
3258
3259	! ----------------------------------------------------------------------
3260	! SET-UP VEGETATION PARAMETERS
3261	! ----------------------------------------------------------------------
3262	NROOT = NROOT_DATA(VEGTYP)
3263	SNUP = SNUPX(VEGTYP)
3264	RSMIN = RSMTBL(VEGTYP)
3265	RGL = RGLTBL(VEGTYP)
3266	HS = HSTBL(VEGTYP)
3267	Z0 = Z0_DATA(VEGTYP)
3268	LAI = LAI_DATA(VEGTYP)
3269	IF (VEGTYP .EQ. BARE) SHDFAC = 0.0
3270
3271	IF (NROOT .GT. NSOIL) THEN
3272	WRITE(wrf_err_message,*) 'Warning: too many root layers'
3273	CALL wrf_error_fatal(wrf_err_message)
3274	! STOP 333
3275	ENDIF
3276
3277	! ----------------------------------------------------------------------
3278	! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
3279	! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
3280	! ----------------------------------------------------------------------
3281	DO I = 1,NROOT
3282	RTDIS(I) = -SLDPTH(I)/ZSOIL(NROOT)
3283	END DO
3284
3285	! ----------------------------------------------------------------------
3286	! SET-UP SLOPE PARAMETER
3287	! ----------------------------------------------------------------------
3288	SLOPE = SLOPE_DATA(SLOPETYP)
3289
3290	! ----------------------------------------------------------------------
3291	! END SUBROUTINE REDPRM
3292	! ----------------------------------------------------------------------
3293	END SUBROUTINE REDPRM
3294
3295	SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3296
3297	IMPLICIT NONE
3298
3299	! ----------------------------------------------------------------------
3300	! SUBROUTINE ROSR12
3301	! ----------------------------------------------------------------------
3302	! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3303	! ### ### ### ### ### ###
3304	! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3305	! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3306	! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3307	! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3308	! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3309	! # . . # # . # = # . #
3310	! # . . # # . # # . #
3311	! # . . # # . # # . #
3312	! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3313	! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3314	! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3315	! ### ### ### ### ### ###
3316	! ----------------------------------------------------------------------
3317	INTEGER K
3318	INTEGER KK
3319	INTEGER NSOIL
3320
3321	REAL A(NSOIL)
3322	REAL B(NSOIL)
3323	REAL C(NSOIL)
3324	REAL D(NSOIL)
3325	REAL DELTA(NSOIL)
3326	REAL P(NSOIL)
3327
3328	! ----------------------------------------------------------------------
3329	! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3330	! ----------------------------------------------------------------------
3331	C(NSOIL) = 0.0
3332
3333	! ----------------------------------------------------------------------
3334	! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
3335	! ----------------------------------------------------------------------
3336	P(1) = -C(1) / B(1)
3337	DELTA(1) = D(1) / B(1)
3338
3339	! ----------------------------------------------------------------------
3340	! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
3341	! ----------------------------------------------------------------------
3342	DO K = 2,NSOIL
3343	P(K) = -C(K) * ( 1.0 / (B(K) + A (K) * P(K-1)) )
3344	DELTA(K) = (D(K)-A(K)DELTA(K-1))(1.0/(B(K)+A(K)*P(K-1)))
3345	END DO
3346
3347	! ----------------------------------------------------------------------
3348	! SET P TO DELTA FOR LOWEST SOIL LAYER
3349	! ----------------------------------------------------------------------
3350	P(NSOIL) = DELTA(NSOIL)
3351
3352	! ----------------------------------------------------------------------
3353	! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3354	! ----------------------------------------------------------------------
3355	DO K = 2,NSOIL
3356	KK = NSOIL - K + 1
3357	P(KK) = P(KK) * P(KK+1) + DELTA(KK)
3358	END DO
3359
3360	! ----------------------------------------------------------------------
3361	! END SUBROUTINE ROSR12
3362	! ----------------------------------------------------------------------
3363	END SUBROUTINE ROSR12
3364
3365	SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3366	& TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3367	& QUARTZ,CSOIL)
3368
3369	IMPLICIT NONE
3370
3371	! ----------------------------------------------------------------------
3372	! SUBROUTINE SHFLX
3373	! ----------------------------------------------------------------------
3374	! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3375	! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3376	! ON THE TEMPERATURE.
3377	! ----------------------------------------------------------------------
3378	INTEGER NSOLD
3379	PARAMETER(NSOLD = 20)
3380
3381	INTEGER I
3382	INTEGER ICE
3383	INTEGER IFRZ
3384	INTEGER NSOIL
3385
3386	REAL AI(NSOLD)
3387	REAL BI(NSOLD)
3388	REAL CI(NSOLD)
3389
3390	REAL BEXP
3391	REAL CSOIL
3392	REAL DF1
3393	REAL DT
3394	REAL F1
3395	REAL PSISAT
3396	REAL QUARTZ
3397	REAL RHSTS(NSOLD)
3398	REAL SSOIL
3399	REAL SH2O(NSOIL)
3400	REAL SMC(NSOIL)
3401	REAL SMCMAX
3402	REAL SMCWLT
3403	REAL STC(NSOIL)
3404	REAL STCF(NSOLD)
3405	REAL T0
3406	REAL T1
3407	REAL TBOT
3408	REAL YY
3409	REAL ZBOT
3410	REAL ZSOIL(NSOIL)
3411	REAL ZZ1
3412
3413	PARAMETER(T0 = 273.15)
3414
3415	! ----------------------------------------------------------------------
3416	! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3417	! ----------------------------------------------------------------------
3418	IF (ICE.EQ.1) THEN
3419
3420	! ----------------------------------------------------------------------
3421	! SEA-ICE CASE
3422	! ----------------------------------------------------------------------
3423	CALL HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
3424
3425	CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3426
3427	ELSE
3428
3429	! ----------------------------------------------------------------------
3430	! LAND-MASS CASE
3431	! ----------------------------------------------------------------------
3432	CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3433	& ZBOT,PSISAT,SH2O,DT, &
3434	& BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
3435
3436	CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3437
3438	ENDIF
3439
3440	DO I = 1,NSOIL
3441	STC(I) = STCF(I)
3442	END DO
3443
3444	! ----------------------------------------------------------------------
3445	! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
3446	! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
3447	! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
3448	! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
3449	! DIFFERENTLY IN ROUTINE SNOPAC)
3450	! ----------------------------------------------------------------------
3451	T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1
3452
3453	! ----------------------------------------------------------------------
3454	! CALCULATE SURFACE SOIL HEAT FLUX
3455	! ----------------------------------------------------------------------
3456	SSOIL = DF1 * (STC(1) - T1) / (0.5 * ZSOIL(1))
3457
3458	! ----------------------------------------------------------------------
3459	! END SUBROUTINE SHFLX
3460	! ----------------------------------------------------------------------
3461	END SUBROUTINE SHFLX
3462
3463	SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3464	& SH2O,SLOPE,KDT,FRZFACT, &
3465	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3466	& SHDFAC,CMCMAX, &
3467	& RUNOFF1,RUNOFF2,RUNOFF3, &
3468	& EDIR1,EC1,ET1, &
3469	& DRIP)
3470
3471	IMPLICIT NONE
3472
3473	! ----------------------------------------------------------------------
3474	! SUBROUTINE SMFLX
3475	! ----------------------------------------------------------------------
3476	! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
3477	! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
3478	! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3479	! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
3480	! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
3481	! ----------------------------------------------------------------------
3482	INTEGER NSOLD
3483	PARAMETER(NSOLD = 20)
3484
3485	INTEGER I
3486	INTEGER K
3487	INTEGER NSOIL
3488
3489	REAL AI(NSOLD)
3490	REAL BI(NSOLD)
3491	REAL CI(NSOLD)
3492
3493	REAL BEXP
3494	REAL CMC
3495	REAL CMCMAX
3496	REAL DKSAT
3497	REAL DRIP
3498	REAL DT
3499	REAL DUMMY
3500	REAL DWSAT
3501	REAL EC1
3502	REAL EDIR1
3503	REAL ET1(NSOIL)
3504	REAL EXCESS
3505	REAL FRZFACT
3506	REAL KDT
3507	REAL PCPDRP
3508	REAL PRCP1
3509	REAL RHSCT
3510	REAL RHSTT(NSOLD)
3511	REAL RUNOFF1
3512	REAL RUNOFF2
3513	REAL RUNOFF3
3514	REAL SHDFAC
3515	REAL SMC(NSOIL)
3516	REAL SH2O(NSOIL)
3517	REAL SICE(NSOLD)
3518	REAL SH2OA(NSOLD)
3519	REAL SH2OFG(NSOLD)
3520	REAL SLOPE
3521	REAL SMCMAX
3522	REAL SMCWLT
3523	REAL TRHSCT
3524	REAL ZSOIL(NSOIL)
3525
3526	! ----------------------------------------------------------------------
3527	! EXECUTABLE CODE BEGINS HERE.
3528	! ----------------------------------------------------------------------
3529	DUMMY = 0.
3530
3531	! ----------------------------------------------------------------------
3532	! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
3533	! ----------------------------------------------------------------------
3534	RHSCT = SHDFAC * PRCP1 - EC1
3535
3536	! ----------------------------------------------------------------------
3537	! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
3538	! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
3539	! FALL TO THE GRND.
3540	! ----------------------------------------------------------------------
3541	DRIP = 0.
3542	TRHSCT = DT * RHSCT
3543	EXCESS = CMC + TRHSCT
3544	IF (EXCESS .GT. CMCMAX) DRIP = EXCESS - CMCMAX
3545
3546	! ----------------------------------------------------------------------
3547	! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
3548	! SOIL
3549	! ----------------------------------------------------------------------
3550	PCPDRP = (1. - SHDFAC) * PRCP1 + DRIP / DT
3551
3552	! ----------------------------------------------------------------------
3553	! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT & SSTEP
3554	! ----------------------------------------------------------------------
3555	DO I = 1,NSOIL
3556	SICE(I) = SMC(I) - SH2O(I)
3557	END DO
3558
3559	! ----------------------------------------------------------------------
3560	! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3561	! TENDENCY EQUATIONS.
3562	!
3563	! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
3564	! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
3565	! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
3566	! THE FIRST SOIL LAYER)
3567	! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
3568	! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
3569	! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
3570	! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
3571	! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
3572	! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
3573	! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
3574	! SOIL MOISTURE STATE
3575	! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
3576	! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
3577	! OF SECTION 2 OF KALNAY AND KANAMITSU
3578	! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
3579	! ----------------------------------------------------------------------
3580	! IF ( PCPDRP .GT. 0.0 ) THEN
3581	IF ( (PCPDRPDT) .GT. (0.0011000.0(-ZSOIL(1))SMCMAX) ) THEN
3582
3583	! ----------------------------------------------------------------------
3584	! FROZEN GROUND VERSION:
3585	! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
3586	! INCLUDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
3587	! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
3588	! ----------------------------------------------------------------------
3589	CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3590	& DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3591	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3592
3593	CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3594	& CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3595
3596	DO K = 1,NSOIL
3597	SH2OA(K) = (SH2O(K) + SH2OFG(K)) * 0.5
3598	END DO
3599
3600	CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
3601	& DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3602	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3603
3604	CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3605	& CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3606
3607	ELSE
3608
3609	CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3610	& DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3611	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3612
3613	CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3614	& CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3615
3616	ENDIF
3617
3618	! RUNOF = RUNOFF
3619
3620	! ----------------------------------------------------------------------
3621	! END SUBROUTINE SMFLX
3622	! ----------------------------------------------------------------------
3623	END SUBROUTINE SMFLX
3624
3625	SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
3626
3627	IMPLICIT NONE
3628
3629	! ----------------------------------------------------------------------
3630	! SUBROUTINE SNFRAC
3631	! ----------------------------------------------------------------------
3632	! CALCULATE SNOW FRACTION (0 -> 1)
3633	! SNEQV SNOW WATER EQUIVALENT (M)
3634	! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
3635	! SALP TUNING PARAMETER
3636	! SNCOVR FRACTIONAL SNOW COVER
3637	! ----------------------------------------------------------------------
3638	REAL SNEQV, SNUP, SALP, SNCOVR, RSNOW, Z0N, SNOWH
3639
3640	! ----------------------------------------------------------------------
3641	! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
3642	! REDPRM) ABOVE WHICH SNOCVR=1.
3643	! ----------------------------------------------------------------------
3644	IF (SNEQV .LT. SNUP) THEN
3645	RSNOW = SNEQV/SNUP
3646	SNCOVR = 1. - ( EXP(-SALPRSNOW) - RSNOWEXP(-SALP))
3647	ELSE
3648	SNCOVR = 1.0
3649	ENDIF
3650
3651	Z0N=0.035
3652	! FORMULATION OF DICKINSON ET AL. 1986
3653
3654	! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
3655
3656	! FORMULATION OF MARSHALL ET AL. 1994
3657	! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
3658
3659	! ----------------------------------------------------------------------
3660	! END SUBROUTINE SNFRAC
3661	! ----------------------------------------------------------------------
3662	END SUBROUTINE SNFRAC
3663
3664	SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, &
3665	& SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
3666	& SBETA,DF1, &
3667	& Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
3668	& SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
3669	& SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT,SNUP, &
3670	& ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
3671	& RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
3672	& ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
3673	& BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW)
3674
3675	IMPLICIT NONE
3676
3677	! ----------------------------------------------------------------------
3678	! SUBROUTINE SNOPAC
3679	! ----------------------------------------------------------------------
3680	! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
3681	! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
3682	! PRESENT.
3683	! ----------------------------------------------------------------------
3684	INTEGER ICE
3685	INTEGER NROOT
3686	INTEGER NSOIL
3687
3688	LOGICAL SNOWNG
3689
3690	REAL BEXP
3691	REAL BETA
3692	REAL CFACTR
3693	REAL CMC
3694	REAL CMCMAX
3695	REAL CP
3696	REAL CPH2O
3697	REAL CPICE
3698	REAL CSOIL
3699	REAL DENOM
3700	REAL DEW
3701	REAL DF1
3702	REAL DKSAT
3703	REAL DRIP
3704	REAL DSOIL
3705	REAL DTOT
3706	REAL DT
3707	REAL DWSAT
3708	REAL EC
3709	REAL EDIR
3710	REAL EPSCA
3711	REAL ESD
3712	REAL ESDMIN
3713	REAL EXPSNO
3714	REAL EXPSOI
3715	REAL ETA
3716	REAL ETA1
3717	REAL ETP
3718	REAL ETP1
3719	REAL ETP2
3720	REAL ET(NSOIL)
3721	REAL ETT
3722	REAL EX
3723	REAL EXPFAC
3724	REAL FDOWN
3725	REAL FXEXP
3726	REAL FLX1
3727	REAL FLX2
3728	REAL FLX3
3729	REAL F1
3730	REAL KDT
3731	REAL LSUBF
3732	REAL LSUBC
3733	REAL LSUBS
3734	REAL PC
3735	REAL PRCP
3736	REAL PRCP1
3737	REAL Q2
3738	REAL RCH
3739	REAL RR
3740	REAL RTDIS(NSOIL)
3741	REAL SSOIL
3742	REAL SBETA
3743	REAL SSOIL1
3744	REAL SFCTMP
3745	REAL SHDFAC
3746	REAL SIGMA
3747	REAL SMC(NSOIL)
3748	REAL SH2O(NSOIL)
3749	REAL SMCDRY
3750	REAL SMCMAX
3751	REAL SMCREF
3752	REAL SMCWLT
3753	REAL SNOMLT
3754	REAL SNOWH
3755	REAL STC(NSOIL)
3756	REAL T1
3757	REAL T11
3758	REAL T12
3759	REAL T12A
3760	REAL T12B
3761	REAL T24
3762	REAL TBOT
3763	REAL ZBOT
3764	REAL TH2
3765	REAL YY
3766	REAL ZSOIL(NSOIL)
3767	REAL ZZ1
3768	REAL TFREEZ
3769	REAL SALP
3770	REAL SFCPRS
3771	REAL SLOPE
3772	REAL FRZFACT
3773	REAL PSISAT
3774	REAL SNUP
3775	REAL RUNOFF1
3776	REAL RUNOFF2
3777	REAL RUNOFF3
3778	REAL QUARTZ
3779	REAL SNDENS
3780	REAL SNCOND
3781	REAL RSNOW
3782	REAL SNCOVR
3783	REAL QSAT
3784	REAL ETP3
3785	REAL SEH
3786	REAL T14
3787	! REAL CSNOW
3788
3789	REAL EC1
3790	REAL EDIR1
3791	REAL ET1(NSOIL)
3792	REAL ETT1
3793
3794	REAL ETNS
3795	REAL ETNS1
3796	REAL ESNOW
3797	REAL ESNOW1
3798	REAL ESNOW2
3799	REAL ETANRG
3800
3801	INTEGER K
3802
3803	REAL SNOEXP
3804
3805	PARAMETER(CP = 1004.5)
3806	PARAMETER(CPH2O = 4.218E+3)
3807	PARAMETER(CPICE = 2.106E+3)
3808	PARAMETER(ESDMIN = 1.E-6)
3809	PARAMETER(LSUBF = 3.335E+5)
3810	PARAMETER(LSUBC = 2.501000E+6)
3811	PARAMETER(LSUBS = 2.83E+6)
3812	PARAMETER(SIGMA = 5.67E-8)
3813	PARAMETER(TFREEZ = 273.15)
3814
3815	! DATA SNOEXP /1.0/
3816	DATA SNOEXP /2.0/
3817
3818	! ----------------------------------------------------------------------
3819	! EXECUTABLE CODE BEGINS HERE:
3820	! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO M S-1 AND THEN TO AN
3821	! AMOUNT (M) GIVEN TIMESTEP (DT) AND CALL IT AN EFFECTIVE SNOWPACK
3822	! REDUCTION AMOUNT, ETP2 (M). THIS IS THE AMOUNT THE SNOWPACK WOULD BE
3823	! REDUCED DUE TO EVAPORATION FROM THE SNOW SFC DURING THE TIMESTEP.
3824	! EVAPORATION WILL PROCEED AT THE POTENTIAL RATE UNLESS THE SNOW DEPTH
3825	! IS LESS THAN THE EXPECTED SNOWPACK REDUCTION.
3826	! IF SEAICE (ICE=1), BETA REMAINS=1.
3827	! ----------------------------------------------------------------------
3828	PRCP1 = PRCP1*0.001
3829
3830	! ETP2 = ETP * 0.001 * DT
3831	! BETA = 1.0
3832	! IF (ICE .NE. 1) THEN
3833	! IF (ESD .LT. ETP2) THEN
3834	! BETA = ESD / ETP2
3835	! ENDIF
3836	! ENDIF
3837
3838	! ----------------------------------------------------------------------
3839	EDIR = 0.0
3840	EDIR1 = 0.0
3841	EC = 0.0
3842	EC1 = 0.0
3843	DO K = 1,NSOIL
3844	ET(K) = 0.0
3845	ET1(K) = 0.0
3846	ENDDO
3847	ETT = 0.0
3848	ETT1 = 0.0
3849	ETNS = 0.0
3850	ETNS1 = 0.0
3851	ESNOW = 0.0
3852	ESNOW1 = 0.0
3853	ESNOW2 = 0.0
3854	! ----------------------------------------------------------------------
3855
3856	! ----------------------------------------------------------------------
3857	! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
3858	! ----------------------------------------------------------------------
3859	DEW = 0.0
3860	ETP1 = ETP*0.001
3861	IF (ETP .LT. 0.0) THEN
3862	! DEW = -ETP * 0.001
3863	DEW = -ETP1
3864	! ESNOW2 = ETP * 0.001 * DT
3865	ESNOW2 = ETP1 * DT
3866	ETANRG = ETP((1.-SNCOVR)LSUBC + SNCOVR*LSUBS)
3867	! ENDIF
3868	ELSE
3869	! ----------------------------------------------------------------------
3870	! ETP1 = 0.0
3871	! ETP1 = ETP*0.001
3872	IF (ICE .NE. 1) THEN
3873	IF (SNCOVR .LT. 1.) THEN
3874	! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
3875	CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3876	& SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3877	& SMCREF,SHDFAC,CMCMAX, &
3878	& SMCDRY,CFACTR, &
3879	& EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
3880	! ENDIF
3881	! ----------------------------------------------------------------------
3882	EDIR1 = EDIR1*(1.-SNCOVR)
3883	EC1 = EC1*(1.-SNCOVR)
3884	DO K = 1,NSOIL
3885	ET1(K) = ET1(K)*(1.-SNCOVR)
3886	END DO
3887	ETT1 = ETT1*(1.-SNCOVR)
3888	ETNS1 = ETNS1*(1.-SNCOVR)
3889	! ----------------------------------------------------------------------
3890	EDIR = EDIR1 * 1000.0
3891	EC = EC1 * 1000.0
3892	DO K = 1,NSOIL
3893	ET(K) = ET1(K) * 1000.0
3894	END DO
3895	ETT = ETT1 * 1000.0
3896	ETNS = ETNS1 * 1000.0
3897	! ----------------------------------------------------------------------
3898	ENDIF
3899	! ESNOW = ETP*SNCOVR
3900	! ESNOW1 = ETP*0.001
3901	! ESNOW1 = ESNOW*0.001
3902	! ESNOW2 = ESNOW1*DT
3903	! ETANRG = ESNOWLSUBS + ETNSLSUBC
3904	ENDIF
3905	ESNOW = ETP*SNCOVR
3906	ESNOW1 = ESNOW*0.001
3907	ESNOW2 = ESNOW1*DT
3908	ETANRG = ESNOWLSUBS + ETNSLSUBC
3909	ENDIF
3910
3911	! ----------------------------------------------------------------------
3912	! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3913	! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3914	! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3915	! SNOWFALL STRIKING THE GOUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3916	! ----------------------------------------------------------------------
3917	FLX1 = 0.0
3918	IF (SNOWNG) THEN
3919	FLX1 = CPICE * PRCP * (T1 - SFCTMP)
3920	ELSE
3921	IF (PRCP .GT. 0.0) FLX1 = CPH2O * PRCP * (T1 - SFCTMP)
3922	ENDIF
3923
3924	! ----------------------------------------------------------------------
3925	! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3926	! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3927	! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3928	! FLUXES.
3929	! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
3930	! PENMAN.
3931	! ----------------------------------------------------------------------
3932	DSOIL = -(0.5 * ZSOIL(1))
3933	DTOT = SNOWH + DSOIL
3934	DENOM = 1.0 + DF1 / (DTOT * RR * RCH)
3935	! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3936	! & + TH2 - SFCTMP - BETA*EPSCA ) / RR
3937	! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3938	! M.Ek, 24Nov04, add snow emissivity
3939	T12A = ((FDOWN-FLX1-FLX2 &
3940	& -(0.95SNCOVR+(1.0-SNCOVR))SIGMA*T24)/RCH &
3941	& + TH2 - SFCTMP - ETANRG/RCH ) / RR
3942	T12B = DF1 * STC(1) / (DTOT * RR * RCH)
3943	T12 = (SFCTMP + T12A + T12B) / DENOM
3944
3945	! ----------------------------------------------------------------------
3946	! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
3947	! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
3948	! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
3949	! DEPENDING ON SIGN OF ETP.
3950	! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
3951	! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
3952	! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
3953	! TO ZERO.
3954	! ----------------------------------------------------------------------
3955	IF (T12 .LE. TFREEZ) THEN
3956	T1 = T12
3957	SSOIL = DF1 * (T1 - STC(1)) / DTOT
3958	! ESD = MAX(0.0, ESD-ETP2)
3959	ESD = MAX(0.0, ESD-ESNOW2)
3960	FLX3 = 0.0
3961	EX = 0.0
3962	SNOMLT = 0.0
3963
3964	ELSE
3965	! ----------------------------------------------------------------------
3966	! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3967	! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
3968	! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
3969	! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
3970	! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
3971	! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
3972	! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
3973	! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
3974	! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
3975	! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
3976	! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
3977	! ----------------------------------------------------------------------
3978	! T1 = TFREEZ * SNCOVR + T12 * (1.0 - SNCOVR)
3979	! mek Feb2004
3980	! non-linear weighting of snow vs non-snow covered portions of gridbox
3981	! so with SNOEXP = 2.0 (>1), surface skin temperature is higher than for
3982	! the linear case (SNOEXP = 1).
3983	T1 = TFREEZ * SNCOVR*SNOEXP + T12 (1.0 - SNCOVR**SNOEXP)
3984	! QSAT = (0.6226.11E2)/(SFCPRS-0.3786.11E2)
3985	! ETP = RCH*(QSAT-Q2)/CP
3986	! ETP2 = ETP0.001DT
3987	! BETA = 1.0
3988	SSOIL = DF1 * (T1 - STC(1)) / DTOT
3989
3990	! ----------------------------------------------------------------------
3991	! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
3992	! BETA<1
3993	! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
3994	! ----------------------------------------------------------------------
3995	! IF (ESD .LE. ETP2) THEN
3996	! IF (ESD .LE. ESNOW2) THEN
3997	IF (ESD-ESNOW2 .LE. ESDMIN) THEN
3998	! BETA = ESD / ETP2
3999	ESD = 0.0
4000	EX = 0.0
4001	SNOMLT = 0.0
4002
4003	ELSE
4004	! ----------------------------------------------------------------------
4005	! POTENTIAL EVAP (SUBLIMATION) LESS THAN DEPTH OF SNOWPACK, RETAIN
4006	! BETA=1.
4007	! SNOWPACK (ESD) REDUCED BY POTENTIAL EVAP RATE
4008	! ETP3 (CONVERT TO FLUX)
4009	! ----------------------------------------------------------------------
4010	! ESD = ESD-ETP2
4011	ESD = ESD-ESNOW2
4012	! ETP3 = ETP*LSUBC
4013	SEH = RCH*(T1-TH2)
4014	T14 = T1*T1
4015	T14 = T14*T14
4016	! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETP3
4017	! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETANRG
4018	! M.Ek, 24Nov04, add snow emissivity
4019	FLX3 = FDOWN - FLX1 - FLX2 - &
4020	& (0.95SNCOVR+(1.0-SNCOVR))SIGMA*T14 - SSOIL - SEH - ETANRG
4021	IF (FLX3 .LE. 0.0) FLX3 = 0.0
4022	EX = FLX3*0.001/LSUBF
4023
4024	! ----------------------------------------------------------------------
4025	! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
4026	! IF SNOW COVER LESS THAN 5% NO SNOWMELT REDUCTION
4027	! ***NOTE: DOES 'IF' BELOW FAIL TO MATCH THE MELT WATER WITH THE MELT
4028	! ENERGY?
4029	! ----------------------------------------------------------------------
4030	! IF (SNCOVR .GT. 0.05) EX = EX * SNCOVR
4031	SNOMLT = EX * DT
4032
4033	! ----------------------------------------------------------------------
4034	! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
4035	! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
4036	! ----------------------------------------------------------------------
4037	IF (ESD-SNOMLT .GE. ESDMIN) THEN
4038	ESD = ESD - SNOMLT
4039
4040	ELSE
4041	! ----------------------------------------------------------------------
4042	! SNOWMELT EXCEEDS SNOW DEPTH
4043	! ----------------------------------------------------------------------
4044	EX = ESD/DT
4045	FLX3 = EX1000.0LSUBF
4046	SNOMLT = ESD
4047	ESD = 0.0
4048
4049	ENDIF
4050	! ----------------------------------------------------------------------
4051	! END OF 'ESD .LE. ETP2' IF-BLOCK
4052	! ----------------------------------------------------------------------
4053	ENDIF
4054
4055	PRCP1 = PRCP1 + EX
4056
4057	! ----------------------------------------------------------------------
4058	! END OF 'T12 .LE. TFREEZ' IF-BLOCK
4059	! ----------------------------------------------------------------------
4060	ENDIF
4061
4062	! ----------------------------------------------------------------------
4063	! FINAL BETA NOW IN HAND, SO COMPUTE EVAPORATION. EVAP EQUALS ETP
4064	! UNLESS BETA<1.
4065	! ----------------------------------------------------------------------
4066	! ETA = BETA*ETP
4067
4068	! ----------------------------------------------------------------------
4069	! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
4070	! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
4071	! (BELOW).
4072	! IF SEAICE (ICE=1) SKIP CALL TO SMFLX.
4073	! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES. IN THIS, THE SNOW PACK
4074	! CASE, ETA1 IS NOT USED IN CALCULATION OF EVAP.
4075	! ----------------------------------------------------------------------
4076	! ETP1 = 0.0
4077	IF (ICE .NE. 1) THEN
4078	! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
4079	! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
4080	! & SMCREF,SHDFAC,CMCMAX,
4081	! & SMCDRY,CFACTR,
4082	! & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
4083	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
4084	& SH2O,SLOPE,KDT,FRZFACT, &
4085	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
4086	& SHDFAC,CMCMAX, &
4087	& RUNOFF1,RUNOFF2,RUNOFF3, &
4088	& EDIR1,EC1,ET1, &
4089	& DRIP)
4090
4091	ENDIF
4092
4093	! ----------------------------------------------------------------------
4094	! EDIR = EDIR1 * 1000.0
4095	! EC = EC1 * 1000.0
4096	! ETT = ETT1 * 1000.0
4097	! ET(1) = ET1(1) * 1000.0
4098	! ET(2) = ET1(2) * 1000.0
4099	! ET(3) = ET1(3) * 1000.0
4100	! ET(4) = ET1(4) * 1000.0
4101	! ----------------------------------------------------------------------
4102
4103	! ----------------------------------------------------------------------
4104	! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
4105	! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
4106	! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
4107	! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
4108	! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
4109	! SKIN TEMP VALUE AS REVISED BY SHFLX.
4110	! ----------------------------------------------------------------------
4111	ZZ1 = 1.0
4112	YY = STC(1)-0.5SSOILZSOIL(1)*ZZ1/DF1
4113	T11 = T1
4114
4115	! ----------------------------------------------------------------------
4116	! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
4117	! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
4118	! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
4119	! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
4120	! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
4121	! ----------------------------------------------------------------------
4122	CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
4123	& TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
4124	& QUARTZ,CSOIL)
4125
4126	! ----------------------------------------------------------------------
4127	! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
4128	! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
4129	! ----------------------------------------------------------------------
4130	IF (ESD .GT. 0.) THEN
4131	CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
4132	ELSE
4133	ESD = 0.
4134	SNOWH = 0.
4135	SNDENS = 0.
4136	SNCOND = 1.
4137	SNCOVR = 0.
4138	ENDIF
4139
4140	! ----------------------------------------------------------------------
4141	! END SUBROUTINE SNOPAC
4142	! ----------------------------------------------------------------------
4143	END SUBROUTINE SNOPAC
4144
4145	SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
4146
4147	IMPLICIT NONE
4148
4149	! ----------------------------------------------------------------------
4150	! SUBROUTINE SNOWPACK
4151	! ----------------------------------------------------------------------
4152	! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
4153	! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
4154	! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
4155	! KOREN, 03/25/95.
4156	! ----------------------------------------------------------------------
4157	! ESD WATER EQUIVALENT OF SNOW (M)
4158	! DTSEC TIME STEP (SEC)
4159	! SNOWH SNOW DEPTH (M)
4160	! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4161	! TSNOW SNOW SURFACE TEMPERATURE (K)
4162	! TSOIL SOIL SURFACE TEMPERATURE (K)
4163	!
4164	! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
4165	! ----------------------------------------------------------------------
4166	INTEGER IPOL, J
4167
4168	REAL BFAC,C1,C2,SNDENS,DSX,DTHR,DTSEC,DW,SNOWHC,SNOWH,PEXP,TAVGC, &
4169	& TSNOW,TSNOWC,TSOIL,TSOILC,ESD,ESDC,ESDCX,G,KN
4170
4171	PARAMETER(C1 = 0.01, C2=21.0, G=9.81, KN=4000.0)
4172
4173	! ----------------------------------------------------------------------
4174	! CONVERSION INTO SIMULATION UNITS
4175	! ----------------------------------------------------------------------
4176	SNOWHC = SNOWH*100.
4177	ESDC = ESD*100.
4178	DTHR = DTSEC/3600.
4179	TSNOWC = TSNOW-273.15
4180	TSOILC = TSOIL-273.15
4181
4182	! ----------------------------------------------------------------------
4183	! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
4184	! ----------------------------------------------------------------------
4185	TAVGC = 0.5*(TSNOWC+TSOILC)
4186
4187	! ----------------------------------------------------------------------
4188	! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
4189	! SNDENS=DS0(EXP(BFACESD)-1.)/(BFAC*ESD)
4190	! BFAC=DTHRC1EXP(0.08TAVGC-C2DS0)
4191	! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
4192	! NUMERICALLY BELOW:
4193	! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
4194	! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
4195	! ----------------------------------------------------------------------
4196	IF (ESDC .GT. 1.E-2) THEN
4197	ESDCX = ESDC
4198	ELSE
4199	ESDCX = 1.E-2
4200	ENDIF
4201	BFAC = DTHRC1EXP(0.08TAVGC-C2SNDENS)
4202
4203	! DSX = SNDENS((DEXP(BFACESDC)-1.)/(BFAC*ESDC))
4204	! ----------------------------------------------------------------------
4205	! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
4206	! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
4207	! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
4208	! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
4209	! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
4210	! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
4211	! POLYNOMIAL EXPANSION.
4212	!
4213	! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
4214	! IS GOVERNED BY ITERATION LIMIT "IPOL".
4215	! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
4216	! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
4217	! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
4218	! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
4219	! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
4220	! ----------------------------------------------------------------------
4221	IPOL = 4
4222	PEXP = 0.
4223	DO J = IPOL,1,-1
4224	! PEXP = (1. + PEXP)BFACESDC/REAL(J+1)
4225	PEXP = (1. + PEXP)BFACESDCX/REAL(J+1)
4226	END DO
4227	PEXP = PEXP + 1.
4228
4229	DSX = SNDENS*(PEXP)
4230	! ----------------------------------------------------------------------
4231	! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
4232	! ----------------------------------------------------------------------
4233	! END OF KOREAN FORMULATION
4234
4235	! BASE FORMULATION (COGLEY ET AL., 1990)
4236	! CONVERT DENSITY FROM G/CM3 TO KG/M3
4237	! DSM=SNDENS*1000.0
4238
4239	! DSX=DSM+DTSEC0.5DSMGESD/
4240	! & (1E7EXP(-0.02DSM+KN/(TAVGC+273.16)-14.643))
4241
4242	! CONVERT DENSITY FROM KG/M3 TO G/CM3
4243	! DSX=DSX/1000.0
4244
4245	! END OF COGLEY ET AL. FORMULATION
4246
4247	! ----------------------------------------------------------------------
4248	! SET UPPER/LOWER LIMIT ON SNOW DENSITY
4249	! ----------------------------------------------------------------------
4250	IF (DSX .GT. 0.40) DSX = 0.40
4251	IF (DSX .LT. 0.05) DSX = 0.05
4252	SNDENS = DSX
4253	! ----------------------------------------------------------------------
4254	! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
4255	! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
4256	! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
4257	! ----------------------------------------------------------------------
4258	IF (TSNOWC .GE. 0.) THEN
4259	DW = 0.13*DTHR/24.
4260	SNDENS = SNDENS*(1.-DW)+DW
4261	IF (SNDENS .GT. 0.40) SNDENS = 0.40
4262	ENDIF
4263
4264	! ----------------------------------------------------------------------
4265	! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
4266	! CHANGE SNOW DEPTH UNITS TO METERS
4267	! ----------------------------------------------------------------------
4268	SNOWHC = ESDC/SNDENS
4269	SNOWH = SNOWHC*0.01
4270
4271	! ----------------------------------------------------------------------
4272	! END SUBROUTINE SNOWPACK
4273	! ----------------------------------------------------------------------
4274	END SUBROUTINE SNOWPACK
4275
4276	SUBROUTINE SNOWZ0 (SNCOVR,Z0)
4277
4278	IMPLICIT NONE
4279
4280	! ----------------------------------------------------------------------
4281	! SUBROUTINE SNOWZ0
4282	! ----------------------------------------------------------------------
4283	! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
4284	! SNCOVR FRACTIONAL SNOW COVER
4285	! Z0 ROUGHNESS LENGTH (m)
4286	! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
4287	! ----------------------------------------------------------------------
4288	REAL SNCOVR, Z0, Z0S
4289	! PARAMETER (Z0S=0.001)
4290
4291	! CURRENT NOAH LSM CONDITION - MBEK, 09-OCT-2001
4292	Z0S = Z0
4293	!
4294	Z0 = (1-SNCOVR)Z0 + SNCOVRZ0S
4295	! ----------------------------------------------------------------------
4296	! END SUBROUTINE SNOWZ0
4297	! ----------------------------------------------------------------------
4298	END SUBROUTINE SNOWZ0
4299
4300	SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
4301
4302	IMPLICIT NONE
4303
4304	! ----------------------------------------------------------------------
4305	! SUBROUTINE SNOW_NEW
4306	! ----------------------------------------------------------------------
4307	! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
4308	! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
4309	!
4310	! TEMP AIR TEMPERATURE (K)
4311	! NEWSN NEW SNOWFALL (M)
4312	! SNOWH SNOW DEPTH (M)
4313	! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4314	! ----------------------------------------------------------------------
4315	REAL SNDENS
4316	REAL DSNEW
4317	REAL SNOWHC
4318	REAL HNEWC
4319	REAL SNOWH
4320	REAL NEWSN
4321	REAL NEWSNC
4322	REAL TEMP
4323	REAL TEMPC
4324
4325	! ----------------------------------------------------------------------
4326	! CONVERSION INTO SIMULATION UNITS
4327	! ----------------------------------------------------------------------
4328	SNOWHC = SNOWH*100.
4329	NEWSNC = NEWSN*100.
4330	TEMPC = TEMP-273.15
4331
4332	! ----------------------------------------------------------------------
4333	! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
4334	! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
4335	! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
4336	! VEMADOLEN, SWEDEN, 1980, 172-177PP.
4337	!-----------------------------------------------------------------------
4338	IF (TEMPC .LE. -15.) THEN
4339	DSNEW = 0.05
4340	ELSE
4341	DSNEW = 0.05+0.0017(TEMPC+15.)*1.5
4342	ENDIF
4343
4344	! ----------------------------------------------------------------------
4345	! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
4346	! ----------------------------------------------------------------------
4347	HNEWC = NEWSNC/DSNEW
4348	SNDENS = (SNOWHCSNDENS+HNEWCDSNEW)/(SNOWHC+HNEWC)
4349	SNOWHC = SNOWHC+HNEWC
4350	SNOWH = SNOWHC*0.01
4351
4352	! ----------------------------------------------------------------------
4353	! END SUBROUTINE SNOW_NEW
4354	! ----------------------------------------------------------------------
4355	END SUBROUTINE SNOW_NEW
4356
4357	SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
4358	& ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
4359	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
4360
4361	IMPLICIT NONE
4362
4363	! ----------------------------------------------------------------------
4364	! SUBROUTINE SRT
4365	! ----------------------------------------------------------------------
4366	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
4367	! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
4368	! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
4369	! ----------------------------------------------------------------------
4370	INTEGER NSOLD
4371	PARAMETER(NSOLD = 20)
4372
4373	INTEGER CVFRZ
4374	INTEGER IALP1
4375	INTEGER IOHINF
4376	INTEGER J
4377	INTEGER JJ
4378	INTEGER K
4379	INTEGER KS
4380	INTEGER NSOIL
4381
4382	REAL ACRT
4383	REAL AI(NSOLD)
4384	REAL BEXP
4385	REAL BI(NSOLD)
4386	REAL CI(NSOLD)
4387	REAL DD
4388	REAL DDT
4389	REAL DDZ
4390	REAL DDZ2
4391	REAL DENOM
4392	REAL DENOM2
4393	REAL DICE
4394	REAL DKSAT
4395	REAL DMAX(NSOLD)
4396	REAL DSMDZ
4397	REAL DSMDZ2
4398	REAL DT
4399	REAL DT1
4400	REAL DWSAT
4401	REAL EDIR
4402	REAL ET(NSOIL)
4403	REAL FCR
4404	REAL FRZX
4405	REAL INFMAX
4406	REAL KDT
4407	REAL MXSMC
4408	REAL MXSMC2
4409	REAL NUMER
4410	REAL PCPDRP
4411	REAL PDDUM
4412	REAL PX
4413	REAL RHSTT(NSOIL)
4414	REAL RUNOFF1
4415	REAL RUNOFF2
4416	REAL SH2O(NSOIL)
4417	REAL SH2OA(NSOIL)
4418	REAL SICE(NSOIL)
4419	REAL SICEMAX
4420	REAL SLOPE
4421	REAL SLOPX
4422	REAL SMCAV
4423	REAL SMCMAX
4424	REAL SMCWLT
4425	REAL SSTT
4426	REAL SUM
4427	REAL VAL
4428	REAL WCND
4429	REAL WCND2
4430	REAL WDF
4431	REAL WDF2
4432	REAL ZSOIL(NSOIL)
4433
4434	! ----------------------------------------------------------------------
4435	! FROZEN GROUND VERSION:
4436	! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4437	! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4438	! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
4439	! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
4440	! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
4441	! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
4442	! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
4443	! ----------------------------------------------------------------------
4444	PARAMETER(CVFRZ = 3)
4445
4446	! ----------------------------------------------------------------------
4447	! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
4448	! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
4449	! MODIFIED BY Q DUAN
4450	! ----------------------------------------------------------------------
4451	IOHINF=1
4452
4453	! ----------------------------------------------------------------------
4454	! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
4455	! LAYERS.
4456	! ----------------------------------------------------------------------
4457	SICEMAX = 0.0
4458	DO KS=1,NSOIL
4459	IF (SICE(KS) .GT. SICEMAX) SICEMAX = SICE(KS)
4460	END DO
4461
4462	! ----------------------------------------------------------------------
4463	! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
4464	! ----------------------------------------------------------------------
4465	PDDUM = PCPDRP
4466	RUNOFF1 = 0.0
4467	IF (PCPDRP .NE. 0.0) THEN
4468
4469	! ----------------------------------------------------------------------
4470	! MODIFIED BY Q. DUAN, 5/16/94
4471	! ----------------------------------------------------------------------
4472	! IF (IOHINF .EQ. 1) THEN
4473
4474	DT1 = DT/86400.
4475	SMCAV = SMCMAX - SMCWLT
4476	DMAX(1)=-ZSOIL(1)*SMCAV
4477
4478	! ----------------------------------------------------------------------
4479	! FROZEN GROUND VERSION:
4480	! ----------------------------------------------------------------------
4481	DICE = -ZSOIL(1) * SICE(1)
4482
4483	DMAX(1)=DMAX(1)*(1.0 - (SH2OA(1)+SICE(1)-SMCWLT)/SMCAV)
4484	DD=DMAX(1)
4485
4486	DO KS=2,NSOIL
4487
4488	! ----------------------------------------------------------------------
4489	! FROZEN GROUND VERSION:
4490	! ----------------------------------------------------------------------
4491	DICE = DICE + ( ZSOIL(KS-1) - ZSOIL(KS) ) * SICE(KS)
4492
4493	DMAX(KS) = (ZSOIL(KS-1)-ZSOIL(KS))*SMCAV
4494	DMAX(KS) = DMAX(KS)*(1.0 - (SH2OA(KS)+SICE(KS)-SMCWLT)/SMCAV)
4495	DD = DD+DMAX(KS)
4496	END DO
4497
4498	! ----------------------------------------------------------------------
4499	! VAL = (1.-EXP(-KDT*SQRT(DT1)))
4500	! IN BELOW, REMOVE THE SQRT IN ABOVE
4501	! ----------------------------------------------------------------------
4502	VAL = (1.-EXP(-KDT*DT1))
4503	DDT = DD*VAL
4504	PX = PCPDRP*DT
4505	IF (PX .LT. 0.0) PX = 0.0
4506	INFMAX = (PX*(DDT/(PX+DDT)))/DT
4507
4508	! ----------------------------------------------------------------------
4509	! FROZEN GROUND VERSION:
4510	! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4511	! ----------------------------------------------------------------------
4512	FCR = 1.
4513	IF (DICE .GT. 1.E-2) THEN
4514	ACRT = CVFRZ * FRZX / DICE
4515	SUM = 1.
4516	IALP1 = CVFRZ - 1
4517	DO J = 1,IALP1
4518	K = 1
4519	DO JJ = J+1,IALP1
4520	K = K * JJ
4521	END DO
4522	SUM = SUM + (ACRT ** ( CVFRZ-J)) / FLOAT (K)
4523	END DO
4524	FCR = 1. - EXP(-ACRT) * SUM
4525	ENDIF
4526	INFMAX = INFMAX * FCR
4527
4528	! ----------------------------------------------------------------------
4529	! CORRECTION OF INFILTRATION LIMITATION:
4530	! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
4531	! HYDROLIC CONDUCTIVITY
4532	! ----------------------------------------------------------------------
4533	! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
4534	MXSMC = SH2OA(1)
4535
4536	CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4537	& SICEMAX)
4538
4539	INFMAX = MAX(INFMAX,WCND)
4540	INFMAX = MIN(INFMAX,PX)
4541
4542	IF (PCPDRP .GT. INFMAX) THEN
4543	RUNOFF1 = PCPDRP - INFMAX
4544	PDDUM = INFMAX
4545	ENDIF
4546
4547	ENDIF
4548
4549	! ----------------------------------------------------------------------
4550	! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
4551	! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4552	! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
4553	! ----------------------------------------------------------------------
4554	MXSMC = SH2OA(1)
4555
4556	CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4557	& SICEMAX)
4558
4559	! ----------------------------------------------------------------------
4560	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
4561	! ----------------------------------------------------------------------
4562	DDZ = 1. / ( -.5 * ZSOIL(2) )
4563	AI(1) = 0.0
4564	BI(1) = WDF * DDZ / ( -ZSOIL(1) )
4565	CI(1) = -BI(1)
4566
4567	! ----------------------------------------------------------------------
4568	! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
4569	! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
4570	! ----------------------------------------------------------------------
4571	DSMDZ = ( SH2O(1) - SH2O(2) ) / ( -.5 * ZSOIL(2) )
4572	RHSTT(1) = (WDF * DSMDZ + WCND - PDDUM + EDIR + ET(1))/ZSOIL(1)
4573	SSTT = WDF * DSMDZ + WCND + EDIR + ET(1)
4574
4575	! ----------------------------------------------------------------------
4576	! INITIALIZE DDZ2
4577	! ----------------------------------------------------------------------
4578	DDZ2 = 0.0
4579
4580	! ----------------------------------------------------------------------
4581	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
4582	! ----------------------------------------------------------------------
4583	DO K = 2,NSOIL
4584	DENOM2 = (ZSOIL(K-1) - ZSOIL(K))
4585	IF (K .NE. NSOIL) THEN
4586	SLOPX = 1.
4587
4588	! ----------------------------------------------------------------------
4589	! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
4590	! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4591	! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
4592	! ----------------------------------------------------------------------
4593	MXSMC2 = SH2OA(K)
4594
4595	CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
4596	& SICEMAX)
4597
4598	! ----------------------------------------------------------------------
4599	! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
4600	! ----------------------------------------------------------------------
4601	DENOM = (ZSOIL(K-1) - ZSOIL(K+1))
4602	DSMDZ2 = (SH2O(K) - SH2O(K+1)) / (DENOM * 0.5)
4603
4604	! ----------------------------------------------------------------------
4605	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
4606	! ----------------------------------------------------------------------
4607	DDZ2 = 2.0 / DENOM
4608	CI(K) = -WDF2 * DDZ2 / DENOM2
4609	ELSE
4610
4611	! ----------------------------------------------------------------------
4612	! SLOPE OF BOTTOM LAYER IS INTRODUCED
4613	! ----------------------------------------------------------------------
4614	SLOPX = SLOPE
4615
4616	! ----------------------------------------------------------------------
4617	! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
4618	! THIS LAYER
4619	! ----------------------------------------------------------------------
4620	CALL WDFCND (WDF2,WCND2,SH2OA(NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
4621	& SICEMAX)
4622
4623	! ----------------------------------------------------------------------
4624	! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
4625	! ----------------------------------------------------------------------
4626	DSMDZ2 = 0.0
4627
4628	! ----------------------------------------------------------------------
4629	! SET MATRIX COEF CI TO ZERO
4630	! ----------------------------------------------------------------------
4631	CI(K) = 0.0
4632	ENDIF
4633
4634	! ----------------------------------------------------------------------
4635	! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
4636	! ----------------------------------------------------------------------
4637	NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2 - (WDF * DSMDZ) &
4638	& - WCND + ET(K)
4639	RHSTT(K) = NUMER / (-DENOM2)
4640
4641	! ----------------------------------------------------------------------
4642	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
4643	! ----------------------------------------------------------------------
4644	AI(K) = -WDF * DDZ / DENOM2
4645	BI(K) = -( AI(K) + CI(K) )
4646
4647	! ----------------------------------------------------------------------
4648	! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
4649	! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
4650	! ----------------------------------------------------------------------
4651	IF (K .EQ. NSOIL) THEN
4652	RUNOFF2 = SLOPX * WCND2
4653	ENDIF
4654
4655	IF (K .NE. NSOIL) THEN
4656	WDF = WDF2
4657	WCND = WCND2
4658	DSMDZ = DSMDZ2
4659	DDZ = DDZ2
4660	ENDIF
4661	END DO
4662
4663	! ----------------------------------------------------------------------
4664	! END SUBROUTINE SRT
4665	! ----------------------------------------------------------------------
4666	END SUBROUTINE SRT
4667
4668	SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
4669	& NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
4670	& AI,BI,CI)
4671
4672	IMPLICIT NONE
4673
4674	! ----------------------------------------------------------------------
4675	! SUBROUTINE SSTEP
4676	! ----------------------------------------------------------------------
4677	! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
4678	! CONTENT VALUES.
4679	! ----------------------------------------------------------------------
4680	INTEGER NSOLD
4681	PARAMETER(NSOLD = 20)
4682
4683	INTEGER I
4684	INTEGER K
4685	INTEGER KK11
4686	INTEGER NSOIL
4687
4688	REAL AI(NSOLD)
4689	REAL BI(NSOLD)
4690	REAL CI(NSOLD)
4691	REAL CIin(NSOLD)
4692	REAL CMC
4693	REAL CMCMAX
4694	REAL DDZ
4695	REAL DT
4696	REAL RHSCT
4697	REAL RHSTT(NSOIL)
4698	REAL RHSTTin(NSOIL)
4699	REAL RUNOFF3
4700	REAL SH2OIN(NSOIL)
4701	REAL SH2OOUT(NSOIL)
4702	REAL SICE(NSOIL)
4703	REAL SMC(NSOIL)
4704	REAL SMCMAX
4705	REAL STOT
4706	REAL WPLUS
4707	REAL ZSOIL(NSOIL)
4708
4709	! ----------------------------------------------------------------------
4710	! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
4711	! TRI-DIAGONAL MATRIX ROUTINE.
4712	! ----------------------------------------------------------------------
4713	DO K = 1,NSOIL
4714	RHSTT(K) = RHSTT(K) * DT
4715	AI(K) = AI(K) * DT
4716	BI(K) = 1. + BI(K) * DT
4717	CI(K) = CI(K) * DT
4718	END DO
4719
4720	! ----------------------------------------------------------------------
4721	! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4722	! ----------------------------------------------------------------------
4723	DO K = 1,NSOIL
4724	RHSTTin(K) = RHSTT(K)
4725	END DO
4726	DO K = 1,NSOIL
4727	CIin(K) = CI(K)
4728	END DO
4729
4730	! ----------------------------------------------------------------------
4731	! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
4732	! ----------------------------------------------------------------------
4733	CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
4734
4735	! ----------------------------------------------------------------------
4736	! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
4737	! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
4738	! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
4739	! ----------------------------------------------------------------------
4740	WPLUS = 0.0
4741	RUNOFF3 = 0.
4742	DDZ = -ZSOIL(1)
4743
4744	DO K = 1,NSOIL
4745	IF (K .NE. 1) DDZ = ZSOIL(K - 1) - ZSOIL(K)
4746	SH2OOUT(K) = SH2OIN(K) + CI(K) + WPLUS / DDZ
4747
4748	STOT = SH2OOUT(K) + SICE(K)
4749	IF (STOT .GT. SMCMAX) THEN
4750	IF (K .EQ. 1) THEN
4751	DDZ = -ZSOIL(1)
4752	ELSE
4753	KK11 = K - 1
4754	DDZ = -ZSOIL(K) + ZSOIL(KK11)
4755	ENDIF
4756	WPLUS = (STOT-SMCMAX) * DDZ
4757	ELSE
4758	WPLUS = 0.
4759	ENDIF
4760	SMC(K) = MAX ( MIN(STOT,SMCMAX),0.02 )
4761	SH2OOUT(K) = MAX((SMC(K)-SICE(K)),0.0)
4762	END DO
4763
4764	RUNOFF3 = WPLUS
4765
4766	! ----------------------------------------------------------------------
4767	! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
4768	! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
4769	! ----------------------------------------------------------------------
4770	CMC = CMC + DT * RHSCT
4771	IF (CMC .LT. 1.E-20) CMC=0.0
4772	CMC = MIN(CMC,CMCMAX)
4773
4774	! ----------------------------------------------------------------------
4775	! END SUBROUTINE SSTEP
4776	! ----------------------------------------------------------------------
4777	END SUBROUTINE SSTEP
4778
4779	SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4780
4781	IMPLICIT NONE
4782
4783	! ----------------------------------------------------------------------
4784	! SUBROUTINE TBND
4785	! ----------------------------------------------------------------------
4786	! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4787	! THE MIDDLE LAYER TEMPERATURES
4788	! ----------------------------------------------------------------------
4789	INTEGER NSOIL
4790	INTEGER K
4791
4792	REAL TBND1
4793	REAL T0
4794	REAL TU
4795	REAL TB
4796	REAL ZB
4797	REAL ZBOT
4798	REAL ZUP
4799	REAL ZSOIL (NSOIL)
4800
4801	PARAMETER(T0 = 273.15)
4802
4803	! ----------------------------------------------------------------------
4804	! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4805	! ----------------------------------------------------------------------
4806	IF (K .EQ. 1) THEN
4807	ZUP = 0.
4808	ELSE
4809	ZUP = ZSOIL(K-1)
4810	ENDIF
4811
4812	! ----------------------------------------------------------------------
4813	! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4814	! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4815	! ----------------------------------------------------------------------
4816	IF (K .EQ. NSOIL) THEN
4817	ZB = 2.*ZBOT-ZSOIL(K)
4818	ELSE
4819	ZB = ZSOIL(K+1)
4820	ENDIF
4821
4822	! ----------------------------------------------------------------------
4823	! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4824	! ----------------------------------------------------------------------
4825	TBND1 = TU+(TB-TU)*(ZUP-ZSOIL(K))/(ZUP-ZB)
4826
4827	! ----------------------------------------------------------------------
4828	! END SUBROUTINE TBND
4829	! ----------------------------------------------------------------------
4830	END SUBROUTINE TBND
4831
4832	SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
4833
4834	IMPLICIT NONE
4835
4836	! ----------------------------------------------------------------------
4837	! SUBROUTINE TDFCND
4838	! ----------------------------------------------------------------------
4839	! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
4840	! POINT AND TIME.
4841	! ----------------------------------------------------------------------
4842	! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4843	! June 2001 CHANGES: FROZEN SOIL CONDITION.
4844	! ----------------------------------------------------------------------
4845	REAL DF
4846	REAL GAMMD
4847	REAL THKDRY
4848	REAL AKE
4849	REAL THKICE
4850	REAL THKO
4851	REAL THKQTZ
4852	REAL THKSAT
4853	REAL THKS
4854	REAL THKW
4855	REAL QZ
4856	REAL SATRATIO
4857	REAL SH2O
4858	REAL SMC
4859	REAL SMCMAX
4860	REAL XU
4861	REAL XUNFROZ
4862
4863	! ----------------------------------------------------------------------
4864	! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
4865	! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
4866	! & 0.35, 0.60, 0.40, 0.82/
4867	! ----------------------------------------------------------------------
4868	! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4869	! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4870	! ----------------------------------------------------------------------
4871	! THKW ......WATER THERMAL CONDUCTIVITY
4872	! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4873	! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4874	! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4875	! THKICE ....ICE THERMAL CONDUCTIVITY
4876	! SMCMAX ....POROSITY (= SMCMAX)
4877	! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4878	! ----------------------------------------------------------------------
4879	! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4880	!
4881	! PABLO GRUNMANN, 08/17/98
4882	! REFS.:
4883	! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4884	! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4885	! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4886	! UNIVERSITY OF TRONDHEIM,
4887	! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4888	! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4889	! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4890	! VOL. 55, PP. 1209-1224.
4891	! ----------------------------------------------------------------------
4892	! NEEDS PARAMETERS
4893	! POROSITY(SOIL TYPE):
4894	! POROS = SMCMAX
4895	! SATURATION RATIO:
4896	SATRATIO = SMC/SMCMAX
4897
4898	! PARAMETERS W/(M.K)
4899	THKICE = 2.2
4900	THKW = 0.57
4901	THKO = 2.0
4902	! IF (QZ .LE. 0.2) THKO = 3.0
4903	THKQTZ = 7.7
4904	! SOLIDS' CONDUCTIVITY
4905	THKS = (THKQTZ*QZ)(THKO**(1.- QZ))
4906
4907	! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
4908	XUNFROZ = (SH2O + 1.E-9) / (SMC + 1.E-9)
4909
4910	! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
4911	XU=XUNFROZ*SMCMAX
4912	! SATURATED THERMAL CONDUCTIVITY
4913	THKSAT = THKS*(1.-SMCMAX)THKICE*(SMCMAX-XU)THKW**(XU)
4914
4915	! DRY DENSITY IN KG/M3
4916	GAMMD = (1. - SMCMAX)*2700.
4917
4918	! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4919	THKDRY = (0.135GAMMD + 64.7)/(2700. - 0.947GAMMD)
4920
4921	IF ( (SH2O + 0.0005) .LT. SMC ) THEN
4922	! FROZEN
4923	AKE = SATRATIO
4924	ELSE
4925	! UNFROZEN
4926	! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
4927	IF ( SATRATIO .GT. 0.1 ) THEN
4928
4929	! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4930	! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4931	! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4932
4933	AKE = LOG10(SATRATIO) + 1.0
4934
4935	ELSE
4936
4937	! USE K = KDRY
4938	AKE = 0.0
4939
4940	ENDIF
4941	ENDIF
4942
4943	! THERMAL CONDUCTIVITY
4944
4945	DF = AKE*(THKSAT - THKDRY) + THKDRY
4946
4947	! ----------------------------------------------------------------------
4948	! END SUBROUTINE TDFCND
4949	! ----------------------------------------------------------------------
4950	END SUBROUTINE TDFCND
4951
4952	SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
4953
4954	IMPLICIT NONE
4955
4956	! ----------------------------------------------------------------------
4957	! SUBROUTINE TMPAVG
4958	! ----------------------------------------------------------------------
4959	! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
4960	! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
4961	! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
4962	! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
4963	! ----------------------------------------------------------------------
4964	INTEGER K
4965	INTEGER NSOIL
4966
4967	REAL DZ
4968	REAL DZH
4969	REAL T0
4970	REAL TAVG
4971	REAL TDN
4972	REAL TM
4973	REAL TUP
4974	REAL X0
4975	REAL XDN
4976	REAL XUP
4977	REAL ZSOIL (NSOIL)
4978
4979	PARAMETER(T0 = 2.7315E2)
4980
4981	! ----------------------------------------------------------------------
4982	IF (K .EQ. 1) THEN
4983	DZ = -ZSOIL(1)
4984	ELSE
4985	DZ = ZSOIL(K-1)-ZSOIL(K)
4986	ENDIF
4987
4988	DZH=DZ*0.5
4989
4990	IF (TUP .LT. T0) THEN
4991	IF (TM .LT. T0) THEN
4992	IF (TDN .LT. T0) THEN
4993	! ----------------------------------------------------------------------
4994	! TUP, TM, TDN < T0
4995	! ----------------------------------------------------------------------
4996	TAVG = (TUP + 2.0*TM + TDN)/ 4.0
4997	ELSE
4998	! ----------------------------------------------------------------------
4999	! TUP & TM < T0, TDN >= T0
5000	! ----------------------------------------------------------------------
5001	X0 = (T0 - TM) * DZH / (TDN - TM)
5002	TAVG = 0.5 * (TUPDZH+TM(DZH+X0)+T0(2.DZH-X0)) / DZ
5003	ENDIF
5004	ELSE
5005	IF (TDN .LT. T0) THEN
5006	! ----------------------------------------------------------------------
5007	! TUP < T0, TM >= T0, TDN < T0
5008	! ----------------------------------------------------------------------
5009	XUP = (T0-TUP) * DZH / (TM-TUP)
5010	XDN = DZH - (T0-TM) * DZH / (TDN-TM)
5011	TAVG = 0.5 * (TUPXUP+T0(2.DZ-XUP-XDN)+TDNXDN) / DZ
5012	ELSE
5013	! ----------------------------------------------------------------------
5014	! TUP < T0, TM >= T0, TDN >= T0
5015	! ----------------------------------------------------------------------
5016	XUP = (T0-TUP) * DZH / (TM-TUP)
5017	TAVG = 0.5 * (TUPXUP+T0(2.*DZ-XUP)) / DZ
5018	ENDIF
5019	ENDIF
5020	ELSE
5021	IF (TM .LT. T0) THEN
5022	IF (TDN .LT. T0) THEN
5023	! ----------------------------------------------------------------------
5024	! TUP >= T0, TM < T0, TDN < T0
5025	! ----------------------------------------------------------------------
5026	XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5027	TAVG = 0.5 * (T0(DZ-XUP)+TM(DZH+XUP)+TDN*DZH) / DZ
5028	ELSE
5029	! ----------------------------------------------------------------------
5030	! TUP >= T0, TM < T0, TDN >= T0
5031	! ----------------------------------------------------------------------
5032	XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5033	XDN = (T0-TM) * DZH / (TDN-TM)
5034	TAVG = 0.5 * (T0(2.DZ-XUP-XDN)+TM*(XUP+XDN)) / DZ
5035	ENDIF
5036	ELSE
5037	IF (TDN .LT. T0) THEN
5038	! ----------------------------------------------------------------------
5039	! TUP >= T0, TM >= T0, TDN < T0
5040	! ----------------------------------------------------------------------
5041	XDN = DZH - (T0-TM) * DZH / (TDN-TM)
5042	TAVG = (T0(DZ-XDN)+0.5(T0+TDN)*XDN) / DZ
5043	ELSE
5044	! ----------------------------------------------------------------------
5045	! TUP >= T0, TM >= T0, TDN >= T0
5046	! ----------------------------------------------------------------------
5047	TAVG = (TUP + 2.0*TM + TDN) / 4.0
5048	ENDIF
5049	ENDIF
5050	ENDIF
5051	! ----------------------------------------------------------------------
5052	! END SUBROUTINE TMPAVG
5053	! ----------------------------------------------------------------------
5054	END SUBROUTINE TMPAVG
5055
5056	SUBROUTINE TRANSP (ET1,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
5057	& CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
5058
5059	IMPLICIT NONE
5060
5061	! ----------------------------------------------------------------------
5062	! SUBROUTINE TRANSP
5063	! ----------------------------------------------------------------------
5064	! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
5065	! ----------------------------------------------------------------------
5066	INTEGER I
5067	INTEGER K
5068	INTEGER NSOIL
5069	INTEGER NROOT
5070
5071	REAL CFACTR
5072	REAL CMC
5073	REAL CMCMAX
5074	REAL DENOM
5075	REAL ET1(NSOIL)
5076	REAL ETP1
5077	REAL ETP1A
5078	REAL GX (7)
5079	!.....REAL PART(NSOIL)
5080	REAL PC
5081	REAL Q2
5082	REAL RTDIS(NSOIL)
5083	REAL RTX
5084	REAL SFCTMP
5085	REAL SGX
5086	REAL SHDFAC
5087	REAL SMC(NSOIL)
5088	REAL SMCREF
5089	REAL SMCWLT
5090	REAL ZSOIL(NSOIL)
5091
5092	! ----------------------------------------------------------------------
5093	! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
5094	! ----------------------------------------------------------------------
5095	DO K = 1,NSOIL
5096	ET1(K) = 0.
5097	END DO
5098
5099	! ----------------------------------------------------------------------
5100	! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
5101	! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5102	! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
5103	! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
5104	! TOTAL ETP1A.
5105	! ----------------------------------------------------------------------
5106	IF (CMC .NE. 0.0) THEN
5107	ETP1A = SHDFAC * PC * ETP1 * (1.0 - (CMC /CMCMAX) ** CFACTR)
5108	ELSE
5109	ETP1A = SHDFAC * PC * ETP1
5110	ENDIF
5111
5112	SGX = 0.0
5113	DO I = 1,NROOT
5114	GX(I) = ( SMC(I) - SMCWLT ) / ( SMCREF - SMCWLT )
5115	GX(I) = MAX ( MIN ( GX(I), 1. ), 0. )
5116	SGX = SGX + GX (I)
5117	END DO
5118	SGX = SGX / NROOT
5119
5120	DENOM = 0.
5121	DO I = 1,NROOT
5122	RTX = RTDIS(I) + GX(I) - SGX
5123	GX(I) = GX(I) * MAX ( RTX, 0. )
5124	DENOM = DENOM + GX(I)
5125	END DO
5126	IF (DENOM .LE. 0.0) DENOM = 1.
5127
5128	DO I = 1,NROOT
5129	ET1(I) = ETP1A * GX(I) / DENOM
5130	END DO
5131
5132	! ----------------------------------------------------------------------
5133	! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
5134	! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
5135	! ----------------------------------------------------------------------
5136	! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
5137	! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
5138	! ----------------------------------------------------------------------
5139	! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5140	! ----------------------------------------------------------------------
5141	! ET(1) = RTDIS(1) * ETP1A
5142	! ET(1) = ETP1A * PART(1)
5143	! ----------------------------------------------------------------------
5144	! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
5145	! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
5146	! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
5147	! ----------------------------------------------------------------------
5148	! DO K = 2,NROOT
5149	! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
5150	! GX = MAX ( MIN ( GX, 1. ), 0. )
5151	! TEST CANOPY RESISTANCE
5152	! GX = 1.0
5153	! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))GXETP1A
5154	! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
5155	! ----------------------------------------------------------------------
5156	! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5157	! ----------------------------------------------------------------------
5158	! ET(K) = RTDIS(K) * ETP1A
5159	! ET(K) = ETP1A*PART(K)
5160	! END DO
5161	! ----------------------------------------------------------------------
5162	! END SUBROUTINE TRANSP
5163	! ----------------------------------------------------------------------
5164	END SUBROUTINE TRANSP
5165
5166	SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
5167	& SICEMAX)
5168
5169	IMPLICIT NONE
5170
5171	! ----------------------------------------------------------------------
5172	! SUBROUTINE WDFCND
5173	! ----------------------------------------------------------------------
5174	! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
5175	! ----------------------------------------------------------------------
5176	REAL BEXP
5177	REAL DKSAT
5178	REAL DWSAT
5179	REAL EXPON
5180	REAL FACTR1
5181	REAL FACTR2
5182	REAL SICEMAX
5183	REAL SMC
5184	REAL SMCMAX
5185	REAL VKwgt
5186	REAL WCND
5187	REAL WDF
5188
5189	! ----------------------------------------------------------------------
5190	! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
5191	! ----------------------------------------------------------------------
5192	SMC = SMC
5193	SMCMAX = SMCMAX
5194	FACTR1 = 0.2 / SMCMAX
5195	FACTR2 = SMC / SMCMAX
5196
5197	! ----------------------------------------------------------------------
5198	! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
5199	! ----------------------------------------------------------------------
5200	EXPON = BEXP + 2.0
5201	WDF = DWSAT * FACTR2 ** EXPON
5202
5203	! ----------------------------------------------------------------------
5204	! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
5205	! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
5206	! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
5207	! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
5208	! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
5209	! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
5210	! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
5211	! --
5212	! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
5213	! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
5214	! ----------------------------------------------------------------------
5215	IF (SICEMAX .GT. 0.0) THEN
5216	VKWGT = 1./(1.+(500.SICEMAX)*3.)
5217	WDF = VKWGTWDF + (1.- VKWGT)DWSATFACTR1*EXPON
5218	ENDIF
5219
5220	! ----------------------------------------------------------------------
5221	! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
5222	! ----------------------------------------------------------------------
5223	EXPON = (2.0 * BEXP) + 3.0
5224	WCND = DKSAT * FACTR2 ** EXPON
5225
5226	! ----------------------------------------------------------------------
5227	! END SUBROUTINE WDFCND
5228	! ----------------------------------------------------------------------
5229	END SUBROUTINE WDFCND
5230
5231	SUBROUTINE nmmlsminit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
5232	SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
5233	ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, & ! STEMP
5234	TMN, &
5235	num_soil_layers, &
5236	allowed_to_read, &
5237	ids,ide, jds,jde, kds,kde, &
5238	ims,ime, jms,jme, kms,kme, &
5239	its,ite, jts,jte, kts,kte )
5240
5241	IMPLICIT NONE
5242
5243	! Arguments
5244	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
5245	ims,ime, jms,jme, kms,kme, &
5246	its,ite, jts,jte, kts,kte
5247
5248	INTEGER, INTENT(IN) :: num_soil_layers
5249
5250	REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
5251
5252	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
5253	INTENT(INOUT) :: SMOIS, &
5254	TSLB !STEMP
5255
5256	REAL, DIMENSION( ims:ime, jms:jme ) , &
5257	INTENT(INOUT) :: SNOW, &
5258	SNOWC, &
5259	CANWAT, &
5260	SMSTAV, &
5261	SMSTOT, &
5262	SFCRUNOFF, &
5263	UDRUNOFF, &
5264	SFCEVP, &
5265	GRDFLX, &
5266	ACSNOW, &
5267	XICE, &
5268	VEGFRA, &
5269	TMN, &
5270	ACSNOM
5271
5272	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
5273	INTENT(INOUT) :: IVGTYP, &
5274	ISLTYP
5275
5276	!
5277
5278	INTEGER, INTENT(IN) :: isn
5279	LOGICAL, INTENT(IN) :: allowed_to_read
5280	! Local
5281	INTEGER :: iseason
5282	INTEGER :: icm,jcm,itf,jtf
5283	INTEGER :: I,J,L
5284
5285
5286	itf=min0(ite,ide-1)
5287	jtf=min0(jte,jde-1)
5288
5289	icm = ide/2
5290	jcm = jde/2
5291
5292	iseason=isn
5293
5294	DO J=jts,jtf
5295	DO I=its,itf
5296	! SNOW(i,j)=0.
5297	SNOWC(i,j)=0.
5298	! SMSTAV(i,j)=
5299	! SMSTOT(i,j)=
5300	! SFCRUNOFF(i,j)=
5301	! UDRUNOFF(i,j)=
5302	! GRDFLX(i,j)=
5303	! ACSNOW(i,j)=
5304	! ACSNOM(i,j)=
5305	ENDDO
5306	ENDDO
5307
5308	END SUBROUTINE nmmlsminit
5309
5310	FUNCTION CSNOW (DSNOW)
5311
5312	IMPLICIT NONE
5313
5314	! ----------------------------------------------------------------------
5315	! FUNCTION CSNOW
5316	! ----------------------------------------------------------------------
5317	! CALCULATE SNOW TERMAL CONDUCTIVITY
5318	! ----------------------------------------------------------------------
5319	REAL C
5320	REAL DSNOW
5321	REAL CSNOW
5322	REAL UNIT
5323
5324	PARAMETER(UNIT = 0.11631)
5325
5326	! ----------------------------------------------------------------------
5327	! CSNOW IN UNITS OF CAL/(CMHRC), RETURNED IN W/(M*C)
5328	! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
5329	! ----------------------------------------------------------------------
5330	C=0.32810(2.25DSNOW)
5331	! CSNOW=UNIT*C
5332	! MEK JAN 2006, DOUBLE SNOW THERMAL CONDUCTIVITY
5333	CSNOW=2.0UNITC
5334
5335	! ----------------------------------------------------------------------
5336	! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
5337	! ----------------------------------------------------------------------
5338	! CSNOW=0.0293(1.+100.DSNOW**2)
5339
5340	! ----------------------------------------------------------------------
5341	! E. ANDERSEN FROM FLERCHINGER
5342	! ----------------------------------------------------------------------
5343	! CSNOW=0.021+2.51DSNOW*2
5344
5345	! ----------------------------------------------------------------------
5346	! END FUNCTION CSNOW
5347	! ----------------------------------------------------------------------
5348	END FUNCTION CSNOW
5349
5350	FUNCTION FRH2O (TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
5351
5352	IMPLICIT NONE
5353
5354	! ----------------------------------------------------------------------
5355	! FUNCTION FRH2O
5356	! ----------------------------------------------------------------------
5357	! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
5358	! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
5359	! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
5360	! (1999, JGR, VOL 104(D16), 19569-19585).
5361	! ----------------------------------------------------------------------
5362	! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
5363	! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
5364	! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
5365	! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
5366	! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
5367	! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
5368	! LIMIT OF FREEZING POINT TEMPERATURE T0.
5369	! ----------------------------------------------------------------------
5370	! INPUT:
5371	!
5372	! TKELV.........TEMPERATURE (Kelvin)
5373	! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
5374	! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
5375	! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
5376	! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
5377	! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
5378	!
5379	! OUTPUT:
5380	! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
5381	! ----------------------------------------------------------------------
5382	REAL BEXP
5383	REAL BLIM
5384	REAL BX
5385	REAL CK
5386	REAL DENOM
5387	REAL DF
5388	REAL DH2O
5389	REAL DICE
5390	REAL DSWL
5391	REAL ERROR
5392	REAL FK
5393	REAL FRH2O
5394	REAL GS
5395	REAL HLICE
5396	REAL PSIS
5397	REAL SH2O
5398	REAL SMC
5399	REAL SMCMAX
5400	REAL SWL
5401	REAL SWLK
5402	REAL TKELV
5403	REAL T0
5404
5405	INTEGER NLOG
5406	INTEGER KCOUNT
5407
5408	PARAMETER(CK = 8.0)
5409	! PARAMETER(CK = 0.0)
5410	PARAMETER(BLIM = 5.5)
5411	PARAMETER(ERROR = 0.005)
5412
5413	PARAMETER(HLICE = 3.335E5)
5414	PARAMETER(GS = 9.81)
5415	PARAMETER(DICE = 920.0)
5416	PARAMETER(DH2O = 1000.0)
5417	PARAMETER(T0 = 273.15)
5418
5419	! ----------------------------------------------------------------------
5420	! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
5421	! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
5422	! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
5423	! ----------------------------------------------------------------------
5424	BX = BEXP
5425	IF (BEXP .GT. BLIM) BX = BLIM
5426
5427	! ----------------------------------------------------------------------
5428	! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
5429	! ----------------------------------------------------------------------
5430	NLOG=0
5431	KCOUNT=0
5432
5433	! ----------------------------------------------------------------------
5434	! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
5435	! ----------------------------------------------------------------------
5436	IF (TKELV .GT. (T0 - 1.E-3)) THEN
5437	FRH2O = SMC
5438	ELSE
5439	IF (CK .NE. 0.0) THEN
5440
5441	! ----------------------------------------------------------------------
5442	! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
5443	! IN KOREN ET AL, JGR, 1999, EQN 17
5444	! ----------------------------------------------------------------------
5445	! INITIAL GUESS FOR SWL (frozen content)
5446	! ----------------------------------------------------------------------
5447	SWL = SMC-SH2O
5448
5449	! ----------------------------------------------------------------------
5450	! KEEP WITHIN BOUNDS.
5451	! ----------------------------------------------------------------------
5452	IF (SWL .GT. (SMC-0.02)) SWL = SMC-0.02
5453	IF (SWL .LT. 0.) SWL = 0.
5454
5455	! ----------------------------------------------------------------------
5456	! START OF ITERATIONS
5457	! ----------------------------------------------------------------------
5458	DO WHILE ( (NLOG .LT. 10) .AND. (KCOUNT .EQ. 0) )
5459	NLOG = NLOG+1
5460	DF = ALOG(( PSISGS/HLICE ) ( ( 1.+CKSWL )2. ) &
5461	& ( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
5462	DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
5463	SWLK = SWL - DF/DENOM
5464	! ----------------------------------------------------------------------
5465	! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
5466	! ----------------------------------------------------------------------
5467	IF (SWLK .GT. (SMC-0.02)) SWLK = SMC - 0.02
5468	IF (SWLK .LT. 0.) SWLK = 0.
5469
5470	! ----------------------------------------------------------------------
5471	! MATHEMATICAL SOLUTION BOUNDS APPLIED.
5472	! ----------------------------------------------------------------------
5473	DSWL = ABS(SWLK-SWL)
5474	SWL = SWLK
5475
5476	! ----------------------------------------------------------------------
5477	! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
5478	! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
5479	! ----------------------------------------------------------------------
5480	IF ( DSWL .LE. ERROR ) THEN
5481	KCOUNT = KCOUNT+1
5482	ENDIF
5483	END DO
5484
5485	! ----------------------------------------------------------------------
5486	! END OF ITERATIONS
5487	! ----------------------------------------------------------------------
5488	! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
5489	! ----------------------------------------------------------------------
5490	FRH2O = SMC - SWL
5491
5492	! ----------------------------------------------------------------------
5493	! END OPTION 1
5494	! ----------------------------------------------------------------------
5495	ENDIF
5496
5497	! ----------------------------------------------------------------------
5498	! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
5499	! IN KOREN ET AL., JGR, 1999, EQN 17
5500	! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
5501	! ----------------------------------------------------------------------
5502	IF (KCOUNT .EQ. 0) THEN
5503	Print*,'Flerchinger used in NEW version. Iterations=',NLOG
5504	FK = (((HLICE/(GS(-PSIS))) &
5505	& ((TKELV-T0)/TKELV))*(-1/BX))SMCMAX
5506	IF (FK .LT. 0.02) FK = 0.02
5507	FRH2O = MIN (FK, SMC)
5508	! ----------------------------------------------------------------------
5509	! END OPTION 2
5510	! ----------------------------------------------------------------------
5511	ENDIF
5512
5513	ENDIF
5514
5515	! ----------------------------------------------------------------------
5516	! END FUNCTION FRH2O
5517	! ----------------------------------------------------------------------
5518	END FUNCTION FRH2O
5519
5520	FUNCTION SNKSRC (TAVG,SMC,SH2O,ZSOIL,NSOIL, &
5521	& SMCMAX,PSISAT,BEXP,DT,K,QTOT)
5522
5523	IMPLICIT NONE
5524
5525	! ----------------------------------------------------------------------
5526	! FUNCTION SNKSRC
5527	! ----------------------------------------------------------------------
5528	! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
5529	! AVAILABLE LIQUED WATER.
5530	! ----------------------------------------------------------------------
5531	INTEGER K
5532	INTEGER NSOIL
5533
5534	REAL BEXP
5535	REAL DF
5536	REAL DH2O
5537	REAL DT
5538	REAL DZ
5539	REAL DZH
5540	REAL FREE
5541	! REAL FRH2O
5542	REAL HLICE
5543	REAL PSISAT
5544	REAL QTOT
5545	REAL SH2O
5546	REAL SMC
5547	REAL SMCMAX
5548	REAL SNKSRC
5549	REAL T0
5550	REAL TAVG
5551	REAL TDN
5552	REAL TM
5553	REAL TUP
5554	REAL TZ
5555	REAL X0
5556	REAL XDN
5557	REAL XH2O
5558	REAL XUP
5559	REAL ZSOIL (NSOIL)
5560
5561	PARAMETER(DH2O = 1.0000E3)
5562	PARAMETER(HLICE = 3.3350E5)
5563	PARAMETER(T0 = 2.7315E2)
5564
5565	IF (K .EQ. 1) THEN
5566	DZ = -ZSOIL(1)
5567	ELSE
5568	DZ = ZSOIL(K-1)-ZSOIL(K)
5569	ENDIF
5570
5571	! ----------------------------------------------------------------------
5572	! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
5573	! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
5574	! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
5575	! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
5576	! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
5577	! ----------------------------------------------------------------------
5578	FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
5579
5580	! ----------------------------------------------------------------------
5581	! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
5582	! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
5583	! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
5584	! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
5585	! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
5586	! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
5587	! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
5588	! ----------------------------------------------------------------------
5589	XH2O = SH2O + QTOTDT/(DH2OHLICE*DZ)
5590
5591	! ----------------------------------------------------------------------
5592	! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
5593	! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
5594	! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
5595	! ----------------------------------------------------------------------
5596	IF ( XH2O .LT. SH2O .AND. XH2O .LT. FREE) THEN
5597	IF ( FREE .GT. SH2O ) THEN
5598	XH2O = SH2O
5599	ELSE
5600	XH2O = FREE
5601	ENDIF
5602	ENDIF
5603
5604	! ----------------------------------------------------------------------
5605	! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
5606	! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
5607	! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
5608	! ----------------------------------------------------------------------
5609	IF ( XH2O .GT. SH2O .AND. XH2O .GT. FREE ) THEN
5610	IF ( FREE .LT. SH2O ) THEN
5611	XH2O = SH2O
5612	ELSE
5613	XH2O = FREE
5614	ENDIF
5615	ENDIF
5616
5617	IF (XH2O .LT. 0.) XH2O = 0.
5618	IF (XH2O .GT. SMC) XH2O = SMC
5619
5620	! ----------------------------------------------------------------------
5621	! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
5622	! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
5623	! ----------------------------------------------------------------------
5624	SNKSRC = -DH2OHLICEDZ*(XH2O-SH2O)/DT
5625	SH2O = XH2O
5626
5627	! ----------------------------------------------------------------------
5628	! END FUNCTION SNKSRC
5629	! ----------------------------------------------------------------------
5630	END FUNCTION SNKSRC
5631
5632	END MODULE module_sf_lsm_nmm
5633

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

Download in other formats:

Original Format