Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_noahlsm.F @ 1544

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 193.6 KB

Line
1	MODULE module_sf_noahlsm
2	USE module_model_constants
3
4	! REAL, PARAMETER :: CP = 1004.5
5	REAL, PARAMETER :: RD = 287.04, SIGMA = 5.67E-8, &
6	CPH2O = 4.218E+3,CPICE = 2.106E+3, &
7	LSUBF = 3.335E+5, &
8	EMISSI_S = 0.95
9
10	! VEGETATION PARAMETERS
11	INTEGER :: LUCATS , BARE
12	INTEGER :: NATURAL
13	integer, PARAMETER :: NLUS=50
14	CHARACTER(LEN=256) LUTYPE
15	INTEGER, DIMENSION(1:NLUS) :: NROTBL
16	real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, &
17	SHDTBL, MAXALB, &
18	EMISSMINTBL, EMISSMAXTBL, &
19	LAIMINTBL, LAIMAXTBL, &
20	Z0MINTBL, Z0MAXTBL, &
21	ALBEDOMINTBL, ALBEDOMAXTBL
22	REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
23
24	! SOIL PARAMETERS
25	INTEGER :: SLCATS
26	INTEGER, PARAMETER :: NSLTYPE=30
27	CHARACTER(LEN=256) SLTYPE
28	REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,F11, &
29	MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
30
31	! LSM GENERAL PARAMETERS
32	INTEGER :: SLPCATS
33	INTEGER, PARAMETER :: NSLOPE=30
34	REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
35	REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
36	REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
37	CZIL_DATA
38	REAL :: LVCOEF_DATA
39
40	CHARACTER*256 :: err_message
41
42	!
43	CONTAINS
44	!
45
46	SUBROUTINE SFLX (FFROZP,ICE,ISURBAN,DT,ZLVL,NSOIL,SLDPTH, & !C
47	LOCAL, & !L
48	LLANDUSE, LSOIL, & !CL
49	LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2,SFCSPD, & !F
50	COSZ,PRCPRAIN, SOLARDIRECT, & !F
51	TH2,Q2SAT,DQSDT2, & !I
52	VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHDMIN,SHDMAX, & !I
53	ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, EMBRD, & !S
54	CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM, & !H
55	! ----------------------------------------------------------------------
56	! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
57	! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
58	! MODEL). OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
59	! ----------------------------------------------------------------------
60	ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
61	EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
62	BETA,ETP,SSOIL, & !O
63	FLX1,FLX2,FLX3, & !O
64	SNOMLT,SNCOVR, & !O
65	RUNOFF1,RUNOFF2,RUNOFF3, & !O
66	RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
67	SOILW,SOILM,Q1,SMAV, & !D
68	RDLAI2D,USEMONALB, &
69	SNOTIME1, &
70	RIBB, &
71	SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT) !P
72	! ----------------------------------------------------------------------
73	! SUBROUTINE SFLX - UNIFIED NOAHLSM VERSION 1.0 JULY 2007
74	! ----------------------------------------------------------------------
75	! SUB-DRIVER FOR "Noah LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
76	! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
77	! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
78	! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
79	! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
80	! RADIATION AND PRECIP)
81	! ----------------------------------------------------------------------
82	! SFLX ARGUMENT LIST KEY:
83	! ----------------------------------------------------------------------
84	! C CONFIGURATION INFORMATION
85	! L LOGICAL
86	! CL 4-string character bearing logical meaning
87	! F FORCING DATA
88	! I OTHER (INPUT) FORCING DATA
89	! S SURFACE CHARACTERISTICS
90	! H HISTORY (STATE) VARIABLES
91	! O OUTPUT VARIABLES
92	! D DIAGNOSTIC OUTPUT
93	! P Parameters
94	! Msic Miscellaneous terms passed from gridded driver
95	! ----------------------------------------------------------------------
96	! 1. CONFIGURATION INFORMATION (C):
97	! ----------------------------------------------------------------------
98	! ICE SEA-ICE FLAG (=1: SEA-ICE, =0: LAND (NO ICE), --1 LAND-ICE).
99	! DT TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
100	! 1800 SECS OR LESS)
101	! ZLVL HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
102	! NSOIL NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
103	! PARAMETER NSOLD SET BELOW)
104	! SLDPTH THE THICKNESS OF EACH SOIL LAYER (M)
105	! ----------------------------------------------------------------------
106	! 2. LOGICAL:
107	! ----------------------------------------------------------------------
108	! LCH Exchange coefficient (Ch) calculation flag (false: using
109	! ch-routine SFCDIF; true: Ch is brought in)
110	! LOCAL Flag for local-site simulation (where there is no
111	! maps for albedo, veg fraction, and roughness
112	! true: all LSM parameters (inluding albedo, veg fraction and
113	! roughness length) will be defined by three tables
114	! LLANDUSE (=USGS, using USGS landuse classification)
115	! LSOIL (=STAS, using FAO/STATSGO soil texture classification)
116	! ----------------------------------------------------------------------
117	! 3. FORCING DATA (F):
118	! ----------------------------------------------------------------------
119	! LWDN LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
120	! SOLDN SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
121	! SOLNET NET DOWNWARD SOLAR RADIATION ((W M-2; POSITIVE)
122	! SFCPRS PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
123	! PRCP PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
124	! SFCTMP AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
125	! TH2 AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
126	! Q2 MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
127	! COSZ Solar zenith angle (not used for now)
128	! PRCPRAIN Liquid-precipitation rate (KG M-2 S-1) (not used)
129	! SOLARDIRECT Direct component of downward solar radiation (W M-2) (not used)
130	! FFROZP FRACTION OF FROZEN PRECIPITATION
131	! ----------------------------------------------------------------------
132	! 4. OTHER FORCING (INPUT) DATA (I):
133	! ----------------------------------------------------------------------
134	! SFCSPD WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
135	! Q2SAT SAT SPECIFIC HUMIDITY AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
136	! DQSDT2 SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
137	! (KG KG-1 K-1)
138	! ----------------------------------------------------------------------
139	! 5. CANOPY/SOIL CHARACTERISTICS (S):
140	! ----------------------------------------------------------------------
141	! VEGTYP VEGETATION TYPE (INTEGER INDEX)
142	! SOILTYP SOIL TYPE (INTEGER INDEX)
143	! SLOPETYP CLASS OF SFC SLOPE (INTEGER INDEX)
144	! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
145	! (FRACTION= 0.0-1.0)
146	! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
147	! (FRACTION= 0.0-1.0) <= SHDFAC
148	! PTU PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
149	! (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
150	! VEG PARMS)
151	! ALB BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
152	! DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
153	! MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
154	! INCLUDE DIURNAL SUN ANGLE EFFECT)
155	! SNOALB UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
156	! ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
157	! TBOT BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
158	! TEMPERATURE)
159	! Z0BRD Background fixed roughness length (M)
160	! Z0 Time varying roughness length (M) as function of snow depth
161	!
162	! EMBRD Background surface emissivity (between 0 and 1)
163	! EMISSI Surface emissivity (between 0 and 1)
164	! ----------------------------------------------------------------------
165	! 6. HISTORY (STATE) VARIABLES (H):
166	! ----------------------------------------------------------------------
167	! CMC CANOPY MOISTURE CONTENT (M)
168	! T1 GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
169	! STC(NSOIL) SOIL TEMP (K)
170	! SMC(NSOIL) TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
171	! SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
172	! NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
173	! SNOWH ACTUAL SNOW DEPTH (M)
174	! SNEQV LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
175	! NOTE: SNOW DENSITY = SNEQV/SNOWH
176	! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
177	! =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
178	! =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
179	! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
180	! (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
181	! IT HAS BEEN MULTIPLIED BY WIND SPEED.
182	! CM SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
183	! CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
184	! MULTIPLIED BY WIND SPEED.
185	! ----------------------------------------------------------------------
186	! 7. OUTPUT (O):
187	! ----------------------------------------------------------------------
188	! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
189	! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL. FOR THIS APPLICATION,
190	! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
191	! NECESSARY. OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
192	! ETA ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
193	! SURFACE)
194	! ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
195	! SHEAT SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
196	! SURFACE)
197	! FDOWN Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
198	! ----------------------------------------------------------------------
199	! EC CANOPY WATER EVAPORATION (W m-2)
200	! EDIR DIRECT SOIL EVAPORATION (W m-2)
201	! ET(NSOIL) PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
202	! (W m-2)
203	! ETT TOTAL PLANT TRANSPIRATION (W m-2)
204	! ESNOW SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
205	! (W m-2)
206	! DRIP THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
207	! WATER-HOLDING CAPACITY (M)
208	! DEW DEWFALL (OR FROSTFALL FOR T<273.15) (M)
209	! ----------------------------------------------------------------------
210	! BETA RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
211	! ETP POTENTIAL EVAPORATION (W m-2)
212	! SSOIL SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
213	! ----------------------------------------------------------------------
214	! FLX1 PRECIP-SNOW SFC (W M-2)
215	! FLX2 FREEZING RAIN LATENT HEAT FLUX (W M-2)
216	! FLX3 PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
217	! ----------------------------------------------------------------------
218	! SNOMLT SNOW MELT (M) (WATER EQUIVALENT)
219	! SNCOVR FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
220	! ----------------------------------------------------------------------
221	! RUNOFF1 SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
222	! RUNOFF2 SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
223	! SOIL LAYER (BASEFLOW)
224	! RUNOFF3 NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
225	! FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP (M S-1).
226	! Note: the above RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
227	! ----------------------------------------------------------------------
228	! RC CANOPY RESISTANCE (S M-1)
229	! PC PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
230	! = ACTUAL TRANSP
231	! XLAI LEAF AREA INDEX (DIMENSIONLESS)
232	! RSMIN MINIMUM CANOPY RESISTANCE (S M-1)
233	! RCS INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
234	! RCT AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
235	! RCQ ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
236	! RCSOIL SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
237	! ----------------------------------------------------------------------
238	! 8. DIAGNOSTIC OUTPUT (D):
239	! ----------------------------------------------------------------------
240	! SOILW AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
241	! BETWEEN SMCWLT AND SMCMAX)
242	! SOILM TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
243	! Q1 Effective mixing ratio at surface (kg kg-1), used for
244	! diagnosing the mixing ratio at 2 meter for coupled model
245	! SMAV Soil Moisture Availability for each layer, as a fraction
246	! between SMCWLT and SMCMAX.
247	! Documentation for SNOTIME1 and SNOABL2 ?????
248	! What categories of arguments do these variables fall into ????
249	! Documentation for RIBB ?????
250	! What category of argument does RIBB fall into ?????
251	! ----------------------------------------------------------------------
252	! 9. PARAMETERS (P):
253	! ----------------------------------------------------------------------
254	! SMCWLT WILTING POINT (VOLUMETRIC)
255	! SMCDRY DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
256	! LAYER ENDS (VOLUMETRIC)
257	! SMCREF SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
258	! STRESS (VOLUMETRIC)
259	! SMCMAX POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
260	! (VOLUMETRIC)
261	! NROOT NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
262	! IN SUBROUTINE REDPRM.
263	! ----------------------------------------------------------------------
264
265
266	IMPLICIT NONE
267	! ----------------------------------------------------------------------
268
269	! DECLARATIONS - LOGICAL AND CHARACTERS
270	! ----------------------------------------------------------------------
271	LOGICAL, INTENT(IN):: LOCAL
272	LOGICAL :: FRZGRA, SNOWNG
273	CHARACTER (LEN=256), INTENT(IN):: LLANDUSE, LSOIL
274
275	! ----------------------------------------------------------------------
276	! 1. CONFIGURATION INFORMATION (C):
277	! ----------------------------------------------------------------------
278	INTEGER,INTENT(IN) :: ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
279	INTEGER, INTENT(IN) :: ISURBAN
280	INTEGER,INTENT(OUT):: NROOT
281	INTEGER KZ, K, iout
282
283	! ----------------------------------------------------------------------
284	! 2. LOGICAL:
285	! ----------------------------------------------------------------------
286	LOGICAL, INTENT(IN) :: RDLAI2D
287	LOGICAL, INTENT(IN) :: USEMONALB
288
289	REAL, INTENT(IN) :: SHDMIN,SHDMAX,DT,DQSDT2,LWDN,PRCP,PRCPRAIN, &
290	Q2,Q2SAT,SFCPRS,SFCSPD,SFCTMP, SNOALB, &
291	SOLDN,SOLNET,TBOT,TH2,ZLVL, &
292	FFROZP
293	REAL, INTENT(OUT) :: EMBRD
294	REAL, INTENT(OUT) :: ALBEDO
295	REAL, INTENT(INOUT):: COSZ, SOLARDIRECT,CH,CM, &
296	CMC,SNEQV,SNCOVR,SNOWH,T1,XLAI,SHDFAC,Z0BRD, &
297	EMISSI, ALB
298	REAL, INTENT(INOUT):: SNOTIME1
299	REAL, INTENT(INOUT):: RIBB
300	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SLDPTH
301	REAL, DIMENSION(1:NSOIL), INTENT(OUT):: ET
302	REAL, DIMENSION(1:NSOIL), INTENT(OUT):: SMAV
303	REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O, SMC, STC
304	REAL,DIMENSION(1:NSOIL):: RTDIS, ZSOIL
305
306	REAL,INTENT(OUT) :: ETA_KINEMATIC,BETA,DEW,DRIP,EC,EDIR,ESNOW,ETA, &
307	ETP,FLX1,FLX2,FLX3,SHEAT,PC,RUNOFF1,RUNOFF2, &
308	RUNOFF3,RC,RSMIN,RCQ,RCS,RCSOIL,RCT,SSOIL, &
309	SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT, SOILM, &
310	SOILW,FDOWN,Q1
311	REAL :: BEXP,CFACTR,CMCMAX,CSOIL,CZIL,DF1,DF1H,DF1A,DKSAT,DWSAT, &
312	DSOIL,DTOT,ETT,FRCSNO,FRCSOI,EPSCA,F1,FXEXP,FRZX,HS, &
313	KDT,LVH2O,PRCP1,PSISAT,QUARTZ,R,RCH,REFKDT,RR,RGL, &
314	RSMAX, &
315	RSNOW,SNDENS,SNCOND,SBETA,SN_NEW,SLOPE,SNUP,SALP,SOILWM, &
316	SOILWW,T1V,T24,T2V,TH2V,TOPT,TFREEZ,TSNOW,ZBOT,Z0,PRCPF, &
317	ETNS,PTU,LSUBS
318	REAL :: LVCOEF
319	REAL :: INTERP_FRACTION
320	REAL :: LAIMIN, LAIMAX
321	REAL :: ALBEDOMIN, ALBEDOMAX
322	REAL :: EMISSMIN, EMISSMAX
323	REAL :: Z0MIN, Z0MAX
324
325	! ----------------------------------------------------------------------
326	! DECLARATIONS - PARAMETERS
327	! ----------------------------------------------------------------------
328	PARAMETER (TFREEZ = 273.15)
329	PARAMETER (LVH2O = 2.501E+6)
330	PARAMETER (LSUBS = 2.83E+6)
331	PARAMETER (R = 287.04)
332	! ----------------------------------------------------------------------
333	! INITIALIZATION
334	! ----------------------------------------------------------------------
335	RUNOFF1 = 0.0
336	RUNOFF2 = 0.0
337	RUNOFF3 = 0.0
338	SNOMLT = 0.0
339
340	! ----------------------------------------------------------------------
341	! THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE / LAND-ICE / ICE-FREE LAND
342	! SEA-ICE CASE, ICE = 1
343	! NON-GLACIAL LAND, ICE = 0
344	! GLACIAL-ICE LAND, ICE = -1
345	IF (ICE /= 0) SHDFAC = 0.0
346	! ----------------------------------------------------------------------
347	! SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
348	! ----------------------------------------------------------------------
349	IF (ICE == 1) THEN
350	DO KZ = 1,NSOIL
351	ZSOIL (KZ) = -3.* FLOAT (KZ)/ FLOAT (NSOIL)
352	END DO
353
354	! ----------------------------------------------------------------------
355	! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
356	! EACH SOIL LAYER. NOTE: SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
357	! GROUND)
358	! ----------------------------------------------------------------------
359	ELSE
360	ZSOIL (1) = - SLDPTH (1)
361	DO KZ = 2,NSOIL
362	ZSOIL (KZ) = - SLDPTH (KZ) + ZSOIL (KZ -1)
363	END DO
364	END IF
365	! ----------------------------------------------------------------------
366	! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
367	! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
368	! ----------------------------------------------------------------------
369	CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,TOPT, &
370	REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
371	PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
372	SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
373	RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL,CZIL, &
374	LAIMIN, LAIMAX, EMISSMIN, EMISSMAX, ALBEDOMIN, &
375	ALBEDOMAX, Z0MIN, Z0MAX, CSOIL, PTU, LLANDUSE, &
376	LSOIL,LOCAL,LVCOEF)
377
378	!urban
379	IF(VEGTYP==ISURBAN)THEN
380	SHDFAC=0.05
381	RSMIN=400.0
382	SMCMAX = 0.45
383	SMCREF = 0.42
384	SMCWLT = 0.40
385	SMCDRY = 0.40
386	ENDIF
387
388	IF ( SHDFAC >= SHDMAX ) THEN
389	EMBRD = EMISSMAX
390	IF (.NOT. RDLAI2D) THEN
391	XLAI = LAIMAX
392	ENDIF
393	IF (.NOT. USEMONALB) THEN
394	ALB = ALBEDOMIN
395	ENDIF
396	Z0BRD = Z0MAX
397	ELSE IF ( SHDFAC <= SHDMIN ) THEN
398	EMBRD = EMISSMIN
399	IF(.NOT. RDLAI2D) THEN
400	XLAI = LAIMIN
401	ENDIF
402	IF(.NOT. USEMONALB) then
403	ALB = ALBEDOMAX
404	ENDIF
405	Z0BRD = Z0MIN
406	ELSE
407
408	IF ( SHDMAX > SHDMIN ) THEN
409
410	INTERP_FRACTION = ( SHDFAC - SHDMIN ) / ( SHDMAX - SHDMIN )
411	! Bound INTERP_FRACTION between 0 and 1
412	INTERP_FRACTION = MIN ( INTERP_FRACTION, 1.0 )
413	INTERP_FRACTION = MAX ( INTERP_FRACTION, 0.0 )
414	! Scale Emissivity and LAI between EMISSMIN and EMISSMAX by INTERP_FRACTION
415	EMBRD = ( ( 1.0 - INTERP_FRACTION ) * EMISSMIN ) + ( INTERP_FRACTION * EMISSMAX )
416	IF (.NOT. RDLAI2D) THEN
417	XLAI = ( ( 1.0 - INTERP_FRACTION ) * LAIMIN ) + ( INTERP_FRACTION * LAIMAX )
418	ENDIF
419	if (.not. USEMONALB) then
420	ALB = ( ( 1.0 - INTERP_FRACTION ) * ALBEDOMAX ) + ( INTERP_FRACTION * ALBEDOMIN )
421	endif
422	Z0BRD = ( ( 1.0 - INTERP_FRACTION ) * Z0MIN ) + ( INTERP_FRACTION * Z0MAX )
423
424	ELSE
425
426	EMBRD = 0.5 * EMISSMIN + 0.5 * EMISSMAX
427	IF (.NOT. RDLAI2D) THEN
428	XLAI = 0.5 * LAIMIN + 0.5 * LAIMAX
429	ENDIF
430	if (.not. USEMONALB) then
431	ALB = 0.5 * ALBEDOMIN + 0.5 * ALBEDOMAX
432	endif
433	Z0BRD = 0.5 * Z0MIN + 0.5 * Z0MAX
434
435	ENDIF
436
437	ENDIF
438	! ----------------------------------------------------------------------
439	! INITIALIZE PRECIPITATION LOGICALS.
440	! ----------------------------------------------------------------------
441	SNOWNG = .FALSE.
442	FRZGRA = .FALSE.
443
444	! ----------------------------------------------------------------------
445	! OVER SEA-ICE OR GLACIAL-ICE, IF S.W.E. (SNEQV) BELOW THRESHOLD LOWER
446	! BOUND (0.01 M FOR SEA-ICE, 0.10 M FOR GLACIAL-ICE), THEN SET AT LOWER
447	! BOUND
448	! ----------------------------------------------------------------------
449	! IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP
450	! ----------------------------------------------------------------------
451	IF (ICE == 1) THEN
452	! Sea-ice case
453	IF ( SNEQV < 0.01 ) THEN
454	SNEQV = 0.01
455	SNOWH = 0.05
456	ENDIF
457	ELSE IF ( ICE == -1 ) THEN
458	! Land-ice case
459	IF ( SNEQV < 0.10 ) THEN
460	SNEQV = 0.10
461	SNOWH = 0.50
462	ENDIF
463	END IF
464	! ----------------------------------------------------------------------
465	! FOR SEA-ICE AND GLACIAL-ICE CASES, SET SMC AND SH20 VALUES = 1.0
466	! ----------------------------------------------------------------------
467	IF ( ICE /= 0 ) THEN
468	DO KZ = 1,NSOIL
469	SMC(KZ) = 1.0
470	SH2O(KZ) = 1.0
471	END DO
472	ENDIF
473	! ----------------------------------------------------------------------
474	! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
475	! SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
476	! SUBROUTINE)
477	! ----------------------------------------------------------------------
478	IF ( SNEQV <= 1.E-7 ) THEN ! safer IF kmh (2008/03/25)
479	SNEQV = 0.0
480	SNDENS = 0.0
481	SNOWH = 0.0
482	SNCOND = 1.0
483	ELSE
484	SNDENS = SNEQV / SNOWH
485	IF(SNDENS > 1.0) THEN
486	CALL wrf_error_fatal ( 'Physical snow depth is less than snow water equiv.' )
487	ENDIF
488	CALL CSNOW (SNCOND,SNDENS)
489	END IF
490	! ----------------------------------------------------------------------
491	! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
492	! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
493	! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
494	! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
495	! ----------------------------------------------------------------------
496	IF (PRCP > 0.0) THEN
497	! snow defined when fraction of frozen precip (FFROZP) > 0.5,
498	! passed in from model microphysics.
499	IF (FFROZP .GT. 0.5) THEN
500	SNOWNG = .TRUE.
501	ELSE
502	IF (T1 <= TFREEZ) FRZGRA = .TRUE.
503	END IF
504	END IF
505	! ----------------------------------------------------------------------
506	! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
507	! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
508	! IT TO THE EXISTING SNOWPACK.
509	! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
510	! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
511	! ----------------------------------------------------------------------
512	IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
513	SN_NEW = PRCP * DT * 0.001
514	SNEQV = SNEQV + SN_NEW
515	PRCPF = 0.0
516
517	! ----------------------------------------------------------------------
518	! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
519	! UPDATE SNOW THERMAL CONDUCTIVITY
520	! ----------------------------------------------------------------------
521	CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
522	!
523	! kmh 09/04/2006 set Snow Density at 0.2 g/cm**3
524	! for "cold permanent ice" or new "dry" snow
525	!
526	IF ( (ICE /= 0) .and. SNCOVR .GT. 0.99 ) THEN
527	! if soil temperature less than 268.15 K, treat as typical Antarctic/Greenland snow firn
528	IF ( STC(1) .LT. (TFREEZ - 5.) ) SNDENS = 0.2
529	IF ( SNOWNG .AND. (T1.LT.273.) .AND. (SFCTMP.LT.273.) ) SNDENS=0.2
530	ENDIF
531	!
532	CALL CSNOW (SNCOND,SNDENS)
533
534	! ----------------------------------------------------------------------
535	! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
536	! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
537	! ANY CANOPY "DRIP" ADDED TO THIS LATER)
538	! ----------------------------------------------------------------------
539	ELSE
540	PRCPF = PRCP
541	END IF
542	! ----------------------------------------------------------------------
543	! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
544	! ----------------------------------------------------------------------
545	! ----------------------------------------------------------------------
546	! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
547	! ----------------------------------------------------------------------
548	IF (ICE == 0 .OR. ICE == -1) THEN
549	IF (SNEQV == 0.0) THEN
550	SNCOVR = 0.0
551	ALBEDO = ALB
552	EMISSI = EMBRD
553	ELSE
554	! ----------------------------------------------------------------------
555	! DETERMINE SNOW FRACTIONAL COVERAGE.
556	! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
557	! ----------------------------------------------------------------------
558	CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
559	! Don't limit snow cover fraction over permanent ice kmh 2008/03/25
560	if ( ICE == 0 ) then
561	SNCOVR = MIN(SNCOVR,0.98)
562	endif
563	CALL ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,T1,ALBEDO,EMISSI, &
564	DT,SNOWNG,SNOTIME1,LVCOEF)
565	END IF
566	! ----------------------------------------------------------------------
567	! SNOW COVER, ALBEDO OVER SEA-ICE, GLACIAL ICE
568	! ----------------------------------------------------------------------
569	ELSE
570	SNCOVR = 1.0
571	!
572	! Albedo of sea ice
573	!
574	! This value should vary seasonally. 0.65 may be good for Arctic Ocean summer bare ice
575	! value could be as low as 0.4 for Arctic bare ice and melt pond combo (Perovich data)
576	! 0.82 may be good for Arctic spring/fall sea ice (Perovich data)
577	! 0.81 may be good for Antarctic sea ice (Wendler et al. December cruise data)
578	!
579	ALBEDO = 0.80
580	!
581	EMISSI = 0.98
582	END IF
583	! ----------------------------------------------------------------------
584	! THERMAL CONDUCTIVITY FOR SEA-ICE CASE, GLACIAL-ICE CASE
585	! ----------------------------------------------------------------------
586	IF ( (ICE == 1) .or. (ICE == -1) ) THEN
587	DF1 = 2.2
588	!
589	! kmh 09/03/2006
590	! kmh 03/25/2008 change SNCOVR threshold to 0.97
591	!
592	! only apply (small) DF1 conductivity for permanent land ice
593	IF ( (ICE == -1) ) THEN
594	IF ( SNCOVR .GT. 0.97 ) THEN
595	DF1 = SNCOND
596	ENDIF
597	ENDIF
598	!
599	ELSE
600	! ----------------------------------------------------------------------
601	! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
602	! CALCULATION OF THE THERMAL DIFFUSIVITY. TREATMENT OF THE
603	! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
604	! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
605	! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
606	! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
607	! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
608	! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
609	! LIMIT OF VERY THIN SNOWPACK. THIS TREATMENT ALSO ELIMINATES
610	! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
611	! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
612	! ----------------------------------------------------------------------
613	! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
614	! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
615	! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
616	! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
617	! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
618	! ----------------------------------------------------------------------
619	! ----------------------------------------------------------------------
620	! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
621	! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
622	! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
623	! ----------------------------------------------------------------------
624	CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
625
626	!urban
627	IF ( VEGTYP == ISURBAN ) DF1=3.24
628
629	DF1 = DF1 * EXP (SBETA * SHDFAC)
630	!
631	! kmh 09/03/2006
632	! kmh 03/25/2008 change SNCOVR threshold to 0.97
633	!
634	IF ( SNCOVR .GT. 0.97 ) THEN
635	DF1 = SNCOND
636	ENDIF
637	!
638	! ----------------------------------------------------------------------
639	! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
640	! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
641	! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
642	! ----------------------------------------------------------------------
643	END IF
644
645	DSOIL = - (0.5 * ZSOIL (1))
646	IF (SNEQV == 0.) THEN
647	SSOIL = DF1 * (T1- STC (1) ) / DSOIL
648	ELSE
649	DTOT = SNOWH + DSOIL
650	FRCSNO = SNOWH / DTOT
651
652	! 1. HARMONIC MEAN (SERIES FLOW)
653	! DF1 = (SNCONDDF1)/(FRCSOISNCOND+FRCSNO*DF1)
654	FRCSOI = DSOIL / DTOT
655	! 2. ARITHMETIC MEAN (PARALLEL FLOW)
656	! DF1 = FRCSNOSNCOND + FRCSOIDF1
657	DF1H = (SNCOND * DF1)/ (FRCSOI * SNCOND+ FRCSNO * DF1)
658
659	! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
660	! DF1 = (SNCOND*FRCSNO)(DF1**FRCSOI)
661	! weigh DF by snow fraction
662	! DF1 = DF1HSNCOVR + DF1A(1.0-SNCOVR)
663	! DF1 = DF1HSNCOVR + DF1(1.0-SNCOVR)
664	DF1A = FRCSNO * SNCOND+ FRCSOI * DF1
665
666	! ----------------------------------------------------------------------
667	! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
668	! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
669	! MID-LAYER SOIL TEMPERATURE
670	! ----------------------------------------------------------------------
671	DF1 = DF1A * SNCOVR + DF1* (1.0- SNCOVR)
672	IF ( ICE /= 0 ) then
673	! kmh 12/15/2005 correct for too deep snow layer
674	! kmh 09/03/2006 adjust DTOT
675	IF ( DTOT .GT. 2.*DSOIL ) then
676	DTOT = 2.*DSOIL
677	ENDIF
678	ENDIF
679	SSOIL = DF1 * (T1- STC (1) ) / DTOT
680	END IF
681	! ----------------------------------------------------------------------
682	! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
683	! THE PREVIOUS TIMESTEP.
684	! ----------------------------------------------------------------------
685	IF (SNCOVR > 0. ) THEN
686	CALL SNOWZ0 (SNCOVR,Z0,Z0BRD,SNOWH)
687	ELSE
688	Z0=Z0BRD
689	END IF
690	! ----------------------------------------------------------------------
691	! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
692	! HEAT AND MOISTURE.
693
694	! NOTE !!!
695	! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
696	! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
697	! (CZIL) ARE SET THERE VIA NAMELIST I/O.
698
699	! NOTE !!!
700	! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
701	! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. HENCE THE CH
702	! RETURNED FROM SFCDIF HAS UNITS OF M/S. THE IMPORTANT COMPANION
703	! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
704	! AIR DENSITY AND PARAMETER "CP". "RCH" IS COMPUTED IN "CALL PENMAN".
705	! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
706
707	! NOTE !!!
708	! ----------------------------------------------------------------------
709	! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
710	! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT. Needed as a state variable
711	! for iterative/implicit solution of CH in SFCDIF
712	! ----------------------------------------------------------------------
713	! IF(.NOT.LCH) THEN
714	! T1V = T1 * (1.0+ 0.61 * Q2)
715	! TH2V = TH2 * (1.0+ 0.61 * Q2)
716	! CALL SFCDIF_off (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
717	! ENDIF
718
719	! ----------------------------------------------------------------------
720	! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
721	! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
722	! CALCULATIONS.
723	! ----------------------------------------------------------------------
724	! ----------------------------------------------------------------------
725	! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
726	! PENMAN EP SUBROUTINE THAT FOLLOWS
727	! ----------------------------------------------------------------------
728	! FDOWN = SOLDN * (1.0- ALBEDO) + LWDN
729	FDOWN = SOLNET + LWDN
730	! ----------------------------------------------------------------------
731	! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
732	! PENMAN.
733	T2V = SFCTMP * (1.0+ 0.61 * Q2 )
734
735	iout=0
736	if(iout.eq.1) then
737	print*,'before penman'
738	print*,' SFCTMP',SFCTMP,'SFCPRS',SFCPRS,'CH',CH,'T2V',T2V, &
739	'TH2',TH2,'PRCP',PRCP,'FDOWN',FDOWN,'T24',T24,'SSOIL',SSOIL, &
740	'Q2',Q2,'Q2SAT',Q2SAT,'ETP',ETP,'RCH',RCH, &
741	'EPSCA',EPSCA,'RR',RR ,'SNOWNG',SNOWNG,'FRZGRA',FRZGRA, &
742	'DQSDT2',DQSDT2,'FLX2',FLX2,'SNOWH',SNOWH,'SNEQV',SNEQV, &
743	' DSOIL',DSOIL,' FRCSNO',FRCSNO,' SNCOVR',SNCOVR,' DTOT',DTOT, &
744	' ZSOIL (1)',ZSOIL(1),' DF1',DF1,'T1',T1,' STC1',STC(1), &
745	'ALBEDO',ALBEDO,'SMC',SMC,'STC',STC,'SH2O',SH2O
746	endif
747
748	CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
749	Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
750	!
751	! kmh 01/09/2007 add T1,ICE,SNCOVR to call
752	!
753	DQSDT2,FLX2,EMISSI,SNEQV,T1,ICE,SNCOVR)
754	!
755	! ----------------------------------------------------------------------
756	! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
757	! IF NONZERO GREENNESS FRACTION
758	! ----------------------------------------------------------------------
759
760	! ----------------------------------------------------------------------
761	! FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
762	! BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
763	! ----------------------------------------------------------------------
764	IF (SHDFAC > 0.) THEN
765	CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, &
766	SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
767	TOPT,RSMAX,RGL,HS,XLAI, &
768	RCS,RCT,RCQ,RCSOIL,EMISSI)
769	ELSE
770	RC = 0.0
771	END IF
772	! ----------------------------------------------------------------------
773	! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
774	! EXISTS OR NOT:
775	! ----------------------------------------------------------------------
776	ESNOW = 0.0
777	IF (SNEQV == 0.0) THEN
778	CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
779	SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
780	SHDFAC, &
781	SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
782	SSOIL, &
783	STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
784	SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL, &
785	DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
786	RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
787	QUARTZ,FXEXP,CSOIL, &
788	BETA,DRIP,DEW,FLX1,FLX3,VEGTYP,ISURBAN)
789	ETA_KINEMATIC = ETA
790	ELSE
791	CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
792	SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
793	SBETA,DF1, &
794	Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
795	SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS,&
796	SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT, &
797	ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
798	RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
799	ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
800	BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI, &
801	RIBB,SOLDN, &
802	ISURBAN, &
803	VEGTYP)
804	ETA_KINEMATIC = ESNOW + ETNS
805	END IF
806
807	! Calculate effective mixing ratio at grnd level (skin)
808	!
809	! Q1=Q2+ETA*CP/RCH
810	Q1=Q2+ETA_KINEMATIC*CP/RCH
811	!
812	! ----------------------------------------------------------------------
813	! DETERMINE SENSIBLE HEAT (H) IN ENERGY UNITS (W M-2)
814	! ----------------------------------------------------------------------
815	SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )
816
817	! ----------------------------------------------------------------------
818	! CONVERT EVAP TERMS FROM KINEMATIC (KG M-2 S-1) TO ENERGY UNITS (W M-2)
819	! ----------------------------------------------------------------------
820	EDIR = EDIR * LVH2O
821	EC = EC * LVH2O
822	DO K=1,4
823	ET(K) = ET(K) * LVH2O
824	ENDDO
825	ETT = ETT * LVH2O
826	ESNOW = ESNOW * LSUBS
827	ETP = ETP((1.-SNCOVR)LVH2O + SNCOVR*LSUBS)
828	IF (ETP .GT. 0.) THEN
829	ETA = EDIR + EC + ETT + ESNOW
830	ELSE
831	ETA = ETP
832	ENDIF
833	! ----------------------------------------------------------------------
834	! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
835	! ----------------------------------------------------------------------
836	IF (ETP == 0.0) THEN
837	BETA = 0.0
838	ELSE
839	BETA = ETA/ETP
840	ENDIF
841
842	! ----------------------------------------------------------------------
843	! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
844	! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
845	! SSOIL<0: COOL THE SURFACE (DAY TIME)
846	! ----------------------------------------------------------------------
847	SSOIL = -1.0* SSOIL
848
849	! ----------------------------------------------------------------------
850	! FOR THE CASE OF LAND (BUT NOT GLACIAL-ICE):
851	! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
852	! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW. RUNOFF2 IS ALREADY
853	! A RATE AT THIS POINT
854	! ----------------------------------------------------------------------
855	IF (ICE == 0) THEN
856	RUNOFF3 = RUNOFF3/ DT
857	RUNOFF2 = RUNOFF2+ RUNOFF3
858	SOILM = -1.0* SMC (1)* ZSOIL (1)
859	DO K = 2,NSOIL
860	SOILM = SOILM + SMC (K)* (ZSOIL (K -1) - ZSOIL (K))
861	END DO
862	SOILWM = -1.0* (SMCMAX - SMCWLT)* ZSOIL (1)
863	SOILWW = -1.0* (SMC (1) - SMCWLT)* ZSOIL (1)
864	!
865	DO K = 1,NSOIL
866	SMAV(K)=(SMC(K) - SMCWLT)/(SMCMAX - SMCWLT)
867	END DO
868
869	IF (NROOT >= 2) THEN
870	DO K = 2,NROOT
871	SOILWM = SOILWM + (SMCMAX - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
872	SOILWW = SOILWW + (SMC(K) - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
873	END DO
874	END IF
875	IF (SOILWM .LT. 1.E-6) THEN
876	SOILWM = 0.0
877	SOILW = 0.0
878	SOILM = 0.0
879	ELSE
880	SOILW = SOILWW / SOILWM
881	END IF
882	ELSE
883	! ----------------------------------------------------------------------
884	! FOR THE CASE OF SEA-ICE (ICE=1) OR GLACIAL-ICE (ICE=-1), ADD ANY
885	! SNOWMELT DIRECTLY TO SURFACE RUNOFF (RUNOFF1) SINCE THERE IS NO
886	! SOIL MEDIUM, AND THUS NO CALL TO SUBROUTINE SMFLX (FOR SOIL MOISTURE
887	! TENDENCY).
888	! ----------------------------------------------------------------------
889	RUNOFF1 = SNOMLT/DT
890	SOILWM = 0.0
891	SOILW = 0.0
892	SOILM = 0.0
893	DO K = 1,NSOIL
894	SMAV(K)= 1.0
895	END DO
896	END IF
897
898	! ----------------------------------------------------------------------
899	END SUBROUTINE SFLX
900	! ----------------------------------------------------------------------
901
902	SUBROUTINE ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO,EMISSI, &
903	DT,SNOWNG,SNOTIME1,LVCOEF)
904
905	! ----------------------------------------------------------------------
906	! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
907	! ALB SNOWFREE ALBEDO
908	! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
909	! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
910	! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
911	! SNCOVR FRACTIONAL SNOW COVER
912	! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
913	! TSNOW SNOW SURFACE TEMPERATURE (K)
914	! ----------------------------------------------------------------------
915	IMPLICIT NONE
916
917	! ----------------------------------------------------------------------
918	! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
919	! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
920	! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
921	! (1985, JCAM, VOL 24, 402-411)
922	! ----------------------------------------------------------------------
923	REAL, INTENT(IN) :: ALB, SNOALB, EMBRD, SHDFAC, SHDMIN, SNCOVR, TSNOW
924	REAL, INTENT(IN) :: DT
925	LOGICAL, INTENT(IN) :: SNOWNG
926	REAL, INTENT(INOUT):: SNOTIME1
927	REAL, INTENT(OUT) :: ALBEDO, EMISSI
928	REAL :: SNOALB2
929	REAL :: TM,SNOALB1
930	REAL, INTENT(IN) :: LVCOEF
931	REAL, PARAMETER :: SNACCA=0.94,SNACCB=0.58,SNTHWA=0.82,SNTHWB=0.46
932	! turn of vegetation effect
933	! ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
934	! ALBEDO = (1.0-SNCOVR)ALB + SNCOVRSNOALB !this is equivalent to below
935	ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
936	EMISSI = EMBRD + SNCOVR*(EMISSI_S - EMBRD)
937
938	! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
939	! IF (TSNOW.LE.263.16) THEN
940	! ALBEDO=SNOALB
941	! ELSE
942	! IF (TSNOW.LT.273.16) THEN
943	! TM=0.1*(TSNOW-263.16)
944	! SNOALB1=0.5((0.9-0.2(TM*3))+(0.8-0.16(TM**3)))
945	! ELSE
946	! SNOALB1=0.67
947	! IF(SNCOVR.GT.0.95) SNOALB1= 0.6
948	! SNOALB1 = ALB + SNCOVR*(SNOALB-ALB)
949	! ENDIF
950	! ENDIF
951	! ALBEDO = ALB + SNCOVR*(SNOALB1-ALB)
952
953	! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
954	! SNOALB1 = SNOALB+COEF*(0.85-SNOALB)
955	! SNOALB2=SNOALB1
956	!!m LSTSNW=LSTSNW+1
957	! SNOTIME1 = SNOTIME1 + DT
958	! IF (SNOWNG) THEN
959	! SNOALB2=SNOALB
960	!!m LSTSNW=0
961	! SNOTIME1 = 0.0
962	! ELSE
963	! IF (TSNOW.LT.273.16) THEN
964	!! SNOALB2=SNOALB-0.008LSTSNWDT/86400
965	!!m SNOALB2=SNOALB-0.008*SNOTIME1/86400
966	! SNOALB2=(SNOALB2-0.65)EXP(-0.05DT/3600)+0.65
967	!! SNOALB2=(ALBEDO-0.65)EXP(-0.01DT/3600)+0.65
968	! ELSE
969	! SNOALB2=(SNOALB2-0.5)EXP(-0.0005DT/3600)+0.5
970	!! SNOALB2=(SNOALB-0.5)EXP(-0.24LSTSNW*DT/86400)+0.5
971	!!m SNOALB2=(SNOALB-0.5)EXP(-0.24SNOTIME1/86400)+0.5
972	! ENDIF
973	! ENDIF
974	!
975	!! print*,'SNOALB2',SNOALB2,'ALBEDO',ALBEDO,'DT',DT
976	! ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
977	! IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
978	!!m LSTSNW1=LSTSNW
979	!! SNOTIME = SNOTIME1
980
981	! formulation by Livneh
982	! ----------------------------------------------------------------------
983	! SNOALB IS CONSIDERED AS THE MAXIMUM SNOW ALBEDO FOR NEW SNOW, AT
984	! A VALUE OF 85%. SNOW ALBEDO CURVE DEFAULTS ARE FROM BRAS P.263. SHOULD
985	! NOT BE CHANGED EXCEPT FOR SERIOUS PROBLEMS WITH SNOW MELT.
986	! TO IMPLEMENT ACCUMULATIN PARAMETERS, SNACCA AND SNACCB, ASSERT THAT IT
987	! IS INDEED ACCUMULATION SEASON. I.E. THAT SNOW SURFACE TEMP IS BELOW
988	! ZERO AND THE DATE FALLS BETWEEN OCTOBER AND FEBRUARY
989	! ----------------------------------------------------------------------
990	SNOALB1 = SNOALB+LVCOEF*(0.85-SNOALB)
991	SNOALB2=SNOALB1
992	! ---------------- Initial LSTSNW --------------------------------------
993	IF (SNOWNG) THEN
994	SNOTIME1 = 0.
995	ELSE
996	SNOTIME1=SNOTIME1+DT
997	! IF (TSNOW.LT.273.16) THEN
998	SNOALB2=SNOALB1(SNACCA((SNOTIME1/86400.0)*SNACCB))
999	! ELSE
1000	! SNOALB2 =SNOALB1(SNTHWA((SNOTIME1/86400.0)*SNTHWB))
1001	! ENDIF
1002	ENDIF
1003	!
1004	SNOALB2 = MAX ( SNOALB2, ALB )
1005	ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
1006	IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
1007
1008	! IF (TSNOW.LT.273.16) THEN
1009	! ALBEDO=SNOALB-0.008*DT/86400
1010	! ELSE
1011	! ALBEDO=(SNOALB-0.5)EXP(-0.24DT/86400)+0.5
1012	! ENDIF
1013
1014	! IF (ALBEDO > SNOALB) ALBEDO = SNOALB
1015
1016	! ----------------------------------------------------------------------
1017	END SUBROUTINE ALCALC
1018	! ----------------------------------------------------------------------
1019
1020	SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
1021	SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1022	TOPT,RSMAX,RGL,HS,XLAI, &
1023	RCS,RCT,RCQ,RCSOIL,EMISSI)
1024
1025	! ----------------------------------------------------------------------
1026	! SUBROUTINE CANRES
1027	! ----------------------------------------------------------------------
1028	! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
1029	! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
1030	! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
1031	! MOISTURE RATHER THAN TOTAL)
1032	! ----------------------------------------------------------------------
1033	! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
1034	! NOILHAN (1990, BLM)
1035	! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
1036	! AND TABLE 2 OF SEC. 3.1.2
1037	! ----------------------------------------------------------------------
1038	! INPUT:
1039	! SOLAR INCOMING SOLAR RADIATION
1040	! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1041	! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
1042	! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1043	! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1044	! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
1045	! SFCPRS SURFACE PRESSURE
1046	! SMC VOLUMETRIC SOIL MOISTURE
1047	! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
1048	! NSOIL NO. OF SOIL LAYERS
1049	! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
1050	! XLAI LEAF AREA INDEX
1051	! SMCWLT WILTING POINT
1052	! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
1053	! SETS IN)
1054	! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
1055	! SURBOUTINE REDPRM
1056	! OUTPUT:
1057	! PC PLANT COEFFICIENT
1058	! RC CANOPY RESISTANCE
1059	! ----------------------------------------------------------------------
1060
1061	IMPLICIT NONE
1062	INTEGER, INTENT(IN) :: NROOT,NSOIL
1063	INTEGER K
1064	REAL, INTENT(IN) :: CH,DQSDT2,HS,Q2,Q2SAT,RSMIN,RGL,RSMAX, &
1065	SFCPRS,SFCTMP,SMCREF,SMCWLT, SOLAR,TOPT,XLAI, &
1066	EMISSI
1067	REAL,DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
1068	REAL, INTENT(OUT):: PC,RC,RCQ,RCS,RCSOIL,RCT
1069	REAL :: DELTA,FF,GX,P,RR
1070	REAL, DIMENSION(1:NSOIL) :: PART
1071	REAL, PARAMETER :: SLV = 2.501000E6
1072
1073
1074	! ----------------------------------------------------------------------
1075	! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
1076	! ----------------------------------------------------------------------
1077	RCS = 0.0
1078	RCT = 0.0
1079	RCQ = 0.0
1080	RCSOIL = 0.0
1081
1082	! ----------------------------------------------------------------------
1083	! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
1084	! ----------------------------------------------------------------------
1085	RC = 0.0
1086	FF = 0.552.0 SOLAR / (RGL * XLAI)
1087	RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
1088
1089	! ----------------------------------------------------------------------
1090	! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
1091	! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
1092	! ----------------------------------------------------------------------
1093	RCS = MAX (RCS,0.0001)
1094	RCT = 1.0- 0.0016* ( (TOPT - SFCTMP)**2.0)
1095
1096	! ----------------------------------------------------------------------
1097	! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
1098	! RCQ EXPRESSION FROM SSIB
1099	! ----------------------------------------------------------------------
1100	RCT = MAX (RCT,0.0001)
1101	RCQ = 1.0/ (1.0+ HS * (Q2SAT - Q2))
1102
1103	! ----------------------------------------------------------------------
1104	! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
1105	! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
1106	! ----------------------------------------------------------------------
1107	RCQ = MAX (RCQ,0.01)
1108	GX = (SMC (1) - SMCWLT) / (SMCREF - SMCWLT)
1109	IF (GX > 1.) GX = 1.
1110	IF (GX < 0.) GX = 0.
1111
1112	! ----------------------------------------------------------------------
1113	! USE SOIL DEPTH AS WEIGHTING FACTOR
1114	! ----------------------------------------------------------------------
1115	! ----------------------------------------------------------------------
1116	! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1117	! PART(1) = RTDIS(1) * GX
1118	! ----------------------------------------------------------------------
1119	PART (1) = (ZSOIL (1)/ ZSOIL (NROOT)) * GX
1120	DO K = 2,NROOT
1121	GX = (SMC (K) - SMCWLT) / (SMCREF - SMCWLT)
1122	IF (GX > 1.) GX = 1.
1123	IF (GX < 0.) GX = 0.
1124	! ----------------------------------------------------------------------
1125	! USE SOIL DEPTH AS WEIGHTING FACTOR
1126	! ----------------------------------------------------------------------
1127	! ----------------------------------------------------------------------
1128	! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1129	! PART(K) = RTDIS(K) * GX
1130	! ----------------------------------------------------------------------
1131	PART (K) = ( (ZSOIL (K) - ZSOIL (K -1))/ ZSOIL (NROOT)) * GX
1132	END DO
1133	DO K = 1,NROOT
1134	RCSOIL = RCSOIL + PART (K)
1135	END DO
1136
1137	! ----------------------------------------------------------------------
1138	! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
1139	! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
1140	! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
1141	! PC * LINERIZED PENMAN POTENTIAL EVAP =
1142	! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
1143	! ----------------------------------------------------------------------
1144	RCSOIL = MAX (RCSOIL,0.0001)
1145
1146	RC = RSMIN / (XLAI * RCS * RCT * RCQ * RCSOIL)
1147	! RR = (4.* SIGMA * RD / CP)* (SFCTMP *4.)/ (SFCPRS CH) + 1.0
1148	RR = (4.* EMISSI SIGMA RD / CP)* (SFCTMP *4.)/ (SFCPRS CH) &
1149	+ 1.0
1150
1151	DELTA = (SLV / CP)* DQSDT2
1152
1153	PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
1154
1155	! ----------------------------------------------------------------------
1156	END SUBROUTINE CANRES
1157	! ----------------------------------------------------------------------
1158
1159	SUBROUTINE CSNOW (SNCOND,DSNOW)
1160
1161	! ----------------------------------------------------------------------
1162	! SUBROUTINE CSNOW
1163	! FUNCTION CSNOW
1164	! ----------------------------------------------------------------------
1165	! CALCULATE SNOW TERMAL CONDUCTIVITY
1166	! ----------------------------------------------------------------------
1167	IMPLICIT NONE
1168	REAL, INTENT(IN) :: DSNOW
1169	REAL, INTENT(OUT):: SNCOND
1170	REAL :: C
1171	REAL, PARAMETER :: UNIT = 0.11631
1172
1173	! ----------------------------------------------------------------------
1174	! SNCOND IN UNITS OF CAL/(CMHRC), RETURNED IN W/(M*C)
1175	! CSNOW IN UNITS OF CAL/(CMHRC), RETURNED IN W/(M*C)
1176	! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
1177	! ----------------------------------------------------------------------
1178	C = 0.32810* (2.25* DSNOW)
1179	! CSNOW=UNIT*C
1180
1181	! ----------------------------------------------------------------------
1182	! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
1183	! ----------------------------------------------------------------------
1184	! SNCOND=0.0293(1.+100.DSNOW**2)
1185	! CSNOW=0.0293(1.+100.DSNOW**2)
1186
1187	! ----------------------------------------------------------------------
1188	! E. ANDERSEN FROM FLERCHINGER
1189	! ----------------------------------------------------------------------
1190	! SNCOND=0.021+2.51DSNOW*2
1191	! CSNOW=0.021+2.51DSNOW*2
1192
1193	! SNCOND = UNIT * C
1194	! double snow thermal conductivity
1195	SNCOND = 2.0 * UNIT * C
1196
1197	! ----------------------------------------------------------------------
1198	END SUBROUTINE CSNOW
1199	! ----------------------------------------------------------------------
1200
1201	SUBROUTINE DEVAP (EDIR,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
1202	DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1203
1204	! ----------------------------------------------------------------------
1205	! SUBROUTINE DEVAP
1206	! FUNCTION DEVAP
1207	! ----------------------------------------------------------------------
1208	! CALCULATE DIRECT SOIL EVAPORATION
1209	! ----------------------------------------------------------------------
1210	IMPLICIT NONE
1211	REAL, INTENT(IN) :: ETP1,SMC,BEXP,DKSAT,DWSAT,FXEXP, &
1212	SHDFAC,SMCDRY,SMCMAX,ZSOIL,SMCREF,SMCWLT
1213	REAL, INTENT(OUT):: EDIR
1214	REAL :: FX, SRATIO
1215
1216
1217	! ----------------------------------------------------------------------
1218	! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
1219	! WHEN FXEXP=1.
1220	! ----------------------------------------------------------------------
1221	! ----------------------------------------------------------------------
1222	! FX > 1 REPRESENTS DEMAND CONTROL
1223	! FX < 1 REPRESENTS FLUX CONTROL
1224	! ----------------------------------------------------------------------
1225
1226	SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
1227	IF (SRATIO > 0.) THEN
1228	FX = SRATIO**FXEXP
1229	FX = MAX ( MIN ( FX, 1. ) ,0. )
1230	ELSE
1231	FX = 0.
1232	ENDIF
1233
1234	! ----------------------------------------------------------------------
1235	! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
1236	! ----------------------------------------------------------------------
1237	EDIR = FX * ( 1.0- SHDFAC ) * ETP1
1238
1239	! ----------------------------------------------------------------------
1240	END SUBROUTINE DEVAP
1241	! ----------------------------------------------------------------------
1242
1243	SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
1244	SH2O, &
1245	SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
1246	SMCREF,SHDFAC,CMCMAX, &
1247	SMCDRY,CFACTR, &
1248	EDIR,EC,ET,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1249
1250	! ----------------------------------------------------------------------
1251	! SUBROUTINE EVAPO
1252	! ----------------------------------------------------------------------
1253	! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
1254	! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
1255	! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
1256	! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
1257	! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
1258	! ----------------------------------------------------------------------
1259	IMPLICIT NONE
1260	INTEGER, INTENT(IN) :: NSOIL, NROOT
1261	INTEGER :: I,K
1262	REAL, INTENT(IN) :: BEXP, CFACTR,CMC,CMCMAX,DKSAT, &
1263	DT,DWSAT,ETP1,FXEXP,PC,Q2,SFCTMP, &
1264	SHDFAC,SMCDRY,SMCMAX,SMCREF,SMCWLT
1265	REAL, INTENT(OUT) :: EC,EDIR,ETA1,ETT
1266	REAL :: CMC2MS
1267	REAL,DIMENSION(1:NSOIL), INTENT(IN) :: RTDIS, SMC, SH2O, ZSOIL
1268	REAL,DIMENSION(1:NSOIL), INTENT(OUT) :: ET
1269
1270	! ----------------------------------------------------------------------
1271	! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
1272	! GREATER THAN ZERO.
1273	! ----------------------------------------------------------------------
1274	EDIR = 0.
1275	EC = 0.
1276	ETT = 0.
1277	DO K = 1,NSOIL
1278	ET (K) = 0.
1279	END DO
1280
1281	! ----------------------------------------------------------------------
1282	! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
1283	! ONLY IF VEG COVER NOT COMPLETE.
1284	! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
1285	! ----------------------------------------------------------------------
1286	IF (ETP1 > 0.0) THEN
1287	IF (SHDFAC < 1.) THEN
1288	CALL DEVAP (EDIR,ETP1,SMC (1),ZSOIL (1),SHDFAC,SMCMAX, &
1289	BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1290	END IF
1291	! ----------------------------------------------------------------------
1292	! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
1293	! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
1294	! ----------------------------------------------------------------------
1295
1296	IF (SHDFAC > 0.0) THEN
1297	CALL TRANSP (ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
1298	CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
1299	DO K = 1,NSOIL
1300	ETT = ETT + ET ( K )
1301	END DO
1302	! ----------------------------------------------------------------------
1303	! CALCULATE CANOPY EVAPORATION.
1304	! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
1305	! ----------------------------------------------------------------------
1306	IF (CMC > 0.0) THEN
1307	EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
1308	ELSE
1309	EC = 0.0
1310	END IF
1311	! ----------------------------------------------------------------------
1312	! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
1313	! CANOPY. -F.CHEN, 18-OCT-1994
1314	! ----------------------------------------------------------------------
1315	CMC2MS = CMC / DT
1316	EC = MIN ( CMC2MS, EC )
1317	END IF
1318	END IF
1319	! ----------------------------------------------------------------------
1320	! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
1321	! ----------------------------------------------------------------------
1322	ETA1 = EDIR + ETT + EC
1323
1324	! ----------------------------------------------------------------------
1325	END SUBROUTINE EVAPO
1326	! ----------------------------------------------------------------------
1327
1328	SUBROUTINE FAC2MIT(SMCMAX,FLIMIT)
1329	IMPLICIT NONE
1330	REAL, INTENT(IN) :: SMCMAX
1331	REAL, INTENT(OUT) :: FLIMIT
1332
1333	FLIMIT = 0.90
1334
1335	IF ( SMCMAX == 0.395 ) THEN
1336	FLIMIT = 0.59
1337	ELSE IF ( ( SMCMAX == 0.434 ) .OR. ( SMCMAX == 0.404 ) ) THEN
1338	FLIMIT = 0.85
1339	ELSE IF ( ( SMCMAX == 0.465 ) .OR. ( SMCMAX == 0.406 ) ) THEN
1340	FLIMIT = 0.86
1341	ELSE IF ( ( SMCMAX == 0.476 ) .OR. ( SMCMAX == 0.439 ) ) THEN
1342	FLIMIT = 0.74
1343	ELSE IF ( ( SMCMAX == 0.200 ) .OR. ( SMCMAX == 0.464 ) ) THEN
1344	FLIMIT = 0.80
1345	ENDIF
1346
1347	! ----------------------------------------------------------------------
1348	END SUBROUTINE FAC2MIT
1349	! ----------------------------------------------------------------------
1350
1351	SUBROUTINE FRH2O (FREE,TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
1352
1353	! ----------------------------------------------------------------------
1354	! SUBROUTINE FRH2O
1355	! ----------------------------------------------------------------------
1356	! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
1357	! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
1358	! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
1359	! (1999, JGR, VOL 104(D16), 19569-19585).
1360	! ----------------------------------------------------------------------
1361	! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
1362	! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
1363	! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
1364	! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
1365	! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
1366	! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
1367	! LIMIT OF FREEZING POINT TEMPERATURE T0.
1368	! ----------------------------------------------------------------------
1369	! INPUT:
1370
1371	! TKELV.........TEMPERATURE (Kelvin)
1372	! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
1373	! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
1374	! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
1375	! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
1376	! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
1377
1378	! OUTPUT:
1379	! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
1380	! FREE..........SUPERCOOLED LIQUID WATER CONTENT
1381	! ----------------------------------------------------------------------
1382	IMPLICIT NONE
1383	REAL, INTENT(IN) :: BEXP,PSIS,SH2O,SMC,SMCMAX,TKELV
1384	REAL, INTENT(OUT) :: FREE
1385	REAL :: BX,DENOM,DF,DSWL,FK,SWL,SWLK
1386	INTEGER :: NLOG,KCOUNT
1387	! PARAMETER(CK = 0.0)
1388	REAL, PARAMETER :: CK = 8.0, BLIM = 5.5, ERROR = 0.005, &
1389	HLICE = 3.335E5, GS = 9.81,DICE = 920.0, &
1390	DH2O = 1000.0, T0 = 273.15
1391
1392	! ----------------------------------------------------------------------
1393	! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
1394	! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
1395	! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
1396	! ----------------------------------------------------------------------
1397	BX = BEXP
1398
1399	! ----------------------------------------------------------------------
1400	! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
1401	! ----------------------------------------------------------------------
1402	IF (BEXP > BLIM) BX = BLIM
1403	NLOG = 0
1404
1405	! ----------------------------------------------------------------------
1406	! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
1407	! ----------------------------------------------------------------------
1408	KCOUNT = 0
1409	! FRH2O = SMC
1410	IF (TKELV > (T0- 1.E-3)) THEN
1411	FREE = SMC
1412	ELSE
1413
1414	! ----------------------------------------------------------------------
1415	! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
1416	! IN KOREN ET AL, JGR, 1999, EQN 17
1417	! ----------------------------------------------------------------------
1418	! INITIAL GUESS FOR SWL (frozen content)
1419	! ----------------------------------------------------------------------
1420	IF (CK /= 0.0) THEN
1421	SWL = SMC - SH2O
1422	! ----------------------------------------------------------------------
1423	! KEEP WITHIN BOUNDS.
1424	! ----------------------------------------------------------------------
1425	IF (SWL > (SMC -0.02)) SWL = SMC -0.02
1426
1427	! ----------------------------------------------------------------------
1428	! START OF ITERATIONS
1429	! ----------------------------------------------------------------------
1430	IF (SWL < 0.) SWL = 0.
1431	1001 Continue
1432	IF (.NOT.( (NLOG < 10) .AND. (KCOUNT == 0))) goto 1002
1433	NLOG = NLOG +1
1434	DF = ALOG ( ( PSIS * GS / HLICE ) * ( ( 1. + CK * SWL )*2.) &
1435	( SMCMAX / (SMC - SWL) )** BX) - ALOG ( - ( &
1436	TKELV - T0)/ TKELV)
1437	DENOM = 2. * CK / ( 1. + CK * SWL ) + BX / ( SMC - SWL )
1438	SWLK = SWL - DF / DENOM
1439	! ----------------------------------------------------------------------
1440	! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
1441	! ----------------------------------------------------------------------
1442	IF (SWLK > (SMC -0.02)) SWLK = SMC - 0.02
1443	IF (SWLK < 0.) SWLK = 0.
1444
1445	! ----------------------------------------------------------------------
1446	! MATHEMATICAL SOLUTION BOUNDS APPLIED.
1447	! ----------------------------------------------------------------------
1448	DSWL = ABS (SWLK - SWL)
1449
1450	! ----------------------------------------------------------------------
1451	! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
1452	! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
1453	! ----------------------------------------------------------------------
1454	SWL = SWLK
1455	IF ( DSWL <= ERROR ) THEN
1456	KCOUNT = KCOUNT +1
1457	END IF
1458	! ----------------------------------------------------------------------
1459	! END OF ITERATIONS
1460	! ----------------------------------------------------------------------
1461	! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
1462	! ----------------------------------------------------------------------
1463	! FRH2O = SMC - SWL
1464	goto 1001
1465	1002 continue
1466	FREE = SMC - SWL
1467	END IF
1468	! ----------------------------------------------------------------------
1469	! END OPTION 1
1470	! ----------------------------------------------------------------------
1471	! ----------------------------------------------------------------------
1472	! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
1473	! IN KOREN ET AL., JGR, 1999, EQN 17
1474	! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
1475	! ----------------------------------------------------------------------
1476	IF (KCOUNT == 0) THEN
1477	PRINT *,'Flerchinger USEd in NEW version. Iterations=',NLOG
1478	FK = ( ( (HLICE / (GS * ( - PSIS)))* &
1479	( (TKELV - T0)/ TKELV))** ( -1/ BX))* SMCMAX
1480	! FRH2O = MIN (FK, SMC)
1481	IF (FK < 0.02) FK = 0.02
1482	FREE = MIN (FK, SMC)
1483	! ----------------------------------------------------------------------
1484	! END OPTION 2
1485	! ----------------------------------------------------------------------
1486	END IF
1487	END IF
1488	! ----------------------------------------------------------------------
1489	END SUBROUTINE FRH2O
1490	! ----------------------------------------------------------------------
1491
1492	SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
1493	TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
1494	F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP,ISURBAN)
1495
1496	! ----------------------------------------------------------------------
1497	! SUBROUTINE HRT
1498	! ----------------------------------------------------------------------
1499	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1500	! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
1501	! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1502	! ----------------------------------------------------------------------
1503	IMPLICIT NONE
1504	LOGICAL :: ITAVG
1505	INTEGER, INTENT(IN) :: NSOIL, VEGTYP
1506	INTEGER, INTENT(IN) :: ISURBAN
1507	INTEGER :: I, K
1508
1509	REAL, INTENT(IN) :: BEXP, CSOIL, DF1, DT,F1,PSISAT,QUARTZ, &
1510	SMCMAX ,TBOT,YY,ZZ1, ZBOT
1511	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,STC,ZSOIL
1512	REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SH2O
1513	REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTS
1514	REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI,CI
1515	REAL :: DDZ, DDZ2, DENOM, DF1N, DF1K, DTSDZ, &
1516	DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK, &
1517	TBK1,TSNSR,TSURF,CSOIL_LOC
1518	REAL, PARAMETER :: T0 = 273.15, CAIR = 1004.0, CICE = 2.106E6,&
1519	CH2O = 4.2E6
1520
1521
1522	!urban
1523	IF( VEGTYP == ISURBAN ) then
1524	CSOIL_LOC=3.0E6
1525	ELSE
1526	CSOIL_LOC=CSOIL
1527	ENDIF
1528
1529	! ----------------------------------------------------------------------
1530	! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
1531	! ----------------------------------------------------------------------
1532	ITAVG = .TRUE.
1533	! ----------------------------------------------------------------------
1534	! BEGIN SECTION FOR TOP SOIL LAYER
1535	! ----------------------------------------------------------------------
1536	! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
1537	! ----------------------------------------------------------------------
1538	HCPCT = SH2O (1)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC (1))&
1539	* CAIR &
1540	+ ( SMC (1) - SH2O (1) )* CICE
1541
1542	! ----------------------------------------------------------------------
1543	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1544	! ----------------------------------------------------------------------
1545	DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
1546	AI (1) = 0.0
1547	CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
1548
1549	! ----------------------------------------------------------------------
1550	! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
1551	! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
1552	! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
1553	! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
1554	! ----------------------------------------------------------------------
1555	BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT * &
1556	ZZ1)
1557	DTSDZ = (STC (1) - STC (2)) / ( -0.5 * ZSOIL (2))
1558	SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
1559	! RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
1560	DENOM = (ZSOIL (1) * HCPCT)
1561
1562	! ----------------------------------------------------------------------
1563	! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
1564	! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
1565	! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
1566	! ----------------------------------------------------------------------
1567	! QTOT = SSOIL - DF1*DTSDZ
1568	RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
1569
1570	! ----------------------------------------------------------------------
1571	! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
1572	! ----------------------------------------------------------------------
1573	QTOT = -1.0* RHSTS (1)* DENOM
1574
1575	! ----------------------------------------------------------------------
1576	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
1577	! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
1578	! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
1579	! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
1580	! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
1581	! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
1582	! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
1583	! LATER IN FUNCTION SUBROUTINE SNKSRC
1584	! ----------------------------------------------------------------------
1585	SICE = SMC (1) - SH2O (1)
1586	IF (ITAVG) THEN
1587	TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
1588	! ----------------------------------------------------------------------
1589	! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
1590	! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
1591	! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
1592	! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
1593	! ----------------------------------------------------------------------
1594	CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
1595	IF ( (SICE > 0.) .OR. (STC (1) < T0) .OR. &
1596	(TSURF < T0) .OR. (TBK < T0) ) THEN
1597	! TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1),
1598	CALL TMPAVG (TAVG,TSURF,STC (1),TBK,ZSOIL,NSOIL,1)
1599	CALL SNKSRC (TSNSR,TAVG,SMC (1),SH2O (1), &
1600	ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1601	! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1602	RHSTS (1) = RHSTS (1) - TSNSR / DENOM
1603	END IF
1604	ELSE
1605	! TSNSR = SNKSRC (STC(1),SMC(1),SH2O(1),
1606	IF ( (SICE > 0.) .OR. (STC (1) < T0) ) THEN
1607	CALL SNKSRC (TSNSR,STC (1),SMC (1),SH2O (1), &
1608	ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1609	! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1610	RHSTS (1) = RHSTS (1) - TSNSR / DENOM
1611	END IF
1612	! ----------------------------------------------------------------------
1613	! THIS ENDS SECTION FOR TOP SOIL LAYER.
1614	! ----------------------------------------------------------------------
1615	END IF
1616
1617	! INITIALIZE DDZ2
1618	! ----------------------------------------------------------------------
1619
1620	DDZ2 = 0.0
1621	DF1K = DF1
1622
1623	! ----------------------------------------------------------------------
1624	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1625	! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
1626	! ----------------------------------------------------------------------
1627	! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
1628	! ----------------------------------------------------------------------
1629	DO K = 2,NSOIL
1630	HCPCT = SH2O (K)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC ( &
1631	K))* CAIR + ( SMC (K) - SH2O (K) )* CICE
1632	! ----------------------------------------------------------------------
1633	! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
1634	! ----------------------------------------------------------------------
1635	! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
1636	! ----------------------------------------------------------------------
1637	IF (K /= NSOIL) THEN
1638
1639	! ----------------------------------------------------------------------
1640	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
1641	! ----------------------------------------------------------------------
1642	CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
1643
1644	!urban
1645	IF ( VEGTYP == ISURBAN ) DF1N = 3.24
1646
1647	DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
1648
1649	! ----------------------------------------------------------------------
1650	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
1651	! ----------------------------------------------------------------------
1652	DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
1653	DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
1654
1655	! ----------------------------------------------------------------------
1656	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
1657	! TEMP AT BOTTOM OF LAYER.
1658	! ----------------------------------------------------------------------
1659	CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) * &
1660	HCPCT)
1661	IF (ITAVG) THEN
1662	CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
1663	END IF
1664
1665	ELSE
1666	! ----------------------------------------------------------------------
1667	! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
1668	! BOTTOM LAYER.
1669	! ----------------------------------------------------------------------
1670
1671	! ----------------------------------------------------------------------
1672	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
1673	! ----------------------------------------------------------------------
1674	CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
1675
1676
1677	!urban
1678	IF ( VEGTYP == ISURBAN ) DF1N = 3.24
1679
1680	DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
1681
1682	! ----------------------------------------------------------------------
1683	! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
1684	! ----------------------------------------------------------------------
1685	DTSDZ2 = (STC (K) - TBOT) / DENOM
1686
1687	! ----------------------------------------------------------------------
1688	! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
1689	! TEMP AT BOTTOM OF LAST LAYER.
1690	! ----------------------------------------------------------------------
1691	CI (K) = 0.
1692	IF (ITAVG) THEN
1693	CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
1694	END IF
1695	! ----------------------------------------------------------------------
1696	! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
1697	END IF
1698	! ----------------------------------------------------------------------
1699	! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
1700	! ----------------------------------------------------------------------
1701	DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
1702	RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
1703	QTOT = -1.0* DENOM * RHSTS (K)
1704
1705	SICE = SMC (K) - SH2O (K)
1706	IF (ITAVG) THEN
1707	CALL TMPAVG (TAVG,TBK,STC (K),TBK1,ZSOIL,NSOIL,K)
1708	! TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL,
1709	IF ( (SICE > 0.) .OR. (STC (K) < T0) .OR. &
1710	(TBK .lt. T0) .OR. (TBK1 .lt. T0) ) THEN
1711	CALL SNKSRC (TSNSR,TAVG,SMC (K),SH2O (K),ZSOIL,NSOIL, &
1712	SMCMAX,PSISAT,BEXP,DT,K,QTOT)
1713	RHSTS (K) = RHSTS (K) - TSNSR / DENOM
1714	END IF
1715	ELSE
1716	! TSNSR = SNKSRC(STC(K),SMC(K),SH2O(K),ZSOIL,NSOIL,
1717	IF ( (SICE > 0.) .OR. (STC (K) < T0) ) THEN
1718	CALL SNKSRC (TSNSR,STC (K),SMC (K),SH2O (K),ZSOIL,NSOIL, &
1719	SMCMAX,PSISAT,BEXP,DT,K,QTOT)
1720	RHSTS (K) = RHSTS (K) - TSNSR / DENOM
1721	END IF
1722	END IF
1723
1724	! ----------------------------------------------------------------------
1725	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
1726	! ----------------------------------------------------------------------
1727	AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
1728
1729	! ----------------------------------------------------------------------
1730	! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
1731	! ----------------------------------------------------------------------
1732	BI (K) = - (AI (K) + CI (K))
1733	TBK = TBK1
1734	DF1K = DF1N
1735	DTSDZ = DTSDZ2
1736	DDZ = DDZ2
1737	END DO
1738	! ----------------------------------------------------------------------
1739	END SUBROUTINE HRT
1740	! ----------------------------------------------------------------------
1741
1742	SUBROUTINE HRTICE (RHSTS,STC,TBOT,ICE,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
1743
1744	! ----------------------------------------------------------------------
1745	! SUBROUTINE HRTICE
1746	! ----------------------------------------------------------------------
1747	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1748	! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE (ICE=1) OR GLACIAL
1749	! ICE (ICE=-1). COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE
1750	! TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1751	!
1752	! (NOTE: THIS SUBROUTINE ONLY CALLED FOR SEA-ICE OR GLACIAL ICE, BUT
1753	! NOT FOR NON-GLACIAL LAND (ICE = 0).
1754	! ----------------------------------------------------------------------
1755	IMPLICIT NONE
1756
1757
1758	INTEGER, INTENT(IN) :: NSOIL
1759	INTEGER :: K
1760
1761	REAL, INTENT(IN) :: DF1,YY,ZZ1
1762	REAL, DIMENSION(1:NSOIL), INTENT(OUT):: AI, BI,CI
1763	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: STC, ZSOIL
1764	REAL, DIMENSION(1:NSOIL), INTENT(OUT):: RHSTS
1765	REAL, INTENT(IN) :: TBOT
1766	INTEGER, INTENT(IN) :: ICE
1767	REAL :: DDZ,DDZ2,DENOM,DTSDZ,DTSDZ2,SSOIL, &
1768	ZBOT
1769	REAL :: HCPCT
1770	REAL :: DF1K
1771	REAL :: DF1N
1772	REAL :: ZMD
1773
1774	! ----------------------------------------------------------------------
1775	! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
1776	! HCPCT = 1880.0*917.0.
1777	! ----------------------------------------------------------------------
1778	IF ( ICE == 1 ) THEN
1779	! Sea-ice values
1780	HCPCT = 1.72396E+6
1781	ELSEIF (ICE == -1) THEN
1782	! SET A NOMINAL UNIVERSAL VALUE OF GLACIAL-ICE SPECIFIC HEAT CAPACITY,
1783	! HCPCT = 2100.0*900.0 = 1.89000E+6 (SOURCE: BOB GRUMBINE, 2005)
1784	! TBOT PASSED IN AS ARGUMENT, VALUE FROM GLOBAL DATA SET
1785	!
1786	! A least-squares fit for the four points provided by
1787	! Keith Hines for the Yen (1981) values for Antarctic
1788	! snow firn.
1789	!
1790	HCPCT = 1.E6 * (0.8194 - 0.13090.5ZSOIL(1))
1791	DF1K = DF1
1792	ENDIF
1793
1794	! ----------------------------------------------------------------------
1795	! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
1796	! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
1797	! ----------------------------------------------------------------------
1798	! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
1799	! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
1800	! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
1801	! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
1802	! ----------------------------------------------------------------------
1803	! ----------------------------------------------------------------------
1804	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1805	! ----------------------------------------------------------------------
1806	IF (ICE == 1) THEN
1807	ZBOT = ZSOIL (NSOIL)
1808	ELSE IF (ICE == -1) THEN
1809	ZBOT = -25.0
1810	ENDIF
1811	DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
1812	AI (1) = 0.0
1813	CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
1814
1815	! ----------------------------------------------------------------------
1816	! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
1817	! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
1818	! RHSTS FOR THE TOP SOIL LAYER.
1819	! ----------------------------------------------------------------------
1820	BI (1) = - CI (1) + DF1/ (0.5 * ZSOIL (1) * ZSOIL (1) * HCPCT * &
1821	ZZ1)
1822	DTSDZ = ( STC (1) - STC (2) ) / ( -0.5 * ZSOIL (2) )
1823	SSOIL = DF1 * ( STC (1) - YY ) / ( 0.5 * ZSOIL (1) * ZZ1 )
1824
1825	! ----------------------------------------------------------------------
1826	! INITIALIZE DDZ2
1827	! ----------------------------------------------------------------------
1828	RHSTS (1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL (1) * HCPCT )
1829
1830	! ----------------------------------------------------------------------
1831	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1832	! ----------------------------------------------------------------------
1833	DDZ2 = 0.0
1834	DF1K = DF1
1835	DF1N = DF1
1836	DO K = 2,NSOIL
1837
1838	IF ( ICE == -1 ) THEN
1839	ZMD = 0.5 * (ZSOIL(K)+ZSOIL(K-1))
1840	! For the land-ice case
1841	! kmh 09/03/2006 use Yen (1981)'s values for Antarctic snow firn
1842	! IF ( K .eq. 2 ) HCPCT = 0.855108E6
1843	! IF ( K .eq. 3 ) HCPCT = 0.922906E6
1844	! IF ( K .eq. 4 ) HCPCT = 1.009986E6
1845
1846	! Least squares fit to the four points supplied by Keith Hines
1847	! from Yen (1981) for Antarctic snow firn. Not optimal, but
1848	! probably better than just a constant.
1849	HCPCT = 1.E6 * ( 0.8194 - 0.1309*ZMD )
1850
1851	! IF ( K .eq. 2 ) DF1N = 0.345356
1852	! IF ( K .eq. 3 ) DF1N = 0.398777
1853	! IF ( K .eq. 4 ) DF1N = 0.472653
1854
1855	! Least squares fit to the three points supplied by Keith Hines
1856	! from Yen (1981) for Antarctic snow firn. Not optimal, but
1857	! probably better than just a constant.
1858	DF1N = 0.32333 - ( 0.10073 * ZMD )
1859	ENDIF
1860	! ----------------------------------------------------------------------
1861	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
1862	! ----------------------------------------------------------------------
1863	IF (K /= NSOIL) THEN
1864	DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
1865
1866	! ----------------------------------------------------------------------
1867	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
1868	! ----------------------------------------------------------------------
1869	DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
1870	DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
1871	CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K))*HCPCT)
1872
1873	! ----------------------------------------------------------------------
1874	! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
1875	! ----------------------------------------------------------------------
1876	ELSE
1877
1878	! ----------------------------------------------------------------------
1879	! SET MATRIX COEF, CI TO ZERO.
1880	! ----------------------------------------------------------------------
1881	DTSDZ2 = (STC (K) - TBOT)/ (.5 * (ZSOIL (K -1) + ZSOIL (K)) &
1882	- ZBOT)
1883	CI (K) = 0.
1884	! ----------------------------------------------------------------------
1885	! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
1886	! ----------------------------------------------------------------------
1887	END IF
1888	DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
1889
1890	! ----------------------------------------------------------------------
1891	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
1892	! ----------------------------------------------------------------------
1893	RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
1894	AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
1895
1896	! ----------------------------------------------------------------------
1897	! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
1898	! ----------------------------------------------------------------------
1899	BI (K) = - (AI (K) + CI (K))
1900	DF1K = DF1N
1901	DTSDZ = DTSDZ2
1902	DDZ = DDZ2
1903	END DO
1904	! ----------------------------------------------------------------------
1905	END SUBROUTINE HRTICE
1906	! ----------------------------------------------------------------------
1907
1908	SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
1909
1910	! ----------------------------------------------------------------------
1911	! SUBROUTINE HSTEP
1912	! ----------------------------------------------------------------------
1913	! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
1914	! ----------------------------------------------------------------------
1915	IMPLICIT NONE
1916	INTEGER, INTENT(IN) :: NSOIL
1917	INTEGER :: K
1918
1919	REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
1920	REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
1921	REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
1922	REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
1923	REAL, DIMENSION(1:NSOIL) :: RHSTSin
1924	REAL, DIMENSION(1:NSOIL) :: CIin
1925	REAL :: DT
1926
1927	! ----------------------------------------------------------------------
1928	! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
1929	! ----------------------------------------------------------------------
1930	DO K = 1,NSOIL
1931	RHSTS (K) = RHSTS (K) * DT
1932	AI (K) = AI (K) * DT
1933	BI (K) = 1. + BI (K) * DT
1934	CI (K) = CI (K) * DT
1935	END DO
1936	! ----------------------------------------------------------------------
1937	! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
1938	! ----------------------------------------------------------------------
1939	DO K = 1,NSOIL
1940	RHSTSin (K) = RHSTS (K)
1941	END DO
1942	DO K = 1,NSOIL
1943	CIin (K) = CI (K)
1944	END DO
1945	! ----------------------------------------------------------------------
1946	! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
1947	! ----------------------------------------------------------------------
1948	CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
1949	! ----------------------------------------------------------------------
1950	! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
1951	! ----------------------------------------------------------------------
1952	DO K = 1,NSOIL
1953	STCOUT (K) = STCIN (K) + CI (K)
1954	END DO
1955	! ----------------------------------------------------------------------
1956	END SUBROUTINE HSTEP
1957	! ----------------------------------------------------------------------
1958
1959	SUBROUTINE NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
1960	SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
1961	SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
1962	SSOIL, &
1963	STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
1964	SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
1965	DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
1966	RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
1967	QUARTZ,FXEXP,CSOIL, &
1968	BETA,DRIP,DEW,FLX1,FLX3,VEGTYP,ISURBAN)
1969
1970	! ----------------------------------------------------------------------
1971	! SUBROUTINE NOPAC
1972	! ----------------------------------------------------------------------
1973	! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
1974	! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
1975	! PRESENT.
1976	! ----------------------------------------------------------------------
1977	IMPLICIT NONE
1978
1979	INTEGER, INTENT(IN) :: ICE, NROOT,NSOIL,VEGTYP
1980	INTEGER, INTENT(IN) :: ISURBAN
1981	INTEGER :: K
1982
1983	REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DKSAT,DT,DWSAT, &
1984	EPSCA,ETP,FDOWN,F1,FXEXP,FRZFACT,KDT,PC, &
1985	PRCP,PSISAT,Q2,QUARTZ,RCH,RR,SBETA,SFCTMP,&
1986	SHDFAC,SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, &
1987	T24,TBOT,TH2,ZBOT,EMISSI
1988	REAL, INTENT(INOUT) :: CMC,BETA,T1
1989	REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR,ETA,ETT,FLX1,FLX3, &
1990	RUNOFF1,RUNOFF2,RUNOFF3,SSOIL
1991	REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
1992	REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
1993	REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
1994	REAL, DIMENSION(1:NSOIL) :: ET1
1995	REAL :: EC1,EDIR1,ETT1,DF1,ETA1,ETP1,PRCP1,YY, &
1996	YYNUM,ZZ1
1997
1998	! ----------------------------------------------------------------------
1999	! EXECUTABLE CODE BEGINS HERE:
2000	! CONVERT ETP Fnd PRCP FROM KG M-2 S-1 TO M S-1 AND INITIALIZE DEW.
2001	! ----------------------------------------------------------------------
2002	PRCP1 = PRCP * 0.001
2003	ETP1 = ETP * 0.001
2004	DEW = 0.0
2005	! ----------------------------------------------------------------------
2006	! INITIALIZE EVAP TERMS.
2007	! ----------------------------------------------------------------------
2008	EDIR = 0.
2009	EDIR1 = 0.
2010	EC1 = 0.
2011	EC = 0.
2012	DO K = 1,NSOIL
2013	ET(K) = 0.
2014	ET1(K) = 0.
2015	END DO
2016	ETT = 0.
2017	ETT1 = 0.
2018
2019	IF (ETP > 0.0) THEN
2020	CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2021	SH2O, &
2022	SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2023	SMCREF,SHDFAC,CMCMAX, &
2024	SMCDRY,CFACTR, &
2025	EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2026	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2027	SH2O,SLOPE,KDT,FRZFACT, &
2028	SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2029	SHDFAC,CMCMAX, &
2030	RUNOFF1,RUNOFF2,RUNOFF3, &
2031	EDIR1,EC1,ET1, &
2032	DRIP)
2033
2034	! ----------------------------------------------------------------------
2035	! CONVERT MODELED EVAPOTRANSPIRATION FROM M S-1 TO KG M-2 S-1.
2036	! ----------------------------------------------------------------------
2037
2038	ETA = ETA1 * 1000.0
2039
2040	! ----------------------------------------------------------------------
2041	! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
2042	! ETP1 TO ZERO).
2043	! ----------------------------------------------------------------------
2044	ELSE
2045	DEW = - ETP1
2046
2047	! ----------------------------------------------------------------------
2048	! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
2049	! ----------------------------------------------------------------------
2050
2051	PRCP1 = PRCP1+ DEW
2052	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2053	SH2O,SLOPE,KDT,FRZFACT, &
2054	SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2055	SHDFAC,CMCMAX, &
2056	RUNOFF1,RUNOFF2,RUNOFF3, &
2057	EDIR1,EC1,ET1, &
2058	DRIP)
2059
2060	! ----------------------------------------------------------------------
2061	! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
2062	! ----------------------------------------------------------------------
2063	! ETA = ETA1 * 1000.0
2064	END IF
2065
2066	! ----------------------------------------------------------------------
2067	! BASED ON ETP AND E VALUES, DETERMINE BETA
2068	! ----------------------------------------------------------------------
2069
2070	IF ( ETP <= 0.0 ) THEN
2071	BETA = 0.0
2072	ETA = ETP
2073	IF ( ETP < 0.0 ) THEN
2074	BETA = 1.0
2075	END IF
2076	ELSE
2077	BETA = ETA / ETP
2078	END IF
2079
2080	! ----------------------------------------------------------------------
2081	! CONVERT MODELED EVAPOTRANSPIRATION COMPONENTS 'M S-1' TO 'KG M-2 S-1'.
2082	! ----------------------------------------------------------------------
2083	EDIR = EDIR1*1000.
2084	EC = EC1*1000.
2085	DO K = 1,NSOIL
2086	ET(K) = ET1(K)*1000.
2087	END DO
2088	ETT = ETT1*1000.
2089
2090	! ----------------------------------------------------------------------
2091	! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
2092	! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
2093	! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
2094	! ----------------------------------------------------------------------
2095
2096	CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
2097
2098	!urban
2099	IF ( VEGTYP == ISURBAN ) DF1=3.24
2100	!
2101
2102	! ----------------------------------------------------------------------
2103	! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
2104	! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
2105	! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
2106	! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
2107	! ROUTINE SFLX)
2108	! ----------------------------------------------------------------------
2109	DF1 = DF1 * EXP (SBETA * SHDFAC)
2110	! ----------------------------------------------------------------------
2111	! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
2112	! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
2113	! ----------------------------------------------------------------------
2114	YYNUM = FDOWN - EMISSISIGMA T24
2115	YY = SFCTMP + (YYNUM / RCH + TH2- SFCTMP - BETA * EPSCA) / RR
2116
2117	ZZ1 = DF1 / ( -0.5 * ZSOIL (1) * RCH * RR ) + 1.0
2118
2119	!urban
2120	CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2121	TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2122	QUARTZ,CSOIL,VEGTYP,ISURBAN)
2123
2124	! ----------------------------------------------------------------------
2125	! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
2126	! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
2127	! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
2128	! ----------------------------------------------------------------------
2129	FLX1 = CPH2O * PRCP * (T1- SFCTMP)
2130	FLX3 = 0.0
2131
2132	! ----------------------------------------------------------------------
2133	END SUBROUTINE NOPAC
2134	! ----------------------------------------------------------------------
2135
2136	SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
2137	& Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
2138	& DQSDT2,FLX2,EMISSI_IN,SNEQV,T1,ICE,SNCOVR)
2139
2140	! ----------------------------------------------------------------------
2141	! SUBROUTINE PENMAN
2142	! ----------------------------------------------------------------------
2143	! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
2144	! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
2145	! CALLING ROUTINE FOR LATER USE.
2146	! ----------------------------------------------------------------------
2147	IMPLICIT NONE
2148	LOGICAL, INTENT(IN) :: SNOWNG, FRZGRA
2149	REAL, INTENT(IN) :: CH, DQSDT2,FDOWN,PRCP, &
2150	Q2, Q2SAT,SSOIL, SFCPRS, SFCTMP, &
2151	T2V, TH2,EMISSI_IN,SNEQV
2152	REAL, INTENT(IN) :: T1 , SNCOVR
2153	!
2154	! kmh 09/13/2006
2155	INTEGER, INTENT(IN) :: ICE
2156	! kmh 09/03/2006
2157	!
2158	REAL, INTENT(OUT) :: EPSCA,ETP,FLX2,RCH,RR,T24
2159	REAL :: A, DELTA, FNET,RAD,RHO,EMISSI,ELCP1,LVS
2160
2161	REAL, PARAMETER :: ELCP = 2.4888E+3, LSUBC = 2.501000E+6,CP = 1004.6
2162	REAL, PARAMETER :: LSUBS = 2.83E+6
2163
2164	! ----------------------------------------------------------------------
2165	! EXECUTABLE CODE BEGINS HERE:
2166	! ----------------------------------------------------------------------
2167	! ----------------------------------------------------------------------
2168	! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
2169	! ----------------------------------------------------------------------
2170	EMISSI=EMISSI_IN
2171	IF (ICE==0) THEN
2172	ELCP1 = (1.0-SNCOVR)ELCP + SNCOVRELCP*LSUBS/LSUBC
2173	LVS = (1.0-SNCOVR)LSUBC + SNCOVRLSUBS
2174	ELSE
2175	IF ( T1 > 273.15 ) THEN
2176	ELCP1=ELCP
2177	LVS=LSUBC
2178	ELSE
2179	ELCP1 = ELCP*LSUBS/LSUBC
2180	LVS = LSUBS
2181	ENDIF
2182	ENDIF
2183
2184	FLX2 = 0.0
2185	! DELTA = ELCP * DQSDT2
2186	DELTA = ELCP1 * DQSDT2
2187	T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
2188	! RR = T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
2189	RR = EMISSIT24 6.48E-8 / (SFCPRS * CH) + 1.0
2190	RHO = SFCPRS / (RD * T2V)
2191
2192	! ----------------------------------------------------------------------
2193	! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
2194	! EFFECTS CAUSED BY FALLING PRECIPITATION.
2195	! ----------------------------------------------------------------------
2196	RCH = RHO * CP * CH
2197	IF (.NOT. SNOWNG) THEN
2198	IF (PRCP > 0.0) RR = RR + CPH2O * PRCP / RCH
2199	ELSE
2200	RR = RR + CPICE * PRCP / RCH
2201	END IF
2202
2203	! ----------------------------------------------------------------------
2204	! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
2205	! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
2206	! ----------------------------------------------------------------------
2207	! FNET = FDOWN - SIGMA * T24- SSOIL
2208	FNET = FDOWN - EMISSISIGMA T24- SSOIL
2209	IF (FRZGRA) THEN
2210	FLX2 = - LSUBF * PRCP
2211	FNET = FNET - FLX2
2212	! ----------------------------------------------------------------------
2213	! FINISH PENMAN EQUATION CALCULATIONS.
2214	! ----------------------------------------------------------------------
2215	END IF
2216	RAD = FNET / RCH + TH2- SFCTMP
2217	! A = ELCP * (Q2SAT - Q2)
2218	A = ELCP1 * (Q2SAT - Q2)
2219	EPSCA = (A * RR + RAD * DELTA) / (DELTA + RR)
2220	! ETP = EPSCA * RCH / LSUBC
2221	ETP = EPSCA * RCH / LVS
2222
2223	! ----------------------------------------------------------------------
2224	END SUBROUTINE PENMAN
2225	! ----------------------------------------------------------------------
2226
2227	SUBROUTINE REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX, &
2228	TOPT, &
2229	REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
2230	PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
2231	SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
2232	RTDIS,SLDPTH,ZSOIL, NROOT,NSOIL,CZIL, &
2233	LAIMIN, LAIMAX, EMISSMIN, EMISSMAX, ALBEDOMIN, &
2234	ALBEDOMAX, Z0MIN, Z0MAX, CSOIL, PTU, LLANDUSE, &
2235	LSOIL, LOCAL,LVCOEF)
2236
2237	IMPLICIT NONE
2238	! ----------------------------------------------------------------------
2239	! Internally set (default valuess)
2240	! all soil and vegetation parameters required for the execusion oF
2241	! the Noah lsm are defined in VEGPARM.TBL, SOILPARM.TB, and GENPARM.TBL.
2242	! ----------------------------------------------------------------------
2243	! Vegetation parameters:
2244	! ALBBRD: SFC background snow-free albedo
2245	! CMXTBL: MAX CNPY Capacity
2246	! Z0BRD: Background roughness length
2247	! SHDFAC: Green vegetation fraction
2248	! NROOT: Rooting depth
2249	! RSMIN: Mimimum stomatal resistance
2250	! RSMAX: Max. stomatal resistance
2251	! RGL: Parameters used in radiation stress function
2252	! HS: Parameter used in vapor pressure deficit functio
2253	! TOPT: Optimum transpiration air temperature.
2254	! CMCMAX: Maximum canopy water capacity
2255	! CFACTR: Parameter used in the canopy inteception calculation
2256	! SNUP: Threshold snow depth (in water equivalent m) that
2257	! implies 100 percent snow cover
2258	! LAI: Leaf area index
2259	!
2260	! ----------------------------------------------------------------------
2261	! Soil parameters:
2262	! SMCMAX: MAX soil moisture content (porosity)
2263	! SMCREF: Reference soil moisture (field capacity)
2264	! SMCWLT: Wilting point soil moisture
2265	! SMCWLT: Air dry soil moist content limits
2266	! SSATPSI: SAT (saturation) soil potential
2267	! DKSAT: SAT soil conductivity
2268	! BEXP: B parameter
2269	! SSATDW: SAT soil diffusivity
2270	! F1: Soil thermal diffusivity/conductivity coef.
2271	! QUARTZ: Soil quartz content
2272	! Modified by F. Chen (12/22/97) to use the STATSGO soil map
2273	! Modified By F. Chen (01/22/00) to include PLaya, Lava, and White San
2274	! Modified By F. Chen (08/05/02) to include additional parameters for the Noah
2275	! NOTE: SATDW = BBSATDK(SATPSI/MAXSMC)
2276	! F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0
2277	! REFSMC1=MAXSMC(5.79E-9/SATDK)(1/(2BB+3)) 5.79E-9 m/s= 0.5 mm
2278	! REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1)
2279	! WLTSMC1=MAXSMC(200./SATPSI)*(-1./BB) (Wetzel and Chang, 198
2280	! WLTSMC=WLTSMC1-0.5*WLTSMC1
2281	! Note: the values for playa is set for it to have a thermal conductivit
2282	! as sand and to have a hydrulic conductivity as clay
2283	!
2284	! ----------------------------------------------------------------------
2285	! Class parameter 'SLOPETYP' was included to estimate linear reservoir
2286	! coefficient 'SLOPE' to the baseflow runoff out of the bottom layer.
2287	! lowest class (slopetyp=0) means highest slope parameter = 1.
2288	! definition of slopetyp from 'zobler' slope type:
2289	! slope class percent slope
2290	! 1 0-8
2291	! 2 8-30
2292	! 3 > 30
2293	! 4 0-30
2294	! 5 0-8 & > 30
2295	! 6 8-30 & > 30
2296	! 7 0-8, 8-30, > 30
2297	! 9 GLACIAL ICE
2298	! BLANK OCEAN/SEA
2299	! SLOPE_DATA: linear reservoir coefficient
2300	! SBETA_DATA: parameter used to caluculate vegetation effect on soil heat
2301	! FXEXP_DAT: soil evaporation exponent used in DEVAP
2302	! CSOIL_DATA: soil heat capacity [J M-3 K-1]
2303	! SALP_DATA: shape parameter of distribution function of snow cover
2304	! REFDK_DATA and REFKDT_DATA: parameters in the surface runoff parameteriz
2305	! FRZK_DATA: frozen ground parameter
2306	! ZBOT_DATA: depth[M] of lower boundary soil temperature
2307	! CZIL_DATA: calculate roughness length of heat
2308	! SMLOW_DATA and MHIGH_DATA: two soil moisture wilt, soil moisture referen
2309	! parameters
2310	! Set maximum number of soil-, veg-, and slopetyp in data statement.
2311	! ----------------------------------------------------------------------
2312	INTEGER, PARAMETER :: MAX_SLOPETYP=30,MAX_SOILTYP=30,MAX_VEGTYP=30
2313	LOGICAL :: LOCAL
2314	CHARACTER (LEN=256), INTENT(IN):: LLANDUSE, LSOIL
2315
2316	! Veg parameters
2317	INTEGER, INTENT(IN) :: VEGTYP
2318	INTEGER, INTENT(OUT) :: NROOT
2319	REAL, INTENT(INOUT) :: SHDFAC
2320	REAL, INTENT(OUT) :: HS,RSMIN,RGL,SNUP, &
2321	CMCMAX,RSMAX,TOPT, &
2322	EMISSMIN, EMISSMAX, &
2323	LAIMIN, LAIMAX, &
2324	Z0MIN, Z0MAX, &
2325	ALBEDOMIN, ALBEDOMAX
2326	! Soil parameters
2327	INTEGER, INTENT(IN) :: SOILTYP
2328	REAL, INTENT(OUT) :: BEXP,DKSAT,DWSAT,F1,QUARTZ,SMCDRY, &
2329	SMCMAX,SMCREF,SMCWLT,PSISAT
2330	! General parameters
2331	INTEGER, INTENT(IN) :: SLOPETYP,NSOIL
2332	INTEGER :: I
2333
2334	REAL, INTENT(OUT) :: SLOPE,CZIL,SBETA,FXEXP, &
2335	CSOIL,SALP,FRZX,KDT,CFACTR, &
2336	ZBOT,REFKDT,PTU
2337	REAL, INTENT(OUT) :: LVCOEF
2338	REAL,DIMENSION(1:NSOIL),INTENT(IN) :: SLDPTH,ZSOIL
2339	REAL,DIMENSION(1:NSOIL),INTENT(OUT):: RTDIS
2340	REAL :: FRZFACT,FRZK,REFDK
2341
2342	! SAVE
2343	! ----------------------------------------------------------------------
2344	!
2345	IF (SOILTYP .gt. SLCATS) THEN
2346	CALL wrf_error_fatal ( 'Warning: too many input soil types' )
2347	END IF
2348	IF (VEGTYP .gt. LUCATS) THEN
2349	CALL wrf_error_fatal ( 'Warning: too many input landuse types' )
2350	END IF
2351	IF (SLOPETYP .gt. SLPCATS) THEN
2352	CALL wrf_error_fatal ( 'Warning: too many input slope types' )
2353	END IF
2354
2355	! ----------------------------------------------------------------------
2356	! SET-UP SOIL PARAMETERS
2357	! ----------------------------------------------------------------------
2358	CSOIL = CSOIL_DATA
2359	BEXP = BB (SOILTYP)
2360	DKSAT = SATDK (SOILTYP)
2361	DWSAT = SATDW (SOILTYP)
2362	F1 = F11 (SOILTYP)
2363	PSISAT = SATPSI (SOILTYP)
2364	QUARTZ = QTZ (SOILTYP)
2365	SMCDRY = DRYSMC (SOILTYP)
2366	SMCMAX = MAXSMC (SOILTYP)
2367	SMCREF = REFSMC (SOILTYP)
2368	SMCWLT = WLTSMC (SOILTYP)
2369	! ----------------------------------------------------------------------
2370	! Set-up universal parameters (not dependent on SOILTYP, VEGTYP or
2371	! SLOPETYP)
2372	! ----------------------------------------------------------------------
2373	ZBOT = ZBOT_DATA
2374	SALP = SALP_DATA
2375	SBETA = SBETA_DATA
2376	REFDK = REFDK_DATA
2377	FRZK = FRZK_DATA
2378	FXEXP = FXEXP_DATA
2379	REFKDT = REFKDT_DATA
2380	PTU = 0. ! (not used yet) to satisify intent(out)
2381	KDT = REFKDT * DKSAT / REFDK
2382	CZIL = CZIL_DATA
2383	SLOPE = SLOPE_DATA (SLOPETYP)
2384	LVCOEF = LVCOEF_DATA
2385
2386	! ----------------------------------------------------------------------
2387	! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
2388	! ----------------------------------------------------------------------
2389	FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
2390	FRZX = FRZK * FRZFACT
2391
2392	! ----------------------------------------------------------------------
2393	! SET-UP VEGETATION PARAMETERS
2394	! ----------------------------------------------------------------------
2395	TOPT = TOPT_DATA
2396	CMCMAX = CMCMAX_DATA
2397	CFACTR = CFACTR_DATA
2398	RSMAX = RSMAX_DATA
2399	NROOT = NROTBL (VEGTYP)
2400	SNUP = SNUPTBL (VEGTYP)
2401	RSMIN = RSTBL (VEGTYP)
2402	RGL = RGLTBL (VEGTYP)
2403	HS = HSTBL (VEGTYP)
2404	EMISSMIN = EMISSMINTBL (VEGTYP)
2405	EMISSMAX = EMISSMAXTBL (VEGTYP)
2406	LAIMIN = LAIMINTBL (VEGTYP)
2407	LAIMAX = LAIMAXTBL (VEGTYP)
2408	Z0MIN = Z0MINTBL (VEGTYP)
2409	Z0MAX = Z0MAXTBL (VEGTYP)
2410	ALBEDOMIN = ALBEDOMINTBL (VEGTYP)
2411	ALBEDOMAX = ALBEDOMAXTBL (VEGTYP)
2412
2413	IF (VEGTYP .eq. BARE) SHDFAC = 0.0
2414	IF (NROOT .gt. NSOIL) THEN
2415	WRITE (err_message,*) 'Error: too many root layers ', &
2416	NSOIL,NROOT
2417	CALL wrf_error_fatal ( err_message )
2418	! ----------------------------------------------------------------------
2419	! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
2420	! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
2421	! ----------------------------------------------------------------------
2422	END IF
2423	DO I = 1,NROOT
2424	RTDIS (I) = - SLDPTH (I)/ ZSOIL (NROOT)
2425	! ----------------------------------------------------------------------
2426	! SET-UP SLOPE PARAMETER
2427	! ----------------------------------------------------------------------
2428	END DO
2429
2430	! print*,'end of PRMRED'
2431	! print*,'VEGTYP',VEGTYP,'SOILTYP',SOILTYP,'SLOPETYP',SLOPETYP, &
2432	! & 'CFACTR',CFACTR,'CMCMAX',CMCMAX,'RSMAX',RSMAX,'TOPT',TOPT, &
2433	! & 'REFKDT',REFKDT,'KDT',KDT,'SBETA',SBETA, 'SHDFAC',SHDFAC, &
2434	! & 'RSMIN',RSMIN,'RGL',RGL,'HS',HS,'ZBOT',ZBOT,'FRZX',FRZX, &
2435	! & 'PSISAT',PSISAT,'SLOPE',SLOPE,'SNUP',SNUP,'SALP',SALP,'BEXP', &
2436	! & BEXP, &
2437	! & 'DKSAT',DKSAT,'DWSAT',DWSAT, &
2438	! & 'SMCMAX',SMCMAX,'SMCWLT',SMCWLT,'SMCREF',SMCREF,'SMCDRY',SMCDRY, &
2439	! & 'F1',F1,'QUARTZ',QUARTZ,'FXEXP',FXEXP, &
2440	! & 'RTDIS',RTDIS,'SLDPTH',SLDPTH,'ZSOIL',ZSOIL, 'NROOT',NROOT, &
2441	! & 'NSOIL',NSOIL,'Z0',Z0,'CZIL',CZIL,'LAI',LAI, &
2442	! & 'CSOIL',CSOIL,'PTU',PTU, &
2443	! & 'LOCAL', LOCAL
2444
2445	END SUBROUTINE REDPRM
2446
2447	SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
2448
2449	! ----------------------------------------------------------------------
2450	! SUBROUTINE ROSR12
2451	! ----------------------------------------------------------------------
2452	! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
2453	! ### ### ### ### ### ###
2454	! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
2455	! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
2456	! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
2457	! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
2458	! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
2459	! # . . # # . # = # . #
2460	! # . . # # . # # . #
2461	! # . . # # . # # . #
2462	! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
2463	! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
2464	! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
2465	! ### ### ### ### ### ###
2466	! ----------------------------------------------------------------------
2467	IMPLICIT NONE
2468
2469	INTEGER, INTENT(IN) :: NSOIL
2470	INTEGER :: K, KK
2471
2472	REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
2473	REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
2474
2475	! ----------------------------------------------------------------------
2476	! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
2477	! ----------------------------------------------------------------------
2478	C (NSOIL) = 0.0
2479	P (1) = - C (1) / B (1)
2480	! ----------------------------------------------------------------------
2481	! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
2482	! ----------------------------------------------------------------------
2483
2484	! ----------------------------------------------------------------------
2485	! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
2486	! ----------------------------------------------------------------------
2487	DELTA (1) = D (1) / B (1)
2488	DO K = 2,NSOIL
2489	P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
2490	DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
2491	* P (K -1)))
2492	END DO
2493	! ----------------------------------------------------------------------
2494	! SET P TO DELTA FOR LOWEST SOIL LAYER
2495	! ----------------------------------------------------------------------
2496	P (NSOIL) = DELTA (NSOIL)
2497
2498	! ----------------------------------------------------------------------
2499	! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
2500	! ----------------------------------------------------------------------
2501	DO K = 2,NSOIL
2502	KK = NSOIL - K + 1
2503	P (KK) = P (KK) * P (KK +1) + DELTA (KK)
2504	END DO
2505	! ----------------------------------------------------------------------
2506	END SUBROUTINE ROSR12
2507	! ----------------------------------------------------------------------
2508
2509
2510	SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2511	TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2512	QUARTZ,CSOIL,VEGTYP,ISURBAN)
2513
2514	! ----------------------------------------------------------------------
2515	! SUBROUTINE SHFLX
2516	! ----------------------------------------------------------------------
2517	! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
2518	! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
2519	! ON THE TEMPERATURE.
2520	! ----------------------------------------------------------------------
2521	IMPLICIT NONE
2522
2523	INTEGER, INTENT(IN) :: ICE, NSOIL, VEGTYP, ISURBAN
2524	INTEGER :: I
2525
2526	REAL, INTENT(IN) :: BEXP,CSOIL,DF1,DT,F1,PSISAT,QUARTZ, &
2527	SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1
2528	REAL, INTENT(INOUT) :: T1
2529	REAL, INTENT(OUT) :: SSOIL
2530	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
2531	REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O
2532	REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
2533	REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS
2534	REAL, PARAMETER :: T0 = 273.15
2535
2536	! ----------------------------------------------------------------------
2537	! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
2538	! ----------------------------------------------------------------------
2539
2540	! ----------------------------------------------------------------------
2541	! SEA-ICE CASE, GLACIAL ICE CASE
2542	! ----------------------------------------------------------------------
2543	IF ( ICE /= 0 ) THEN
2544
2545	CALL HRTICE (RHSTS,STC,TBOT,ICE,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2546
2547	CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
2548
2549	! ----------------------------------------------------------------------
2550	! LAND-MASS CASE
2551	! ----------------------------------------------------------------------
2552	ELSE
2553	CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
2554	ZBOT,PSISAT,SH2O,DT, &
2555	BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP,ISURBAN)
2556
2557	CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
2558	END IF
2559	DO I = 1,NSOIL
2560	STC (I) = STCF (I)
2561	END DO
2562
2563	! ----------------------------------------------------------------------
2564	! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
2565	! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
2566	! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
2567	! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
2568	! DIFFERENTLY IN ROUTINE SNOPAC)
2569	! ----------------------------------------------------------------------
2570	! ----------------------------------------------------------------------
2571	! CALCULATE SURFACE SOIL HEAT FLUX
2572	! ----------------------------------------------------------------------
2573	T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
2574	SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
2575
2576	! ----------------------------------------------------------------------
2577	END SUBROUTINE SHFLX
2578	! ----------------------------------------------------------------------
2579
2580	SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2581	& SH2O,SLOPE,KDT,FRZFACT, &
2582	& SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2583	& SHDFAC,CMCMAX, &
2584	& RUNOFF1,RUNOFF2,RUNOFF3, &
2585	& EDIR,EC,ET, &
2586	& DRIP)
2587
2588	! ----------------------------------------------------------------------
2589	! SUBROUTINE SMFLX
2590	! ----------------------------------------------------------------------
2591	! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
2592	! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
2593	! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
2594	! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
2595	! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
2596	! ----------------------------------------------------------------------
2597	IMPLICIT NONE
2598
2599	INTEGER, INTENT(IN) :: NSOIL
2600	INTEGER :: I,K
2601
2602	REAL, INTENT(IN) :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR, &
2603	KDT, PRCP1, SHDFAC, SLOPE, SMCMAX, SMCWLT
2604	REAL, INTENT(OUT) :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
2605	REAL, INTENT(INOUT) :: CMC
2606	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET,ZSOIL
2607	REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SMC, SH2O
2608	REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS, RHSTT, &
2609	SICE, SH2OA, SH2OFG
2610	REAL :: DUMMY, EXCESS,FRZFACT,PCPDRP,RHSCT,TRHSCT
2611	REAL :: FAC2
2612	REAL :: FLIMIT
2613
2614	! ----------------------------------------------------------------------
2615	! EXECUTABLE CODE BEGINS HERE.
2616	! ----------------------------------------------------------------------
2617	! ----------------------------------------------------------------------
2618	! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
2619	! ----------------------------------------------------------------------
2620	DUMMY = 0.
2621
2622	! ----------------------------------------------------------------------
2623	! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
2624	! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
2625	! FALL TO THE GRND.
2626	! ----------------------------------------------------------------------
2627	RHSCT = SHDFAC * PRCP1- EC
2628	DRIP = 0.
2629	TRHSCT = DT * RHSCT
2630	EXCESS = CMC + TRHSCT
2631
2632	! ----------------------------------------------------------------------
2633	! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
2634	! SOIL
2635	! ----------------------------------------------------------------------
2636	IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
2637	PCPDRP = (1. - SHDFAC) * PRCP1+ DRIP / DT
2638
2639	! ----------------------------------------------------------------------
2640	! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT and SSTEP
2641	!
2642	DO I = 1,NSOIL
2643	SICE (I) = SMC (I) - SH2O (I)
2644	END DO
2645	! ----------------------------------------------------------------------
2646	! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
2647	! TENDENCY EQUATIONS.
2648	! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
2649	! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
2650	! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
2651	! THE FIRST SOIL LAYER)
2652	! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
2653	! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
2654	! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
2655	! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
2656	! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
2657	! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
2658	! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
2659	! SOIL MOISTURE STATE
2660	! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
2661	! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
2662	! OF SECTION 2 OF KALNAY AND KANAMITSU
2663	! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
2664	! ----------------------------------------------------------------------
2665	! According to Dr. Ken Mitchell's suggestion, add the second contraint
2666	! to remove numerical instability of runoff and soil moisture
2667	! FLIMIT is a limit value for FAC2
2668	FAC2=0.0
2669	DO I=1,NSOIL
2670	FAC2=MAX(FAC2,SH2O(I)/SMCMAX)
2671	ENDDO
2672	CALL FAC2MIT(SMCMAX,FLIMIT)
2673
2674	! ----------------------------------------------------------------------
2675	! FROZEN GROUND VERSION:
2676	! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
2677	! INC&UDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
2678	! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
2679	! ----------------------------------------------------------------------
2680	IF ( ( (PCPDRP * DT) > (0.00011000.0 (- ZSOIL (1))* SMCMAX) ) &
2681	.OR. (FAC2 > FLIMIT) ) THEN
2682	CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
2683	DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
2684	RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2685	CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
2686	CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2687	DO K = 1,NSOIL
2688	SH2OA (K) = (SH2O (K) + SH2OFG (K)) * 0.5
2689	END DO
2690	CALL SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
2691	DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
2692	RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2693	CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
2694	CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2695
2696	ELSE
2697	CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
2698	DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
2699	RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2700	CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
2701	CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2702	! RUNOF = RUNOFF
2703
2704	END IF
2705
2706	! ----------------------------------------------------------------------
2707	END SUBROUTINE SMFLX
2708	! ----------------------------------------------------------------------
2709
2710
2711	SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
2712
2713	! ----------------------------------------------------------------------
2714	! SUBROUTINE SNFRAC
2715	! ----------------------------------------------------------------------
2716	! CALCULATE SNOW FRACTION (0 -> 1)
2717	! SNEQV SNOW WATER EQUIVALENT (M)
2718	! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
2719	! SALP TUNING PARAMETER
2720	! SNCOVR FRACTIONAL SNOW COVER
2721	! ----------------------------------------------------------------------
2722	IMPLICIT NONE
2723
2724	REAL, INTENT(IN) :: SNEQV,SNUP,SALP,SNOWH
2725	REAL, INTENT(OUT) :: SNCOVR
2726	REAL :: RSNOW, Z0N
2727
2728	! ----------------------------------------------------------------------
2729	! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
2730	! REDPRM) ABOVE WHICH SNOCVR=1.
2731	! ----------------------------------------------------------------------
2732	IF (SNEQV < SNUP) THEN
2733	RSNOW = SNEQV / SNUP
2734	SNCOVR = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
2735	ELSE
2736	SNCOVR = 1.0
2737	END IF
2738
2739	! FORMULATION OF DICKINSON ET AL. 1986
2740	! Z0N = 0.035
2741
2742	! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
2743
2744	! FORMULATION OF MARSHALL ET AL. 1994
2745	! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
2746
2747	! ----------------------------------------------------------------------
2748	END SUBROUTINE SNFRAC
2749	! ----------------------------------------------------------------------
2750
2751	SUBROUTINE SNKSRC (TSNSR,TAVG,SMC,SH2O,ZSOIL,NSOIL, &
2752	& SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2753	! ----------------------------------------------------------------------
2754	! SUBROUTINE SNKSRC
2755	! ----------------------------------------------------------------------
2756	! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
2757	! AVAILABLE LIQUED WATER.
2758	! ----------------------------------------------------------------------
2759	IMPLICIT NONE
2760
2761	INTEGER, INTENT(IN) :: K,NSOIL
2762	REAL, INTENT(IN) :: BEXP, DT, PSISAT, QTOT, SMC, SMCMAX, &
2763	TAVG
2764	REAL, INTENT(INOUT) :: SH2O
2765
2766	REAL, DIMENSION(1:NSOIL), INTENT(IN):: ZSOIL
2767
2768	REAL :: DF, DZ, DZH, FREE, TSNSR, &
2769	TDN, TM, TUP, TZ, X0, XDN, XH2O, XUP
2770
2771	REAL, PARAMETER :: DH2O = 1.0000E3, HLICE = 3.3350E5, &
2772	T0 = 2.7315E2
2773
2774	IF (K == 1) THEN
2775	DZ = - ZSOIL (1)
2776	ELSE
2777	DZ = ZSOIL (K -1) - ZSOIL (K)
2778	END IF
2779	! ----------------------------------------------------------------------
2780	! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
2781	! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
2782	! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
2783	! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
2784	! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
2785	! ----------------------------------------------------------------------
2786	! FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
2787
2788	! ----------------------------------------------------------------------
2789	! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
2790	! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
2791	! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
2792	! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
2793	! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
2794	! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
2795	! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
2796	! ----------------------------------------------------------------------
2797	CALL FRH2O (FREE,TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
2798
2799	! ----------------------------------------------------------------------
2800	! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
2801	! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
2802	! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
2803	! ----------------------------------------------------------------------
2804	XH2O = SH2O + QTOT * DT / (DH2O * HLICE * DZ)
2805	IF ( XH2O < SH2O .AND. XH2O < FREE) THEN
2806	IF ( FREE > SH2O ) THEN
2807	XH2O = SH2O
2808	ELSE
2809	XH2O = FREE
2810	END IF
2811	END IF
2812	! ----------------------------------------------------------------------
2813	! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
2814	! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
2815	! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
2816	! ----------------------------------------------------------------------
2817	IF ( XH2O > SH2O .AND. XH2O > FREE ) THEN
2818	IF ( FREE < SH2O ) THEN
2819	XH2O = SH2O
2820	ELSE
2821	XH2O = FREE
2822	END IF
2823	END IF
2824
2825	! ----------------------------------------------------------------------
2826	! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
2827	! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
2828	! ----------------------------------------------------------------------
2829	! SNKSRC = -DH2OHLICEDZ*(XH2O-SH2O)/DT
2830	IF (XH2O < 0.) XH2O = 0.
2831	IF (XH2O > SMC) XH2O = SMC
2832	TSNSR = - DH2O * HLICE * DZ * (XH2O - SH2O)/ DT
2833	SH2O = XH2O
2834
2835	! ----------------------------------------------------------------------
2836	END SUBROUTINE SNKSRC
2837	! ----------------------------------------------------------------------
2838
2839	SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
2840	SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
2841	SBETA,DF1, &
2842	Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,&
2843	SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
2844	SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT, &
2845	ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
2846	RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
2847	ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
2848	BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI,&
2849	RIBB,SOLDN, &
2850	ISURBAN, &
2851
2852	VEGTYP)
2853
2854	! ----------------------------------------------------------------------
2855	! SUBROUTINE SNOPAC
2856	! ----------------------------------------------------------------------
2857	! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
2858	! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
2859	! PRESENT.
2860	! ----------------------------------------------------------------------
2861	IMPLICIT NONE
2862
2863	INTEGER, INTENT(IN) :: ICE, NROOT, NSOIL,VEGTYP
2864	INTEGER, INTENT(IN) :: ISURBAN
2865	INTEGER :: K
2866	!
2867	! kmh 09/03/2006 add IT16 for surface temperature iteration
2868	!
2869	INTEGER :: IT16
2870	LOGICAL, INTENT(IN) :: SNOWNG
2871	REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DF1,DKSAT, &
2872	DT,DWSAT, EPSCA,FDOWN,F1,FXEXP, &
2873	FRZFACT,KDT,PC, PRCP,PSISAT,Q2,QUARTZ, &
2874	RCH,RR,SBETA,SFCPRS, SFCTMP, SHDFAC, &
2875	SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, T24, &
2876	TBOT,TH2,ZBOT,EMISSI,SOLDN
2877	REAL, INTENT(INOUT) :: CMC, BETA, ESD,FLX2,PRCPF,SNOWH,SNCOVR, &
2878	SNDENS, T1, RIBB, ETP
2879	REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR, ETNS, ESNOW,ETT, &
2880	FLX1,FLX3, RUNOFF1,RUNOFF2,RUNOFF3, &
2881	SSOIL,SNOMLT
2882	REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
2883	REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
2884	REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
2885	REAL, DIMENSION(1:NSOIL) :: ET1
2886	REAL :: DENOM,DSOIL,DTOT,EC1,EDIR1,ESDFLX,ETA, &
2887	ETT1, ESNOW1, ESNOW2, ETA1,ETP1,ETP2, &
2888	ETP3, ETNS1, ETANRG, ETAX, EX, FLX3X, &
2889	FRCSNO,FRCSOI, PRCP1, QSAT,RSNOW, SEH, &
2890	SNCOND,SSOIL1, T11,T12, T12A, T12AX, &
2891	T12B, T14, YY, ZZ1
2892	! T12B, T14, YY, ZZ1,EMISSI_S
2893	!
2894	! kmh 01/11/2007 add T15, T16, and DTOT2 for SFC T iteration and snow heat flux
2895	!
2896	REAL :: T15, T16, DTOT2
2897	REAL, PARAMETER :: ESDMIN = 1.E-6, LSUBC = 2.501000E+6, &
2898	LSUBS = 2.83E+6, TFREEZ = 273.15, &
2899	SNOEXP = 2.0
2900
2901	! ----------------------------------------------------------------------
2902	! EXECUTABLE CODE BEGINS HERE:
2903	! ----------------------------------------------------------------------
2904	! IF SEA-ICE (ICE=1) OR GLACIAL-ICE (ICE=-1), SNOWCOVER FRACTION = 1.0,
2905	! AND SUBLIMATION IS AT THE POTENTIAL RATE.
2906	! FOR NON-GLACIAL LAND (ICE=0), IF SNOWCOVER FRACTION < 1.0, TOTAL
2907	! EVAPORATION < POTENTIAL DUE TO NON-POTENTIAL CONTRIBUTION FROM
2908	! NON-SNOW COVERED FRACTION.
2909	! ----------------------------------------------------------------------
2910	! INITIALIZE EVAP TERMS.
2911	! ----------------------------------------------------------------------
2912	! conversions:
2913	! ESNOW [KG M-2 S-1]
2914	! ESDFLX [KG M-2 S-1] .le. ESNOW
2915	! ESNOW1 [M S-1]
2916	! ESNOW2 [M]
2917	! ETP [KG M-2 S-1]
2918	! ETP1 [M S-1]
2919	! ETP2 [M]
2920	! ----------------------------------------------------------------------
2921	DEW = 0.
2922	EDIR = 0.
2923	EDIR1 = 0.
2924	EC1 = 0.
2925	EC = 0.
2926	! EMISSI_S=0.95 ! For snow
2927
2928	DO K = 1,NSOIL
2929	ET (K) = 0.
2930	ET1 (K) = 0.
2931	END DO
2932	ETT = 0.
2933	ETT1 = 0.
2934	ETNS = 0.
2935	ETNS1 = 0.
2936	ESNOW = 0.
2937	ESNOW1 = 0.
2938	ESNOW2 = 0.
2939
2940	! ----------------------------------------------------------------------
2941	! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO ETP1 IN M S-1
2942	! ----------------------------------------------------------------------
2943	PRCP1 = PRCPF *0.001
2944	! ----------------------------------------------------------------------
2945	! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
2946	! ----------------------------------------------------------------------
2947	BETA = 1.0
2948	IF (ETP <= 0.0) THEN
2949	IF ( ( RIBB >= 0.1 ) .AND. ( FDOWN > 150.0 ) ) THEN
2950	ETP=(MIN(ETP(1.0-RIBB),0.)SNCOVR/0.980 + ETP*(0.980-SNCOVR))/0.980
2951	ENDIF
2952	IF(ETP == 0.) BETA = 0.0
2953	ETP1 = ETP * 0.001
2954	DEW = -ETP1
2955	ESNOW2 = ETP1*DT
2956	ETANRG = ETP((1.-SNCOVR)LSUBC + SNCOVR*LSUBS)
2957	ELSE
2958	ETP1 = ETP * 0.001
2959	IF ( ICE /= 0 ) THEN
2960	! SEA-ICE AND GLACIAL-ICE CASE
2961	ESNOW = ETP
2962	ESNOW1 = ESNOW*0.001
2963	ESNOW2 = ESNOW1*DT
2964	ETANRG = ESNOW*LSUBS
2965	ELSE IF ( ICE == 0) THEN
2966	! NON-GLACIAL LAND CASE
2967	IF (SNCOVR < 1.) THEN
2968	CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2969	SH2O, &
2970	SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2971	SMCREF,SHDFAC,CMCMAX, &
2972	SMCDRY,CFACTR, &
2973	EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS, &
2974	FXEXP)
2975	! ----------------------------------------------------------------------------
2976	EDIR1 = EDIR1* (1. - SNCOVR)
2977	EC1 = EC1* (1. - SNCOVR)
2978	DO K = 1,NSOIL
2979	ET1 (K) = ET1 (K)* (1. - SNCOVR)
2980	END DO
2981	ETT1 = ETT1*(1.-SNCOVR)
2982	! ETNS1 = EDIR1+ EC1+ ETT1
2983	ETNS1 = ETNS1*(1.-SNCOVR)
2984	! ----------------------------------------------------------------------------
2985	EDIR = EDIR1*1000.
2986	EC = EC1*1000.
2987	DO K = 1,NSOIL
2988	ET (K) = ET1 (K)*1000.
2989	END DO
2990	ETT = ETT1*1000.
2991	ETNS = ETNS1*1000.
2992	! ----------------------------------------------------------------------
2993
2994	! end IF (SNCOVR .lt. 1.)
2995	END IF
2996	! end IF (ICE .ne. 1)
2997	END IF
2998	ESNOW = ETP*SNCOVR
2999	ESNOW1 = ESNOW*0.001
3000	ESNOW2 = ESNOW1*DT
3001	ETANRG = ESNOWLSUBS + ETNSLSUBC
3002	! end IF (ETP .le. 0.0)
3003	END IF
3004
3005	! ----------------------------------------------------------------------
3006	! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3007	! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3008	! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3009	! SNOWFALL STRIKING THE GROUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3010	! ----------------------------------------------------------------------
3011	FLX1 = 0.0
3012	IF (SNOWNG) THEN
3013	FLX1 = CPICE * PRCP * (T1- SFCTMP)
3014	ELSE
3015	IF (PRCP > 0.0) FLX1 = CPH2O * PRCP * (T1- SFCTMP)
3016	! ----------------------------------------------------------------------
3017	! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3018	! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3019	! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3020	! FLUXES. FLX1 FROM ABOVE, FLX2 BROUGHT IN VIA COMMOM BLOCK RITE.
3021	! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
3022	! PENMAN.
3023	! ----------------------------------------------------------------------
3024	END IF
3025	DSOIL = - (0.5 * ZSOIL (1))
3026	DTOT = SNOWH + DSOIL
3027	DENOM = 1.0+ DF1 / (DTOT * RR * RCH)
3028	! surface emissivity weighted by snow cover fraction
3029	! T12A = ( (FDOWN - FLX1 - FLX2 - &
3030	! & ((SNCOVREMISSI_S)+EMISSI(1.0-SNCOVR))SIGMA T24)/RCH &
3031	! & + TH2 - SFCTMP - ETANRG/RCH ) / RR
3032	T12A = ( (FDOWN - FLX1- FLX2- EMISSI * SIGMA * T24)/ RCH &
3033	+ TH2- SFCTMP - ETANRG / RCH ) / RR
3034
3035	T12B = DF1 * STC (1) / (DTOT * RR * RCH)
3036
3037	! ----------------------------------------------------------------------
3038	! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
3039	! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
3040	! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
3041	! DEPENDING ON SIGN OF ETP.
3042	! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
3043	! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
3044	! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
3045	! TO ZERO.
3046	! ----------------------------------------------------------------------
3047	! SUB-FREEZING BLOCK
3048	! ----------------------------------------------------------------------
3049	T12 = (SFCTMP + T12A + T12B) / DENOM
3050	IF (T12 <= TFREEZ) THEN
3051	T1 = T12
3052	SSOIL = DF1 * (T1- STC (1)) / DTOT
3053	! ESD = MAX (0.0, ESD- ETP2)
3054	ESD = MAX(0.0, ESD-ESNOW2)
3055	FLX3 = 0.0
3056	EX = 0.0
3057
3058	SNOMLT = 0.0
3059	! ----------------------------------------------------------------------
3060	! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3061	! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
3062	! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
3063	! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
3064	! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
3065	! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
3066	! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
3067	! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
3068	! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
3069	! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
3070	! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
3071	! ----------------------------------------------------------------------
3072	! ABOVE FREEZING BLOCK
3073	! ----------------------------------------------------------------------
3074	ELSE
3075	T1 = TFREEZ * SNCOVR ** SNOEXP + T12 * (1.0- SNCOVR ** SNOEXP)
3076	BETA = 1.0
3077
3078	! ----------------------------------------------------------------------
3079	! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
3080	! BETA<1
3081	! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
3082	! ----------------------------------------------------------------------
3083	IF ( ICE /= 0 ) then
3084	! kmh 12/15/2005 modify SSOIL
3085	! kmh 09/03/2006 modify DTOT
3086	IF ( DTOT .GT. 2.0*DSOIL ) THEN
3087	DTOT = 2.0*DSOIL
3088	ENDIF
3089	ENDIF
3090	SSOIL = DF1 * (T1- STC (1)) / DTOT
3091	IF (ESD-ESNOW2 <= ESDMIN) THEN
3092	ESD = 0.0
3093	EX = 0.0
3094	SNOMLT = 0.0
3095	FLX3 = 0.0
3096	! ----------------------------------------------------------------------
3097	! SUBLIMATION LESS THAN DEPTH OF SNOWPACK
3098	! SNOWPACK (ESD) REDUCED BY ESNOW2 (DEPTH OF SUBLIMATED SNOW)
3099	! ----------------------------------------------------------------------
3100	ELSE
3101	ESD = ESD-ESNOW2
3102	ETP3 = ETP * LSUBC
3103	SEH = RCH * (T1- TH2)
3104	T14 = T1* T1
3105	T14 = T14* T14
3106	! FLX3 = FDOWN - FLX1 - FLX2 - &
3107	! ((SNCOVREMISSI_S)+EMISSI(1-SNCOVR))SIGMAT14 - &
3108	! SSOIL - SEH - ETANRG
3109	FLX3 = FDOWN - FLX1- FLX2- EMISSISIGMA T14- SSOIL - SEH - ETANRG
3110	IF (FLX3 <= 0.0) FLX3 = 0.0
3111	! ----------------------------------------------------------------------
3112	! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
3113	! ----------------------------------------------------------------------
3114	EX = FLX3*0.001/ LSUBF
3115
3116	! ----------------------------------------------------------------------
3117	! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
3118	! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
3119	! ----------------------------------------------------------------------
3120	SNOMLT = EX * DT
3121	IF (ESD- SNOMLT >= ESDMIN) THEN
3122	ESD = ESD- SNOMLT
3123	! ----------------------------------------------------------------------
3124	! SNOWMELT EXCEEDS SNOW DEPTH
3125	! ----------------------------------------------------------------------
3126	ELSE
3127	EX = ESD / DT
3128	FLX3 = EX 1000.0 LSUBF
3129	SNOMLT = ESD
3130
3131	ESD = 0.0
3132	! ----------------------------------------------------------------------
3133	! END OF 'ESD .LE. ETP2' IF-BLOCK
3134	! ----------------------------------------------------------------------
3135	END IF
3136	END IF
3137
3138	! ----------------------------------------------------------------------
3139	! END OF 'T12 .LE. TFREEZ' IF-BLOCK
3140	! ----------------------------------------------------------------------
3141	! ----------------------------------------------------------------------
3142	! IF NON-GLACIAL LAND, ADD SNOWMELT RATE (EX) TO PRECIP RATE TO BE USED
3143	! IN SUBROUTINE SMFLX (SOIL MOISTURE EVOLUTION) VIA INFILTRATION.
3144	!
3145	! FOR SEA-ICE AND GLACIAL-ICE, THE SNOWMELT WILL BE ADDED TO SUBSURFACE
3146	! RUNOFF/BASEFLOW LATER NEAR THE END OF SFLX (AFTER RETURN FROM CALL TO
3147	! SUBROUTINE SNOPAC)
3148	! ----------------------------------------------------------------------
3149	IF (ICE == 0) PRCP1 = PRCP1+ EX
3150
3151	! ----------------------------------------------------------------------
3152	! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
3153	! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
3154	! (BELOW).
3155	! IF SEAICE (ICE==1) SKIP CALL TO SMFLX.
3156	! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES FOR NON-GLACIAL LAND.
3157	! IF SEA-ICE (ICE==1) OR GLACIAL-ICE (ICE==-1), SKIP CALL TO SMFLX,
3158	! SINCE NO SOIL MEDIUM FOR SEA-ICE OR GLACIAL-ICE.
3159	! ----------------------------------------------------------------------
3160	END IF
3161	IF (ICE == 0) THEN
3162	CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3163	SH2O,SLOPE,KDT,FRZFACT, &
3164	SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3165	SHDFAC,CMCMAX, &
3166	RUNOFF1,RUNOFF2,RUNOFF3, &
3167	EDIR1,EC1,ET1, &
3168	DRIP)
3169	! ----------------------------------------------------------------------
3170	! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
3171	! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
3172	! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
3173	! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
3174	! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
3175	! SKIN TEMP VALUE AS REVISED BY SHFLX.
3176	! ----------------------------------------------------------------------
3177	END IF
3178	ZZ1 = 1.0
3179	YY = STC (1) -0.5* SSOIL * ZSOIL (1)* ZZ1/ DF1
3180
3181	! ----------------------------------------------------------------------
3182	! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
3183	! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
3184	! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
3185	! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
3186	! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
3187	! ----------------------------------------------------------------------
3188	T11 = T1
3189	CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
3190	TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3191	QUARTZ,CSOIL,VEGTYP,ISURBAN)
3192
3193	! ----------------------------------------------------------------------
3194	! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
3195	! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
3196	! ----------------------------------------------------------------------
3197	IF (ICE == 0) THEN
3198	! NON-GLACIAL LAND
3199	IF (ESD > 0.) THEN
3200	CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3201	ELSE
3202	ESD = 0.
3203	SNOWH = 0.
3204	SNDENS = 0.
3205	SNCOND = 1.
3206	SNCOVR = 0.
3207	END IF
3208	ELSEIF (ICE == 1) THEN
3209	! SEA-ICE
3210	IF (ESD .GE. 0.01) THEN
3211	CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3212	ELSE
3213	ESD = 0.01
3214	SNOWH = 0.05
3215	!KWM???? SNDENS =
3216	!KWM???? SNCOND =
3217	SNCOVR = 1.0
3218	ENDIF
3219	ELSEIF (ICE == -1) THEN
3220	! GLACIAL-ICE
3221	IF (ESD .GE. 0.10) THEN
3222	CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3223	ELSE
3224	ESD = 0.10
3225	SNOWH = 0.50
3226	!KWM???? SNDENS =
3227	!KWM???? SNCOND =
3228	SNCOVR = 1.0
3229	ENDIF
3230	ENDIF
3231	! ----------------------------------------------------------------------
3232	END SUBROUTINE SNOPAC
3233	! ----------------------------------------------------------------------
3234
3235
3236	SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
3237
3238	! ----------------------------------------------------------------------
3239	! SUBROUTINE SNOWPACK
3240	! ----------------------------------------------------------------------
3241	! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
3242	! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
3243	! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
3244	! KOREN, 03/25/95.
3245	! ----------------------------------------------------------------------
3246	! ESD WATER EQUIVALENT OF SNOW (M)
3247	! DTSEC TIME STEP (SEC)
3248	! SNOWH SNOW DEPTH (M)
3249	! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
3250	! TSNOW SNOW SURFACE TEMPERATURE (K)
3251	! TSOIL SOIL SURFACE TEMPERATURE (K)
3252
3253	! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
3254	! ----------------------------------------------------------------------
3255	IMPLICIT NONE
3256
3257	INTEGER :: IPOL, J
3258	REAL, INTENT(IN) :: ESD, DTSEC,TSNOW,TSOIL
3259	REAL, INTENT(INOUT) :: SNOWH, SNDENS
3260	REAL :: BFAC,DSX,DTHR,DW,SNOWHC,PEXP, &
3261	TAVGC,TSNOWC,TSOILC,ESDC,ESDCX
3262	REAL, PARAMETER :: C1 = 0.01, C2 = 21.0, G = 9.81, &
3263	KN = 4000.0
3264	! ----------------------------------------------------------------------
3265	! CONVERSION INTO SIMULATION UNITS
3266	! ----------------------------------------------------------------------
3267	SNOWHC = SNOWH *100.
3268	ESDC = ESD *100.
3269	DTHR = DTSEC /3600.
3270	TSNOWC = TSNOW -273.15
3271	TSOILC = TSOIL -273.15
3272
3273	! ----------------------------------------------------------------------
3274	! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
3275	! ----------------------------------------------------------------------
3276	! ----------------------------------------------------------------------
3277	! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
3278	! SNDENS=DS0(EXP(BFACESD)-1.)/(BFAC*ESD)
3279	! BFAC=DTHRC1EXP(0.08TAVGC-C2DS0)
3280	! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
3281	! NUMERICALLY BELOW:
3282	! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
3283	! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
3284	! ----------------------------------------------------------------------
3285	TAVGC = 0.5* (TSNOWC + TSOILC)
3286	IF (ESDC > 1.E-2) THEN
3287	ESDCX = ESDC
3288	ELSE
3289	ESDCX = 1.E-2
3290	END IF
3291
3292	! DSX = SNDENS((DEXP(BFACESDC)-1.)/(BFAC*ESDC))
3293	! ----------------------------------------------------------------------
3294	! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
3295	! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
3296	! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
3297	! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
3298	! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
3299	! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
3300	! POLYNOMIAL EXPANSION.
3301
3302	! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
3303	! IS GOVERNED BY ITERATION LIMIT "IPOL".
3304	! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
3305	! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
3306	! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
3307	! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
3308	! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
3309	! ----------------------------------------------------------------------
3310	BFAC = DTHR * C1* EXP (0.08* TAVGC - C2* SNDENS)
3311	IPOL = 4
3312	PEXP = 0.
3313	! PEXP = (1. + PEXP)BFACESDC/REAL(J+1)
3314	DO J = IPOL,1, -1
3315	PEXP = (1. + PEXP)* BFAC * ESDCX / REAL (J +1)
3316	END DO
3317
3318	PEXP = PEXP + 1.
3319	! ----------------------------------------------------------------------
3320	! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
3321	! ----------------------------------------------------------------------
3322	! END OF KOREAN FORMULATION
3323
3324	! BASE FORMULATION (COGLEY ET AL., 1990)
3325	! CONVERT DENSITY FROM G/CM3 TO KG/M3
3326	! DSM=SNDENS*1000.0
3327
3328	! DSX=DSM+DTSEC0.5DSMGESD/
3329	! & (1E7EXP(-0.02DSM+KN/(TAVGC+273.16)-14.643))
3330
3331	! & CONVERT DENSITY FROM KG/M3 TO G/CM3
3332	! DSX=DSX/1000.0
3333
3334	! END OF COGLEY ET AL. FORMULATION
3335
3336	! ----------------------------------------------------------------------
3337	! SET UPPER/LOWER LIMIT ON SNOW DENSITY
3338	! ----------------------------------------------------------------------
3339	DSX = SNDENS * (PEXP)
3340	IF (DSX > 0.40) DSX = 0.40
3341	IF (DSX < 0.05) DSX = 0.05
3342	! ----------------------------------------------------------------------
3343	! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
3344	! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
3345	! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
3346	! ----------------------------------------------------------------------
3347	SNDENS = DSX
3348	IF (TSNOWC >= 0.) THEN
3349	DW = 0.13* DTHR /24.
3350	SNDENS = SNDENS * (1. - DW) + DW
3351	IF (SNDENS >= 0.40) SNDENS = 0.40
3352	! ----------------------------------------------------------------------
3353	! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
3354	! CHANGE SNOW DEPTH UNITS TO METERS
3355	! ----------------------------------------------------------------------
3356	END IF
3357	SNOWHC = ESDC / SNDENS
3358	SNOWH = SNOWHC *0.01
3359
3360	! ----------------------------------------------------------------------
3361	END SUBROUTINE SNOWPACK
3362	! ----------------------------------------------------------------------
3363
3364	SUBROUTINE SNOWZ0 (SNCOVR,Z0, Z0BRD, SNOWH)
3365
3366	! ----------------------------------------------------------------------
3367	! SUBROUTINE SNOWZ0
3368	! ----------------------------------------------------------------------
3369	! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
3370	! SNCOVR FRACTIONAL SNOW COVER
3371	! Z0 ROUGHNESS LENGTH (m)
3372	! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
3373	! ----------------------------------------------------------------------
3374	IMPLICIT NONE
3375	REAL, INTENT(IN) :: SNCOVR, Z0BRD
3376	REAL, INTENT(OUT) :: Z0
3377	REAL, PARAMETER :: Z0S=0.001
3378	REAL, INTENT(IN) :: SNOWH
3379	REAL :: BURIAL
3380	REAL :: Z0EFF
3381
3382	!m Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0S
3383	BURIAL = 7.0*Z0BRD - SNOWH
3384	IF(BURIAL.LE.0.0007) THEN
3385	Z0EFF = Z0S
3386	ELSE
3387	Z0EFF = BURIAL/7.0
3388	ENDIF
3389
3390	Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0EFF
3391
3392	! ----------------------------------------------------------------------
3393	END SUBROUTINE SNOWZ0
3394	! ----------------------------------------------------------------------
3395
3396
3397	SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
3398
3399	! ----------------------------------------------------------------------
3400	! SUBROUTINE SNOW_NEW
3401	! ----------------------------------------------------------------------
3402	! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
3403	! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
3404
3405	! TEMP AIR TEMPERATURE (K)
3406	! NEWSN NEW SNOWFALL (M)
3407	! SNOWH SNOW DEPTH (M)
3408	! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
3409	! ----------------------------------------------------------------------
3410	IMPLICIT NONE
3411	REAL, INTENT(IN) :: NEWSN, TEMP
3412	REAL, INTENT(INOUT) :: SNDENS, SNOWH
3413	REAL :: DSNEW, HNEWC, SNOWHC,NEWSNC,TEMPC
3414
3415	! ----------------------------------------------------------------------
3416	! CONVERSION INTO SIMULATION UNITS
3417	! ----------------------------------------------------------------------
3418	SNOWHC = SNOWH *100.
3419	NEWSNC = NEWSN *100.
3420
3421	! ----------------------------------------------------------------------
3422	! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
3423	! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
3424	! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
3425	! VEMADOLEN, SWEDEN, 1980, 172-177PP.
3426	!-----------------------------------------------------------------------
3427	TEMPC = TEMP -273.15
3428	IF (TEMPC <= -15.) THEN
3429	DSNEW = 0.05
3430	ELSE
3431	DSNEW = 0.05+0.0017* (TEMPC +15.)**1.5
3432	END IF
3433	! ----------------------------------------------------------------------
3434	! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
3435	! ----------------------------------------------------------------------
3436	HNEWC = NEWSNC / DSNEW
3437	IF (SNOWHC + HNEWC .LT. 1.0E-3) THEN
3438	SNDENS = MAX(DSNEW,SNDENS)
3439	ELSE
3440	SNDENS = (SNOWHC * SNDENS + HNEWC * DSNEW)/ (SNOWHC + HNEWC)
3441	ENDIF
3442	SNOWHC = SNOWHC + HNEWC
3443	SNOWH = SNOWHC *0.01
3444
3445	! ----------------------------------------------------------------------
3446	END SUBROUTINE SNOW_NEW
3447	! ----------------------------------------------------------------------
3448
3449	SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
3450	ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3451	RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
3452
3453	! ----------------------------------------------------------------------
3454	! SUBROUTINE SRT
3455	! ----------------------------------------------------------------------
3456	! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
3457	! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
3458	! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
3459	! ----------------------------------------------------------------------
3460	IMPLICIT NONE
3461	INTEGER, INTENT(IN) :: NSOIL
3462	INTEGER :: IALP1, IOHINF, J, JJ, K, KS
3463	REAL, INTENT(IN) :: BEXP, DKSAT, DT, DWSAT, EDIR, FRZX, &
3464	KDT, PCPDRP, SLOPE, SMCMAX, SMCWLT
3465	REAL, INTENT(OUT) :: RUNOFF1, RUNOFF2
3466	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET, SH2O, SH2OA, SICE, &
3467	ZSOIL
3468	REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTT
3469	REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI, CI
3470	REAL, DIMENSION(1:NSOIL) :: DMAX
3471	REAL :: ACRT, DD, DDT, DDZ, DDZ2, DENOM, &
3472	DENOM2,DICE, DSMDZ, DSMDZ2, DT1, &
3473	FCR,INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
3474	PX, SICEMAX,SLOPX, SMCAV, SSTT, &
3475	SUM, VAL, WCND, WCND2, WDF, WDF2
3476	INTEGER, PARAMETER :: CVFRZ = 3
3477
3478	! ----------------------------------------------------------------------
3479	! FROZEN GROUND VERSION:
3480	! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
3481	! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
3482	! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
3483	! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
3484	! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
3485	! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
3486	! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
3487	! ----------------------------------------------------------------------
3488
3489	! ----------------------------------------------------------------------
3490	! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
3491	! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
3492	! MODIFIED BY Q DUAN
3493	! ----------------------------------------------------------------------
3494	! ----------------------------------------------------------------------
3495	! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
3496	! LAYERS.
3497	! ----------------------------------------------------------------------
3498	IOHINF = 1
3499	SICEMAX = 0.0
3500	DO KS = 1,NSOIL
3501	IF (SICE (KS) > SICEMAX) SICEMAX = SICE (KS)
3502	! ----------------------------------------------------------------------
3503	! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
3504	! ----------------------------------------------------------------------
3505	END DO
3506	PDDUM = PCPDRP
3507	RUNOFF1 = 0.0
3508
3509	! ----------------------------------------------------------------------
3510	! MODIFIED BY Q. DUAN, 5/16/94
3511	! ----------------------------------------------------------------------
3512	! IF (IOHINF == 1) THEN
3513
3514	IF (PCPDRP /= 0.0) THEN
3515	DT1 = DT /86400.
3516	SMCAV = SMCMAX - SMCWLT
3517
3518	! ----------------------------------------------------------------------
3519	! FROZEN GROUND VERSION:
3520	! ----------------------------------------------------------------------
3521	DMAX (1)= - ZSOIL (1)* SMCAV
3522
3523	DICE = - ZSOIL (1) * SICE (1)
3524	DMAX (1)= DMAX (1)* (1.0- (SH2OA (1) + SICE (1) - SMCWLT)/ &
3525	SMCAV)
3526
3527	DD = DMAX (1)
3528
3529	! ----------------------------------------------------------------------
3530	! FROZEN GROUND VERSION:
3531	! ----------------------------------------------------------------------
3532	DO KS = 2,NSOIL
3533
3534	DICE = DICE+ ( ZSOIL (KS -1) - ZSOIL (KS) ) * SICE (KS)
3535	DMAX (KS) = (ZSOIL (KS -1) - ZSOIL (KS))* SMCAV
3536	DMAX (KS) = DMAX (KS)* (1.0- (SH2OA (KS) + SICE (KS) &
3537	- SMCWLT)/ SMCAV)
3538	DD = DD+ DMAX (KS)
3539	! ----------------------------------------------------------------------
3540	! VAL = (1.-EXP(-KDT*SQRT(DT1)))
3541	! IN BELOW, REMOVE THE SQRT IN ABOVE
3542	! ----------------------------------------------------------------------
3543	END DO
3544	VAL = (1. - EXP ( - KDT * DT1))
3545	DDT = DD * VAL
3546	PX = PCPDRP * DT
3547	IF (PX < 0.0) PX = 0.0
3548
3549	! ----------------------------------------------------------------------
3550	! FROZEN GROUND VERSION:
3551	! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
3552	! ----------------------------------------------------------------------
3553	INFMAX = (PX * (DDT / (PX + DDT)))/ DT
3554	FCR = 1.
3555	IF (DICE > 1.E-2) THEN
3556	ACRT = CVFRZ * FRZX / DICE
3557	SUM = 1.
3558	IALP1 = CVFRZ - 1
3559	DO J = 1,IALP1
3560	K = 1
3561	DO JJ = J +1,IALP1
3562	K = K * JJ
3563	END DO
3564	SUM = SUM + (ACRT ** ( CVFRZ - J)) / FLOAT (K)
3565	END DO
3566	FCR = 1. - EXP ( - ACRT) * SUM
3567	END IF
3568
3569	! ----------------------------------------------------------------------
3570	! CORRECTION OF INFILTRATION LIMITATION:
3571	! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
3572	! HYDROLIC CONDUCTIVITY
3573	! ----------------------------------------------------------------------
3574	! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
3575	INFMAX = INFMAX * FCR
3576
3577	MXSMC = SH2OA (1)
3578	CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
3579	SICEMAX)
3580	INFMAX = MAX (INFMAX,WCND)
3581
3582	INFMAX = MIN (INFMAX,PX/DT)
3583	IF (PCPDRP > INFMAX) THEN
3584	RUNOFF1 = PCPDRP - INFMAX
3585	PDDUM = INFMAX
3586	END IF
3587	! ----------------------------------------------------------------------
3588	! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
3589	! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
3590	! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
3591	! ----------------------------------------------------------------------
3592	END IF
3593
3594	MXSMC = SH2OA (1)
3595	CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
3596	SICEMAX)
3597	! ----------------------------------------------------------------------
3598	! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
3599	! ----------------------------------------------------------------------
3600	DDZ = 1. / ( - .5 * ZSOIL (2) )
3601	AI (1) = 0.0
3602	BI (1) = WDF * DDZ / ( - ZSOIL (1) )
3603
3604	! ----------------------------------------------------------------------
3605	! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
3606	! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
3607	! ----------------------------------------------------------------------
3608	CI (1) = - BI (1)
3609	DSMDZ = ( SH2O (1) - SH2O (2) ) / ( - .5 * ZSOIL (2) )
3610	RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET (1))/ ZSOIL (1)
3611
3612	! ----------------------------------------------------------------------
3613	! INITIALIZE DDZ2
3614	! ----------------------------------------------------------------------
3615	SSTT = WDF * DSMDZ + WCND+ EDIR + ET (1)
3616
3617	! ----------------------------------------------------------------------
3618	! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
3619	! ----------------------------------------------------------------------
3620	DDZ2 = 0.0
3621	DO K = 2,NSOIL
3622	DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
3623	IF (K /= NSOIL) THEN
3624
3625	! ----------------------------------------------------------------------
3626	! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
3627	! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
3628	! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
3629	! ----------------------------------------------------------------------
3630	SLOPX = 1.
3631
3632	MXSMC2 = SH2OA (K)
3633	CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
3634	SICEMAX)
3635	! -----------------------------------------------------------------------
3636	! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
3637	! ----------------------------------------------------------------------
3638	DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
3639
3640	! ----------------------------------------------------------------------
3641	! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
3642	! ----------------------------------------------------------------------
3643	DSMDZ2 = (SH2O (K) - SH2O (K +1)) / (DENOM * 0.5)
3644	DDZ2 = 2.0 / DENOM
3645	CI (K) = - WDF2 * DDZ2 / DENOM2
3646
3647	ELSE
3648	! ----------------------------------------------------------------------
3649	! SLOPE OF BOTTOM LAYER IS INTRODUCED
3650	! ----------------------------------------------------------------------
3651
3652	! ----------------------------------------------------------------------
3653	! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
3654	! THIS LAYER
3655	! ----------------------------------------------------------------------
3656	SLOPX = SLOPE
3657	CALL WDFCND (WDF2,WCND2,SH2OA (NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
3658	SICEMAX)
3659
3660	! ----------------------------------------------------------------------
3661	! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
3662	! ----------------------------------------------------------------------
3663
3664	! ----------------------------------------------------------------------
3665	! SET MATRIX COEF CI TO ZERO
3666	! ----------------------------------------------------------------------
3667	DSMDZ2 = 0.0
3668	CI (K) = 0.0
3669	! ----------------------------------------------------------------------
3670	! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
3671	! ----------------------------------------------------------------------
3672	END IF
3673	NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2- (WDF * DSMDZ) &
3674	- WCND+ ET (K)
3675
3676	! ----------------------------------------------------------------------
3677	! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
3678	! ----------------------------------------------------------------------
3679	RHSTT (K) = NUMER / ( - DENOM2)
3680	AI (K) = - WDF * DDZ / DENOM2
3681
3682	! ----------------------------------------------------------------------
3683	! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
3684	! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
3685	! ----------------------------------------------------------------------
3686	BI (K) = - ( AI (K) + CI (K) )
3687	IF (K .eq. NSOIL) THEN
3688	RUNOFF2 = SLOPX * WCND2
3689	END IF
3690	IF (K .ne. NSOIL) THEN
3691	WDF = WDF2
3692	WCND = WCND2
3693	DSMDZ = DSMDZ2
3694	DDZ = DDZ2
3695	END IF
3696	END DO
3697	! ----------------------------------------------------------------------
3698	END SUBROUTINE SRT
3699	! ----------------------------------------------------------------------
3700
3701	SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
3702	NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
3703	AI,BI,CI)
3704
3705	! ----------------------------------------------------------------------
3706	! SUBROUTINE SSTEP
3707	! ----------------------------------------------------------------------
3708	! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
3709	! CONTENT VALUES.
3710	! ----------------------------------------------------------------------
3711	IMPLICIT NONE
3712	INTEGER, INTENT(IN) :: NSOIL
3713	INTEGER :: I, K, KK11
3714
3715	REAL, INTENT(IN) :: CMCMAX, DT, SMCMAX
3716	REAL, INTENT(OUT) :: RUNOFF3
3717	REAL, INTENT(INOUT) :: CMC
3718	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SH2OIN, SICE, ZSOIL
3719	REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: SH2OOUT
3720	REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: RHSTT, SMC
3721	REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: AI, BI, CI
3722	REAL, DIMENSION(1:NSOIL) :: RHSTTin
3723	REAL, DIMENSION(1:NSOIL) :: CIin
3724	REAL :: DDZ, RHSCT, STOT, WPLUS
3725
3726	! ----------------------------------------------------------------------
3727	! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
3728	! TRI-DIAGONAL MATRIX ROUTINE.
3729	! ----------------------------------------------------------------------
3730	DO K = 1,NSOIL
3731	RHSTT (K) = RHSTT (K) * DT
3732	AI (K) = AI (K) * DT
3733	BI (K) = 1. + BI (K) * DT
3734	CI (K) = CI (K) * DT
3735	END DO
3736	! ----------------------------------------------------------------------
3737	! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
3738	! ----------------------------------------------------------------------
3739	DO K = 1,NSOIL
3740	RHSTTin (K) = RHSTT (K)
3741	END DO
3742	DO K = 1,NSOIL
3743	CIin (K) = CI (K)
3744	END DO
3745	! ----------------------------------------------------------------------
3746	! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
3747	! ----------------------------------------------------------------------
3748	CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
3749	! ----------------------------------------------------------------------
3750	! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
3751	! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
3752	! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
3753	! ----------------------------------------------------------------------
3754	WPLUS = 0.0
3755	RUNOFF3 = 0.
3756
3757	DDZ = - ZSOIL (1)
3758	DO K = 1,NSOIL
3759	IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
3760	SH2OOUT (K) = SH2OIN (K) + CI (K) + WPLUS / DDZ
3761	STOT = SH2OOUT (K) + SICE (K)
3762	IF (STOT > SMCMAX) THEN
3763	IF (K .eq. 1) THEN
3764	DDZ = - ZSOIL (1)
3765	ELSE
3766	KK11 = K - 1
3767	DDZ = - ZSOIL (K) + ZSOIL (KK11)
3768	END IF
3769	WPLUS = (STOT - SMCMAX) * DDZ
3770	ELSE
3771	WPLUS = 0.
3772	END IF
3773	SMC (K) = MAX ( MIN (STOT,SMCMAX),0.02 )
3774	SH2OOUT (K) = MAX ( (SMC (K) - SICE (K)),0.0)
3775	END DO
3776
3777	! ----------------------------------------------------------------------
3778	! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
3779	! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
3780	! ----------------------------------------------------------------------
3781	RUNOFF3 = WPLUS
3782	CMC = CMC + DT * RHSCT
3783	IF (CMC < 1.E-20) CMC = 0.0
3784	CMC = MIN (CMC,CMCMAX)
3785
3786	! ----------------------------------------------------------------------
3787	END SUBROUTINE SSTEP
3788	! ----------------------------------------------------------------------
3789
3790	SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
3791
3792	! ----------------------------------------------------------------------
3793	! SUBROUTINE TBND
3794	! ----------------------------------------------------------------------
3795	! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
3796	! THE MIDDLE LAYER TEMPERATURES
3797	! ----------------------------------------------------------------------
3798	IMPLICIT NONE
3799	INTEGER, INTENT(IN) :: NSOIL
3800	INTEGER :: K
3801	REAL, INTENT(IN) :: TB, TU, ZBOT
3802	REAL, INTENT(OUT) :: TBND1
3803	REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL
3804	REAL :: ZB, ZUP
3805	REAL, PARAMETER :: T0 = 273.15
3806
3807	! ----------------------------------------------------------------------
3808	! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
3809	! ----------------------------------------------------------------------
3810	IF (K == 1) THEN
3811	ZUP = 0.
3812	ELSE
3813	ZUP = ZSOIL (K -1)
3814	END IF
3815	! ----------------------------------------------------------------------
3816	! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
3817	! TEMPERATURE INTO THE LAST LAYER BOUNDARY
3818	! ----------------------------------------------------------------------
3819	IF (K == NSOIL) THEN
3820	ZB = 2.* ZBOT - ZSOIL (K)
3821	ELSE
3822	ZB = ZSOIL (K +1)
3823	END IF
3824	! ----------------------------------------------------------------------
3825	! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
3826	! ----------------------------------------------------------------------
3827
3828	TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
3829	! ----------------------------------------------------------------------
3830	END SUBROUTINE TBND
3831	! ----------------------------------------------------------------------
3832
3833
3834	SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
3835
3836	! ----------------------------------------------------------------------
3837	! SUBROUTINE TDFCND
3838	! ----------------------------------------------------------------------
3839	! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
3840	! POINT AND TIME.
3841	! ----------------------------------------------------------------------
3842	! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
3843	! June 2001 CHANGES: FROZEN SOIL CONDITION.
3844	! ----------------------------------------------------------------------
3845	IMPLICIT NONE
3846	REAL, INTENT(IN) :: QZ, SMC, SMCMAX, SH2O
3847	REAL, INTENT(OUT) :: DF
3848	REAL :: AKE, GAMMD, THKDRY, THKICE, THKO, &
3849	THKQTZ,THKSAT,THKS,THKW,SATRATIO,XU, &
3850	XUNFROZ
3851
3852	! ----------------------------------------------------------------------
3853	! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
3854	! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
3855	! & 0.35, 0.60, 0.40, 0.82/
3856	! ----------------------------------------------------------------------
3857	! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
3858	! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
3859	! ----------------------------------------------------------------------
3860	! THKW ......WATER THERMAL CONDUCTIVITY
3861	! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
3862	! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
3863	! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
3864	! THKICE ....ICE THERMAL CONDUCTIVITY
3865	! SMCMAX ....POROSITY (= SMCMAX)
3866	! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
3867	! ----------------------------------------------------------------------
3868	! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
3869
3870	! PABLO GRUNMANN, 08/17/98
3871	! REFS.:
3872	! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
3873	! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
3874	! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
3875	! UNIVERSITY OF TRONDHEIM,
3876	! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
3877	! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
3878	! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
3879	! VOL. 55, PP. 1209-1224.
3880	! ----------------------------------------------------------------------
3881	! NEEDS PARAMETERS
3882	! POROSITY(SOIL TYPE):
3883	! POROS = SMCMAX
3884	! SATURATION RATIO:
3885	! PARAMETERS W/(M.K)
3886	SATRATIO = SMC / SMCMAX
3887	! ICE CONDUCTIVITY:
3888	THKICE = 2.2
3889	! WATER CONDUCTIVITY:
3890	THKW = 0.57
3891	! THERMAL CONDUCTIVITY OF "OTHER" SOIL COMPONENTS
3892	! IF (QZ .LE. 0.2) THKO = 3.0
3893	THKO = 2.0
3894	! QUARTZ' CONDUCTIVITY
3895	THKQTZ = 7.7
3896	! SOLIDS' CONDUCTIVITY
3897	THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
3898
3899	! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
3900	XUNFROZ = SH2O / SMC
3901	! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
3902	XU = XUNFROZ * SMCMAX
3903
3904	! SATURATED THERMAL CONDUCTIVITY
3905	THKSAT = THKS ** (1. - SMCMAX)* THKICE ** (SMCMAX - XU)* THKW ** &
3906	(XU)
3907
3908	! DRY DENSITY IN KG/M3
3909	GAMMD = (1. - SMCMAX)*2700.
3910
3911	! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
3912	THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
3913	! FROZEN
3914	IF ( (SH2O + 0.0005) < SMC ) THEN
3915	AKE = SATRATIO
3916	! UNFROZEN
3917	! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
3918	ELSE
3919
3920	! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
3921	! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
3922	! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
3923
3924	IF ( SATRATIO > 0.1 ) THEN
3925
3926	AKE = LOG10 (SATRATIO) + 1.0
3927
3928	! USE K = KDRY
3929	ELSE
3930
3931	AKE = 0.0
3932	END IF
3933	! THERMAL CONDUCTIVITY
3934
3935	END IF
3936
3937	DF = AKE * (THKSAT - THKDRY) + THKDRY
3938	! ----------------------------------------------------------------------
3939	END SUBROUTINE TDFCND
3940	! ----------------------------------------------------------------------
3941
3942	SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
3943
3944	! ----------------------------------------------------------------------
3945	! SUBROUTINE TMPAVG
3946	! ----------------------------------------------------------------------
3947	! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
3948	! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
3949	! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
3950	! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
3951	! ----------------------------------------------------------------------
3952	IMPLICIT NONE
3953	INTEGER K
3954
3955	INTEGER NSOIL
3956	REAL DZ
3957	REAL DZH
3958	REAL T0
3959	REAL TAVG
3960	REAL TDN
3961	REAL TM
3962	REAL TUP
3963	REAL X0
3964	REAL XDN
3965	REAL XUP
3966
3967	REAL ZSOIL (NSOIL)
3968
3969	! ----------------------------------------------------------------------
3970	PARAMETER (T0 = 2.7315E2)
3971	IF (K .eq. 1) THEN
3972	DZ = - ZSOIL (1)
3973	ELSE
3974	DZ = ZSOIL (K -1) - ZSOIL (K)
3975	END IF
3976
3977	DZH = DZ *0.5
3978	IF (TUP .lt. T0) THEN
3979	IF (TM .lt. T0) THEN
3980	! ----------------------------------------------------------------------
3981	! TUP, TM, TDN < T0
3982	! ----------------------------------------------------------------------
3983	IF (TDN .lt. T0) THEN
3984	TAVG = (TUP + 2.0* TM + TDN)/ 4.0
3985	! ----------------------------------------------------------------------
3986	! TUP & TM < T0, TDN .ge. T0
3987	! ----------------------------------------------------------------------
3988	ELSE
3989	X0 = (T0- TM) * DZH / (TDN - TM)
3990	TAVG = 0.5 * (TUP * DZH + TM * (DZH + X0) + T0* ( &
3991	& 2.* DZH - X0)) / DZ
3992	END IF
3993	ELSE
3994	! ----------------------------------------------------------------------
3995	! TUP < T0, TM .ge. T0, TDN < T0
3996	! ----------------------------------------------------------------------
3997	IF (TDN .lt. T0) THEN
3998	XUP = (T0- TUP) * DZH / (TM - TUP)
3999	XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4000	TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP - XDN) &
4001	& + TDN * XDN) / DZ
4002	! ----------------------------------------------------------------------
4003	! TUP < T0, TM .ge. T0, TDN .ge. T0
4004	! ----------------------------------------------------------------------
4005	ELSE
4006	XUP = (T0- TUP) * DZH / (TM - TUP)
4007	TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP)) / DZ
4008	END IF
4009	END IF
4010	ELSE
4011	IF (TM .lt. T0) THEN
4012	! ----------------------------------------------------------------------
4013	! TUP .ge. T0, TM < T0, TDN < T0
4014	! ----------------------------------------------------------------------
4015	IF (TDN .lt. T0) THEN
4016	XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4017	TAVG = 0.5 * (T0* (DZ - XUP) + TM * (DZH + XUP) &
4018	& + TDN * DZH) / DZ
4019	! ----------------------------------------------------------------------
4020	! TUP .ge. T0, TM < T0, TDN .ge. T0
4021	! ----------------------------------------------------------------------
4022	ELSE
4023	XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4024	XDN = (T0- TM) * DZH / (TDN - TM)
4025	TAVG = 0.5 * (T0* (2.* DZ - XUP - XDN) + TM * &
4026	& (XUP + XDN)) / DZ
4027	END IF
4028	ELSE
4029	! ----------------------------------------------------------------------
4030	! TUP .ge. T0, TM .ge. T0, TDN < T0
4031	! ----------------------------------------------------------------------
4032	IF (TDN .lt. T0) THEN
4033	XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4034	TAVG = (T0* (DZ - XDN) +0.5* (T0+ TDN)* XDN) / DZ
4035	! ----------------------------------------------------------------------
4036	! TUP .ge. T0, TM .ge. T0, TDN .ge. T0
4037	! ----------------------------------------------------------------------
4038	ELSE
4039	TAVG = (TUP + 2.0* TM + TDN) / 4.0
4040	END IF
4041	END IF
4042	END IF
4043	! ----------------------------------------------------------------------
4044	END SUBROUTINE TMPAVG
4045	! ----------------------------------------------------------------------
4046
4047	SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
4048	& CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT, &
4049	& RTDIS)
4050
4051	! ----------------------------------------------------------------------
4052	! SUBROUTINE TRANSP
4053	! ----------------------------------------------------------------------
4054	! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
4055	! ----------------------------------------------------------------------
4056	IMPLICIT NONE
4057	INTEGER I
4058	INTEGER K
4059	INTEGER NSOIL
4060
4061	INTEGER NROOT
4062	REAL CFACTR
4063	REAL CMC
4064	REAL CMCMAX
4065	REAL DENOM
4066	REAL ET (NSOIL)
4067	REAL ETP1
4068	REAL ETP1A
4069	!.....REAL PART(NSOIL)
4070	REAL GX (NROOT)
4071	REAL PC
4072	REAL Q2
4073	REAL RTDIS (NSOIL)
4074	REAL RTX
4075	REAL SFCTMP
4076	REAL SGX
4077	REAL SHDFAC
4078	REAL SMC (NSOIL)
4079	REAL SMCREF
4080	REAL SMCWLT
4081
4082	! ----------------------------------------------------------------------
4083	! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
4084	! ----------------------------------------------------------------------
4085	REAL ZSOIL (NSOIL)
4086	DO K = 1,NSOIL
4087	ET (K) = 0.
4088	! ----------------------------------------------------------------------
4089	! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
4090	! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
4091	! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
4092	! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
4093	! TOTAL ETP1A.
4094	! ----------------------------------------------------------------------
4095	END DO
4096	IF (CMC .ne. 0.0) THEN
4097	ETP1A = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR)
4098	ELSE
4099	ETP1A = SHDFAC * PC * ETP1
4100	END IF
4101	SGX = 0.0
4102	DO I = 1,NROOT
4103	GX (I) = ( SMC (I) - SMCWLT ) / ( SMCREF - SMCWLT )
4104	GX (I) = MAX ( MIN ( GX (I), 1. ), 0. )
4105	SGX = SGX + GX (I)
4106	END DO
4107
4108	SGX = SGX / NROOT
4109	DENOM = 0.
4110	DO I = 1,NROOT
4111	RTX = RTDIS (I) + GX (I) - SGX
4112	GX (I) = GX (I) * MAX ( RTX, 0. )
4113	DENOM = DENOM + GX (I)
4114	END DO
4115
4116	IF (DENOM .le. 0.0) DENOM = 1.
4117	DO I = 1,NROOT
4118	ET (I) = ETP1A * GX (I) / DENOM
4119	! ----------------------------------------------------------------------
4120	! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
4121	! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
4122	! ----------------------------------------------------------------------
4123	! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
4124	! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
4125	! ----------------------------------------------------------------------
4126	! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
4127	! ----------------------------------------------------------------------
4128	! ET(1) = RTDIS(1) * ETP1A
4129	! ET(1) = ETP1A * PART(1)
4130	! ----------------------------------------------------------------------
4131	! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
4132	! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
4133	! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
4134	! ----------------------------------------------------------------------
4135	! DO K = 2,NROOT
4136	! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
4137	! GX = MAX ( MIN ( GX, 1. ), 0. )
4138	! TEST CANOPY RESISTANCE
4139	! GX = 1.0
4140	! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))GXETP1A
4141	! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
4142	! ----------------------------------------------------------------------
4143	! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
4144	! ----------------------------------------------------------------------
4145	! ET(K) = RTDIS(K) * ETP1A
4146	! ET(K) = ETP1A*PART(K)
4147	! END DO
4148	END DO
4149	! ----------------------------------------------------------------------
4150	END SUBROUTINE TRANSP
4151	! ----------------------------------------------------------------------
4152
4153	SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4154	& SICEMAX)
4155
4156	! ----------------------------------------------------------------------
4157	! SUBROUTINE WDFCND
4158	! ----------------------------------------------------------------------
4159	! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
4160	! ----------------------------------------------------------------------
4161	IMPLICIT NONE
4162	REAL BEXP
4163	REAL DKSAT
4164	REAL DWSAT
4165	REAL EXPON
4166	REAL FACTR1
4167	REAL FACTR2
4168	REAL SICEMAX
4169	REAL SMC
4170	REAL SMCMAX
4171	REAL VKwgt
4172	REAL WCND
4173
4174	! ----------------------------------------------------------------------
4175	! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
4176	! ----------------------------------------------------------------------
4177	REAL WDF
4178	FACTR1 = 0.05 / SMCMAX
4179
4180	! ----------------------------------------------------------------------
4181	! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
4182	! ----------------------------------------------------------------------
4183	FACTR2 = SMC / SMCMAX
4184	FACTR1 = MIN(FACTR1,FACTR2)
4185	EXPON = BEXP + 2.0
4186
4187	! ----------------------------------------------------------------------
4188	! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
4189	! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
4190	! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
4191	! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
4192	! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
4193	! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
4194	! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
4195	! --
4196	! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
4197	! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
4198	! ----------------------------------------------------------------------
4199	WDF = DWSAT * FACTR2 ** EXPON
4200	IF (SICEMAX .gt. 0.0) THEN
4201	VKWGT = 1./ (1. + (500.* SICEMAX)**3.)
4202	WDF = VKWGT * WDF + (1. - VKWGT)* DWSAT * FACTR1** EXPON
4203	! ----------------------------------------------------------------------
4204	! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
4205	! ----------------------------------------------------------------------
4206	END IF
4207	EXPON = (2.0 * BEXP) + 3.0
4208	WCND = DKSAT * FACTR2 ** EXPON
4209
4210	! ----------------------------------------------------------------------
4211	END SUBROUTINE WDFCND
4212	! ----------------------------------------------------------------------
4213
4214	SUBROUTINE SFCDIF_off (ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS)
4215
4216	! ----------------------------------------------------------------------
4217	! SUBROUTINE SFCDIF (renamed SFCDIF_off to avoid clash with Eta PBL)
4218	! ----------------------------------------------------------------------
4219	! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
4220	! SEE CHEN ET AL (1997, BLM)
4221	! ----------------------------------------------------------------------
4222
4223	IMPLICIT NONE
4224	REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
4225	REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, &
4226	& SQVISC
4227	REAL RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU, &
4228	& PSLHS
4229	REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
4230	REAL SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH
4231	REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
4232	REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
4233	!CC ......REAL ZTFC
4234
4235	REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
4236	& RLMA
4237
4238	INTEGER ITRMX, ILECH, ITR
4239	PARAMETER &
4240	& (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40, &
4241	& EXCM = 0.001 &
4242	& ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG &
4243	& ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05, &
4244	& PIHF = 3.14159265/2.)
4245	PARAMETER &
4246	& (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8 &
4247	& ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0 &
4248	& ,SQVISC = 258.2)
4249	PARAMETER &
4250	& (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191 &
4251	& ,RFAC = RIC / (FHNEU * RFC * RFC))
4252
4253	! ----------------------------------------------------------------------
4254	! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
4255	! ----------------------------------------------------------------------
4256	! LECH'S SURFACE FUNCTIONS
4257	! ----------------------------------------------------------------------
4258	PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
4259	PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
4260	PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
4261
4262	! ----------------------------------------------------------------------
4263	! PAULSON'S SURFACE FUNCTIONS
4264	! ----------------------------------------------------------------------
4265	PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
4266	PSPMU (XX)= -2.* log ( (XX +1.)0.5) - log ( (XX XX +1.)*0.5) &
4267	& +2.* ATAN (XX) &
4268	&- PIHF
4269	PSPMS (YY)= 5.* YY
4270	PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
4271
4272	! ----------------------------------------------------------------------
4273	! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
4274	! OVER SOLID SURFACE (LAND, SEA-ICE).
4275	! ----------------------------------------------------------------------
4276	PSPHS (YY)= 5.* YY
4277
4278	! ----------------------------------------------------------------------
4279	! ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1
4280	! C......ZTFC=0.1
4281	! CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
4282	! ----------------------------------------------------------------------
4283	ILECH = 0
4284
4285	! ----------------------------------------------------------------------
4286	ZILFC = - CZIL * VKRM * SQVISC
4287	! C.......ZT=Z0*ZTFC
4288	ZU = Z0
4289	RDZ = 1./ ZLM
4290	CXCH = EXCM * RDZ
4291	DTHV = THLM - THZ0
4292
4293	! ----------------------------------------------------------------------
4294	! BELJARS CORRECTION OF USTAR
4295	! ----------------------------------------------------------------------
4296	DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
4297	!cc If statements to avoid TANGENT LINEAR problems near zero
4298	BTGH = BTG * HPBL
4299	IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
4300	WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
4301	ELSE
4302	WSTAR2 = 0.0
4303	END IF
4304
4305	! ----------------------------------------------------------------------
4306	! ZILITINKEVITCH APPROACH FOR ZT
4307	! ----------------------------------------------------------------------
4308	USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
4309
4310	! ----------------------------------------------------------------------
4311	ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
4312	ZSLU = ZLM + ZU
4313	! PRINT*,'ZSLT=',ZSLT
4314	! PRINT*,'ZLM=',ZLM
4315	! PRINT*,'ZT=',ZT
4316
4317	ZSLT = ZLM + ZT
4318	RLOGU = log (ZSLU / ZU)
4319
4320	RLOGT = log (ZSLT / ZT)
4321	! PRINT*,'RLMO=',RLMO
4322	! PRINT*,'ELFC=',ELFC
4323	! PRINT*,'AKHS=',AKHS
4324	! PRINT*,'DTHV=',DTHV
4325	! PRINT*,'USTAR=',USTAR
4326
4327	RLMO = ELFC * AKHS * DTHV / USTAR **3
4328	! ----------------------------------------------------------------------
4329	! 1./MONIN-OBUKKHOV LENGTH-SCALE
4330	! ----------------------------------------------------------------------
4331	DO ITR = 1,ITRMX
4332	ZETALT = MAX (ZSLT * RLMO,ZTMIN)
4333	RLMO = ZETALT / ZSLT
4334	ZETALU = ZSLU * RLMO
4335	ZETAU = ZU * RLMO
4336
4337	ZETAT = ZT * RLMO
4338	IF (ILECH .eq. 0) THEN
4339	IF (RLMO .lt. 0.)THEN
4340	XLU4 = 1. -16.* ZETALU
4341	XLT4 = 1. -16.* ZETALT
4342	XU4 = 1. -16.* ZETAU
4343
4344	XT4 = 1. -16.* ZETAT
4345	XLU = SQRT (SQRT (XLU4))
4346	XLT = SQRT (SQRT (XLT4))
4347	XU = SQRT (SQRT (XU4))
4348
4349	XT = SQRT (SQRT (XT4))
4350	! PRINT*,'-----------1------------'
4351	! PRINT*,'PSMZ=',PSMZ
4352	! PRINT*,'PSPMU(ZETAU)=',PSPMU(ZETAU)
4353	! PRINT*,'XU=',XU
4354	! PRINT*,'------------------------'
4355	PSMZ = PSPMU (XU)
4356	SIMM = PSPMU (XLU) - PSMZ + RLOGU
4357	PSHZ = PSPHU (XT)
4358	SIMH = PSPHU (XLT) - PSHZ + RLOGT
4359	ELSE
4360	ZETALU = MIN (ZETALU,ZTMAX)
4361	ZETALT = MIN (ZETALT,ZTMAX)
4362	! PRINT*,'-----------2------------'
4363	! PRINT*,'PSMZ=',PSMZ
4364	! PRINT*,'PSPMS(ZETAU)=',PSPMS(ZETAU)
4365	! PRINT*,'ZETAU=',ZETAU
4366	! PRINT*,'------------------------'
4367	PSMZ = PSPMS (ZETAU)
4368	SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
4369	PSHZ = PSPHS (ZETAT)
4370	SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
4371	END IF
4372	! ----------------------------------------------------------------------
4373	! LECH'S FUNCTIONS
4374	! ----------------------------------------------------------------------
4375	ELSE
4376	IF (RLMO .lt. 0.)THEN
4377	! PRINT*,'-----------3------------'
4378	! PRINT*,'PSMZ=',PSMZ
4379	! PRINT*,'PSLMU(ZETAU)=',PSLMU(ZETAU)
4380	! PRINT*,'ZETAU=',ZETAU
4381	! PRINT*,'------------------------'
4382	PSMZ = PSLMU (ZETAU)
4383	SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
4384	PSHZ = PSLHU (ZETAT)
4385	SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
4386	ELSE
4387	ZETALU = MIN (ZETALU,ZTMAX)
4388
4389	ZETALT = MIN (ZETALT,ZTMAX)
4390	! PRINT*,'-----------4------------'
4391	! PRINT*,'PSMZ=',PSMZ
4392	! PRINT*,'PSLMS(ZETAU)=',PSLMS(ZETAU)
4393	! PRINT*,'ZETAU=',ZETAU
4394	! PRINT*,'------------------------'
4395	PSMZ = PSLMS (ZETAU)
4396	SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
4397	PSHZ = PSLHS (ZETAT)
4398	SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
4399	END IF
4400	! ----------------------------------------------------------------------
4401	! BELJAARS CORRECTION FOR USTAR
4402	! ----------------------------------------------------------------------
4403	END IF
4404
4405	! ----------------------------------------------------------------------
4406	! ZILITINKEVITCH FIX FOR ZT
4407	! ----------------------------------------------------------------------
4408	USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
4409
4410	ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
4411	ZSLT = ZLM + ZT
4412	!-----------------------------------------------------------------------
4413	RLOGT = log (ZSLT / ZT)
4414	USTARK = USTAR * VKRM
4415	AKMS = MAX (USTARK / SIMM,CXCH)
4416	!-----------------------------------------------------------------------
4417	! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
4418	!-----------------------------------------------------------------------
4419	AKHS = MAX (USTARK / SIMH,CXCH)
4420	IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
4421	WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
4422	ELSE
4423	WSTAR2 = 0.0
4424	END IF
4425	!-----------------------------------------------------------------------
4426	RLMN = ELFC * AKHS * DTHV / USTAR **3
4427	!-----------------------------------------------------------------------
4428	! IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110
4429	!-----------------------------------------------------------------------
4430	RLMA = RLMO * WOLD+ RLMN * WNEW
4431	!-----------------------------------------------------------------------
4432	RLMO = RLMA
4433	! PRINT*,'----------------------------'
4434	! PRINT*,'SFCDIF OUTPUT ! ! ! ! ! ! ! ! ! ! ! !'
4435
4436	! PRINT*,'ZLM=',ZLM
4437	! PRINT*,'Z0=',Z0
4438	! PRINT*,'THZ0=',THZ0
4439	! PRINT*,'THLM=',THLM
4440	! PRINT*,'SFCSPD=',SFCSPD
4441	! PRINT*,'CZIL=',CZIL
4442	! PRINT*,'AKMS=',AKMS
4443	! PRINT*,'AKHS=',AKHS
4444	! PRINT*,'----------------------------'
4445
4446	END DO
4447	! ----------------------------------------------------------------------
4448	END SUBROUTINE SFCDIF_off
4449	! ----------------------------------------------------------------------
4450
4451	END MODULE module_sf_noahlsm

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

Download in other formats:

Original Format