Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_lsm_nmm.F @ 3094

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

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

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

Download in other formats:

Original Format