Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_noahlsm.F @ 3584

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

Rev	Line
[2759]	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