Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_noahlsm.F @ 3431

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

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

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

Download in other formats:

Original Format