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 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 220.5 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	BETA = ETA/ETP
1382	! ----------------------------------------------------------------------
1383
1384	! ----------------------------------------------------------------------
1385	! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
1386	! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
1387	! SSOIL<0: COOL THE SURFACE (DAY TIME)
1388	! ----------------------------------------------------------------------
1389	SSOIL = -1.0*SSOIL
1390
1391	! ----------------------------------------------------------------------
1392	! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
1393	! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW
1394	! ----------------------------------------------------------------------
1395	RUNOFF3 = RUNOFF3/DT
1396	RUNOFF2 = RUNOFF2+RUNOFF3
1397
1398	! ----------------------------------------------------------------------
1399	! TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE
1400	! SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION
1401	! ----------------------------------------------------------------------
1402	SOILM = -1.0SMC(1)ZSOIL(1)
1403	DO K = 2,NSOIL
1404	SOILM = SOILM+SMC(K)*(ZSOIL(K-1)-ZSOIL(K))
1405	END DO
1406	SOILWM = -1.0(SMCMAX-SMCWLT)ZSOIL(1)
1407	SOILWW = -1.0(SMC(1)-SMCWLT)ZSOIL(1)
1408	DO K = 2,NROOT
1409	SOILWM = SOILWM+(SMCMAX-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1410	SOILWW = SOILWW+(SMC(K)-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1411	END DO
1412	SOILW = SOILWW/SOILWM
1413
1414	! ----------------------------------------------------------------------
1415	! END SUBROUTINE SFLX
1416	! ----------------------------------------------------------------------
1417	END SUBROUTINE SFLX
1418
1419	SUBROUTINE ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1420
1421	IMPLICIT NONE
1422
1423	! ----------------------------------------------------------------------
1424	! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
1425	! ALB SNOWFREE ALBEDO
1426	! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
1427	! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1428	! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1429	! SNCOVR FRACTIONAL SNOW COVER
1430	! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
1431	! TSNOW SNOW SURFACE TEMPERATURE (K)
1432	! ----------------------------------------------------------------------
1433	REAL ALB, SNOALB, SHDFAC, SHDMIN, SNCOVR, ALBEDO, TSNOW
1434
1435	! ----------------------------------------------------------------------
1436	! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
1437	! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
1438	! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
1439	! (1985, JCAM, VOL 24, 402-411)
1440	! ----------------------------------------------------------------------
1441	! changed in version 2.6 on June 2nd 2003
1442	! ALBEDO = ALB + (1.0-(SHDFAC-SHDMIN))SNCOVR(SNOALB-ALB)
1443	ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
1444	IF (ALBEDO .GT. SNOALB) ALBEDO=SNOALB
1445
1446	! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
1447	! IF (TSNOW.LE.263.16) THEN
1448	! ALBEDO=SNOALB
1449	! ELSE
1450	! IF (TSNOW.LT.273.16) THEN
1451	! TM=0.1*(TSNOW-263.16)
1452	! ALBEDO=0.5((0.9-0.2(TM*3))+(0.8-0.16(TM**3)))
1453	! ELSE
1454	! ALBEDO=0.67
1455	! ENDIF
1456	! ENDIF
1457
1458	! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
1459	! IF (TSNOW.LT.273.16) THEN
1460	! ALBEDO=SNOALB-0.008*DT/86400
1461	! ELSE
1462	! ALBEDO=(SNOALB-0.5)EXP(-0.24DT/86400)+0.5
1463	! ENDIF
1464
1465	! ----------------------------------------------------------------------
1466	! END SUBROUTINE ALCALC
1467	! ----------------------------------------------------------------------
1468	END SUBROUTINE ALCALC
1469
1470	SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
1471	& SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1472	& TOPT,RSMAX,RGL,HS,XLAI, &
1473	& RCS,RCT,RCQ,RCSOIL)
1474
1475	IMPLICIT NONE
1476
1477	! ----------------------------------------------------------------------
1478	! SUBROUTINE CANRES
1479	! ----------------------------------------------------------------------
1480	! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
1481	! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
1482	! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
1483	! MOISTURE RATHER THAN TOTAL)
1484	! ----------------------------------------------------------------------
1485	! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
1486	! NOILHAN (1990, BLM)
1487	! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
1488	! AND TABLE 2 OF SEC. 3.1.2
1489	! ----------------------------------------------------------------------
1490	! INPUT:
1491	! SOLAR INCOMING SOLAR RADIATION
1492	! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1493	! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
1494	! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1495	! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1496	! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
1497	! SFCPRS SURFACE PRESSURE
1498	! SMC VOLUMETRIC SOIL MOISTURE
1499	! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
1500	! NSOIL NO. OF SOIL LAYERS
1501	! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
1502	! XLAI LEAF AREA INDEX
1503	! SMCWLT WILTING POINT
1504	! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
1505	! SETS IN)
1506	! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
1507	! SURBOUTINE REDPRM
1508	! OUTPUT:
1509	! PC PLANT COEFFICIENT
1510	! RC CANOPY RESISTANCE
1511	! ----------------------------------------------------------------------
1512	INTEGER NSOLD
1513	PARAMETER(NSOLD = 20)
1514
1515	INTEGER K
1516	INTEGER NROOT
1517	INTEGER NSOIL
1518
1519	REAL CH
1520	REAL CP
1521	REAL DELTA
1522	REAL DQSDT2
1523	REAL FF
1524	REAL GX
1525	REAL HS
1526	REAL P
1527	REAL PART(NSOLD)
1528	REAL PC
1529	REAL Q2
1530	REAL Q2SAT
1531	REAL RC
1532	REAL RSMIN
1533	REAL RCQ
1534	REAL RCS
1535	REAL RCSOIL
1536	REAL RCT
1537	REAL RD
1538	REAL RGL
1539	REAL RR
1540	REAL RSMAX
1541	REAL SFCPRS
1542	REAL SFCTMP
1543	REAL SIGMA
1544	REAL SLV
1545	REAL SMC(NSOIL)
1546	REAL SMCREF
1547	REAL SMCWLT
1548	REAL SOLAR
1549	REAL TOPT
1550	REAL SLVCP
1551	REAL ST1
1552	REAL TAIR4
1553	REAL XLAI
1554	REAL ZSOIL(NSOIL)
1555
1556	PARAMETER(CP = 1004.5)
1557	PARAMETER(RD = 287.04)
1558	PARAMETER(SIGMA = 5.67E-8)
1559	PARAMETER(SLV = 2.501000E6)
1560
1561	! ----------------------------------------------------------------------
1562	! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
1563	! ----------------------------------------------------------------------
1564	RCS = 0.0
1565	RCT = 0.0
1566	RCQ = 0.0
1567	RCSOIL = 0.0
1568	RC = 0.0
1569
1570	! ----------------------------------------------------------------------
1571	! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
1572	! ----------------------------------------------------------------------
1573	FF = 0.552.0SOLAR/(RGL*XLAI)
1574	RCS = (FF + RSMIN/RSMAX) / (1.0 + FF)
1575	RCS = MAX(RCS,0.0001)
1576
1577	! ----------------------------------------------------------------------
1578	! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
1579	! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
1580	! ----------------------------------------------------------------------
1581	RCT = 1.0 - 0.0016((TOPT-SFCTMP)*2.0)
1582	RCT = MAX(RCT,0.0001)
1583
1584	! ----------------------------------------------------------------------
1585	! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
1586	! RCQ EXPRESSION FROM SSIB
1587	! ----------------------------------------------------------------------
1588	RCQ = 1.0/(1.0+HS*(Q2SAT-Q2))
1589	RCQ = MAX(RCQ,0.01)
1590
1591	! ----------------------------------------------------------------------
1592	! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
1593	! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
1594	! ----------------------------------------------------------------------
1595	GX = (SMC(1) - SMCWLT) / (SMCREF - SMCWLT)
1596	IF (GX .GT. 1.) GX = 1.
1597	IF (GX .LT. 0.) GX = 0.
1598
1599	! ----------------------------------------------------------------------
1600	! USE SOIL DEPTH AS WEIGHTING FACTOR
1601	! ----------------------------------------------------------------------
1602	PART(1) = (ZSOIL(1)/ZSOIL(NROOT)) * GX
1603	! ----------------------------------------------------------------------
1604	! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1605	! PART(1) = RTDIS(1) * GX
1606	! ----------------------------------------------------------------------
1607	IF (NROOT .GT. 1) THEN
1608	DO K = 2,NROOT
1609	GX = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT)
1610	IF (GX .GT. 1.) GX = 1.
1611	IF (GX .LT. 0.) GX = 0.
1612	! ----------------------------------------------------------------------
1613	! USE SOIL DEPTH AS WEIGHTING FACTOR
1614	! ----------------------------------------------------------------------
1615	PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT)) * GX
1616	! ----------------------------------------------------------------------
1617	! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1618	! PART(K) = RTDIS(K) * GX
1619	! ----------------------------------------------------------------------
1620	END DO
1621	ENDIF
1622
1623	DO K = 1,NROOT
1624	RCSOIL = RCSOIL+PART(K)
1625	END DO
1626	RCSOIL = MAX(RCSOIL,0.0001)
1627
1628	! ----------------------------------------------------------------------
1629	! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
1630	! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
1631	! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
1632	! PC * LINERIZED PENMAN POTENTIAL EVAP =
1633	! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
1634	! ----------------------------------------------------------------------
1635	RC = RSMIN/(XLAIRCSRCTRCQRCSOIL)
1636
1637	! TAIR4 = SFCTMP**4.
1638	! ST1 = (4.SIGMARD)/CP
1639	! SLVCP = SLV/CP
1640	! RR = ST1TAIR4/(SFCPRSCH) + 1.0
1641	RR = (4.SIGMARD/CP)(SFCTMP4.)/(SFCPRSCH) + 1.0
1642	DELTA = (SLV/CP)*DQSDT2
1643
1644	PC = (RR+DELTA)/(RR(1.+RCCH)+DELTA)
1645
1646	! ----------------------------------------------------------------------
1647	! END SUBROUTINE CANRES
1648	! ----------------------------------------------------------------------
1649	END SUBROUTINE CANRES
1650
1651	SUBROUTINE DEVAP (EDIR1,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
1652	& DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1653
1654	IMPLICIT NONE
1655
1656	! ----------------------------------------------------------------------
1657	! SUBROUTINE DEVAP
1658	! ----------------------------------------------------------------------
1659	! CALCULATE DIRECT SOIL EVAPORATION
1660	! ----------------------------------------------------------------------
1661	REAL BEXP
1662	! REAL DEVAP
1663	REAL EDIR1
1664	REAL DKSAT
1665	REAL DWSAT
1666	REAL ETP1
1667	REAL FX
1668	REAL FXEXP
1669	REAL SHDFAC
1670	REAL SMC
1671	REAL SMCDRY
1672	REAL SMCMAX
1673	REAL ZSOIL
1674	REAL SMCREF
1675	REAL SMCWLT
1676	REAL SRATIO
1677
1678	! ----------------------------------------------------------------------
1679	! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
1680	! WHEN FXEXP=1.
1681	! FX > 1 REPRESENTS DEMAND CONTROL
1682	! FX < 1 REPRESENTS FLUX CONTROL
1683	! ----------------------------------------------------------------------
1684	SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
1685	IF (SRATIO .GT. 0.) THEN
1686	FX = SRATIO**FXEXP
1687	FX = MAX ( MIN ( FX, 1. ) ,0. )
1688	ELSE
1689	FX = 0.
1690	ENDIF
1691
1692	! ----------------------------------------------------------------------
1693	! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
1694	! ----------------------------------------------------------------------
1695	! DEVAP = FX * ( 1.0 - SHDFAC ) * ETP1
1696	EDIR1 = FX * ( 1.0 - SHDFAC ) * ETP1
1697
1698	! ----------------------------------------------------------------------
1699	! END SUBROUTINE DEVAP
1700	! ----------------------------------------------------------------------
1701	END SUBROUTINE DEVAP
1702
1703	SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
1704	& SH2O, &
1705	& SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
1706	& SMCREF,SHDFAC,CMCMAX, &
1707	& SMCDRY,CFACTR, &
1708	& EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1709
1710	IMPLICIT NONE
1711
1712	! ----------------------------------------------------------------------
1713	! SUBROUTINE EVAPO
1714	! ----------------------------------------------------------------------
1715	! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
1716	! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
1717	! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
1718	! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
1719	! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
1720	! ----------------------------------------------------------------------
1721	INTEGER NSOLD
1722	PARAMETER(NSOLD = 20)
1723
1724	INTEGER I
1725	INTEGER K
1726	INTEGER NSOIL
1727	INTEGER NROOT
1728
1729	REAL BEXP
1730	REAL CFACTR
1731	REAL CMC
1732	REAL CMC2MS
1733	REAL CMCMAX
1734	! REAL DEVAP
1735	REAL DKSAT
1736	REAL DT
1737	REAL DWSAT
1738	REAL EC1
1739	REAL EDIR1
1740	REAL ET1(NSOIL)
1741	REAL ETA1
1742	REAL ETP1
1743	REAL ETT1
1744	REAL FXEXP
1745	REAL PC
1746	REAL Q2
1747	REAL RTDIS(NSOIL)
1748	REAL SFCTMP
1749	REAL SHDFAC
1750	REAL SMC(NSOIL)
1751	REAL SH2O(NSOIL)
1752	REAL SMCDRY
1753	REAL SMCMAX
1754	REAL SMCREF
1755	REAL SMCWLT
1756	REAL ZSOIL(NSOIL)
1757
1758	! ----------------------------------------------------------------------
1759	! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
1760	! GREATER THAN ZERO.
1761	! ----------------------------------------------------------------------
1762	EDIR1 = 0.
1763	EC1 = 0.
1764	DO K = 1,NSOIL
1765	ET1(K) = 0.
1766	END DO
1767	ETT1 = 0.
1768
1769	IF (ETP1 .GT. 0.0) THEN
1770
1771	! ----------------------------------------------------------------------
1772	! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
1773	! ONLY IF VEG COVER NOT COMPLETE.
1774	! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
1775	! ----------------------------------------------------------------------
1776	IF (SHDFAC .LT. 1.) THEN
1777	CALL DEVAP (EDIR1,ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1778	! EDIR = DEVAP(ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1779	& BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1780	ENDIF
1781
1782	! ----------------------------------------------------------------------
1783	! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
1784	! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
1785	! ----------------------------------------------------------------------
1786	IF (SHDFAC.GT.0.0) THEN
1787
1788	CALL TRANSP (ET1,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
1789	& CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
1790
1791	DO K = 1,NSOIL
1792	ETT1 = ETT1 + ET1(K)
1793	END DO
1794
1795	! ----------------------------------------------------------------------
1796	! CALCULATE CANOPY EVAPORATION.
1797	! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
1798	! ----------------------------------------------------------------------
1799	IF (CMC .GT. 0.0) THEN
1800	EC1 = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
1801	ELSE
1802	EC1 = 0.0
1803	ENDIF
1804
1805	! ----------------------------------------------------------------------
1806	! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
1807	! CANOPY. -F.CHEN, 18-OCT-1994
1808	! ----------------------------------------------------------------------
1809	CMC2MS = CMC / DT
1810	EC1 = MIN ( CMC2MS, EC1 )
1811	ENDIF
1812	ENDIF
1813
1814	! ----------------------------------------------------------------------
1815	! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
1816	! ----------------------------------------------------------------------
1817	ETA1 = EDIR1 + ETT1 + EC1
1818
1819	! ----------------------------------------------------------------------
1820	! END SUBROUTINE EVAPO
1821	! ----------------------------------------------------------------------
1822	END SUBROUTINE EVAPO
1823
1824	SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
1825	& TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
1826	& F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
1827
1828	IMPLICIT NONE
1829
1830	! ----------------------------------------------------------------------
1831	! SUBROUTINE HRT
1832	! ----------------------------------------------------------------------
1833	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1834	! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
1835	! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1836	! ----------------------------------------------------------------------
1837	INTEGER NSOLD
1838	PARAMETER(NSOLD = 20)
1839
1840	LOGICAL ITAVG
1841
1842	INTEGER I
1843	INTEGER K
1844	INTEGER NSOIL
1845
1846	! ----------------------------------------------------------------------
1847	! DECLARE WORK ARRAYS NEEDED IN TRI-DIAGONAL IMPLICIT SOLVER
1848	! ----------------------------------------------------------------------
1849	REAL AI(NSOLD)
1850	REAL BI(NSOLD)
1851	REAL CI(NSOLD)
1852
1853	! ----------------------------------------------------------------------
1854	! DECLARATIONS
1855	! ----------------------------------------------------------------------
1856	REAL BEXP
1857	REAL CAIR
1858	REAL CH2O
1859	REAL CICE
1860	REAL CSOIL
1861	REAL DDZ
1862	REAL DDZ2
1863	REAL DENOM
1864	REAL DF1
1865	REAL DF1N
1866	REAL DF1K
1867	REAL DT
1868	REAL DTSDZ
1869	REAL DTSDZ2
1870	REAL F1
1871	REAL HCPCT
1872	REAL PSISAT
1873	REAL QUARTZ
1874	REAL QTOT
1875	REAL RHSTS(NSOIL)
1876	REAL SSOIL
1877	REAL SICE
1878	REAL SMC(NSOIL)
1879	REAL SH2O(NSOIL)
1880	REAL SMCMAX
1881	! REAL SNKSRC
1882	REAL STC(NSOIL)
1883	REAL T0
1884	REAL TAVG
1885	REAL TBK
1886	REAL TBK1
1887	REAL TBOT
1888	REAL ZBOT
1889	REAL TSNSR
1890	REAL TSURF
1891	REAL YY
1892	REAL ZSOIL(NSOIL)
1893	REAL ZZ1
1894
1895	PARAMETER(T0 = 273.15)
1896
1897	! ----------------------------------------------------------------------
1898	! SET SPECIFIC HEAT CAPACITIES OF AIR, WATER, ICE, SOIL MINERAL
1899	! ----------------------------------------------------------------------
1900	PARAMETER(CAIR = 1004.0)
1901	PARAMETER(CH2O = 4.2E6)
1902	PARAMETER(CICE = 2.106E6)
1903	! NOTE: CSOIL NOW SET IN ROUTINE REDPRM AND PASSED IN
1904	! PARAMETER(CSOIL = 1.26E6)
1905
1906	! ----------------------------------------------------------------------
1907	! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
1908	! ----------------------------------------------------------------------
1909	ITAVG = .TRUE.
1910	! ITAVG = .FALSE.
1911
1912	! ----------------------------------------------------------------------
1913	! BEGIN SECTION FOR TOP SOIL LAYER
1914	! ----------------------------------------------------------------------
1915	! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
1916	! ----------------------------------------------------------------------
1917	HCPCT = SH2O(1)CH2O + (1.0-SMCMAX)CSOIL + (SMCMAX-SMC(1))*CAIR &
1918	& + ( SMC(1) - SH2O(1) )*CICE
1919
1920	! ----------------------------------------------------------------------
1921	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1922	! ----------------------------------------------------------------------
1923	DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
1924	AI(1) = 0.0
1925	CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
1926	BI(1) = -CI(1) + DF1 / (0.5 * ZSOIL(1) * ZSOIL(1)HCPCTZZ1)
1927
1928	! ----------------------------------------------------------------------
1929	! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
1930	! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
1931	! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
1932	! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
1933	! ----------------------------------------------------------------------
1934	DTSDZ = (STC(1) - STC(2)) / (-0.5 * ZSOIL(2))
1935	SSOIL = DF1 * (STC(1) - YY) / (0.5 * ZSOIL(1) * ZZ1)
1936	RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
1937
1938	! ----------------------------------------------------------------------
1939	! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
1940	! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
1941	! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
1942	! ----------------------------------------------------------------------
1943	QTOT = SSOIL - DF1*DTSDZ
1944
1945	! ----------------------------------------------------------------------
1946	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
1947	! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
1948	! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
1949	! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
1950	! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
1951	! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
1952	! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
1953	! LATER IN FUNCTION SUBROUTINE SNKSRC
1954	! ----------------------------------------------------------------------
1955	IF (ITAVG) THEN
1956	TSURF = (YY + (ZZ1-1) * STC(1)) / ZZ1
1957	CALL TBND (STC(1),STC(2),ZSOIL,ZBOT,1,NSOIL,TBK)
1958	ENDIF
1959
1960	! ----------------------------------------------------------------------
1961	! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
1962	! ----------------------------------------------------------------------
1963	SICE = SMC(1) - SH2O(1)
1964
1965	! ----------------------------------------------------------------------
1966	! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
1967	! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
1968	! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
1969	! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
1970	! ----------------------------------------------------------------------
1971	IF ( (SICE .GT. 0.) .OR. (TSURF .LT. T0) .OR. &
1972	& (STC(1) .LT. T0) .OR. (TBK .LT. T0) ) THEN
1973
1974	IF (ITAVG) THEN
1975	CALL TMPAVG(TAVG,TSURF,STC(1),TBK,ZSOIL,NSOIL,1)
1976	ELSE
1977	TAVG = STC(1)
1978	ENDIF
1979	TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1), &
1980	& ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1981
1982	RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1983	ENDIF
1984
1985	! ----------------------------------------------------------------------
1986	! THIS ENDS SECTION FOR TOP SOIL LAYER.
1987	! ----------------------------------------------------------------------
1988	! INITIALIZE DDZ2
1989	! ----------------------------------------------------------------------
1990	DDZ2 = 0.0
1991
1992	! ----------------------------------------------------------------------
1993	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1994	! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
1995	! ----------------------------------------------------------------------
1996	DF1K = DF1
1997	DO K = 2,NSOIL
1998
1999	! ----------------------------------------------------------------------
2000	! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
2001	! ----------------------------------------------------------------------
2002	HCPCT = SH2O(K)CH2O +(1.0-SMCMAX)CSOIL +(SMCMAX-SMC(K))*CAIR &
2003	& + ( SMC(K) - SH2O(K) )*CICE
2004
2005	IF (K .NE. NSOIL) THEN
2006	! ----------------------------------------------------------------------
2007	! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
2008	! ----------------------------------------------------------------------
2009	! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
2010	! ----------------------------------------------------------------------
2011	CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2012
2013	! ----------------------------------------------------------------------
2014	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
2015	! ----------------------------------------------------------------------
2016	DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2017	DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2018
2019	! ----------------------------------------------------------------------
2020	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
2021	! ----------------------------------------------------------------------
2022	DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2023	CI(K) = -DF1N * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2024
2025	! ----------------------------------------------------------------------
2026	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2027	! TEMP AT BOTTOM OF LAYER.
2028	! ----------------------------------------------------------------------
2029	IF (ITAVG) THEN
2030	CALL TBND (STC(K),STC(K+1),ZSOIL,ZBOT,K,NSOIL,TBK1)
2031	ENDIF
2032	ELSE
2033
2034	! ----------------------------------------------------------------------
2035	! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
2036	! BOTTOM LAYER.
2037	! ----------------------------------------------------------------------
2038	CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2039
2040	! ----------------------------------------------------------------------
2041	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
2042	! ----------------------------------------------------------------------
2043	DENOM = .5 * (ZSOIL(K-1) + ZSOIL(K)) - ZBOT
2044	DTSDZ2 = (STC(K)-TBOT) / DENOM
2045
2046	! ----------------------------------------------------------------------
2047	! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
2048	! ----------------------------------------------------------------------
2049	CI(K) = 0.
2050
2051	! ----------------------------------------------------------------------
2052	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2053	! TEMP AT BOTTOM OF LAST LAYER.
2054	! ----------------------------------------------------------------------
2055	IF (ITAVG) THEN
2056	CALL TBND (STC(K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
2057	ENDIF
2058
2059	ENDIF
2060	! ----------------------------------------------------------------------
2061	! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
2062	! ----------------------------------------------------------------------
2063	! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2064	! ----------------------------------------------------------------------
2065	DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2066	RHSTS(K) = ( DF1N * DTSDZ2 - DF1K * DTSDZ ) / DENOM
2067	QTOT = -1.0DENOMRHSTS(K)
2068	SICE = SMC(K) - SH2O(K)
2069
2070	IF ( (SICE .GT. 0.) .OR. (TBK .LT. T0) .OR. &
2071	& (STC(K) .LT. T0) .OR. (TBK1 .LT. T0) ) THEN
2072
2073	IF (ITAVG) THEN
2074	CALL TMPAVG(TAVG,TBK,STC(K),TBK1,ZSOIL,NSOIL,K)
2075	ELSE
2076	TAVG = STC(K)
2077	ENDIF
2078	TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL, &
2079	& SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2080	RHSTS(K) = RHSTS(K) - TSNSR / DENOM
2081	ENDIF
2082
2083	! ----------------------------------------------------------------------
2084	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2085	! ----------------------------------------------------------------------
2086	AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2087	BI(K) = -(AI(K) + CI(K))
2088
2089	! ----------------------------------------------------------------------
2090	! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
2091	! ----------------------------------------------------------------------
2092	TBK = TBK1
2093	DF1K = DF1N
2094	DTSDZ = DTSDZ2
2095	DDZ = DDZ2
2096	END DO
2097
2098	! ----------------------------------------------------------------------
2099	! END SUBROUTINE HRT
2100	! ----------------------------------------------------------------------
2101	END SUBROUTINE HRT
2102
2103	SUBROUTINE HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2104
2105	IMPLICIT NONE
2106
2107	! ----------------------------------------------------------------------
2108	! SUBROUTINE HRTICE
2109	! ----------------------------------------------------------------------
2110	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2111	! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE PACK. ALSO TO
2112	! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
2113	! OF THE IMPLICIT TIME SCHEME.
2114	! ----------------------------------------------------------------------
2115	INTEGER NSOLD
2116	PARAMETER(NSOLD = 20)
2117
2118	INTEGER K
2119	INTEGER NSOIL
2120
2121	REAL AI(NSOLD)
2122	REAL BI(NSOLD)
2123	REAL CI(NSOLD)
2124
2125	REAL DDZ
2126	REAL DDZ2
2127	REAL DENOM
2128	REAL DF1
2129	REAL DTSDZ
2130	REAL DTSDZ2
2131	REAL HCPCT
2132	REAL RHSTS(NSOIL)
2133	REAL SSOIL
2134	REAL STC(NSOIL)
2135	REAL TBOT
2136	REAL YY
2137	REAL ZBOT
2138	REAL ZSOIL(NSOIL)
2139	REAL ZZ1
2140
2141	DATA TBOT /271.16/
2142
2143	! ----------------------------------------------------------------------
2144	! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
2145	! HCPCT = 1880.0*917.0.
2146	! ----------------------------------------------------------------------
2147	PARAMETER(HCPCT = 1.72396E+6)
2148
2149	! ----------------------------------------------------------------------
2150	! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
2151	! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
2152	! ----------------------------------------------------------------------
2153	! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
2154	! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
2155	! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
2156	! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
2157	! ----------------------------------------------------------------------
2158	ZBOT = ZSOIL(NSOIL)
2159
2160	! ----------------------------------------------------------------------
2161	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2162	! ----------------------------------------------------------------------
2163	DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
2164	AI(1) = 0.0
2165	CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
2166	BI(1) = -CI(1) + DF1/(0.5 * ZSOIL(1) * ZSOIL(1) * HCPCT * ZZ1)
2167
2168	! ----------------------------------------------------------------------
2169	! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
2170	! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
2171	! RHSTS FOR THE TOP SOIL LAYER.
2172	! ----------------------------------------------------------------------
2173	DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) )
2174	SSOIL = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 )
2175	RHSTS(1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL(1) * HCPCT )
2176
2177	! ----------------------------------------------------------------------
2178	! INITIALIZE DDZ2
2179	! ----------------------------------------------------------------------
2180	DDZ2 = 0.0
2181
2182	! ----------------------------------------------------------------------
2183	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2184	! ----------------------------------------------------------------------
2185	DO K = 2,NSOIL
2186	IF (K .NE. NSOIL) THEN
2187
2188	! ----------------------------------------------------------------------
2189	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
2190	! ----------------------------------------------------------------------
2191	DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2192	DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2193
2194	! ----------------------------------------------------------------------
2195	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
2196	! ----------------------------------------------------------------------
2197	DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2198	CI(K) = -DF1 * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2199	ELSE
2200
2201	! ----------------------------------------------------------------------
2202	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
2203	! ----------------------------------------------------------------------
2204	DTSDZ2 = (STC(K)-TBOT)/(.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT)
2205
2206	! ----------------------------------------------------------------------
2207	! SET MATRIX COEF, CI TO ZERO.
2208	! ----------------------------------------------------------------------
2209	CI(K) = 0.
2210	ENDIF
2211
2212	! ----------------------------------------------------------------------
2213	! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2214	! ----------------------------------------------------------------------
2215	DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2216	RHSTS(K) = ( DF1 * DTSDZ2 - DF1 * DTSDZ ) / DENOM
2217
2218	! ----------------------------------------------------------------------
2219	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2220	! ----------------------------------------------------------------------
2221	AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2222	BI(K) = -(AI(K) + CI(K))
2223
2224	! ----------------------------------------------------------------------
2225	! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
2226	! ----------------------------------------------------------------------
2227	DTSDZ = DTSDZ2
2228	DDZ = DDZ2
2229
2230	END DO
2231	! ----------------------------------------------------------------------
2232	! END SUBROUTINE HRTICE
2233	! ----------------------------------------------------------------------
2234	END SUBROUTINE HRTICE
2235
2236	SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
2237
2238	IMPLICIT NONE
2239
2240	! ----------------------------------------------------------------------
2241	! SUBROUTINE HSTEP
2242	! ----------------------------------------------------------------------
2243	! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
2244	! ----------------------------------------------------------------------
2245	INTEGER NSOLD
2246	PARAMETER(NSOLD = 20)
2247
2248	INTEGER K
2249	INTEGER NSOIL
2250
2251	REAL AI(NSOLD)
2252	REAL BI(NSOLD)
2253	REAL CI(NSOLD)
2254	REAL CIin(NSOLD)
2255	REAL DT
2256	REAL RHSTS(NSOIL)
2257	REAL RHSTSin(NSOIL)
2258	REAL STCIN(NSOIL)
2259	REAL STCOUT(NSOIL)
2260
2261	! ----------------------------------------------------------------------
2262	! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
2263	! ----------------------------------------------------------------------
2264	DO K = 1,NSOIL
2265	RHSTS(K) = RHSTS(K) * DT
2266	AI(K) = AI(K) * DT
2267	BI(K) = 1. + BI(K) * DT
2268	CI(K) = CI(K) * DT
2269	END DO
2270
2271	! ----------------------------------------------------------------------
2272	! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
2273	! ----------------------------------------------------------------------
2274	DO K = 1,NSOIL
2275	RHSTSin(K) = RHSTS(K)
2276	END DO
2277	DO K = 1,NSOIL
2278	CIin(K) = CI(K)
2279	END DO
2280
2281	! ----------------------------------------------------------------------
2282	! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
2283	! ----------------------------------------------------------------------
2284	CALL ROSR12(CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
2285
2286	! ----------------------------------------------------------------------
2287	! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
2288	! ----------------------------------------------------------------------
2289	DO K = 1,NSOIL
2290	STCOUT(K) = STCIN(K) + CI(K)
2291	END DO
2292
2293	! ----------------------------------------------------------------------
2294	! END SUBROUTINE HSTEP
2295	! ----------------------------------------------------------------------
2296	END SUBROUTINE HSTEP
2297
2298	SUBROUTINE NOPAC(ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
2299	& SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
2300	& SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL, &
2301	& STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
2302	& SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
2303	& DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
2304	& RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
2305	& QUARTZ,FXEXP,CSOIL, &
2306	& BETA,DRIP,DEW,FLX1,FLX2,FLX3)
2307
2308	IMPLICIT NONE
2309
2310	! ----------------------------------------------------------------------
2311	! SUBROUTINE NOPAC
2312	! ----------------------------------------------------------------------
2313	! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
2314	! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
2315	! PRESENT.
2316	! ----------------------------------------------------------------------
2317	INTEGER ICE
2318	INTEGER NROOT
2319	INTEGER NSOIL
2320
2321	REAL BEXP
2322	REAL BETA
2323	REAL CFACTR
2324	REAL CMC
2325	REAL CMCMAX
2326	REAL CP
2327	REAL CSOIL
2328	REAL DEW
2329	REAL DF1
2330	REAL DKSAT
2331	REAL DRIP
2332	REAL DT
2333	REAL DWSAT
2334	REAL EC
2335	REAL EDIR
2336	REAL EPSCA
2337	REAL ETA
2338	REAL ETA1
2339	REAL ETP
2340	REAL ETP1
2341	REAL ET(NSOIL)
2342	REAL ETT
2343	REAL FDOWN
2344	REAL F1
2345	REAL FXEXP
2346	REAL FLX1
2347	REAL FLX2
2348	REAL FLX3
2349	REAL FRZFACT
2350	REAL KDT
2351	REAL PC
2352	REAL PRCP
2353	REAL PRCP1
2354	REAL PSISAT
2355	REAL Q2
2356	REAL QUARTZ
2357	REAL RCH
2358	REAL RR
2359	REAL RTDIS(NSOIL)
2360	REAL RUNOFF1
2361	REAL RUNOFF2
2362	REAL RUNOFF3
2363	REAL SSOIL
2364	REAL SBETA
2365	REAL SFCTMP
2366	REAL SHDFAC
2367	REAL SH2O(NSOIL)
2368	REAL SIGMA
2369	REAL SLOPE
2370	REAL SMC(NSOIL)
2371	REAL SMCDRY
2372	REAL SMCMAX
2373	REAL SMCREF
2374	REAL SMCWLT
2375	REAL STC(NSOIL)
2376	REAL T1
2377	REAL T24
2378	REAL TBOT
2379	REAL TH2
2380	REAL YY
2381	REAL YYNUM
2382	REAL ZBOT
2383	REAL ZSOIL(NSOIL)
2384	REAL ZZ1
2385
2386	REAL EC1
2387	REAL EDIR1
2388	REAL ET1(NSOIL)
2389	REAL ETT1
2390
2391	INTEGER K
2392
2393	PARAMETER(CP = 1004.5)
2394	PARAMETER(SIGMA = 5.67E-8)
2395
2396	! ----------------------------------------------------------------------
2397	! EXECUTABLE CODE BEGINS HERE:
2398	! CONVERT ETP FROM KG M-2 S-1 TO MS-1 AND INITIALIZE DEW.
2399	! ----------------------------------------------------------------------
2400	PRCP1 = PRCP * 0.001
2401	ETP1 = ETP * 0.001
2402	DEW = 0.0
2403
2404	EDIR = 0.
2405	EDIR1 = 0.
2406	EC = 0.
2407	EC1 = 0.
2408	DO K = 1,NSOIL
2409	ET(K) = 0.
2410	ET1(K) = 0.
2411	END DO
2412	ETT = 0.
2413	ETT1 = 0.
2414
2415	IF (ETP .GT. 0.0) THEN
2416
2417	! ----------------------------------------------------------------------
2418	! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1'.
2419	! ----------------------------------------------------------------------
2420	CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2421	& SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2422	& SMCREF,SHDFAC,CMCMAX, &
2423	& SMCDRY,CFACTR, &
2424	& EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2425	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2426	& SH2O,SLOPE,KDT,FRZFACT, &
2427	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2428	& SHDFAC,CMCMAX, &
2429	& RUNOFF1,RUNOFF2,RUNOFF3, &
2430	& EDIR1,EC1,ET1, &
2431	& DRIP)
2432
2433	! ----------------------------------------------------------------------
2434	! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2435	! ----------------------------------------------------------------------
2436	! ETA = ETA1 * 1000.0
2437
2438	! ----------------------------------------------------------------------
2439	! EDIR = EDIR1 * 1000.0
2440	! EC = EC1 * 1000.0
2441	! ETT = ETT1 * 1000.0
2442	! ET(1) = ET1(1) * 1000.0
2443	! ET(2) = ET1(2) * 1000.0
2444	! ET(3) = ET1(3) * 1000.0
2445	! ET(4) = ET1(4) * 1000.0
2446	! ----------------------------------------------------------------------
2447
2448	ELSE
2449
2450	! ----------------------------------------------------------------------
2451	! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
2452	! ETP1 TO ZERO).
2453	! ----------------------------------------------------------------------
2454	DEW = -ETP1
2455	! ETP1 = 0.0
2456
2457	! ----------------------------------------------------------------------
2458	! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
2459	! ----------------------------------------------------------------------
2460	PRCP1 = PRCP1 + DEW
2461	!
2462	! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
2463	! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
2464	! & SMCREF,SHDFAC,CMCMAX,
2465	! & SMCDRY,CFACTR,
2466	! & EDIR1,EC1,ET1,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2467	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2468	& SH2O,SLOPE,KDT,FRZFACT, &
2469	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2470	& SHDFAC,CMCMAX, &
2471	& RUNOFF1,RUNOFF2,RUNOFF3, &
2472	& EDIR1,EC1,ET1, &
2473	& DRIP)
2474
2475	! ----------------------------------------------------------------------
2476	! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
2477	! ----------------------------------------------------------------------
2478	! ETA = ETA1 * 1000.0
2479
2480	! ----------------------------------------------------------------------
2481	! EDIR = EDIR1 * 1000.0
2482	! EC = EC1 * 1000.0
2483	! ETT = ETT1 * 1000.0
2484	! ET(1) = ET1(1) * 1000.0
2485	! ET(2) = ET1(2) * 1000.0
2486	! ET(3) = ET1(3) * 1000.0
2487	! ET(4) = ET1(4) * 1000.0
2488	! ----------------------------------------------------------------------
2489
2490	ENDIF
2491
2492	! ----------------------------------------------------------------------
2493	! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2494	! ----------------------------------------------------------------------
2495	ETA = ETA1 * 1000.0
2496
2497	! ----------------------------------------------------------------------
2498	EDIR = EDIR1 * 1000.0
2499	EC = EC1 * 1000.0
2500	DO K = 1,NSOIL
2501	ET(K) = ET1(K) * 1000.0
2502	! ET(1) = ET1(1) * 1000.0
2503	! ET(2) = ET1(2) * 1000.0
2504	! ET(3) = ET1(3) * 1000.0
2505	! ET(4) = ET1(4) * 1000.0
2506	ENDDO
2507	ETT = ETT1 * 1000.0
2508	! ----------------------------------------------------------------------
2509
2510	! ----------------------------------------------------------------------
2511	! BASED ON ETP AND E VALUES, DETERMINE BETA
2512	! ----------------------------------------------------------------------
2513	IF ( ETP .LE. 0.0 ) THEN
2514	BETA = 0.0
2515	IF ( ETP .LT. 0.0 ) THEN
2516	BETA = 1.0
2517	! ETA = ETP
2518	ENDIF
2519	ELSE
2520	BETA = ETA / ETP
2521	ENDIF
2522
2523	! ----------------------------------------------------------------------
2524	! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
2525	! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
2526	! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
2527	! ----------------------------------------------------------------------
2528	CALL TDFCND (DF1,SMC(1),QUARTZ,SMCMAX,SH2O(1))
2529
2530	! ----------------------------------------------------------------------
2531	! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
2532	! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
2533	! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
2534	! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
2535	! ROUTINE SFLX)
2536	! ----------------------------------------------------------------------
2537	DF1 = DF1 * EXP(SBETA*SHDFAC)
2538
2539	! ----------------------------------------------------------------------
2540	! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
2541	! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
2542	! ----------------------------------------------------------------------
2543	YYNUM = FDOWN - SIGMA * T24
2544	YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-BETA*EPSCA) / RR
2545	ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0
2546
2547	CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2548	& TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2549	& QUARTZ,CSOIL)
2550
2551	! ----------------------------------------------------------------------
2552	! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
2553	! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
2554	! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
2555	! ----------------------------------------------------------------------
2556	FLX1 = 0.0
2557	FLX3 = 0.0
2558
2559	! ----------------------------------------------------------------------
2560	! END SUBROUTINE NOPAC
2561	! ----------------------------------------------------------------------
2562	END SUBROUTINE NOPAC
2563
2564	SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
2565	& Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
2566	& DQSDT2,FLX2)
2567
2568	IMPLICIT NONE
2569
2570	! ----------------------------------------------------------------------
2571	! SUBROUTINE PENMAN
2572	! ----------------------------------------------------------------------
2573	! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
2574	! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
2575	! CALLING ROUTINE FOR LATER USE.
2576	! ----------------------------------------------------------------------
2577	LOGICAL SNOWNG
2578	LOGICAL FRZGRA
2579
2580	REAL A
2581	REAL BETA
2582	REAL CH
2583	REAL CP
2584	REAL CPH2O
2585	REAL CPICE
2586	REAL DELTA
2587	REAL DQSDT2
2588	REAL ELCP
2589	REAL EPSCA
2590	REAL ETP
2591	REAL FDOWN
2592	REAL FLX2
2593	REAL FNET
2594	REAL LSUBC
2595	REAL LSUBF
2596	REAL PRCP
2597	REAL Q2
2598	REAL Q2SAT
2599	REAL R
2600	REAL RAD
2601	REAL RCH
2602	REAL RHO
2603	REAL RR
2604	REAL SSOIL
2605	REAL SFCPRS
2606	REAL SFCTMP
2607	REAL SIGMA
2608	REAL T24
2609	REAL T2V
2610	REAL TH2
2611
2612	PARAMETER(CP = 1004.6)
2613	PARAMETER(CPH2O = 4.218E+3)
2614	PARAMETER(CPICE = 2.106E+3)
2615	PARAMETER(R = 287.04)
2616	PARAMETER(ELCP = 2.4888E+3)
2617	PARAMETER(LSUBF = 3.335E+5)
2618	PARAMETER(LSUBC = 2.501000E+6)
2619	PARAMETER(SIGMA = 5.67E-8)
2620
2621	! ----------------------------------------------------------------------
2622	! EXECUTABLE CODE BEGINS HERE:
2623	! ----------------------------------------------------------------------
2624	FLX2 = 0.0
2625
2626	! ----------------------------------------------------------------------
2627	! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
2628	! ----------------------------------------------------------------------
2629	DELTA = ELCP * DQSDT2
2630	T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
2631	RR = T24 * 6.48E-8 /(SFCPRS * CH) + 1.0
2632	RHO = SFCPRS / (R * T2V)
2633	RCH = RHO * CP * CH
2634
2635	! ----------------------------------------------------------------------
2636	! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
2637	! EFFECTS CAUSED BY FALLING PRECIPITATION.
2638	! ----------------------------------------------------------------------
2639	IF (.NOT. SNOWNG) THEN
2640	IF (PRCP .GT. 0.0) RR = RR + CPH2O*PRCP/RCH
2641	ELSE
2642	RR = RR + CPICE*PRCP/RCH
2643	ENDIF
2644
2645	FNET = FDOWN - SIGMA*T24 - SSOIL
2646
2647	! ----------------------------------------------------------------------
2648	! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
2649	! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
2650	! ----------------------------------------------------------------------
2651	IF (FRZGRA) THEN
2652	FLX2 = -LSUBF * PRCP
2653	FNET = FNET - FLX2
2654	ENDIF
2655
2656	! ----------------------------------------------------------------------
2657	! FINISH PENMAN EQUATION CALCULATIONS.
2658	! ----------------------------------------------------------------------
2659	RAD = FNET/RCH + TH2 - SFCTMP
2660	A = ELCP * (Q2SAT - Q2)
2661	EPSCA = (ARR + RADDELTA) / (DELTA + RR)
2662	ETP = EPSCA * RCH / LSUBC
2663
2664	! ----------------------------------------------------------------------
2665	! END SUBROUTINE PENMAN
2666	! ----------------------------------------------------------------------
2667	END SUBROUTINE PENMAN
2668
2669	SUBROUTINE REDPRM ( &
2670	& VEGTYP,SOILTYP,SLOPETYP, &
2671	& CFACTR,CMCMAX,RSMAX,TOPT,REFKDT,KDT,SBETA, &
2672	& SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,PSISAT,SLOPE, &
2673	& SNUP,SALP,BEXP,DKSAT,DWSAT, &
2674	& SMCMAX,SMCWLT,SMCREF, &
2675	& SMCDRY,F1,QUARTZ,FXEXP,RTDIS,SLDPTH,ZSOIL, &
2676	& NROOT,NSOIL,Z0,CZIL,LAI,CSOIL,PTU)
2677
2678	IMPLICIT NONE
2679	! ----------------------------------------------------------------------
2680	! SUBROUTINE REDPRM
2681	! ----------------------------------------------------------------------
2682	! INTERNALLY SET (DEFAULT VALUESS), OR OPTIONALLY READ-IN VIA NAMELIST
2683	! I/O, ALL SOIL AND VEGETATION PARAMETERS REQUIRED FOR THE EXECUSION OF
2684	! THE NOAH LSM.
2685	!
2686	! OPTIONAL NON-DEFAULT PARAMETERS CAN BE READ IN, ACCOMMODATING UP TO 30
2687	! SOIL, VEG, OR SLOPE CLASSES, IF THE DEFAULT MAX NUMBER OF SOIL, VEG,
2688	! AND/OR SLOPE TYPES IS RESET.
2689	!
2690	! FUTURE UPGRADES OF ROUTINE REDPRM MUST EXPAND TO INCORPORATE SOME OF
2691	! THE EMPIRICAL PARAMETERS OF THE FROZEN SOIL AND SNOWPACK PHYSICS (SUCH
2692	! AS IN ROUTINES FRH2O, SNOWPACK, AND SNOW_NEW) NOT YET SET IN THIS
2693	! REDPRM ROUTINE, BUT RATHER SET IN LOWER LEVEL SUBROUTINES.
2694	!
2695	! SET MAXIMUM NUMBER OF SOIL-, VEG-, AND SLOPETYP IN DATA STATEMENT.
2696	! ----------------------------------------------------------------------
2697	INTEGER MAX_SLOPETYP
2698	INTEGER MAX_SOILTYP
2699	INTEGER MAX_VEGTYP
2700
2701	PARAMETER(MAX_SLOPETYP = 30)
2702	PARAMETER(MAX_SOILTYP = 30)
2703	PARAMETER(MAX_VEGTYP = 30)
2704
2705	! ----------------------------------------------------------------------
2706	! NUMBER OF DEFINED SOIL-, VEG-, AND SLOPETYPS USED.
2707	! ----------------------------------------------------------------------
2708	INTEGER DEFINED_VEG
2709	INTEGER DEFINED_SOIL
2710	INTEGER DEFINED_SLOPE
2711
2712	DATA DEFINED_VEG/27/
2713	DATA DEFINED_SOIL/19/
2714	DATA DEFINED_SLOPE/9/
2715
2716	! ----------------------------------------------------------------------
2717	! SET-UP SOIL PARAMETERS FOR GIVEN SOIL TYPE
2718	! INPUT: SOLTYP: SOIL TYPE (INTEGER INDEX)
2719	! OUTPUT: SOIL PARAMETERS:
2720	! MAXSMC: MAX SOIL MOISTURE CONTENT (POROSITY)
2721	! REFSMC: REFERENCE SOIL MOISTURE (ONSET OF SOIL MOISTURE
2722	! STRESS IN TRANSPIRATION)
2723	! WLTSMC: WILTING PT SOIL MOISTURE CONTENTS
2724	! DRYSMC: AIR DRY SOIL MOIST CONTENT LIMITS
2725	! SATPSI: SATURATED SOIL POTENTIAL
2726	! SATDK: SATURATED SOIL HYDRAULIC CONDUCTIVITY
2727	! BB: THE 'B' PARAMETER
2728	! SATDW: SATURATED SOIL DIFFUSIVITY
2729	! F11: USED TO COMPUTE SOIL DIFFUSIVITY/CONDUCTIVITY
2730	! QUARTZ: SOIL QUARTZ CONTENT
2731	! ----------------------------------------------------------------------
2732	! SOIL STATSGO
2733	! TYPE CLASS
2734	! ---- -------
2735	! 1 SAND
2736	! 2 LOAMY SAND
2737	! 3 SANDY LOAM
2738	! 4 SILT LOAM
2739	! 5 SILT
2740	! 6 LOAM
2741	! 7 SANDY CLAY LOAM
2742	! 8 SILTY CLAY LOAM
2743	! 9 CLAY LOAM
2744	! 10 SANDY CLAY
2745	! 11 SILTY CLAY
2746	! 12 CLAY
2747	! 13 ORGANIC MATERIAL
2748	! 14 WATER
2749	! 15 BEDROCK
2750	! 16 OTHER(land-ice)
2751	! 17 PLAYA
2752	! 18 LAVA
2753	! 19 WHITE SAND
2754	! ----------------------------------------------------------------------
2755
2756	REAL BB(MAX_SOILTYP)
2757	REAL DRYSMC(MAX_SOILTYP)
2758	REAL F11(MAX_SOILTYP)
2759	REAL MAXSMC(MAX_SOILTYP)
2760	REAL REFSMC(MAX_SOILTYP)
2761	REAL SATPSI(MAX_SOILTYP)
2762	REAL SATDK(MAX_SOILTYP)
2763	REAL SATDW(MAX_SOILTYP)
2764	REAL WLTSMC(MAX_SOILTYP)
2765	REAL QTZ(MAX_SOILTYP)
2766
2767	REAL BEXP
2768	REAL DKSAT
2769	REAL DWSAT
2770	REAL F1
2771	REAL PTU
2772	REAL QUARTZ
2773	REAL REFSMC1
2774	REAL SMCDRY
2775	REAL SMCMAX
2776	REAL SMCREF
2777	REAL SMCWLT
2778	REAL WLTSMC1
2779
2780	! ----------------------------------------------------------------------
2781	! SOIL TEXTURE-RELATED ARRAYS.
2782	! ----------------------------------------------------------------------
2783	DATA MAXSMC/0.395, 0.421, 0.434, 0.476, 0.476, 0.439, &
2784	& 0.404, 0.464, 0.465, 0.406, 0.468, 0.457, &
2785	& 0.464, 0.464, 0.200, 0.421, 0.457, 0.200, &
2786	& 0.395, 0.000, 0.000, 0.000, 0.000, 0.000, &
2787	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2788	! ----------------------------------------------------------------------
2789	DATA SATPSI/0.0350, 0.0363, 0.1413, 0.7586, 0.7586, 0.3548, &
2790	& 0.1349, 0.6166, 0.2630, 0.0977, 0.3236, 0.4677, &
2791	& 0.3548, 0.3548, 0.0350, 0.0363, 0.4677, 0.0350, &
2792	& 0.0350, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, &
2793	& 0.000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000/
2794	! ----------------------------------------------------------------------
2795	DATA SATDK /1.7600E-4, 1.4078E-5, 5.2304E-6, 2.8089E-6, 2.8089E-6,&
2796	& 3.3770E-6, 4.4518E-6, 2.0348E-6, 2.4464E-6, 7.2199E-6,&
2797	& 1.3444E-6, 9.7394E-7, 3.3770E-6, 3.3770E-6, 1.4078E-5,&
2798	& 1.4078E-5, 9.7394E-7, 1.4078E-5, 1.7600E-4, 0.0,&
2799	& 0.0, 0.0, 0.0, 0.0, 0.0,&
2800	& 0.0, 0.0, 0.0, 0.0, 0.0/
2801	! ----------------------------------------------------------------------
2802	DATA BB /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
2803	& 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
2804	& 5.25, 5.25, 4.05, 4.26, 11.55, 4.05, &
2805	& 4.05, 0.00, 0.00, 0.00, 0.00, 0.00, &
2806	& 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2807	! ----------------------------------------------------------------------
2808	DATA QTZ /0.92, 0.82, 0.60, 0.25, 0.10, 0.40, &
2809	& 0.60, 0.10, 0.35, 0.52, 0.10, 0.25, &
2810	& 0.05, 0.05, 0.07, 0.25, 0.60, 0.52, &
2811	& 0.92, 0.00, 0.00, 0.00, 0.00, 0.00, &
2812	& 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2813
2814	! ----------------------------------------------------------------------
2815	! THE FOLLOWING 5 PARAMETERS ARE DERIVED LATER IN REDPRM.F FROM THE SOIL
2816	! DATA, AND ARE JUST GIVEN HERE FOR REFERENCE AND TO FORCE STATIC
2817	! STORAGE ALLOCATION. -DAG LOHMANN, FEB. 2001
2818	! ----------------------------------------------------------------------
2819	! !!!!!!!!!!!!!! The following values in the table are NOT used
2820	! !!!!!!!!!!!!!! and are just given for reference
2821	DATA REFSMC/0.196, 0.248, 0.282, 0.332, 0.332, 0.301, &
2822	& 0.293, 0.368, 0.361, 0.320, 0.388, 0.389, &
2823	& 0.319, 0.000, 0.116, 0.248, 0.389, 0.116, &
2824	& 0.196, 0.000, 0.000, 0.000, 0.000, 0.000, &
2825	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2826	! !!!!!!!!!!!!!! The following values in the table are NOT used
2827	! !!!!!!!!!!!!!! and are just given for reference
2828	DATA WLTSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2829	& 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2830	& 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2831	& 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2832	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2833	! !!!!!!!!!!!!!! The following values in the table are NOT used
2834	! !!!!!!!!!!!!!! and are just given for reference
2835	DATA DRYSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2836	& 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2837	& 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2838	& 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2839	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2840
2841	! !!!!!!!!!!!!!! The following values in the table are NOT used
2842	! !!!!!!!!!!!!!! and are just given for reference
2843	DATA SATDW /0.632E-4, 0.517E-5, 0.807E-5, 0.239E-4, 0.239E-4, &
2844	& 0.143E-4, 0.101E-4, 0.236E-4, 0.113E-4, 0.186E-4, &
2845	& 0.966E-5, 0.115E-4, 0.136E-4, 0.0, 0.998E-5, &
2846	& 0.517E-5, 0.115E-4, 0.998E-5, 0.632E-4, 0.0, &
2847	& 0.0, 0.0, 0.0, 0.0, 0.0, &
2848	& 0.0, 0.0, 0.0, 0.0, 0.0/
2849	! !!!!!!!!!!!!!! The following values in the table are NOT used
2850	! !!!!!!!!!!!!!! and are just given for reference
2851	DATA F11 /-1.090, -1.041, -0.568, 0.162, 0.162, -0.327, &
2852	& -1.535, -1.118, -1.297, -3.211, -1.916, -2.258, &
2853	& -0.201, 0.000, -2.287, -1.041, -2.258, -2.287, &
2854	& -1.090, 0.000, 0.000, 0.000, 0.000, 0.000, &
2855	& 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2856
2857	! ----------------------------------------------------------------------
2858	! SET-UP VEGETATION PARAMETERS FOR A GIVEN VEGETAION TYPE:
2859	! INPUT: VEGTYP = VEGETATION TYPE (INTEGER INDEX)
2860	! OUPUT: VEGETATION PARAMETERS
2861	! SHDFAC: VEGETATION GREENNESS FRACTION
2862	! RSMIN: MIMIMUM STOMATAL RESISTANCE
2863	! RGL: PARAMETER USED IN SOLAR RAD TERM OF
2864	! CANOPY RESISTANCE FUNCTION
2865	! HS: PARAMETER USED IN VAPOR PRESSURE DEFICIT TERM OF
2866	! CANOPY RESISTANCE FUNCTION
2867	! SNUP: THRESHOLD SNOW DEPTH (IN WATER EQUIVALENT M) THAT
2868	! IMPLIES 100% SNOW COVER
2869	! ----------------------------------------------------------------------
2870	! CLASS USGS-WRF VEGETATION/SURFACE TYPE
2871	! 1 Urban and Built-Up Land
2872	! 2 Dryland Cropland and Pasture
2873	! 3 Irrigated Cropland and Pasture
2874	! 4 Mixed Dryland/Irrigated Cropland and Pasture
2875	! 5 Cropland/Grassland Mosaic
2876	! 6 Cropland/Woodland Mosaic
2877	! 7 Grassland
2878	! 8 Shrubland
2879	! 9 Mixed Shrubland/Grassland
2880	! 10 Savanna
2881	! 11 Deciduous Broadleaf Forest
2882	! 12 Deciduous Needleleaf Forest
2883	! 13 Evergreen Broadleaf Forest
2884	! 14 Evergreen Needleleaf Forest
2885	! 15 Mixed Forest
2886	! 16 Water Bodies
2887	! 17 Herbaceous Wetland
2888	! 18 Wooded Wetland
2889	! 19 Barren or Sparsely Vegetated
2890	! 20 Herbaceous Tundra
2891	! 21 Wooded Tundra
2892	! 22 Mixed Tundra
2893	! 23 Bare Ground Tundra
2894	! 24 Snow or Ice
2895	! 25 Playa
2896	! 26 Lava
2897	! 27 White Sand
2898	! ----------------------------------------------------------------------
2899
2900	INTEGER NROOT
2901	INTEGER NROOT_DATA(MAX_VEGTYP)
2902
2903	REAL FRZFACT
2904	REAL HS
2905	REAL HSTBL(MAX_VEGTYP)
2906	REAL LAI
2907	REAL LAI_DATA(MAX_VEGTYP)
2908	REAL PSISAT
2909	REAL RSMIN
2910	REAL RGL
2911	REAL RGLTBL(MAX_VEGTYP)
2912	REAL RSMTBL(MAX_VEGTYP)
2913	REAL SHDFAC
2914	REAL SNUP
2915	REAL SNUPX(MAX_VEGTYP)
2916	REAL Z0
2917	REAL Z0_DATA(MAX_VEGTYP)
2918
2919	! ----------------------------------------------------------------------
2920	! VEGETATION CLASS-RELATED ARRAYS
2921	! ----------------------------------------------------------------------
2922	! DATA NROOT_DATA /2,3,3,3,3,3,3,3,3,3,
2923	! & 4,4,4,4,4,2,2,2,2,3,
2924	! & 3,3,2,2,2,2,2,0,0,0/
2925	DATA NROOT_DATA /1,3,3,3,3,3,3,3,3,3, &
2926	& 4,4,4,4,4,0,2,2,1,3, &
2927	& 3,3,2,1,1,1,1,0,0,0/
2928	DATA RSMTBL /200.0, 70.0, 70.0, 70.0, 70.0, 70.0, &
2929	& 70.0, 300.0, 170.0, 70.0, 100.0, 150.0, &
2930	& 150.0, 125.0, 125.0, 100.0, 40.0, 100.0, &
2931	& 300.0, 150.0, 150.0, 150.0, 200.0, 200.0, &
2932	& 40.0, 100.0, 300.0, 0.0, 0.0, 0.0/
2933	DATA RGLTBL /100.0, 100.0, 100.0, 100.0, 100.0, 65.0, &
2934	& 100.0, 100.0, 100.0, 65.0, 30.0, 30.0, &
2935	& 30.0, 30.0, 30.0, 30.0, 100.0, 30.0, &
2936	& 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, &
2937	& 100.0, 100.0, 100.0, 0.0, 0.0, 0.0/
2938	DATA HSTBL /42.00, 36.25, 36.25, 36.25, 36.25, 44.14, &
2939	& 36.35, 42.00, 39.18, 54.53, 54.53, 47.35, &
2940	& 41.69, 47.35, 51.93, 51.75, 60.00, 51.93, &
2941	& 42.00, 42.00, 42.00, 42.00, 42.00, 42.00, &
2942	& 36.25, 42.00, 42.00, 0.00, 0.00, 0.00/
2943	DATA SNUPX /0.020, 0.020, 0.020, 0.020, 0.020, 0.020, &
2944	& 0.020, 0.020, 0.020, 0.040, 0.040, 0.040, &
2945	& 0.040, 0.040, 0.040, 0.010, 0.013, 0.020, &
2946	& 0.013, 0.020, 0.020, 0.020, 0.020, 0.013, &
2947	& 0.013, 0.013, 0.013, 0.000, 0.000, 0.000/
2948	DATA Z0_DATA / 1.00, 0.07, 0.07, 0.07, 0.07, 0.15, &
2949	& 0.08, 0.03, 0.05, 0.86, 0.80, 0.85, &
2950	& 2.65, 1.09, 0.80, 0.001, 0.04, 0.05, &
2951	& 0.01, 0.04, 0.06, 0.05, 0.03, 0.001, &
2952	& 0.01, 0.15, 0.01, 0.00, 0.00, 0.00/
2953	DATA LAI_DATA /4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2954	& 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2955	& 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2956	& 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2957	& 4.0, 4.0, 4.0, 0.0, 0.0, 0.0/
2958
2959	! ----------------------------------------------------------------------
2960	! CLASS PARAMETER 'SLOPETYP' WAS INCLUDED TO ESTIMATE LINEAR RESERVOIR
2961	! COEFFICIENT 'SLOPE' TO THE BASEFLOW RUNOFF OUT OF THE BOTTOM LAYER.
2962	! LOWEST CLASS (SLOPETYP=0) MEANS HIGHEST SLOPE PARAMETER = 1.
2963	! DEFINITION OF SLOPETYP FROM 'ZOBLER' SLOPE TYPE:
2964	! SLOPE CLASS PERCENT SLOPE
2965	! 1 0-8
2966	! 2 8-30
2967	! 3 > 30
2968	! 4 0-30
2969	! 5 0-8 & > 30
2970	! 6 8-30 & > 30
2971	! 7 0-8, 8-30, > 30
2972	! 9 GLACIAL ICE
2973	! BLANK OCEAN/SEA
2974	! ----------------------------------------------------------------------
2975	! NOTE:
2976	! CLASS 9 FROM 'ZOBLER' FILE SHOULD BE REPLACED BY 8 AND 'BLANK' 9
2977	! ----------------------------------------------------------------------
2978	REAL SLOPE
2979	REAL SLOPE_DATA(MAX_SLOPETYP)
2980
2981	DATA SLOPE_DATA /0.1, 0.6, 1.0, 0.35, 0.55, 0.8, &
2982	& 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, &
2983	& 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2984	& 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2985	& 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/
2986
2987	! ----------------------------------------------------------------------
2988	! SET NAMELIST FILE NAME
2989	! ----------------------------------------------------------------------
2990	CHARACTER*50 NAMELIST_NAME
2991
2992	! ----------------------------------------------------------------------
2993	! SET UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOIL, VEG, SLOPE TYPE)
2994	! ----------------------------------------------------------------------
2995	INTEGER I
2996	INTEGER NSOIL
2997	INTEGER SLOPETYP
2998	INTEGER SOILTYP
2999	INTEGER VEGTYP
3000
3001	INTEGER BARE
3002	! DATA BARE /11/
3003	DATA BARE /19/
3004
3005	LOGICAL LPARAM
3006	DATA LPARAM /.TRUE./
3007
3008	LOGICAL LFIRST
3009	DATA LFIRST /.TRUE./
3010
3011	! ----------------------------------------------------------------------
3012	! PARAMETER USED TO CALCULATE ROUGHNESS LENGTH OF HEAT.
3013	! ----------------------------------------------------------------------
3014	REAL CZIL
3015	REAL CZIL_DATA
3016	! changed in version 2.6 June 2nd 2003
3017	! DATA CZIL_DATA /0.2/
3018	DATA CZIL_DATA /0.1/
3019
3020	! ----------------------------------------------------------------------
3021	! PARAMETER USED TO CALUCULATE VEGETATION EFFECT ON SOIL HEAT FLUX.
3022	! ----------------------------------------------------------------------
3023	REAL SBETA
3024	REAL SBETA_DATA
3025	DATA SBETA_DATA /-2.0/
3026
3027	! ----------------------------------------------------------------------
3028	! BARE SOIL EVAPORATION EXPONENT USED IN DEVAP.
3029	! ----------------------------------------------------------------------
3030	REAL FXEXP
3031	REAL FXEXP_DATA
3032	DATA FXEXP_DATA /2.0/
3033
3034	! ----------------------------------------------------------------------
3035	! SOIL HEAT CAPACITY [J M-3 K-1]
3036	! ----------------------------------------------------------------------
3037	REAL CSOIL
3038	REAL CSOIL_DATA
3039	! DATA CSOIL_DATA /1.26E+6/
3040	DATA CSOIL_DATA /2.00E+6/
3041
3042	! ----------------------------------------------------------------------
3043	! SPECIFY SNOW DISTRIBUTION SHAPE PARAMETER SALP - SHAPE PARAMETER OF
3044	! DISTRIBUTION FUNCTION OF SNOW COVER. FROM ANDERSON'S DATA (HYDRO-17)
3045	! BEST FIT IS WHEN SALP = 2.6
3046	! ----------------------------------------------------------------------
3047	REAL SALP
3048	REAL SALP_DATA
3049	! changed for version 2.6 June 2nd 2003
3050	! DATA SALP_DATA /2.6/
3051	DATA SALP_DATA /4.0/
3052
3053	! ----------------------------------------------------------------------
3054	! KDT IS DEFINED BY REFERENCE REFKDT AND DKSAT; REFDK=2.E-6 IS THE SAT.
3055	! DK. VALUE FOR THE SOIL TYPE 2
3056	! ----------------------------------------------------------------------
3057	REAL REFDK
3058	REAL REFDK_DATA
3059	DATA REFDK_DATA /2.0E-6/
3060
3061	REAL REFKDT
3062	REAL REFKDT_DATA
3063	DATA REFKDT_DATA /3.0/
3064
3065	REAL FRZX
3066	REAL KDT
3067
3068	! ----------------------------------------------------------------------
3069	! FROZEN GROUND PARAMETER, FRZK, DEFINITION: ICE CONTENT THRESHOLD ABOVE
3070	! WHICH FROZEN SOIL IS IMPERMEABLE REFERENCE VALUE OF THIS PARAMETER FOR
3071	! THE LIGHT CLAY SOIL (TYPE=3) FRZK = 0.15 M.
3072	! ----------------------------------------------------------------------
3073	REAL FRZK
3074	REAL FRZK_DATA
3075	DATA FRZK_DATA /0.15/
3076
3077	REAL RTDIS(NSOIL)
3078	REAL SLDPTH(NSOIL)
3079	REAL ZSOIL(NSOIL)
3080
3081	! ----------------------------------------------------------------------
3082	! SET TWO CANOPY WATER PARAMETERS.
3083	! ----------------------------------------------------------------------
3084	REAL CFACTR
3085	REAL CFACTR_DATA
3086	DATA CFACTR_DATA /0.5/
3087
3088	REAL CMCMAX
3089	REAL CMCMAX_DATA
3090	DATA CMCMAX_DATA /0.5E-3/
3091
3092	! ----------------------------------------------------------------------
3093	! SET MAX. STOMATAL RESISTANCE.
3094	! ----------------------------------------------------------------------
3095	REAL RSMAX
3096	REAL RSMAX_DATA
3097	DATA RSMAX_DATA /5000.0/
3098
3099	! ----------------------------------------------------------------------
3100	! SET OPTIMUM TRANSPIRATION AIR TEMPERATURE.
3101	! ----------------------------------------------------------------------
3102	REAL TOPT
3103	REAL TOPT_DATA
3104	DATA TOPT_DATA /298.0/
3105
3106	! ----------------------------------------------------------------------
3107	! SPECIFY DEPTH[M] OF LOWER BOUNDARY SOIL TEMPERATURE.
3108	! ----------------------------------------------------------------------
3109	REAL ZBOT
3110	REAL ZBOT_DATA
3111	! changed for version 2.5.2
3112	! DATA ZBOT_DATA /-3.0/
3113	DATA ZBOT_DATA /-8.0/
3114
3115	! ----------------------------------------------------------------------
3116	! SET TWO SOIL MOISTURE WILT, SOIL MOISTURE REFERENCE PARAMETERS
3117	! ----------------------------------------------------------------------
3118	REAL SMLOW
3119	REAL SMLOW_DATA
3120	DATA SMLOW_DATA /0.5/
3121
3122	REAL SMHIGH
3123	REAL SMHIGH_DATA
3124	! changed in 2.6 from 3 to 6 on June 2nd 2003
3125	DATA SMHIGH_DATA /3.0/
3126	! DATA SMHIGH_DATA /6.0/
3127
3128	! ----------------------------------------------------------------------
3129	! NAMELIST DEFINITION:
3130	! ----------------------------------------------------------------------
3131	NAMELIST /SOIL_VEG/ SLOPE_DATA, RSMTBL, RGLTBL, HSTBL, SNUPX, &
3132	& BB, DRYSMC, F11, MAXSMC, REFSMC, SATPSI, SATDK, SATDW, &
3133	& WLTSMC, QTZ, LPARAM, ZBOT_DATA, SALP_DATA, CFACTR_DATA, &
3134	& CMCMAX_DATA, SBETA_DATA, RSMAX_DATA, TOPT_DATA, &
3135	& REFDK_DATA, FRZK_DATA, BARE, DEFINED_VEG, DEFINED_SOIL, &
3136	& DEFINED_SLOPE, FXEXP_DATA, NROOT_DATA, REFKDT_DATA, Z0_DATA, &
3137	& CZIL_DATA, LAI_DATA, CSOIL_DATA
3138
3139	! ----------------------------------------------------------------------
3140	! READ NAMELIST FILE TO OVERRIDE DEFAULT PARAMETERS ONLY ONCE.
3141	! NAMELIST_NAME must be 50 characters or less.
3142	! ----------------------------------------------------------------------
3143	IF (LFIRST) THEN
3144	! WRITE(,) 'READ NAMELIST'
3145	! OPEN(58, FILE = 'namelist_filename.txt')
3146	! READ(58,'(A)') NAMELIST_NAME
3147	! CLOSE(58)
3148	! WRITE(,) 'Namelist Filename is ', NAMELIST_NAME
3149	! OPEN(59, FILE = NAMELIST_NAME)
3150	! 50 CONTINUE
3151	! READ(59, SOIL_VEG, END=100)
3152	! IF (LPARAM) GOTO 50
3153	! 100 CONTINUE
3154	! CLOSE(59)
3155	! WRITE(*,NML=SOIL_VEG)
3156	LFIRST = .FALSE.
3157	IF (DEFINED_SOIL .GT. MAX_SOILTYP) THEN
3158	WRITE(,) 'Warning: DEFINED_SOIL too large in namelist'
3159	STOP 222
3160	ENDIF
3161	IF (DEFINED_VEG .GT. MAX_VEGTYP) THEN
3162	WRITE(,) 'Warning: DEFINED_VEG too large in namelist'
3163	STOP 222
3164	ENDIF
3165	IF (DEFINED_SLOPE .GT. MAX_SLOPETYP) THEN
3166	WRITE(,) 'Warning: DEFINED_SLOPE too large in namelist'
3167	STOP 222
3168	ENDIF
3169
3170	SMLOW = SMLOW_DATA
3171	SMHIGH = SMHIGH_DATA
3172
3173	DO I = 1,DEFINED_SOIL
3174	SATDW(I) = BB(I)SATDK(I)(SATPSI(I)/MAXSMC(I))
3175	F11(I) = ALOG10(SATPSI(I)) + BB(I)*ALOG10(MAXSMC(I)) + 2.0
3176	REFSMC1 = MAXSMC(I)*(5.79E-9/SATDK(I)) &
3177	& *(1.0/(2.0BB(I)+3.0))
3178	REFSMC(I) = REFSMC1 + (MAXSMC(I)-REFSMC1) / SMHIGH
3179	WLTSMC1 = MAXSMC(I) * (200.0/SATPSI(I))**(-1.0/BB(I))
3180	WLTSMC(I) = WLTSMC1 - SMLOW * WLTSMC1
3181
3182	! ----------------------------------------------------------------------
3183	! CURRENT VERSION DRYSMC VALUES THAT EQUATE TO WLTSMC.
3184	! FUTURE VERSION COULD LET DRYSMC BE INDEPENDENTLY SET VIA NAMELIST.
3185	! ----------------------------------------------------------------------
3186	DRYSMC(I) = WLTSMC(I)
3187	END DO
3188
3189	! ----------------------------------------------------------------------
3190	! END LFIRST BLOCK
3191	! ----------------------------------------------------------------------
3192	ENDIF
3193
3194	IF (SOILTYP .GT. DEFINED_SOIL) THEN
3195	WRITE(,) 'Warning: too many soil types'
3196	STOP 333
3197	ENDIF
3198	IF (VEGTYP .GT. DEFINED_VEG) THEN
3199	WRITE(,) 'Warning: too many veg types'
3200	STOP 333
3201	ENDIF
3202	IF (SLOPETYP .GT. DEFINED_SLOPE) THEN
3203	WRITE(,) 'Warning: too many slope types'
3204	STOP 333
3205	ENDIF
3206
3207	! ----------------------------------------------------------------------
3208	! SET-UP UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOILTYP, VEGTYP OR
3209	! SLOPETYP)
3210	! ----------------------------------------------------------------------
3211	ZBOT = ZBOT_DATA
3212	SALP = SALP_DATA
3213	CFACTR = CFACTR_DATA
3214	CMCMAX = CMCMAX_DATA
3215	SBETA = SBETA_DATA
3216	RSMAX = RSMAX_DATA
3217	TOPT = TOPT_DATA
3218	REFDK = REFDK_DATA
3219	FRZK = FRZK_DATA
3220	FXEXP = FXEXP_DATA
3221	REFKDT = REFKDT_DATA
3222	CZIL = CZIL_DATA
3223	CSOIL = CSOIL_DATA
3224
3225	! ----------------------------------------------------------------------
3226	! SET-UP SOIL PARAMETERS
3227	! ----------------------------------------------------------------------
3228	BEXP = BB(SOILTYP)
3229	DKSAT = SATDK(SOILTYP)
3230	DWSAT = SATDW(SOILTYP)
3231	F1 = F11(SOILTYP)
3232	KDT = REFKDT * DKSAT/REFDK
3233	PSISAT = SATPSI(SOILTYP)
3234	QUARTZ = QTZ(SOILTYP)
3235	SMCDRY = DRYSMC(SOILTYP)
3236	SMCMAX = MAXSMC(SOILTYP)
3237	SMCREF = REFSMC(SOILTYP)
3238	SMCWLT = WLTSMC(SOILTYP)
3239	FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
3240
3241	! ----------------------------------------------------------------------
3242	! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
3243	! ----------------------------------------------------------------------
3244	FRZX = FRZK * FRZFACT
3245
3246	! ----------------------------------------------------------------------
3247	! SET-UP VEGETATION PARAMETERS
3248	! ----------------------------------------------------------------------
3249	NROOT = NROOT_DATA(VEGTYP)
3250	SNUP = SNUPX(VEGTYP)
3251	RSMIN = RSMTBL(VEGTYP)
3252	RGL = RGLTBL(VEGTYP)
3253	HS = HSTBL(VEGTYP)
3254	Z0 = Z0_DATA(VEGTYP)
3255	LAI = LAI_DATA(VEGTYP)
3256	IF (VEGTYP .EQ. BARE) SHDFAC = 0.0
3257
3258	IF (NROOT .GT. NSOIL) THEN
3259	WRITE(,) 'Warning: too many root layers'
3260	STOP 333
3261	ENDIF
3262
3263	! ----------------------------------------------------------------------
3264	! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
3265	! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
3266	! ----------------------------------------------------------------------
3267	DO I = 1,NROOT
3268	RTDIS(I) = -SLDPTH(I)/ZSOIL(NROOT)
3269	END DO
3270
3271	! ----------------------------------------------------------------------
3272	! SET-UP SLOPE PARAMETER
3273	! ----------------------------------------------------------------------
3274	SLOPE = SLOPE_DATA(SLOPETYP)
3275
3276	! ----------------------------------------------------------------------
3277	! END SUBROUTINE REDPRM
3278	! ----------------------------------------------------------------------
3279	END SUBROUTINE REDPRM
3280
3281	SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3282
3283	IMPLICIT NONE
3284
3285	! ----------------------------------------------------------------------
3286	! SUBROUTINE ROSR12
3287	! ----------------------------------------------------------------------
3288	! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3289	! ### ### ### ### ### ###
3290	! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3291	! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3292	! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3293	! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3294	! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3295	! # . . # # . # = # . #
3296	! # . . # # . # # . #
3297	! # . . # # . # # . #
3298	! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3299	! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3300	! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3301	! ### ### ### ### ### ###
3302	! ----------------------------------------------------------------------
3303	INTEGER K
3304	INTEGER KK
3305	INTEGER NSOIL
3306
3307	REAL A(NSOIL)
3308	REAL B(NSOIL)
3309	REAL C(NSOIL)
3310	REAL D(NSOIL)
3311	REAL DELTA(NSOIL)
3312	REAL P(NSOIL)
3313
3314	! ----------------------------------------------------------------------
3315	! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3316	! ----------------------------------------------------------------------
3317	C(NSOIL) = 0.0
3318
3319	! ----------------------------------------------------------------------
3320	! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
3321	! ----------------------------------------------------------------------
3322	P(1) = -C(1) / B(1)
3323	DELTA(1) = D(1) / B(1)
3324
3325	! ----------------------------------------------------------------------
3326	! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
3327	! ----------------------------------------------------------------------
3328	DO K = 2,NSOIL
3329	P(K) = -C(K) * ( 1.0 / (B(K) + A (K) * P(K-1)) )
3330	DELTA(K) = (D(K)-A(K)DELTA(K-1))(1.0/(B(K)+A(K)*P(K-1)))
3331	END DO
3332
3333	! ----------------------------------------------------------------------
3334	! SET P TO DELTA FOR LOWEST SOIL LAYER
3335	! ----------------------------------------------------------------------
3336	P(NSOIL) = DELTA(NSOIL)
3337
3338	! ----------------------------------------------------------------------
3339	! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3340	! ----------------------------------------------------------------------
3341	DO K = 2,NSOIL
3342	KK = NSOIL - K + 1
3343	P(KK) = P(KK) * P(KK+1) + DELTA(KK)
3344	END DO
3345
3346	! ----------------------------------------------------------------------
3347	! END SUBROUTINE ROSR12
3348	! ----------------------------------------------------------------------
3349	END SUBROUTINE ROSR12
3350
3351	SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3352	& TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3353	& QUARTZ,CSOIL)
3354
3355	IMPLICIT NONE
3356
3357	! ----------------------------------------------------------------------
3358	! SUBROUTINE SHFLX
3359	! ----------------------------------------------------------------------
3360	! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3361	! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3362	! ON THE TEMPERATURE.
3363	! ----------------------------------------------------------------------
3364	INTEGER NSOLD
3365	PARAMETER(NSOLD = 20)
3366
3367	INTEGER I
3368	INTEGER ICE
3369	INTEGER IFRZ
3370	INTEGER NSOIL
3371
3372	REAL AI(NSOLD)
3373	REAL BI(NSOLD)
3374	REAL CI(NSOLD)
3375
3376	REAL BEXP
3377	REAL CSOIL
3378	REAL DF1
3379	REAL DT
3380	REAL F1
3381	REAL PSISAT
3382	REAL QUARTZ
3383	REAL RHSTS(NSOLD)
3384	REAL SSOIL
3385	REAL SH2O(NSOIL)
3386	REAL SMC(NSOIL)
3387	REAL SMCMAX
3388	REAL SMCWLT
3389	REAL STC(NSOIL)
3390	REAL STCF(NSOLD)
3391	REAL T0
3392	REAL T1
3393	REAL TBOT
3394	REAL YY
3395	REAL ZBOT
3396	REAL ZSOIL(NSOIL)
3397	REAL ZZ1
3398
3399	PARAMETER(T0 = 273.15)
3400
3401	! ----------------------------------------------------------------------
3402	! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3403	! ----------------------------------------------------------------------
3404	IF (ICE.EQ.1) THEN
3405
3406	! ----------------------------------------------------------------------
3407	! SEA-ICE CASE
3408	! ----------------------------------------------------------------------
3409	CALL HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
3410
3411	CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3412
3413	ELSE
3414
3415	! ----------------------------------------------------------------------
3416	! LAND-MASS CASE
3417	! ----------------------------------------------------------------------
3418	CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3419	& ZBOT,PSISAT,SH2O,DT, &
3420	& BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
3421
3422	CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3423
3424	ENDIF
3425
3426	DO I = 1,NSOIL
3427	STC(I) = STCF(I)
3428	END DO
3429
3430	! ----------------------------------------------------------------------
3431	! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
3432	! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
3433	! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
3434	! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
3435	! DIFFERENTLY IN ROUTINE SNOPAC)
3436	! ----------------------------------------------------------------------
3437	T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1
3438
3439	! ----------------------------------------------------------------------
3440	! CALCULATE SURFACE SOIL HEAT FLUX
3441	! ----------------------------------------------------------------------
3442	SSOIL = DF1 * (STC(1) - T1) / (0.5 * ZSOIL(1))
3443
3444	! ----------------------------------------------------------------------
3445	! END SUBROUTINE SHFLX
3446	! ----------------------------------------------------------------------
3447	END SUBROUTINE SHFLX
3448
3449	SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3450	& SH2O,SLOPE,KDT,FRZFACT, &
3451	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3452	& SHDFAC,CMCMAX, &
3453	& RUNOFF1,RUNOFF2,RUNOFF3, &
3454	& EDIR1,EC1,ET1, &
3455	& DRIP)
3456
3457	IMPLICIT NONE
3458
3459	! ----------------------------------------------------------------------
3460	! SUBROUTINE SMFLX
3461	! ----------------------------------------------------------------------
3462	! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
3463	! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
3464	! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3465	! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
3466	! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
3467	! ----------------------------------------------------------------------
3468	INTEGER NSOLD
3469	PARAMETER(NSOLD = 20)
3470
3471	INTEGER I
3472	INTEGER K
3473	INTEGER NSOIL
3474
3475	REAL AI(NSOLD)
3476	REAL BI(NSOLD)
3477	REAL CI(NSOLD)
3478
3479	REAL BEXP
3480	REAL CMC
3481	REAL CMCMAX
3482	REAL DKSAT
3483	REAL DRIP
3484	REAL DT
3485	REAL DUMMY
3486	REAL DWSAT
3487	REAL EC1
3488	REAL EDIR1
3489	REAL ET1(NSOIL)
3490	REAL EXCESS
3491	REAL FRZFACT
3492	REAL KDT
3493	REAL PCPDRP
3494	REAL PRCP1
3495	REAL RHSCT
3496	REAL RHSTT(NSOLD)
3497	REAL RUNOFF1
3498	REAL RUNOFF2
3499	REAL RUNOFF3
3500	REAL SHDFAC
3501	REAL SMC(NSOIL)
3502	REAL SH2O(NSOIL)
3503	REAL SICE(NSOLD)
3504	REAL SH2OA(NSOLD)
3505	REAL SH2OFG(NSOLD)
3506	REAL SLOPE
3507	REAL SMCMAX
3508	REAL SMCWLT
3509	REAL TRHSCT
3510	REAL ZSOIL(NSOIL)
3511
3512	! ----------------------------------------------------------------------
3513	! EXECUTABLE CODE BEGINS HERE.
3514	! ----------------------------------------------------------------------
3515	DUMMY = 0.
3516
3517	! ----------------------------------------------------------------------
3518	! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
3519	! ----------------------------------------------------------------------
3520	RHSCT = SHDFAC * PRCP1 - EC1
3521
3522	! ----------------------------------------------------------------------
3523	! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
3524	! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
3525	! FALL TO THE GRND.
3526	! ----------------------------------------------------------------------
3527	DRIP = 0.
3528	TRHSCT = DT * RHSCT
3529	EXCESS = CMC + TRHSCT
3530	IF (EXCESS .GT. CMCMAX) DRIP = EXCESS - CMCMAX
3531
3532	! ----------------------------------------------------------------------
3533	! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
3534	! SOIL
3535	! ----------------------------------------------------------------------
3536	PCPDRP = (1. - SHDFAC) * PRCP1 + DRIP / DT
3537
3538	! ----------------------------------------------------------------------
3539	! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT & SSTEP
3540	! ----------------------------------------------------------------------
3541	DO I = 1,NSOIL
3542	SICE(I) = SMC(I) - SH2O(I)
3543	END DO
3544
3545	! ----------------------------------------------------------------------
3546	! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3547	! TENDENCY EQUATIONS.
3548	!
3549	! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
3550	! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
3551	! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
3552	! THE FIRST SOIL LAYER)
3553	! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
3554	! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
3555	! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
3556	! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
3557	! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
3558	! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
3559	! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
3560	! SOIL MOISTURE STATE
3561	! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
3562	! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
3563	! OF SECTION 2 OF KALNAY AND KANAMITSU
3564	! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
3565	! ----------------------------------------------------------------------
3566	! IF ( PCPDRP .GT. 0.0 ) THEN
3567	IF ( (PCPDRPDT) .GT. (0.0011000.0(-ZSOIL(1))SMCMAX) ) THEN
3568
3569	! ----------------------------------------------------------------------
3570	! FROZEN GROUND VERSION:
3571	! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
3572	! INCLUDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
3573	! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
3574	! ----------------------------------------------------------------------
3575	CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3576	& DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3577	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3578
3579	CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3580	& CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3581
3582	DO K = 1,NSOIL
3583	SH2OA(K) = (SH2O(K) + SH2OFG(K)) * 0.5
3584	END DO
3585
3586	CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
3587	& DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3588	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3589
3590	CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3591	& CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3592
3593	ELSE
3594
3595	CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3596	& DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3597	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3598
3599	CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3600	& CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3601
3602	ENDIF
3603
3604	! RUNOF = RUNOFF
3605
3606	! ----------------------------------------------------------------------
3607	! END SUBROUTINE SMFLX
3608	! ----------------------------------------------------------------------
3609	END SUBROUTINE SMFLX
3610
3611	SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
3612
3613	IMPLICIT NONE
3614
3615	! ----------------------------------------------------------------------
3616	! SUBROUTINE SNFRAC
3617	! ----------------------------------------------------------------------
3618	! CALCULATE SNOW FRACTION (0 -> 1)
3619	! SNEQV SNOW WATER EQUIVALENT (M)
3620	! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
3621	! SALP TUNING PARAMETER
3622	! SNCOVR FRACTIONAL SNOW COVER
3623	! ----------------------------------------------------------------------
3624	REAL SNEQV, SNUP, SALP, SNCOVR, RSNOW, Z0N, SNOWH
3625
3626	! ----------------------------------------------------------------------
3627	! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
3628	! REDPRM) ABOVE WHICH SNOCVR=1.
3629	! ----------------------------------------------------------------------
3630	IF (SNEQV .LT. SNUP) THEN
3631	RSNOW = SNEQV/SNUP
3632	SNCOVR = 1. - ( EXP(-SALPRSNOW) - RSNOWEXP(-SALP))
3633	ELSE
3634	SNCOVR = 1.0
3635	ENDIF
3636
3637	Z0N=0.035
3638	! FORMULATION OF DICKINSON ET AL. 1986
3639
3640	! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
3641
3642	! FORMULATION OF MARSHALL ET AL. 1994
3643	! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
3644
3645	! ----------------------------------------------------------------------
3646	! END SUBROUTINE SNFRAC
3647	! ----------------------------------------------------------------------
3648	END SUBROUTINE SNFRAC
3649
3650	SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, &
3651	& SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
3652	& SBETA,DF1, &
3653	& Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
3654	& SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
3655	& SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT,SNUP, &
3656	& ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
3657	& RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
3658	& ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
3659	& BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW)
3660
3661	IMPLICIT NONE
3662
3663	! ----------------------------------------------------------------------
3664	! SUBROUTINE SNOPAC
3665	! ----------------------------------------------------------------------
3666	! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
3667	! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
3668	! PRESENT.
3669	! ----------------------------------------------------------------------
3670	INTEGER ICE
3671	INTEGER NROOT
3672	INTEGER NSOIL
3673
3674	LOGICAL SNOWNG
3675
3676	REAL BEXP
3677	REAL BETA
3678	REAL CFACTR
3679	REAL CMC
3680	REAL CMCMAX
3681	REAL CP
3682	REAL CPH2O
3683	REAL CPICE
3684	REAL CSOIL
3685	REAL DENOM
3686	REAL DEW
3687	REAL DF1
3688	REAL DKSAT
3689	REAL DRIP
3690	REAL DSOIL
3691	REAL DTOT
3692	REAL DT
3693	REAL DWSAT
3694	REAL EC
3695	REAL EDIR
3696	REAL EPSCA
3697	REAL ESD
3698	REAL ESDMIN
3699	REAL EXPSNO
3700	REAL EXPSOI
3701	REAL ETA
3702	REAL ETA1
3703	REAL ETP
3704	REAL ETP1
3705	REAL ETP2
3706	REAL ET(NSOIL)
3707	REAL ETT
3708	REAL EX
3709	REAL EXPFAC
3710	REAL FDOWN
3711	REAL FXEXP
3712	REAL FLX1
3713	REAL FLX2
3714	REAL FLX3
3715	REAL F1
3716	REAL KDT
3717	REAL LSUBF
3718	REAL LSUBC
3719	REAL LSUBS
3720	REAL PC
3721	REAL PRCP
3722	REAL PRCP1
3723	REAL Q2
3724	REAL RCH
3725	REAL RR
3726	REAL RTDIS(NSOIL)
3727	REAL SSOIL
3728	REAL SBETA
3729	REAL SSOIL1
3730	REAL SFCTMP
3731	REAL SHDFAC
3732	REAL SIGMA
3733	REAL SMC(NSOIL)
3734	REAL SH2O(NSOIL)
3735	REAL SMCDRY
3736	REAL SMCMAX
3737	REAL SMCREF
3738	REAL SMCWLT
3739	REAL SNOMLT
3740	REAL SNOWH
3741	REAL STC(NSOIL)
3742	REAL T1
3743	REAL T11
3744	REAL T12
3745	REAL T12A
3746	REAL T12B
3747	REAL T24
3748	REAL TBOT
3749	REAL ZBOT
3750	REAL TH2
3751	REAL YY
3752	REAL ZSOIL(NSOIL)
3753	REAL ZZ1
3754	REAL TFREEZ
3755	REAL SALP
3756	REAL SFCPRS
3757	REAL SLOPE
3758	REAL FRZFACT
3759	REAL PSISAT
3760	REAL SNUP
3761	REAL RUNOFF1
3762	REAL RUNOFF2
3763	REAL RUNOFF3
3764	REAL QUARTZ
3765	REAL SNDENS
3766	REAL SNCOND
3767	REAL RSNOW
3768	REAL SNCOVR
3769	REAL QSAT
3770	REAL ETP3
3771	REAL SEH
3772	REAL T14
3773	! REAL CSNOW
3774
3775	REAL EC1
3776	REAL EDIR1
3777	REAL ET1(NSOIL)
3778	REAL ETT1
3779
3780	REAL ETNS
3781	REAL ETNS1
3782	REAL ESNOW
3783	REAL ESNOW1
3784	REAL ESNOW2
3785	REAL ETANRG
3786
3787	INTEGER K
3788
3789	REAL SNOEXP
3790
3791	PARAMETER(CP = 1004.5)
3792	PARAMETER(CPH2O = 4.218E+3)
3793	PARAMETER(CPICE = 2.106E+3)
3794	PARAMETER(ESDMIN = 1.E-6)
3795	PARAMETER(LSUBF = 3.335E+5)
3796	PARAMETER(LSUBC = 2.501000E+6)
3797	PARAMETER(LSUBS = 2.83E+6)
3798	PARAMETER(SIGMA = 5.67E-8)
3799	PARAMETER(TFREEZ = 273.15)
3800
3801	! DATA SNOEXP /1.0/
3802	DATA SNOEXP /2.0/
3803
3804	! ----------------------------------------------------------------------
3805	! EXECUTABLE CODE BEGINS HERE:
3806	! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO M S-1 AND THEN TO AN
3807	! AMOUNT (M) GIVEN TIMESTEP (DT) AND CALL IT AN EFFECTIVE SNOWPACK
3808	! REDUCTION AMOUNT, ETP2 (M). THIS IS THE AMOUNT THE SNOWPACK WOULD BE
3809	! REDUCED DUE TO EVAPORATION FROM THE SNOW SFC DURING THE TIMESTEP.
3810	! EVAPORATION WILL PROCEED AT THE POTENTIAL RATE UNLESS THE SNOW DEPTH
3811	! IS LESS THAN THE EXPECTED SNOWPACK REDUCTION.
3812	! IF SEAICE (ICE=1), BETA REMAINS=1.
3813	! ----------------------------------------------------------------------
3814	PRCP1 = PRCP1*0.001
3815
3816	! ETP2 = ETP * 0.001 * DT
3817	BETA = 1.0
3818	IF (ICE .NE. 1) THEN
3819	IF (ESD .LT. ETP2) THEN
3820	! BETA = ESD / ETP2
3821	ENDIF
3822	ENDIF
3823
3824	! ----------------------------------------------------------------------
3825	EDIR = 0.0
3826	EDIR1 = 0.0
3827	EC = 0.0
3828	EC1 = 0.0
3829	DO K = 1,NSOIL
3830	ET(K) = 0.0
3831	ET1(K) = 0.0
3832	ENDDO
3833	ETT = 0.0
3834	ETT1 = 0.0
3835	ETNS = 0.0
3836	ETNS1 = 0.0
3837	ESNOW = 0.0
3838	ESNOW1 = 0.0
3839	ESNOW2 = 0.0
3840	! ----------------------------------------------------------------------
3841
3842	! ----------------------------------------------------------------------
3843	! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
3844	! ----------------------------------------------------------------------
3845	DEW = 0.0
3846	ETP1 = ETP*0.001
3847	IF (ETP .LT. 0.0) THEN
3848	! DEW = -ETP * 0.001
3849	DEW = -ETP1
3850	! ESNOW2 = ETP * 0.001 * DT
3851	ESNOW2 = ETP1 * DT
3852	ETANRG = ETP((1.-SNCOVR)LSUBC + SNCOVR*LSUBS)
3853	! ENDIF
3854	ELSE
3855	! ----------------------------------------------------------------------
3856	! ETP1 = 0.0
3857	! ETP1 = ETP*0.001
3858	IF (ICE .NE. 1) THEN
3859	IF (SNCOVR .LT. 1.) THEN
3860	! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
3861	CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3862	& SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3863	& SMCREF,SHDFAC,CMCMAX, &
3864	& SMCDRY,CFACTR, &
3865	& EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
3866	! ENDIF
3867	! ----------------------------------------------------------------------
3868	EDIR1 = EDIR1*(1.-SNCOVR)
3869	EC1 = EC1*(1.-SNCOVR)
3870	DO K = 1,NSOIL
3871	ET1(K) = ET1(K)*(1.-SNCOVR)
3872	END DO
3873	ETT1 = ETT1*(1.-SNCOVR)
3874	ETNS1 = ETNS1*(1.-SNCOVR)
3875	! ----------------------------------------------------------------------
3876	EDIR = EDIR1 * 1000.0
3877	EC = EC1 * 1000.0
3878	DO K = 1,NSOIL
3879	ET(K) = ET1(K) * 1000.0
3880	END DO
3881	ETT = ETT1 * 1000.0
3882	ETNS = ETNS1 * 1000.0
3883	! ----------------------------------------------------------------------
3884	ENDIF
3885	ESNOW = ETP*SNCOVR
3886	! ESNOW1 = ETP*0.001
3887	ESNOW1 = ESNOW*0.001
3888	ESNOW2 = ESNOW1*DT
3889	ETANRG = ESNOWLSUBS + ETNSLSUBC
3890	ENDIF
3891	ENDIF
3892
3893	! ----------------------------------------------------------------------
3894	! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3895	! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3896	! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3897	! SNOWFALL STRIKING THE GOUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3898	! ----------------------------------------------------------------------
3899	FLX1 = 0.0
3900	IF (SNOWNG) THEN
3901	FLX1 = CPICE * PRCP * (T1 - SFCTMP)
3902	ELSE
3903	IF (PRCP .GT. 0.0) FLX1 = CPH2O * PRCP * (T1 - SFCTMP)
3904	ENDIF
3905
3906	! ----------------------------------------------------------------------
3907	! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3908	! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3909	! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3910	! FLUXES.
3911	! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
3912	! PENMAN.
3913	! ----------------------------------------------------------------------
3914	DSOIL = -(0.5 * ZSOIL(1))
3915	DTOT = SNOWH + DSOIL
3916	DENOM = 1.0 + DF1 / (DTOT * RR * RCH)
3917	! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3918	! & + TH2 - SFCTMP - BETA*EPSCA ) / RR
3919	! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3920	! M.Ek, 24Nov04, add snow emissivity
3921	T12A = ((FDOWN-FLX1-FLX2 &
3922	& -(0.95SNCOVR+(1.0-SNCOVR))SIGMA*T24)/RCH &
3923	& + TH2 - SFCTMP - ETANRG/RCH ) / RR
3924	T12B = DF1 * STC(1) / (DTOT * RR * RCH)
3925	T12 = (SFCTMP + T12A + T12B) / DENOM
3926
3927	! ----------------------------------------------------------------------
3928	! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
3929	! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
3930	! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
3931	! DEPENDING ON SIGN OF ETP.
3932	! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
3933	! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
3934	! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
3935	! TO ZERO.
3936	! ----------------------------------------------------------------------
3937	IF (T12 .LE. TFREEZ) THEN
3938	T1 = T12
3939	SSOIL = DF1 * (T1 - STC(1)) / DTOT
3940	! ESD = MAX(0.0, ESD-ETP2)
3941	ESD = MAX(0.0, ESD-ESNOW2)
3942	FLX3 = 0.0
3943	EX = 0.0
3944	SNOMLT = 0.0
3945
3946	ELSE
3947	! ----------------------------------------------------------------------
3948	! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3949	! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
3950	! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
3951	! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
3952	! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
3953	! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
3954	! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
3955	! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
3956	! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
3957	! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
3958	! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
3959	! ----------------------------------------------------------------------
3960	! T1 = TFREEZ * SNCOVR + T12 * (1.0 - SNCOVR)
3961	! mek Feb2004
3962	! non-linear weighting of snow vs non-snow covered portions of gridbox
3963	! so with SNOEXP = 2.0 (>1), surface skin temperature is higher than for
3964	! the linear case (SNOEXP = 1).
3965	T1 = TFREEZ * SNCOVR*SNOEXP + T12 (1.0 - SNCOVR**SNOEXP)
3966	! QSAT = (0.6226.11E2)/(SFCPRS-0.3786.11E2)
3967	! ETP = RCH*(QSAT-Q2)/CP
3968	! ETP2 = ETP0.001DT
3969	BETA = 1.0
3970	SSOIL = DF1 * (T1 - STC(1)) / DTOT
3971
3972	! ----------------------------------------------------------------------
3973	! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
3974	! BETA<1
3975	! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
3976	! ----------------------------------------------------------------------
3977	! IF (ESD .LE. ETP2) THEN
3978	! IF (ESD .LE. ESNOW2) THEN
3979	IF (ESD-ESNOW2 .LE. ESDMIN) THEN
3980	! BETA = ESD / ETP2
3981	ESD = 0.0
3982	EX = 0.0
3983	SNOMLT = 0.0
3984
3985	ELSE
3986	! ----------------------------------------------------------------------
3987	! POTENTIAL EVAP (SUBLIMATION) LESS THAN DEPTH OF SNOWPACK, RETAIN
3988	! BETA=1.
3989	! SNOWPACK (ESD) REDUCED BY POTENTIAL EVAP RATE
3990	! ETP3 (CONVERT TO FLUX)
3991	! ----------------------------------------------------------------------
3992	! ESD = ESD-ETP2
3993	ESD = ESD-ESNOW2
3994	! ETP3 = ETP*LSUBC
3995	SEH = RCH*(T1-TH2)
3996	T14 = T1*T1
3997	T14 = T14*T14
3998	! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETP3
3999	! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETANRG
4000	! M.Ek, 24Nov04, add snow emissivity
4001	FLX3 = FDOWN - FLX1 - FLX2 - &
4002	& (0.95SNCOVR+(1.0-SNCOVR))SIGMA*T14 - SSOIL - SEH - ETANRG
4003	IF (FLX3 .LE. 0.0) FLX3 = 0.0
4004	EX = FLX3*0.001/LSUBF
4005
4006	! ----------------------------------------------------------------------
4007	! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
4008	! IF SNOW COVER LESS THAN 5% NO SNOWMELT REDUCTION
4009	! ***NOTE: DOES 'IF' BELOW FAIL TO MATCH THE MELT WATER WITH THE MELT
4010	! ENERGY?
4011	! ----------------------------------------------------------------------
4012	! IF (SNCOVR .GT. 0.05) EX = EX * SNCOVR
4013	SNOMLT = EX * DT
4014
4015	! ----------------------------------------------------------------------
4016	! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
4017	! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
4018	! ----------------------------------------------------------------------
4019	IF (ESD-SNOMLT .GE. ESDMIN) THEN
4020	ESD = ESD - SNOMLT
4021
4022	ELSE
4023	! ----------------------------------------------------------------------
4024	! SNOWMELT EXCEEDS SNOW DEPTH
4025	! ----------------------------------------------------------------------
4026	EX = ESD/DT
4027	FLX3 = EX1000.0LSUBF
4028	SNOMLT = ESD
4029	ESD = 0.0
4030
4031	ENDIF
4032	! ----------------------------------------------------------------------
4033	! END OF 'ESD .LE. ETP2' IF-BLOCK
4034	! ----------------------------------------------------------------------
4035	ENDIF
4036
4037	PRCP1 = PRCP1 + EX
4038
4039	! ----------------------------------------------------------------------
4040	! END OF 'T12 .LE. TFREEZ' IF-BLOCK
4041	! ----------------------------------------------------------------------
4042	ENDIF
4043
4044	! ----------------------------------------------------------------------
4045	! FINAL BETA NOW IN HAND, SO COMPUTE EVAPORATION. EVAP EQUALS ETP
4046	! UNLESS BETA<1.
4047	! ----------------------------------------------------------------------
4048	! ETA = BETA*ETP
4049
4050	! ----------------------------------------------------------------------
4051	! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
4052	! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
4053	! (BELOW).
4054	! IF SEAICE (ICE=1) SKIP CALL TO SMFLX.
4055	! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES. IN THIS, THE SNOW PACK
4056	! CASE, ETA1 IS NOT USED IN CALCULATION OF EVAP.
4057	! ----------------------------------------------------------------------
4058	! ETP1 = 0.0
4059	IF (ICE .NE. 1) THEN
4060	! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
4061	! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
4062	! & SMCREF,SHDFAC,CMCMAX,
4063	! & SMCDRY,CFACTR,
4064	! & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
4065	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
4066	& SH2O,SLOPE,KDT,FRZFACT, &
4067	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
4068	& SHDFAC,CMCMAX, &
4069	& RUNOFF1,RUNOFF2,RUNOFF3, &
4070	& EDIR1,EC1,ET1, &
4071	& DRIP)
4072
4073	ENDIF
4074
4075	! ----------------------------------------------------------------------
4076	! EDIR = EDIR1 * 1000.0
4077	! EC = EC1 * 1000.0
4078	! ETT = ETT1 * 1000.0
4079	! ET(1) = ET1(1) * 1000.0
4080	! ET(2) = ET1(2) * 1000.0
4081	! ET(3) = ET1(3) * 1000.0
4082	! ET(4) = ET1(4) * 1000.0
4083	! ----------------------------------------------------------------------
4084
4085	! ----------------------------------------------------------------------
4086	! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
4087	! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
4088	! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
4089	! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
4090	! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
4091	! SKIN TEMP VALUE AS REVISED BY SHFLX.
4092	! ----------------------------------------------------------------------
4093	ZZ1 = 1.0
4094	YY = STC(1)-0.5SSOILZSOIL(1)*ZZ1/DF1
4095	T11 = T1
4096
4097	! ----------------------------------------------------------------------
4098	! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
4099	! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
4100	! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
4101	! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
4102	! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
4103	! ----------------------------------------------------------------------
4104	CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
4105	& TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
4106	& QUARTZ,CSOIL)
4107
4108	! ----------------------------------------------------------------------
4109	! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
4110	! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
4111	! ----------------------------------------------------------------------
4112	IF (ESD .GT. 0.) THEN
4113	CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
4114	ELSE
4115	ESD = 0.
4116	SNOWH = 0.
4117	SNDENS = 0.
4118	SNCOND = 1.
4119	SNCOVR = 0.
4120	ENDIF
4121
4122	! ----------------------------------------------------------------------
4123	! END SUBROUTINE SNOPAC
4124	! ----------------------------------------------------------------------
4125	END SUBROUTINE SNOPAC
4126
4127	SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
4128
4129	IMPLICIT NONE
4130
4131	! ----------------------------------------------------------------------
4132	! SUBROUTINE SNOWPACK
4133	! ----------------------------------------------------------------------
4134	! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
4135	! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
4136	! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
4137	! KOREN, 03/25/95.
4138	! ----------------------------------------------------------------------
4139	! ESD WATER EQUIVALENT OF SNOW (M)
4140	! DTSEC TIME STEP (SEC)
4141	! SNOWH SNOW DEPTH (M)
4142	! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4143	! TSNOW SNOW SURFACE TEMPERATURE (K)
4144	! TSOIL SOIL SURFACE TEMPERATURE (K)
4145	!
4146	! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
4147	! ----------------------------------------------------------------------
4148	INTEGER IPOL, J
4149
4150	REAL BFAC,C1,C2,SNDENS,DSX,DTHR,DTSEC,DW,SNOWHC,SNOWH,PEXP,TAVGC, &
4151	& TSNOW,TSNOWC,TSOIL,TSOILC,ESD,ESDC,ESDCX,G,KN
4152
4153	PARAMETER(C1 = 0.01, C2=21.0, G=9.81, KN=4000.0)
4154
4155	! ----------------------------------------------------------------------
4156	! CONVERSION INTO SIMULATION UNITS
4157	! ----------------------------------------------------------------------
4158	SNOWHC = SNOWH*100.
4159	ESDC = ESD*100.
4160	DTHR = DTSEC/3600.
4161	TSNOWC = TSNOW-273.15
4162	TSOILC = TSOIL-273.15
4163
4164	! ----------------------------------------------------------------------
4165	! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
4166	! ----------------------------------------------------------------------
4167	TAVGC = 0.5*(TSNOWC+TSOILC)
4168
4169	! ----------------------------------------------------------------------
4170	! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
4171	! SNDENS=DS0(EXP(BFACESD)-1.)/(BFAC*ESD)
4172	! BFAC=DTHRC1EXP(0.08TAVGC-C2DS0)
4173	! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
4174	! NUMERICALLY BELOW:
4175	! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
4176	! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
4177	! ----------------------------------------------------------------------
4178	IF (ESDC .GT. 1.E-2) THEN
4179	ESDCX = ESDC
4180	ELSE
4181	ESDCX = 1.E-2
4182	ENDIF
4183	BFAC = DTHRC1EXP(0.08TAVGC-C2SNDENS)
4184
4185	! DSX = SNDENS((DEXP(BFACESDC)-1.)/(BFAC*ESDC))
4186	! ----------------------------------------------------------------------
4187	! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
4188	! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
4189	! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
4190	! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
4191	! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
4192	! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
4193	! POLYNOMIAL EXPANSION.
4194	!
4195	! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
4196	! IS GOVERNED BY ITERATION LIMIT "IPOL".
4197	! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
4198	! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
4199	! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
4200	! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
4201	! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
4202	! ----------------------------------------------------------------------
4203	IPOL = 4
4204	PEXP = 0.
4205	DO J = IPOL,1,-1
4206	! PEXP = (1. + PEXP)BFACESDC/REAL(J+1)
4207	PEXP = (1. + PEXP)BFACESDCX/REAL(J+1)
4208	END DO
4209	PEXP = PEXP + 1.
4210
4211	DSX = SNDENS*(PEXP)
4212	! ----------------------------------------------------------------------
4213	! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
4214	! ----------------------------------------------------------------------
4215	! END OF KOREAN FORMULATION
4216
4217	! BASE FORMULATION (COGLEY ET AL., 1990)
4218	! CONVERT DENSITY FROM G/CM3 TO KG/M3
4219	! DSM=SNDENS*1000.0
4220
4221	! DSX=DSM+DTSEC0.5DSMGESD/
4222	! & (1E7EXP(-0.02DSM+KN/(TAVGC+273.16)-14.643))
4223
4224	! CONVERT DENSITY FROM KG/M3 TO G/CM3
4225	! DSX=DSX/1000.0
4226
4227	! END OF COGLEY ET AL. FORMULATION
4228
4229	! ----------------------------------------------------------------------
4230	! SET UPPER/LOWER LIMIT ON SNOW DENSITY
4231	! ----------------------------------------------------------------------
4232	IF (DSX .GT. 0.40) DSX = 0.40
4233	IF (DSX .LT. 0.05) DSX = 0.05
4234	SNDENS = DSX
4235	! ----------------------------------------------------------------------
4236	! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
4237	! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
4238	! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
4239	! ----------------------------------------------------------------------
4240	IF (TSNOWC .GE. 0.) THEN
4241	DW = 0.13*DTHR/24.
4242	SNDENS = SNDENS*(1.-DW)+DW
4243	IF (SNDENS .GT. 0.40) SNDENS = 0.40
4244	ENDIF
4245
4246	! ----------------------------------------------------------------------
4247	! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
4248	! CHANGE SNOW DEPTH UNITS TO METERS
4249	! ----------------------------------------------------------------------
4250	SNOWHC = ESDC/SNDENS
4251	SNOWH = SNOWHC*0.01
4252
4253	! ----------------------------------------------------------------------
4254	! END SUBROUTINE SNOWPACK
4255	! ----------------------------------------------------------------------
4256	END SUBROUTINE SNOWPACK
4257
4258	SUBROUTINE SNOWZ0 (SNCOVR,Z0)
4259
4260	IMPLICIT NONE
4261
4262	! ----------------------------------------------------------------------
4263	! SUBROUTINE SNOWZ0
4264	! ----------------------------------------------------------------------
4265	! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
4266	! SNCOVR FRACTIONAL SNOW COVER
4267	! Z0 ROUGHNESS LENGTH (m)
4268	! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
4269	! ----------------------------------------------------------------------
4270	REAL SNCOVR, Z0, Z0S
4271	! PARAMETER (Z0S=0.001)
4272
4273	! CURRENT NOAH LSM CONDITION - MBEK, 09-OCT-2001
4274	Z0S = Z0
4275	!
4276	Z0 = (1-SNCOVR)Z0 + SNCOVRZ0S
4277	! ----------------------------------------------------------------------
4278	! END SUBROUTINE SNOWZ0
4279	! ----------------------------------------------------------------------
4280	END SUBROUTINE SNOWZ0
4281
4282	SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
4283
4284	IMPLICIT NONE
4285
4286	! ----------------------------------------------------------------------
4287	! SUBROUTINE SNOW_NEW
4288	! ----------------------------------------------------------------------
4289	! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
4290	! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
4291	!
4292	! TEMP AIR TEMPERATURE (K)
4293	! NEWSN NEW SNOWFALL (M)
4294	! SNOWH SNOW DEPTH (M)
4295	! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4296	! ----------------------------------------------------------------------
4297	REAL SNDENS
4298	REAL DSNEW
4299	REAL SNOWHC
4300	REAL HNEWC
4301	REAL SNOWH
4302	REAL NEWSN
4303	REAL NEWSNC
4304	REAL TEMP
4305	REAL TEMPC
4306
4307	! ----------------------------------------------------------------------
4308	! CONVERSION INTO SIMULATION UNITS
4309	! ----------------------------------------------------------------------
4310	SNOWHC = SNOWH*100.
4311	NEWSNC = NEWSN*100.
4312	TEMPC = TEMP-273.15
4313
4314	! ----------------------------------------------------------------------
4315	! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
4316	! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
4317	! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
4318	! VEMADOLEN, SWEDEN, 1980, 172-177PP.
4319	!-----------------------------------------------------------------------
4320	IF (TEMPC .LE. -15.) THEN
4321	DSNEW = 0.05
4322	ELSE
4323	DSNEW = 0.05+0.0017(TEMPC+15.)*1.5
4324	ENDIF
4325
4326	! ----------------------------------------------------------------------
4327	! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
4328	! ----------------------------------------------------------------------
4329	HNEWC = NEWSNC/DSNEW
4330	SNDENS = (SNOWHCSNDENS+HNEWCDSNEW)/(SNOWHC+HNEWC)
4331	SNOWHC = SNOWHC+HNEWC
4332	SNOWH = SNOWHC*0.01
4333
4334	! ----------------------------------------------------------------------
4335	! END SUBROUTINE SNOW_NEW
4336	! ----------------------------------------------------------------------
4337	END SUBROUTINE SNOW_NEW
4338
4339	SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
4340	& ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
4341	& RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
4342
4343	IMPLICIT NONE
4344
4345	! ----------------------------------------------------------------------
4346	! SUBROUTINE SRT
4347	! ----------------------------------------------------------------------
4348	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
4349	! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
4350	! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
4351	! ----------------------------------------------------------------------
4352	INTEGER NSOLD
4353	PARAMETER(NSOLD = 20)
4354
4355	INTEGER CVFRZ
4356	INTEGER IALP1
4357	INTEGER IOHINF
4358	INTEGER J
4359	INTEGER JJ
4360	INTEGER K
4361	INTEGER KS
4362	INTEGER NSOIL
4363
4364	REAL ACRT
4365	REAL AI(NSOLD)
4366	REAL BEXP
4367	REAL BI(NSOLD)
4368	REAL CI(NSOLD)
4369	REAL DD
4370	REAL DDT
4371	REAL DDZ
4372	REAL DDZ2
4373	REAL DENOM
4374	REAL DENOM2
4375	REAL DICE
4376	REAL DKSAT
4377	REAL DMAX(NSOLD)
4378	REAL DSMDZ
4379	REAL DSMDZ2
4380	REAL DT
4381	REAL DT1
4382	REAL DWSAT
4383	REAL EDIR
4384	REAL ET(NSOIL)
4385	REAL FCR
4386	REAL FRZX
4387	REAL INFMAX
4388	REAL KDT
4389	REAL MXSMC
4390	REAL MXSMC2
4391	REAL NUMER
4392	REAL PCPDRP
4393	REAL PDDUM
4394	REAL PX
4395	REAL RHSTT(NSOIL)
4396	REAL RUNOFF1
4397	REAL RUNOFF2
4398	REAL SH2O(NSOIL)
4399	REAL SH2OA(NSOIL)
4400	REAL SICE(NSOIL)
4401	REAL SICEMAX
4402	REAL SLOPE
4403	REAL SLOPX
4404	REAL SMCAV
4405	REAL SMCMAX
4406	REAL SMCWLT
4407	REAL SSTT
4408	REAL SUM
4409	REAL VAL
4410	REAL WCND
4411	REAL WCND2
4412	REAL WDF
4413	REAL WDF2
4414	REAL ZSOIL(NSOIL)
4415
4416	! ----------------------------------------------------------------------
4417	! FROZEN GROUND VERSION:
4418	! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4419	! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4420	! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
4421	! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
4422	! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
4423	! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
4424	! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
4425	! ----------------------------------------------------------------------
4426	PARAMETER(CVFRZ = 3)
4427
4428	! ----------------------------------------------------------------------
4429	! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
4430	! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
4431	! MODIFIED BY Q DUAN
4432	! ----------------------------------------------------------------------
4433	IOHINF=1
4434
4435	! ----------------------------------------------------------------------
4436	! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
4437	! LAYERS.
4438	! ----------------------------------------------------------------------
4439	SICEMAX = 0.0
4440	DO KS=1,NSOIL
4441	IF (SICE(KS) .GT. SICEMAX) SICEMAX = SICE(KS)
4442	END DO
4443
4444	! ----------------------------------------------------------------------
4445	! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
4446	! ----------------------------------------------------------------------
4447	PDDUM = PCPDRP
4448	RUNOFF1 = 0.0
4449	IF (PCPDRP .NE. 0.0) THEN
4450
4451	! ----------------------------------------------------------------------
4452	! MODIFIED BY Q. DUAN, 5/16/94
4453	! ----------------------------------------------------------------------
4454	! IF (IOHINF .EQ. 1) THEN
4455
4456	DT1 = DT/86400.
4457	SMCAV = SMCMAX - SMCWLT
4458	DMAX(1)=-ZSOIL(1)*SMCAV
4459
4460	! ----------------------------------------------------------------------
4461	! FROZEN GROUND VERSION:
4462	! ----------------------------------------------------------------------
4463	DICE = -ZSOIL(1) * SICE(1)
4464
4465	DMAX(1)=DMAX(1)*(1.0 - (SH2OA(1)+SICE(1)-SMCWLT)/SMCAV)
4466	DD=DMAX(1)
4467
4468	DO KS=2,NSOIL
4469
4470	! ----------------------------------------------------------------------
4471	! FROZEN GROUND VERSION:
4472	! ----------------------------------------------------------------------
4473	DICE = DICE + ( ZSOIL(KS-1) - ZSOIL(KS) ) * SICE(KS)
4474
4475	DMAX(KS) = (ZSOIL(KS-1)-ZSOIL(KS))*SMCAV
4476	DMAX(KS) = DMAX(KS)*(1.0 - (SH2OA(KS)+SICE(KS)-SMCWLT)/SMCAV)
4477	DD = DD+DMAX(KS)
4478	END DO
4479
4480	! ----------------------------------------------------------------------
4481	! VAL = (1.-EXP(-KDT*SQRT(DT1)))
4482	! IN BELOW, REMOVE THE SQRT IN ABOVE
4483	! ----------------------------------------------------------------------
4484	VAL = (1.-EXP(-KDT*DT1))
4485	DDT = DD*VAL
4486	PX = PCPDRP*DT
4487	IF (PX .LT. 0.0) PX = 0.0
4488	INFMAX = (PX*(DDT/(PX+DDT)))/DT
4489
4490	! ----------------------------------------------------------------------
4491	! FROZEN GROUND VERSION:
4492	! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4493	! ----------------------------------------------------------------------
4494	FCR = 1.
4495	IF (DICE .GT. 1.E-2) THEN
4496	ACRT = CVFRZ * FRZX / DICE
4497	SUM = 1.
4498	IALP1 = CVFRZ - 1
4499	DO J = 1,IALP1
4500	K = 1
4501	DO JJ = J+1,IALP1
4502	K = K * JJ
4503	END DO
4504	SUM = SUM + (ACRT ** ( CVFRZ-J)) / FLOAT (K)
4505	END DO
4506	FCR = 1. - EXP(-ACRT) * SUM
4507	ENDIF
4508	INFMAX = INFMAX * FCR
4509
4510	! ----------------------------------------------------------------------
4511	! CORRECTION OF INFILTRATION LIMITATION:
4512	! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
4513	! HYDROLIC CONDUCTIVITY
4514	! ----------------------------------------------------------------------
4515	! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
4516	MXSMC = SH2OA(1)
4517
4518	CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4519	& SICEMAX)
4520
4521	INFMAX = MAX(INFMAX,WCND)
4522	INFMAX = MIN(INFMAX,PX)
4523
4524	IF (PCPDRP .GT. INFMAX) THEN
4525	RUNOFF1 = PCPDRP - INFMAX
4526	PDDUM = INFMAX
4527	ENDIF
4528
4529	ENDIF
4530
4531	! ----------------------------------------------------------------------
4532	! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
4533	! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4534	! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
4535	! ----------------------------------------------------------------------
4536	MXSMC = SH2OA(1)
4537
4538	CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4539	& SICEMAX)
4540
4541	! ----------------------------------------------------------------------
4542	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
4543	! ----------------------------------------------------------------------
4544	DDZ = 1. / ( -.5 * ZSOIL(2) )
4545	AI(1) = 0.0
4546	BI(1) = WDF * DDZ / ( -ZSOIL(1) )
4547	CI(1) = -BI(1)
4548
4549	! ----------------------------------------------------------------------
4550	! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
4551	! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
4552	! ----------------------------------------------------------------------
4553	DSMDZ = ( SH2O(1) - SH2O(2) ) / ( -.5 * ZSOIL(2) )
4554	RHSTT(1) = (WDF * DSMDZ + WCND - PDDUM + EDIR + ET(1))/ZSOIL(1)
4555	SSTT = WDF * DSMDZ + WCND + EDIR + ET(1)
4556
4557	! ----------------------------------------------------------------------
4558	! INITIALIZE DDZ2
4559	! ----------------------------------------------------------------------
4560	DDZ2 = 0.0
4561
4562	! ----------------------------------------------------------------------
4563	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
4564	! ----------------------------------------------------------------------
4565	DO K = 2,NSOIL
4566	DENOM2 = (ZSOIL(K-1) - ZSOIL(K))
4567	IF (K .NE. NSOIL) THEN
4568	SLOPX = 1.
4569
4570	! ----------------------------------------------------------------------
4571	! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
4572	! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4573	! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
4574	! ----------------------------------------------------------------------
4575	MXSMC2 = SH2OA(K)
4576
4577	CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
4578	& SICEMAX)
4579
4580	! ----------------------------------------------------------------------
4581	! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
4582	! ----------------------------------------------------------------------
4583	DENOM = (ZSOIL(K-1) - ZSOIL(K+1))
4584	DSMDZ2 = (SH2O(K) - SH2O(K+1)) / (DENOM * 0.5)
4585
4586	! ----------------------------------------------------------------------
4587	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
4588	! ----------------------------------------------------------------------
4589	DDZ2 = 2.0 / DENOM
4590	CI(K) = -WDF2 * DDZ2 / DENOM2
4591	ELSE
4592
4593	! ----------------------------------------------------------------------
4594	! SLOPE OF BOTTOM LAYER IS INTRODUCED
4595	! ----------------------------------------------------------------------
4596	SLOPX = SLOPE
4597
4598	! ----------------------------------------------------------------------
4599	! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
4600	! THIS LAYER
4601	! ----------------------------------------------------------------------
4602	CALL WDFCND (WDF2,WCND2,SH2OA(NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
4603	& SICEMAX)
4604
4605	! ----------------------------------------------------------------------
4606	! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
4607	! ----------------------------------------------------------------------
4608	DSMDZ2 = 0.0
4609
4610	! ----------------------------------------------------------------------
4611	! SET MATRIX COEF CI TO ZERO
4612	! ----------------------------------------------------------------------
4613	CI(K) = 0.0
4614	ENDIF
4615
4616	! ----------------------------------------------------------------------
4617	! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
4618	! ----------------------------------------------------------------------
4619	NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2 - (WDF * DSMDZ) &
4620	& - WCND + ET(K)
4621	RHSTT(K) = NUMER / (-DENOM2)
4622
4623	! ----------------------------------------------------------------------
4624	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
4625	! ----------------------------------------------------------------------
4626	AI(K) = -WDF * DDZ / DENOM2
4627	BI(K) = -( AI(K) + CI(K) )
4628
4629	! ----------------------------------------------------------------------
4630	! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
4631	! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
4632	! ----------------------------------------------------------------------
4633	IF (K .EQ. NSOIL) THEN
4634	RUNOFF2 = SLOPX * WCND2
4635	ENDIF
4636
4637	IF (K .NE. NSOIL) THEN
4638	WDF = WDF2
4639	WCND = WCND2
4640	DSMDZ = DSMDZ2
4641	DDZ = DDZ2
4642	ENDIF
4643	END DO
4644
4645	! ----------------------------------------------------------------------
4646	! END SUBROUTINE SRT
4647	! ----------------------------------------------------------------------
4648	END SUBROUTINE SRT
4649
4650	SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
4651	& NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
4652	& AI,BI,CI)
4653
4654	IMPLICIT NONE
4655
4656	! ----------------------------------------------------------------------
4657	! SUBROUTINE SSTEP
4658	! ----------------------------------------------------------------------
4659	! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
4660	! CONTENT VALUES.
4661	! ----------------------------------------------------------------------
4662	INTEGER NSOLD
4663	PARAMETER(NSOLD = 20)
4664
4665	INTEGER I
4666	INTEGER K
4667	INTEGER KK11
4668	INTEGER NSOIL
4669
4670	REAL AI(NSOLD)
4671	REAL BI(NSOLD)
4672	REAL CI(NSOLD)
4673	REAL CIin(NSOLD)
4674	REAL CMC
4675	REAL CMCMAX
4676	REAL DDZ
4677	REAL DT
4678	REAL RHSCT
4679	REAL RHSTT(NSOIL)
4680	REAL RHSTTin(NSOIL)
4681	REAL RUNOFF3
4682	REAL SH2OIN(NSOIL)
4683	REAL SH2OOUT(NSOIL)
4684	REAL SICE(NSOIL)
4685	REAL SMC(NSOIL)
4686	REAL SMCMAX
4687	REAL STOT
4688	REAL WPLUS
4689	REAL ZSOIL(NSOIL)
4690
4691	! ----------------------------------------------------------------------
4692	! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
4693	! TRI-DIAGONAL MATRIX ROUTINE.
4694	! ----------------------------------------------------------------------
4695	DO K = 1,NSOIL
4696	RHSTT(K) = RHSTT(K) * DT
4697	AI(K) = AI(K) * DT
4698	BI(K) = 1. + BI(K) * DT
4699	CI(K) = CI(K) * DT
4700	END DO
4701
4702	! ----------------------------------------------------------------------
4703	! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4704	! ----------------------------------------------------------------------
4705	DO K = 1,NSOIL
4706	RHSTTin(K) = RHSTT(K)
4707	END DO
4708	DO K = 1,NSOIL
4709	CIin(K) = CI(K)
4710	END DO
4711
4712	! ----------------------------------------------------------------------
4713	! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
4714	! ----------------------------------------------------------------------
4715	CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
4716
4717	! ----------------------------------------------------------------------
4718	! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
4719	! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
4720	! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
4721	! ----------------------------------------------------------------------
4722	WPLUS = 0.0
4723	RUNOFF3 = 0.
4724	DDZ = -ZSOIL(1)
4725
4726	DO K = 1,NSOIL
4727	IF (K .NE. 1) DDZ = ZSOIL(K - 1) - ZSOIL(K)
4728	SH2OOUT(K) = SH2OIN(K) + CI(K) + WPLUS / DDZ
4729
4730	STOT = SH2OOUT(K) + SICE(K)
4731	IF (STOT .GT. SMCMAX) THEN
4732	IF (K .EQ. 1) THEN
4733	DDZ = -ZSOIL(1)
4734	ELSE
4735	KK11 = K - 1
4736	DDZ = -ZSOIL(K) + ZSOIL(KK11)
4737	ENDIF
4738	WPLUS = (STOT-SMCMAX) * DDZ
4739	ELSE
4740	WPLUS = 0.
4741	ENDIF
4742	SMC(K) = MAX ( MIN(STOT,SMCMAX),0.02 )
4743	SH2OOUT(K) = MAX((SMC(K)-SICE(K)),0.0)
4744	END DO
4745
4746	RUNOFF3 = WPLUS
4747
4748	! ----------------------------------------------------------------------
4749	! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
4750	! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
4751	! ----------------------------------------------------------------------
4752	CMC = CMC + DT * RHSCT
4753	IF (CMC .LT. 1.E-20) CMC=0.0
4754	CMC = MIN(CMC,CMCMAX)
4755
4756	! ----------------------------------------------------------------------
4757	! END SUBROUTINE SSTEP
4758	! ----------------------------------------------------------------------
4759	END SUBROUTINE SSTEP
4760
4761	SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4762
4763	IMPLICIT NONE
4764
4765	! ----------------------------------------------------------------------
4766	! SUBROUTINE TBND
4767	! ----------------------------------------------------------------------
4768	! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4769	! THE MIDDLE LAYER TEMPERATURES
4770	! ----------------------------------------------------------------------
4771	INTEGER NSOIL
4772	INTEGER K
4773
4774	REAL TBND1
4775	REAL T0
4776	REAL TU
4777	REAL TB
4778	REAL ZB
4779	REAL ZBOT
4780	REAL ZUP
4781	REAL ZSOIL (NSOIL)
4782
4783	PARAMETER(T0 = 273.15)
4784
4785	! ----------------------------------------------------------------------
4786	! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4787	! ----------------------------------------------------------------------
4788	IF (K .EQ. 1) THEN
4789	ZUP = 0.
4790	ELSE
4791	ZUP = ZSOIL(K-1)
4792	ENDIF
4793
4794	! ----------------------------------------------------------------------
4795	! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4796	! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4797	! ----------------------------------------------------------------------
4798	IF (K .EQ. NSOIL) THEN
4799	ZB = 2.*ZBOT-ZSOIL(K)
4800	ELSE
4801	ZB = ZSOIL(K+1)
4802	ENDIF
4803
4804	! ----------------------------------------------------------------------
4805	! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4806	! ----------------------------------------------------------------------
4807	TBND1 = TU+(TB-TU)*(ZUP-ZSOIL(K))/(ZUP-ZB)
4808
4809	! ----------------------------------------------------------------------
4810	! END SUBROUTINE TBND
4811	! ----------------------------------------------------------------------
4812	END SUBROUTINE TBND
4813
4814	SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
4815
4816	IMPLICIT NONE
4817
4818	! ----------------------------------------------------------------------
4819	! SUBROUTINE TDFCND
4820	! ----------------------------------------------------------------------
4821	! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
4822	! POINT AND TIME.
4823	! ----------------------------------------------------------------------
4824	! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4825	! June 2001 CHANGES: FROZEN SOIL CONDITION.
4826	! ----------------------------------------------------------------------
4827	REAL DF
4828	REAL GAMMD
4829	REAL THKDRY
4830	REAL AKE
4831	REAL THKICE
4832	REAL THKO
4833	REAL THKQTZ
4834	REAL THKSAT
4835	REAL THKS
4836	REAL THKW
4837	REAL QZ
4838	REAL SATRATIO
4839	REAL SH2O
4840	REAL SMC
4841	REAL SMCMAX
4842	REAL XU
4843	REAL XUNFROZ
4844
4845	! ----------------------------------------------------------------------
4846	! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
4847	! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
4848	! & 0.35, 0.60, 0.40, 0.82/
4849	! ----------------------------------------------------------------------
4850	! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4851	! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4852	! ----------------------------------------------------------------------
4853	! THKW ......WATER THERMAL CONDUCTIVITY
4854	! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4855	! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4856	! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4857	! THKICE ....ICE THERMAL CONDUCTIVITY
4858	! SMCMAX ....POROSITY (= SMCMAX)
4859	! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4860	! ----------------------------------------------------------------------
4861	! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4862	!
4863	! PABLO GRUNMANN, 08/17/98
4864	! REFS.:
4865	! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4866	! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4867	! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4868	! UNIVERSITY OF TRONDHEIM,
4869	! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4870	! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4871	! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4872	! VOL. 55, PP. 1209-1224.
4873	! ----------------------------------------------------------------------
4874	! NEEDS PARAMETERS
4875	! POROSITY(SOIL TYPE):
4876	! POROS = SMCMAX
4877	! SATURATION RATIO:
4878	SATRATIO = SMC/SMCMAX
4879
4880	! PARAMETERS W/(M.K)
4881	THKICE = 2.2
4882	THKW = 0.57
4883	THKO = 2.0
4884	! IF (QZ .LE. 0.2) THKO = 3.0
4885	THKQTZ = 7.7
4886	! SOLIDS' CONDUCTIVITY
4887	THKS = (THKQTZ*QZ)(THKO**(1.- QZ))
4888
4889	! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
4890	XUNFROZ = (SH2O + 1.E-9) / (SMC + 1.E-9)
4891
4892	! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
4893	XU=XUNFROZ*SMCMAX
4894	! SATURATED THERMAL CONDUCTIVITY
4895	THKSAT = THKS*(1.-SMCMAX)THKICE*(SMCMAX-XU)THKW**(XU)
4896
4897	! DRY DENSITY IN KG/M3
4898	GAMMD = (1. - SMCMAX)*2700.
4899
4900	! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4901	THKDRY = (0.135GAMMD + 64.7)/(2700. - 0.947GAMMD)
4902
4903	IF ( (SH2O + 0.0005) .LT. SMC ) THEN
4904	! FROZEN
4905	AKE = SATRATIO
4906	ELSE
4907	! UNFROZEN
4908	! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
4909	IF ( SATRATIO .GT. 0.1 ) THEN
4910
4911	! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4912	! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4913	! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4914
4915	AKE = LOG10(SATRATIO) + 1.0
4916
4917	ELSE
4918
4919	! USE K = KDRY
4920	AKE = 0.0
4921
4922	ENDIF
4923	ENDIF
4924
4925	! THERMAL CONDUCTIVITY
4926
4927	DF = AKE*(THKSAT - THKDRY) + THKDRY
4928
4929	! ----------------------------------------------------------------------
4930	! END SUBROUTINE TDFCND
4931	! ----------------------------------------------------------------------
4932	END SUBROUTINE TDFCND
4933
4934	SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
4935
4936	IMPLICIT NONE
4937
4938	! ----------------------------------------------------------------------
4939	! SUBROUTINE TMPAVG
4940	! ----------------------------------------------------------------------
4941	! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
4942	! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
4943	! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
4944	! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
4945	! ----------------------------------------------------------------------
4946	INTEGER K
4947	INTEGER NSOIL
4948
4949	REAL DZ
4950	REAL DZH
4951	REAL T0
4952	REAL TAVG
4953	REAL TDN
4954	REAL TM
4955	REAL TUP
4956	REAL X0
4957	REAL XDN
4958	REAL XUP
4959	REAL ZSOIL (NSOIL)
4960
4961	PARAMETER(T0 = 2.7315E2)
4962
4963	! ----------------------------------------------------------------------
4964	IF (K .EQ. 1) THEN
4965	DZ = -ZSOIL(1)
4966	ELSE
4967	DZ = ZSOIL(K-1)-ZSOIL(K)
4968	ENDIF
4969
4970	DZH=DZ*0.5
4971
4972	IF (TUP .LT. T0) THEN
4973	IF (TM .LT. T0) THEN
4974	IF (TDN .LT. T0) THEN
4975	! ----------------------------------------------------------------------
4976	! TUP, TM, TDN < T0
4977	! ----------------------------------------------------------------------
4978	TAVG = (TUP + 2.0*TM + TDN)/ 4.0
4979	ELSE
4980	! ----------------------------------------------------------------------
4981	! TUP & TM < T0, TDN >= T0
4982	! ----------------------------------------------------------------------
4983	X0 = (T0 - TM) * DZH / (TDN - TM)
4984	TAVG = 0.5 * (TUPDZH+TM(DZH+X0)+T0(2.DZH-X0)) / DZ
4985	ENDIF
4986	ELSE
4987	IF (TDN .LT. T0) THEN
4988	! ----------------------------------------------------------------------
4989	! TUP < T0, TM >= T0, TDN < T0
4990	! ----------------------------------------------------------------------
4991	XUP = (T0-TUP) * DZH / (TM-TUP)
4992	XDN = DZH - (T0-TM) * DZH / (TDN-TM)
4993	TAVG = 0.5 * (TUPXUP+T0(2.DZ-XUP-XDN)+TDNXDN) / DZ
4994	ELSE
4995	! ----------------------------------------------------------------------
4996	! TUP < T0, TM >= T0, TDN >= T0
4997	! ----------------------------------------------------------------------
4998	XUP = (T0-TUP) * DZH / (TM-TUP)
4999	TAVG = 0.5 * (TUPXUP+T0(2.*DZ-XUP)) / DZ
5000	ENDIF
5001	ENDIF
5002	ELSE
5003	IF (TM .LT. T0) THEN
5004	IF (TDN .LT. T0) THEN
5005	! ----------------------------------------------------------------------
5006	! TUP >= T0, TM < T0, TDN < T0
5007	! ----------------------------------------------------------------------
5008	XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5009	TAVG = 0.5 * (T0(DZ-XUP)+TM(DZH+XUP)+TDN*DZH) / DZ
5010	ELSE
5011	! ----------------------------------------------------------------------
5012	! TUP >= T0, TM < T0, TDN >= T0
5013	! ----------------------------------------------------------------------
5014	XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5015	XDN = (T0-TM) * DZH / (TDN-TM)
5016	TAVG = 0.5 * (T0(2.DZ-XUP-XDN)+TM*(XUP+XDN)) / DZ
5017	ENDIF
5018	ELSE
5019	IF (TDN .LT. T0) THEN
5020	! ----------------------------------------------------------------------
5021	! TUP >= T0, TM >= T0, TDN < T0
5022	! ----------------------------------------------------------------------
5023	XDN = DZH - (T0-TM) * DZH / (TDN-TM)
5024	TAVG = (T0(DZ-XDN)+0.5(T0+TDN)*XDN) / DZ
5025	ELSE
5026	! ----------------------------------------------------------------------
5027	! TUP >= T0, TM >= T0, TDN >= T0
5028	! ----------------------------------------------------------------------
5029	TAVG = (TUP + 2.0*TM + TDN) / 4.0
5030	ENDIF
5031	ENDIF
5032	ENDIF
5033	! ----------------------------------------------------------------------
5034	! END SUBROUTINE TMPAVG
5035	! ----------------------------------------------------------------------
5036	END SUBROUTINE TMPAVG
5037
5038	SUBROUTINE TRANSP (ET1,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
5039	& CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
5040
5041	IMPLICIT NONE
5042
5043	! ----------------------------------------------------------------------
5044	! SUBROUTINE TRANSP
5045	! ----------------------------------------------------------------------
5046	! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
5047	! ----------------------------------------------------------------------
5048	INTEGER I
5049	INTEGER K
5050	INTEGER NSOIL
5051	INTEGER NROOT
5052
5053	REAL CFACTR
5054	REAL CMC
5055	REAL CMCMAX
5056	REAL DENOM
5057	REAL ET1(NSOIL)
5058	REAL ETP1
5059	REAL ETP1A
5060	REAL GX (7)
5061	!.....REAL PART(NSOIL)
5062	REAL PC
5063	REAL Q2
5064	REAL RTDIS(NSOIL)
5065	REAL RTX
5066	REAL SFCTMP
5067	REAL SGX
5068	REAL SHDFAC
5069	REAL SMC(NSOIL)
5070	REAL SMCREF
5071	REAL SMCWLT
5072	REAL ZSOIL(NSOIL)
5073
5074	! ----------------------------------------------------------------------
5075	! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
5076	! ----------------------------------------------------------------------
5077	DO K = 1,NSOIL
5078	ET1(K) = 0.
5079	END DO
5080
5081	! ----------------------------------------------------------------------
5082	! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
5083	! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5084	! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
5085	! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
5086	! TOTAL ETP1A.
5087	! ----------------------------------------------------------------------
5088	IF (CMC .NE. 0.0) THEN
5089	ETP1A = SHDFAC * PC * ETP1 * (1.0 - (CMC /CMCMAX) ** CFACTR)
5090	ELSE
5091	ETP1A = SHDFAC * PC * ETP1
5092	ENDIF
5093
5094	SGX = 0.0
5095	DO I = 1,NROOT
5096	GX(I) = ( SMC(I) - SMCWLT ) / ( SMCREF - SMCWLT )
5097	GX(I) = MAX ( MIN ( GX(I), 1. ), 0. )
5098	SGX = SGX + GX (I)
5099	END DO
5100	SGX = SGX / NROOT
5101
5102	DENOM = 0.
5103	DO I = 1,NROOT
5104	RTX = RTDIS(I) + GX(I) - SGX
5105	GX(I) = GX(I) * MAX ( RTX, 0. )
5106	DENOM = DENOM + GX(I)
5107	END DO
5108	IF (DENOM .LE. 0.0) DENOM = 1.
5109
5110	DO I = 1,NROOT
5111	ET1(I) = ETP1A * GX(I) / DENOM
5112	END DO
5113
5114	! ----------------------------------------------------------------------
5115	! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
5116	! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
5117	! ----------------------------------------------------------------------
5118	! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
5119	! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
5120	! ----------------------------------------------------------------------
5121	! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5122	! ----------------------------------------------------------------------
5123	! ET(1) = RTDIS(1) * ETP1A
5124	! ET(1) = ETP1A * PART(1)
5125	! ----------------------------------------------------------------------
5126	! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
5127	! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
5128	! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
5129	! ----------------------------------------------------------------------
5130	! DO K = 2,NROOT
5131	! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
5132	! GX = MAX ( MIN ( GX, 1. ), 0. )
5133	! TEST CANOPY RESISTANCE
5134	! GX = 1.0
5135	! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))GXETP1A
5136	! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
5137	! ----------------------------------------------------------------------
5138	! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5139	! ----------------------------------------------------------------------
5140	! ET(K) = RTDIS(K) * ETP1A
5141	! ET(K) = ETP1A*PART(K)
5142	! END DO
5143	! ----------------------------------------------------------------------
5144	! END SUBROUTINE TRANSP
5145	! ----------------------------------------------------------------------
5146	END SUBROUTINE TRANSP
5147
5148	SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
5149	& SICEMAX)
5150
5151	IMPLICIT NONE
5152
5153	! ----------------------------------------------------------------------
5154	! SUBROUTINE WDFCND
5155	! ----------------------------------------------------------------------
5156	! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
5157	! ----------------------------------------------------------------------
5158	REAL BEXP
5159	REAL DKSAT
5160	REAL DWSAT
5161	REAL EXPON
5162	REAL FACTR1
5163	REAL FACTR2
5164	REAL SICEMAX
5165	REAL SMC
5166	REAL SMCMAX
5167	REAL VKwgt
5168	REAL WCND
5169	REAL WDF
5170
5171	! ----------------------------------------------------------------------
5172	! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
5173	! ----------------------------------------------------------------------
5174	SMC = SMC
5175	SMCMAX = SMCMAX
5176	FACTR1 = 0.2 / SMCMAX
5177	FACTR2 = SMC / SMCMAX
5178
5179	! ----------------------------------------------------------------------
5180	! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
5181	! ----------------------------------------------------------------------
5182	EXPON = BEXP + 2.0
5183	WDF = DWSAT * FACTR2 ** EXPON
5184
5185	! ----------------------------------------------------------------------
5186	! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
5187	! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
5188	! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
5189	! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
5190	! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
5191	! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
5192	! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
5193	! --
5194	! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
5195	! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
5196	! ----------------------------------------------------------------------
5197	IF (SICEMAX .GT. 0.0) THEN
5198	VKWGT = 1./(1.+(500.SICEMAX)*3.)
5199	WDF = VKWGTWDF + (1.- VKWGT)DWSATFACTR1*EXPON
5200	ENDIF
5201
5202	! ----------------------------------------------------------------------
5203	! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
5204	! ----------------------------------------------------------------------
5205	EXPON = (2.0 * BEXP) + 3.0
5206	WCND = DKSAT * FACTR2 ** EXPON
5207
5208	! ----------------------------------------------------------------------
5209	! END SUBROUTINE WDFCND
5210	! ----------------------------------------------------------------------
5211	END SUBROUTINE WDFCND
5212
5213	SUBROUTINE nmmlsminit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
5214	SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
5215	ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, & ! STEMP
5216	TMN, &
5217	num_soil_layers, &
5218	allowed_to_read, &
5219	ids,ide, jds,jde, kds,kde, &
5220	ims,ime, jms,jme, kms,kme, &
5221	its,ite, jts,jte, kts,kte )
5222
5223	IMPLICIT NONE
5224
5225	! Arguments
5226	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
5227	ims,ime, jms,jme, kms,kme, &
5228	its,ite, jts,jte, kts,kte
5229
5230	INTEGER, INTENT(IN) :: num_soil_layers
5231
5232	REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
5233
5234	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
5235	INTENT(INOUT) :: SMOIS, &
5236	TSLB !STEMP
5237
5238	REAL, DIMENSION( ims:ime, jms:jme ) , &
5239	INTENT(INOUT) :: SNOW, &
5240	SNOWC, &
5241	CANWAT, &
5242	SMSTAV, &
5243	SMSTOT, &
5244	SFCRUNOFF, &
5245	UDRUNOFF, &
5246	SFCEVP, &
5247	GRDFLX, &
5248	ACSNOW, &
5249	XICE, &
5250	VEGFRA, &
5251	TMN, &
5252	ACSNOM
5253
5254	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
5255	INTENT(INOUT) :: IVGTYP, &
5256	ISLTYP
5257
5258	!
5259
5260	INTEGER, INTENT(IN) :: isn
5261	LOGICAL, INTENT(IN) :: allowed_to_read
5262	! Local
5263	INTEGER :: iseason
5264	INTEGER :: icm,jcm,itf,jtf
5265	INTEGER :: I,J,L
5266
5267
5268	itf=min0(ite,ide-1)
5269	jtf=min0(jte,jde-1)
5270
5271	icm = ide/2
5272	jcm = jde/2
5273
5274	iseason=isn
5275
5276	DO J=jts,jtf
5277	DO I=its,itf
5278	! SNOW(i,j)=0.
5279	SNOWC(i,j)=0.
5280	! SMSTAV(i,j)=
5281	! SMSTOT(i,j)=
5282	! SFCRUNOFF(i,j)=
5283	! UDRUNOFF(i,j)=
5284	! GRDFLX(i,j)=
5285	! ACSNOW(i,j)=
5286	! ACSNOM(i,j)=
5287	ENDDO
5288	ENDDO
5289
5290	END SUBROUTINE nmmlsminit
5291
5292	FUNCTION CSNOW (DSNOW)
5293
5294	IMPLICIT NONE
5295
5296	! ----------------------------------------------------------------------
5297	! FUNCTION CSNOW
5298	! ----------------------------------------------------------------------
5299	! CALCULATE SNOW TERMAL CONDUCTIVITY
5300	! ----------------------------------------------------------------------
5301	REAL C
5302	REAL DSNOW
5303	REAL CSNOW
5304	REAL UNIT
5305
5306	PARAMETER(UNIT = 0.11631)
5307
5308	! ----------------------------------------------------------------------
5309	! CSNOW IN UNITS OF CAL/(CMHRC), RETURNED IN W/(M*C)
5310	! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
5311	! ----------------------------------------------------------------------
5312	C=0.32810(2.25DSNOW)
5313	! CSNOW=UNIT*C
5314	! MEK JAN 2006, DOUBLE SNOW THERMAL CONDUCTIVITY
5315	CSNOW=2.0UNITC
5316
5317	! ----------------------------------------------------------------------
5318	! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
5319	! ----------------------------------------------------------------------
5320	! CSNOW=0.0293(1.+100.DSNOW**2)
5321
5322	! ----------------------------------------------------------------------
5323	! E. ANDERSEN FROM FLERCHINGER
5324	! ----------------------------------------------------------------------
5325	! CSNOW=0.021+2.51DSNOW*2
5326
5327	! ----------------------------------------------------------------------
5328	! END FUNCTION CSNOW
5329	! ----------------------------------------------------------------------
5330	END FUNCTION CSNOW
5331
5332	FUNCTION FRH2O (TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
5333
5334	IMPLICIT NONE
5335
5336	! ----------------------------------------------------------------------
5337	! FUNCTION FRH2O
5338	! ----------------------------------------------------------------------
5339	! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
5340	! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
5341	! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
5342	! (1999, JGR, VOL 104(D16), 19569-19585).
5343	! ----------------------------------------------------------------------
5344	! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
5345	! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
5346	! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
5347	! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
5348	! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
5349	! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
5350	! LIMIT OF FREEZING POINT TEMPERATURE T0.
5351	! ----------------------------------------------------------------------
5352	! INPUT:
5353	!
5354	! TKELV.........TEMPERATURE (Kelvin)
5355	! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
5356	! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
5357	! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
5358	! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
5359	! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
5360	!
5361	! OUTPUT:
5362	! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
5363	! ----------------------------------------------------------------------
5364	REAL BEXP
5365	REAL BLIM
5366	REAL BX
5367	REAL CK
5368	REAL DENOM
5369	REAL DF
5370	REAL DH2O
5371	REAL DICE
5372	REAL DSWL
5373	REAL ERROR
5374	REAL FK
5375	REAL FRH2O
5376	REAL GS
5377	REAL HLICE
5378	REAL PSIS
5379	REAL SH2O
5380	REAL SMC
5381	REAL SMCMAX
5382	REAL SWL
5383	REAL SWLK
5384	REAL TKELV
5385	REAL T0
5386
5387	INTEGER NLOG
5388	INTEGER KCOUNT
5389
5390	PARAMETER(CK = 8.0)
5391	! PARAMETER(CK = 0.0)
5392	PARAMETER(BLIM = 5.5)
5393	PARAMETER(ERROR = 0.005)
5394
5395	PARAMETER(HLICE = 3.335E5)
5396	PARAMETER(GS = 9.81)
5397	PARAMETER(DICE = 920.0)
5398	PARAMETER(DH2O = 1000.0)
5399	PARAMETER(T0 = 273.15)
5400
5401	! ----------------------------------------------------------------------
5402	! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
5403	! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
5404	! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
5405	! ----------------------------------------------------------------------
5406	BX = BEXP
5407	IF (BEXP .GT. BLIM) BX = BLIM
5408
5409	! ----------------------------------------------------------------------
5410	! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
5411	! ----------------------------------------------------------------------
5412	NLOG=0
5413	KCOUNT=0
5414
5415	! ----------------------------------------------------------------------
5416	! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
5417	! ----------------------------------------------------------------------
5418	IF (TKELV .GT. (T0 - 1.E-3)) THEN
5419	FRH2O = SMC
5420	ELSE
5421	IF (CK .NE. 0.0) THEN
5422
5423	! ----------------------------------------------------------------------
5424	! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
5425	! IN KOREN ET AL, JGR, 1999, EQN 17
5426	! ----------------------------------------------------------------------
5427	! INITIAL GUESS FOR SWL (frozen content)
5428	! ----------------------------------------------------------------------
5429	SWL = SMC-SH2O
5430
5431	! ----------------------------------------------------------------------
5432	! KEEP WITHIN BOUNDS.
5433	! ----------------------------------------------------------------------
5434	IF (SWL .GT. (SMC-0.02)) SWL = SMC-0.02
5435	IF (SWL .LT. 0.) SWL = 0.
5436
5437	! ----------------------------------------------------------------------
5438	! START OF ITERATIONS
5439	! ----------------------------------------------------------------------
5440	DO WHILE ( (NLOG .LT. 10) .AND. (KCOUNT .EQ. 0) )
5441	NLOG = NLOG+1
5442	DF = ALOG(( PSISGS/HLICE ) ( ( 1.+CKSWL )2. ) &
5443	& ( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
5444	DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
5445	SWLK = SWL - DF/DENOM
5446	! ----------------------------------------------------------------------
5447	! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
5448	! ----------------------------------------------------------------------
5449	IF (SWLK .GT. (SMC-0.02)) SWLK = SMC - 0.02
5450	IF (SWLK .LT. 0.) SWLK = 0.
5451
5452	! ----------------------------------------------------------------------
5453	! MATHEMATICAL SOLUTION BOUNDS APPLIED.
5454	! ----------------------------------------------------------------------
5455	DSWL = ABS(SWLK-SWL)
5456	SWL = SWLK
5457
5458	! ----------------------------------------------------------------------
5459	! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
5460	! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
5461	! ----------------------------------------------------------------------
5462	IF ( DSWL .LE. ERROR ) THEN
5463	KCOUNT = KCOUNT+1
5464	ENDIF
5465	END DO
5466
5467	! ----------------------------------------------------------------------
5468	! END OF ITERATIONS
5469	! ----------------------------------------------------------------------
5470	! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
5471	! ----------------------------------------------------------------------
5472	FRH2O = SMC - SWL
5473
5474	! ----------------------------------------------------------------------
5475	! END OPTION 1
5476	! ----------------------------------------------------------------------
5477	ENDIF
5478
5479	! ----------------------------------------------------------------------
5480	! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
5481	! IN KOREN ET AL., JGR, 1999, EQN 17
5482	! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
5483	! ----------------------------------------------------------------------
5484	IF (KCOUNT .EQ. 0) THEN
5485	Print*,'Flerchinger used in NEW version. Iterations=',NLOG
5486	FK = (((HLICE/(GS(-PSIS))) &
5487	& ((TKELV-T0)/TKELV))*(-1/BX))SMCMAX
5488	IF (FK .LT. 0.02) FK = 0.02
5489	FRH2O = MIN (FK, SMC)
5490	! ----------------------------------------------------------------------
5491	! END OPTION 2
5492	! ----------------------------------------------------------------------
5493	ENDIF
5494
5495	ENDIF
5496
5497	! ----------------------------------------------------------------------
5498	! END FUNCTION FRH2O
5499	! ----------------------------------------------------------------------
5500	END FUNCTION FRH2O
5501
5502	FUNCTION SNKSRC (TAVG,SMC,SH2O,ZSOIL,NSOIL, &
5503	& SMCMAX,PSISAT,BEXP,DT,K,QTOT)
5504
5505	IMPLICIT NONE
5506
5507	! ----------------------------------------------------------------------
5508	! FUNCTION SNKSRC
5509	! ----------------------------------------------------------------------
5510	! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
5511	! AVAILABLE LIQUED WATER.
5512	! ----------------------------------------------------------------------
5513	INTEGER K
5514	INTEGER NSOIL
5515
5516	REAL BEXP
5517	REAL DF
5518	REAL DH2O
5519	REAL DT
5520	REAL DZ
5521	REAL DZH
5522	REAL FREE
5523	! REAL FRH2O
5524	REAL HLICE
5525	REAL PSISAT
5526	REAL QTOT
5527	REAL SH2O
5528	REAL SMC
5529	REAL SMCMAX
5530	REAL SNKSRC
5531	REAL T0
5532	REAL TAVG
5533	REAL TDN
5534	REAL TM
5535	REAL TUP
5536	REAL TZ
5537	REAL X0
5538	REAL XDN
5539	REAL XH2O
5540	REAL XUP
5541	REAL ZSOIL (NSOIL)
5542
5543	PARAMETER(DH2O = 1.0000E3)
5544	PARAMETER(HLICE = 3.3350E5)
5545	PARAMETER(T0 = 2.7315E2)
5546
5547	IF (K .EQ. 1) THEN
5548	DZ = -ZSOIL(1)
5549	ELSE
5550	DZ = ZSOIL(K-1)-ZSOIL(K)
5551	ENDIF
5552
5553	! ----------------------------------------------------------------------
5554	! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
5555	! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
5556	! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
5557	! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
5558	! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
5559	! ----------------------------------------------------------------------
5560	FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
5561
5562	! ----------------------------------------------------------------------
5563	! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
5564	! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
5565	! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
5566	! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
5567	! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
5568	! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
5569	! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
5570	! ----------------------------------------------------------------------
5571	XH2O = SH2O + QTOTDT/(DH2OHLICE*DZ)
5572
5573	! ----------------------------------------------------------------------
5574	! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
5575	! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
5576	! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
5577	! ----------------------------------------------------------------------
5578	IF ( XH2O .LT. SH2O .AND. XH2O .LT. FREE) THEN
5579	IF ( FREE .GT. SH2O ) THEN
5580	XH2O = SH2O
5581	ELSE
5582	XH2O = FREE
5583	ENDIF
5584	ENDIF
5585
5586	! ----------------------------------------------------------------------
5587	! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
5588	! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
5589	! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
5590	! ----------------------------------------------------------------------
5591	IF ( XH2O .GT. SH2O .AND. XH2O .GT. FREE ) THEN
5592	IF ( FREE .LT. SH2O ) THEN
5593	XH2O = SH2O
5594	ELSE
5595	XH2O = FREE
5596	ENDIF
5597	ENDIF
5598
5599	IF (XH2O .LT. 0.) XH2O = 0.
5600	IF (XH2O .GT. SMC) XH2O = SMC
5601
5602	! ----------------------------------------------------------------------
5603	! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
5604	! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
5605	! ----------------------------------------------------------------------
5606	SNKSRC = -DH2OHLICEDZ*(XH2O-SH2O)/DT
5607	SH2O = XH2O
5608
5609	! ----------------------------------------------------------------------
5610	! END FUNCTION SNKSRC
5611	! ----------------------------------------------------------------------
5612	END FUNCTION SNKSRC
5613
5614	END MODULE module_sf_lsm_nmm
5615

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

Download in other formats:

Original Format