      PROGRAM  MAIN

      IMPLICIT NONE

CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC  NCEP/OH/OSU UNCOUPLED 1-D  LAND-SURFACE MODEL: VERSION 1.0  01 MAR 99
CC
CC  THIS MAIN PROGRAM AND ITS FAMILY OF SUBROUTINES COMPRISE VERSION 1.0 
CC  OF THE  "PUBLIC" RELEASE OF THE UNCOUPLED 1-D COLUMN VERSION OF THE 
CC  NCEP LAND-SURFACE MODEL (LSM).   THE NCEP LSM IS A DESCENDANT OF AN 
CC  EARLY GENERATION OF THE OREGON STATE UNIVERSITY LSM,  BUT INCLUDES 
CC  SUBSTANTIAL PHYSICS EXTENSIONS  AND RECODING ACCOMPLISHED ALONG THE 
CC  WAY BY NCEP, OH, AFGWC, AND AFGL.
CC 
CC  FOR DOCUMENTATION OF THIS CODE AND INSTRUCTIONS ON ITS EXECUTION AND 
CC  INPUT/OUTPUT FILES, SEE DOCUMENTATION FILE README_1.0 IN SAME PUBLIC 
CC  SERVER DIRECTORY AS THIS CODE. 
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC PROGRAM HISTORY LOG
CC
CC VERSION 1.0  --  01 MARCH 99 :  ORIGINAL PUBLIC RELEASE VERSION
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC  COMING SOON:
CC
CC  1 - ARBITRARY NUMBER AND THICKNESS OF SOIL LAYERS
CC
CC  2 - IMPROVED SUBSURFACE HEAT FLUX UNDER SHALLOW SNOWPACK
CC
CC  3 - NAMELIST DIRECTED OPTIONAL INPUT OF NON-DEFAULT PARAMETER 
CC      VALUES IN SUB REDPRM
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC

      CHARACTER*100 CNTRFL, FILENAME

      INTEGER NSOLD
      PARAMETER (NSOLD=4)

      LOGICAL LSATURATED
            
      INTEGER ICE
      INTEGER IDAY
      INTEGER IIDAY
      INTEGER IMONTH
      INTEGER IREC1
      INTEGER IREC3
      INTEGER IREC5
      INTEGER IRECD
      INTEGER NDAILY
      INTEGER NOUT1
      INTEGER NOUT3
      INTEGER NOUT5
      INTEGER NREAD1
      INTEGER NRECL
      INTEGER NROOT
      INTEGER NRUN
      INTEGER NSOIL
      INTEGER ISLOPTYP
      INTEGER ISOILTP
      INTEGER IVEGTYP 
      INTEGER ITIME
      INTEGER jday
      INTEGER ITIME_ctl
      INTEGER jday_ctl
      INTEGER jday0
      INTEGER IJ
      INTEGER INDI
      
      REAL  R
      REAL  CP
      
      PARAMETER (R = 287.04, CP = 1004.5)

      REAL  BETA
      REAL  DRIP
      REAL  EC
      REAL  EDIR
      REAL  ETT
      REAL  F
      REAL  FLX1
      REAL  FLX2
      REAL  FLX3
      REAL  RHO
      REAL  RUNOF
      REAL  DEW
      REAL  RIB
      REAL  RUNOFF3
      REAL  SMSCDIF_O
      REAL  SMSCDIF
C LW      REAL  SIGMA
      REAL  LATITUDE
      REAL  LONGITUDE
C LW      REAL  LWDNCALC
      REAL  RUNOFF2
      REAL  Q1
      REAL  AET
      REAL  ALBEDO
      REAL  CH
      REAL  CHKFF
      REAL  CM
      REAL  CMC
      REAL  DQSDT
      REAL  CZIL
      REAL  REFKDT
      REAL  DQSDT2
      REAL  DT
C LW      REAL  EMISS
      REAL  ESAT
      REAL  ETA
      REAL  ETA_OBS
      REAL  ETP
      REAL  FUP
      REAL  GHF
      REAL  H
      REAL  H_OBS
      REAL  LE
      REAL  LWDN
      REAL  LW_in
      REAL  PRCP
      REAL  Par_in
      REAL  Par_out
      REAL  Q2
      REAL  Q2SAT
      REAL  RES
      REAL  RH
      REAL  RNET
      REAL  RUNOFF1
      REAL  Rg
      REAL  SFCSPD
      REAL  SFCPRS
      REAL  SFCTMP
      REAL  SHDFAC

      REAL  ALBEDOM(13)
      REAL  SHDFACM(13)
      REAL  MLAI   (13)
      REAL  XLAI

      REAL  SKN_IRT
      REAL  SMC(NSOLD)
      REAL  SMCMAX
      REAL  SNMAX
      REAL  SNOALB
      REAL  SOILW
      REAL  SOILM
      REAL  SOLDN
      REAL  STC(NSOLD)
      REAL  S
      REAL  S_OBS
      REAL  T1
      REAL  T1_OBS
      REAL  T14
      REAL  T1V
      REAL  T2V
C LW      REAL  TAK
      REAL  TBOT
      REAL  TH2
      REAL  TH2V
      REAL  T_16cm
      REAL  T_2cm
      REAL  T_32cm
      REAL  T_4cm
      REAL  T_64cm
      REAL  T_8cm
      REAL  Z
      REAL  Z0
      REAL  eddyuw
      REAL  sm_20cm
      REAL  sm_5cm
      REAL  sm_60cm
      REAL  u_bar
      REAL  uprim2
      REAL  vprim2
      REAL  w_dir
      REAL  wet
      REAL  wprim2
C
      REAL  PCPDAY
      REAL  PCPSUM
      REAL  ETADAY
      REAL  ETSUM
      REAL  RUNOFFSUM
      REAL  RUNOFFDAY
      REAL  SMCMM
      REAL  SMCDIF
      REAL  SMCNW
      REAL  SMC_OBS(NSOLD)
      REAL  ETADAY_O
      REAL  ETSUM_O
      REAL  SMCDIF_O
      REAL  SMCMM_O
      REAL  SMCNW_O
C 
      REAL  SH2O(NSOLD)
      REAL  SLDPTH(NSOLD)
      REAL  SNOWH
      REAL  SNEQV
C*************************************************************

CMIC$ TASKCOMMON RITE
      COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOF,
     &  	   DEW,RIB,RUNOFF3

C-------------------------------------------------------------
C     THERE ARE 10 STEPS IN THIS MAIN PROGRAM DRIVER
C
C DRIVER STEP 1 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 1. READ CONTROL FILE <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C >>>>>>>>>>>>>>>>>>>>>>>>>        <<<<<<<<<<<<<<<<<<<<<<<<<<<
	CNTRFL = 'controlfile'

      CALL READCNTL (CNTRFL,NRUN,DT,NSOIL,NROOT,Z,SLDPTH,
     & ISOILTP,IVEGTYP,ISLOPTYP,ALBEDOM,SHDFACM,MLAI,SNOALB,ICE,
     & TBOT,Z0, CZIL,REFKDT,T1,
     & STC,SMC,SH2O,CMC,SNOWH,SNEQV,FILENAME,
     & LATITUDE,
     & LONGITUDE,
     & jday_ctl,
     & ITIME_ctl)

C THIS IS FOR THE DAILY WATER BALANCE ------------------
C Initialize daily accum vars after reading control file
      PCPDAY = 0.0
      PCPSUM = 0.0
      ETADAY = 0.0
      ETSUM = 0.0
      ETADAY_O = 0.0
      ETSUM_O = 0.0
      RUNOFFSUM = 0.0
      RUNOFFDAY = 0.0
      SMCMM = 0.0

C ----- DO NOT INITIALIZE THIS INSIDE TIME LOOP ! -------------
C initialize variables for water balance (daily sm differences)
C initial SMCMM from control file !!!!!!!!!!!!!   
C but only 3 layers !!!!!!!!!!!!!   
          do ij = 1,3
        SMCMM = SLDPTH(ij)*1000.*SMC(ij) + SMCMM
          end do	  
C initial OBS SMCMM from control file 
C but only 3 layers !!!!!!!!!!!!!   
      SMCMM_O = SMCMM 
      SMSCDIF_O = 0.0
      SMCNW_O = 0.0
      SMSCDIF = 0.0
      SMCNW = 0.0
C THIS PART ABOVE IS FOR THE DAILY WATER BALANCE -------
C ----- DO NOT INITIALIZE THE ABOVE INSIDE TIME LOOP ! --------

C DRIVER STEP 2  OPEN INPUT/OUTPUT FILES  <<<<<<<<<<<<<<<<<<<
C 2. OPEN INPUT/OUTPUT FILES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

C >> OPEN FILE CONTAINING ATM. STATION DATA <<<<<
        NREAD1=25
        
        OPEN(NREAD1,FILE=FILENAME,STATUS='OLD')
C
C ------------------------------------------------------------
C
c (NOUT's must be positive for binary output - GrADS)
c (NOUT's must be negative for ascii output intead)
      NOUT1 = -43
      NOUT3 = -45
      NOUT5 = -47
      NDAILY= -49
C DATA WORD LENGTH (RECORD LENGTH)
      NRECL=1
C HYDROLOGY
      IF (NOUT1 .GT. 0) THEN
        IREC1 = 1
        OPEN(UNIT=NOUT1,FILE= 
     & 'HYDRO.GRS',
     & STATUS='NEW'
     & ,ACCESS='DIRECT'
     & ,FORM='UNFORMATTED'
     & ,RECL=NRECL
     &)
      ELSE
        OPEN(UNIT=-NOUT1,FILE= 
     & 'HYDRO.TXT',
     & STATUS='NEW')
      END IF
      
C THERMODYNAMICS
      IF (NOUT3 .GT. 0) THEN
        IREC3 = 1
        OPEN(UNIT=NOUT3,FILE=
     & 'THERMO.GRS',
     & STATUS='NEW'
     &,ACCESS='DIRECT'
     &,FORM='UNFORMATTED'
     &,RECL=NRECL
     &)
      ELSE
        OPEN(UNIT=-NOUT3,FILE=
     & 'THERMO.TXT',
     & STATUS='NEW')
      END IF

C OBSERVATIONAL DATA
      IF (NOUT5 .GT. 0) THEN
        IREC5 = 1
        OPEN(UNIT=NOUT5,
     &  FILE='OBS_DATA.GRS',
     &  STATUS='NEW'
     &, ACCESS='DIRECT'
     &, FORM='UNFORMATTED'
     &,RECL=NRECL
     &)
      ELSE
        OPEN(UNIT=-NOUT5,
     &  FILE='OBS_DATA.TXT',
     &  STATUS='NEW')
      END IF

C DAILY ACCUMULATION VALUES
      IF (NDAILY .GT. 0) THEN
        IRECD = 1
        OPEN(UNIT=NDAILY,
     &  FILE='DAILY.GRS',
     &  STATUS='NEW'
     &, ACCESS='DIRECT'
     &, FORM='UNFORMATTED'
     &,RECL=NRECL
     &)
      ELSE
        OPEN(UNIT=-NDAILY,
     &  FILE='DAILY.TXT',
     &  STATUS='NEW')
      END IF
C 
C >>>>>>>> OPENING OUTPUT FILES FINISHED <<<<<<<<<<<<<<<<<<<<<<<<<<<
C ------------------------------------------------------------
C Initialize Day number
	           jday0 = 0 - 1
		   IIDAY = 0

C READ FIRST RECORD OF FORCING TO CHECK THE START DATE <<<<<<<<<<<<
C 
C >>>>>>>>> READPILPS-->READBND - Read Tilden Meyers BND data <<<<<

         CALL READBND(jday,  ITIME, SFCTMP,   RH, SFCPRS,  Rg,
     * Par_in, Par_out, RNET, LW_in, GHF,  PRCP,   wet,  SKN_IRT,
     * T_2cm,  T_4cm,  T_8cm,  T_16cm, T_32cm, T_64cm, sm_5cm,
     * sm_20cm, sm_60cm,  w_dir, u_bar,     eddyuw, 
     * uprim2,      vprim2,      wprim2,         H,     LE,
     * IMONTH,IDAY,NREAD1)
     
      rewind NREAD1
     
C CHECK DATE and time for simulation on first step
      IF ((ITIME_ctl .NE. ITIME) .or. (jday_ctl .NE. jday)) THEN
       PRINT*,'date/time specified in control file does not match
     & initial read from forcing data file, the program stops'
       PRINT*,' '
       PRINT*,'Julian day in control file= ',jday_ctl
       PRINT*,'Julian day in forcing data file= ',jday
       PRINT*,' '
       PRINT*,'Initial time in control file= ',ITIME_ctl
       PRINT*,'Initial time in forcing data file= ',ITIME
       stop 999
      ENDIF

C DRIVER STEP 3  (time step loop) <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C ---------------------------------------------------------------------|
C ---------------------------------------------------------------------|
C							               |
C     START TIME LOOP					               |
C 							               |C  							               |
      DO INDI = 1,NRUN

C DRIVER STEP 4 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 4. READ FORCING DATA <<<<<<<<<<<<<<<<<<<<<<<<<
C 
C >>>>>>>>> READPILPS-->READBND - Read Tilden Meyers BND data <<<<<

         CALL READBND(jday,  ITIME, SFCTMP,   RH, SFCPRS,  Rg,
     * Par_in, Par_out, RNET, LW_in, GHF,  PRCP,   wet,  SKN_IRT,
     * T_2cm,  T_4cm,  T_8cm,  T_16cm, T_32cm, T_64cm, sm_5cm,
     * sm_20cm, sm_60cm,  w_dir, u_bar,     eddyuw, 
     * uprim2,      vprim2,      wprim2,         H,     LE,
     * IMONTH,IDAY,NREAD1)
     
C SAVE OBSERVATIONS FOR LATER VALIDATION OUTPUT

C Skin Temperature in Kelvin, (T1 IN SFLX)
      T1_OBS = SKN_IRT + 273.16
C
C >>>>>>>> SIGN CONVENTIONS FOR READ FLUXES <<<<<<<<<<<<<<<<<<
C 
C S_OBS: SOIL HEAT FLUX READ FROM DATA (S IN SFLX)
C       (W M-2: NEGATIVE, IF DOWNWARD FROM SURFACE)
      S_OBS =  GHF
C
C H_OBS: SENSIBLE HEAT FLUX READ FROM DATA (H IN SFLX)
C    (W M-2: NEGATIVE, IF UPWARD FROM SURFACE)
      H_OBS =  H
C
C ETA_OBS: LATENT HEAT FLUX READ FROM DATA (ETA IN SFLX)
C    (W M-2: NEGATIVE, IF UPWARD FROM SURFACE)
      ETA_OBS = LE
C
C SMC_OBS: Soil moisture content (volumetric) FROM DATA
C  (just for water balance purposes) - don't have 4th layer
       SMC_OBS(1) = sm_5cm
       SMC_OBS(2) = sm_20cm
       SMC_OBS(3) = sm_60cm
C       SMC_OBS(4) = 0.0
C
C Keep the day number of the simulation 
	IF (jday0 .NE. jday) THEN
	   IIDAY = IIDAY + 1
           jday0 = jday
C
C DRIVER STEP 5 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 5. INTERPOLATE DAILY LAND SURFACE CHARACTERISTICS FROM MONTHLY 
C    VALUES:
C
C........................... monthly, daily ............	   
	   CALL MONTH_D(JDAY,ALBEDOM,ALBEDO)
	   CALL MONTH_D(JDAY,SHDFACM,SHDFAC)
	   CALL MONTH_D(JDAY,MLAI   ,XLAI)
	   
	END IF
C.......................................................

C
C >>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

C  ***********************************************************
C
C                         RADIATION

C The following step (OBTENTION OF LWDN) has been 
C commented out, given that
C the current forcing data file available provides the 
C downward long wave radiative forcing as measured by NOAA's
C SURFRAD site located within 5 Km of the flux/meteorological
C measurements system site.
C 
C OBTENTION OF LWDN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 
C COMPUTATION OF LWDN (INCOMING LW RADIATION) FROM TAIR AND Q:
C   
C...............LWDN = EMISS*SIGMA*(TAK)**4.   
C
C WHERE:   TAK = AIR TEMP IN KELVIN
C        EMISS = 0.7  OR  (IDSO AND JACKSON, 1969): 
C 
C        EMISS = (1 - 0.261 EXP(-7.77*10^(-4)X(273-TAK)^2)
C
C      NEED STEFAN-BOLTZMANN CONSTANT, SIGMA 
C         SIGMA = 5.672 * 10^-8  W M^-2 T^-4
C
C	   SIGMA = 5.672E-8
C	   TAK = SFCTMP
C	   EMISS = 1 - 0.261*EXP((-7.77E-4)*(273-TAK)**2.)
C
C	   LWDN = EMISS*SIGMA*TAK**4.

C  ***********************************************************

C DRIVER STEP 6 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<	  
C 6. ASSIGN DOWNWARD SOLAR AND LONGWAVE RADIATION <<<<<<<<<<<<<

C Using LW_in READ FROM DATA instead of LWDN calculated

      LWDN = LW_in      

      SOLDN = Rg

C DRIVER STEP 7 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 7. CALCULATE A SATURATION MIX RATIO (Q2SAT) <<<<<<<<<<<<<<<<
C          
C....SET SATURATION FLAG (SEE COMMENT IN SFLX)
	  LSATURATED = .FALSE.
	  	  
C NEED Q2 (FROM REL.HUMID.) USE SUBROUTINE QDATAP

          CALL QDATAP(SFCTMP,SFCPRS,RH,Q2,Q2SAT, ESAT)

          IF (Q2 .LT. 0.1E-5) Q2=0.1E-5

          IF (Q2 .GE. Q2SAT) THEN
	  	Q2 = Q2SAT*0.99
	  	LSATURATED = .TRUE.
	  ENDIF
          
C CALCULATE SLOPE OF SAT SPECIFIC HUMIDITY CURVE FOR PENMAN: DQSDT2
          DQSDT2 = DQSDT (SFCTMP, SFCPRS)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  CALC VIRTUAL TEMPS AND POTENTIAL TEMPS AT GRND (SUB 1) AND AT
C  THE 1ST MDL LVL ABV THE GRND (SUB 2). EXPON IS CP DIVD BY R.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      T1V = T1 * (1.0 + 0.61 * Q2 )
      T2V = SFCTMP * ( 1.0 + 0.61 * Q2 )
      TH2 = SFCTMP + ( 0.0098 * Z )
      TH2V = TH2 * (1.0 + 0.61 * Q2 )          

C
C DRIVER STEP 8 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 8. CALCULATE CH (EXCHANGE COEFFICIENT) <<<<<<<<<<<<<<<<<<<<<

C SFCSPD = sqrt(u*u+v*v), if individual (u,v) components provided,
C but Tilden Meyers data provides wind speed as u_bar.

          SFCSPD = u_bar

          CALL SFCDIF ( Z, Z0, T1V, TH2V, SFCSPD,CZIL, CM, CH )

C CALCULATING CHKFF 
C PARAMETERS R = 287.04, AND CP = 1004.5
C
          RHO = SFCPRS/(R * T2V)
          CHKFF = CH * CP * RHO

C DRIVER STEP 9 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 9. CALL LAND-SURFACE PHYSICS  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
	  CALL SFLX(ICE,LSATURATED,DT,Z,NSOIL,NROOT,SLDPTH,
     I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2,TBOT,CH,CHKFF,
     I IVEGTYP,ISOILTP,XLAI,REFKDT,
     I ISLOPTYP, SHDFAC,SNOALB,
     O ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX,
     2 T1,CMC,SMC,
     2 SH2O,SNOWH,
     2 STC,SNEQV,ALBEDO,
     O SOILW,SOILM)

C 
      AET = ETA

C CALCULATE UPWARD LONGWAVE RAD USING UPDATED SKIN TEMPERATURE

      T14 = T1 * T1 * T1 * T1
      FUP = 5.67E-8 * T14

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  CALCULATE RESIDUAL OF ALL SURFACE ENERGY BALANCE EQN TERMS.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      S = -S

      F = SOLDN*(1.0-ALBEDO) + LWDN 

      RES = F - H - S - AET - FUP - 
     &      FLX1 - FLX2 - FLX3
            
CC.....S = S + RES
                
         IF ((INDI .LT. 25) .OR. (MOD(INDI,50) .EQ. 0)) THEN
C Debug      PRINT*,' CALCULATE RESIDUAL OF ALL SURFACE ENERGY'
      PRINT*,'  --------------------------------------'
      PRINT*,' RES=',RES,'     TIMESTEP=',INDI
      PRINT*,' RES/(S+ETA) (in %)=',100*RES/(S+ETA),'%'
      PRINT*,'  --------------------------------------'
         ENDIF
C
C WRITE HYDROLOGICAL VARIABLES IN GRD-FILES
C
C Debug      PRINT*,'  ---------------------------------------'
C Debug      PRINT*,'      IIDAY=',IIDAY,'	TIMESTEP=',INDI
C Debug      PRINT*,'  ---------------------------------------'

C ACCUMULATE DAILY VALUES FOR WATER BALANCE
      PCPSUM = PCPSUM + PRCP*DT
      ETSUM = ETSUM + ETA*DT/2.5E6
      ETSUM_O = ETSUM_O + ETA_OBS*DT/2.5E6
      RUNOFFSUM = RUNOFFSUM + RUNOFF1*DT*1000. + RUNOFF2*DT*1000.
      IF (ITIME .EQ. 0) THEN
C Remember to initialize SMCMM, SMSCDIF, ... 
C after reading controlfile.
        SMCNW = 0.0
        SMCNW_O =0.0
          do ij = 1,3
        SMCNW = SLDPTH(ij)*1000.*SMC(ij) + SMCNW
        SMCNW_O = SLDPTH(ij)*1000.*SMC_OBS(ij) + SMCNW_O
          end do
        SMCDIF = SMCNW - SMCMM
	SMCMM = SMCNW
        SMCDIF_O = SMCNW_O - SMCMM_O
	SMCMM_O = SMCNW_O
	RUNOFFDAY = RUNOFFSUM
	PCPDAY = PCPSUM
        ETADAY = -ETSUM
        ETADAY_O = -ETSUM_O
        PCPSUM = 0.0
        ETSUM = 0.0
        ETSUM_O = 0.0
	RUNOFFSUM = 0.0
C
C DRIVER STEP 10 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 10. Write to model output files

C  EXTENSION .GRS FILES - GrADS-readable format output (unformatted)
C  EXTENSION .TXT FILES - ASCII formatted.

C
C  WRITE DAILY ACUMMULATION VARIABLES, DAILY.GRS, GrADS FORMAT
C                                   or DAILY.TXT, ASCII formatted
C
      CALL PRTDAILY(INDI,NDAILY,NSOIL,jday,
     & IIDAY,ITIME,ETADAY,ETADAY_O,
     & PCPDAY,SMCDIF,SMCDIF_O,RUNOFFDAY,
     & DT,IRECD)

      ENDIF

C
C WRITE VARIABLES IN WATER AMMOUNT UNITS (HYDRO.GRS or .TXT)
C
      CALL PRTHYDF(INDI,NOUT1,NSOIL,jday,
     & IIDAY,ITIME,ETA,
     & ETP,PRCP,SH2O,SMC,ALBEDO,
     & SNEQV, SNMAX, SNOWH, SOILM, SOILW,
     & RUNOFF1, RUNOFF2, RUNOFF3,
     & DT,EDIR,EC,ETT,CMC,IREC1)
C
C  WRITE VARIABLES IN THERMODYNAMIC ENERGY UNITS (THERMO.GRS or .TXT)
C
      CALL PRTHMF(INDI,NOUT3,NSOIL,jday,
     &  IIDAY,ITIME,CH,CM,Z,F,S,CMC,
     &  SMCMAX,SFCTMP,T1,Q1,SFCPRS,SFCSPD,ETA,ETP,PRCP,
     &  STC,DT,RNET,H,AET,RES,IREC3)
C  
C WRITE ATNOSPHERIC INPUT FORCING DATA AND INDEPENDENT VALIDATION DATA 
C IN FILE (BND_DATA.GRS or .TXT)
C
      CALL PRTBND(INDI,NOUT5,jday,
     & IIDAY,  ITIME,	  SFCTMP,   Q2,     SFCPRS,   Rg,
     * LW_in, Par_in,  Par_out,  RNET, S_OBS, PRCP, wet, T1_OBS,
     * T_2cm,	T_4cm,    T_8cm,    T_16cm, T_32cm,   T_64cm, sm_5cm,
     * sm_20cm, sm_60cm,  w_dir,    u_bar,  eddyuw,
     * uprim2,  vprim2,   wprim2,   H_OBS,  ETA_OBS,
     * IMONTH,  IDAY,	  IREC5)
C
C =============================================================
      END DO
C							               |
C   END OF TIME LOOP					               |
C 							               |
C ---------------------------------------------------------------------|
C ---------------------------------------------------------------------|
	
      STOP 0
      END
      SUBROUTINE CANRES(SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL,
     &                  SMCWLT,SMCREF,RCMIN,RC,PC,NROOT,Q2SAT,DQSDT2, 
     &                  TOPT,RSMAX,RGL,HS,XLAI)

      IMPLICIT NONE

C **********************************************************************
C			SUBROUTINE CANRES
C			-----------------
C       THIS ROUTINE CALCULATES THE CANOPY RESISTANCE WHICH DEPENDS
C       ON INCOMING SOLAR RADIATION, AIR TEMPERATURE, ATMOSPHERIC
C       WATER VAPOR PRESSURE DEFICIT AT THE LOWEST MODEL LEVEL, ANS SOIL
C       MOISTURE.
C ----------------------------------------------------------------------
C	SOURCE:  JARVIS (1976), JACQUEMIN AND NOILHAN (1990 BLM)
C ----------------------------------------------------------------------
C ----------------------------------------------------------------------
C	INPUT:   SOLAR: INCOMING SOLAR RADIATION
C                CH: SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
C                SFCTMP: AIR TEMPERATURE AT 1ST LEVEL
C                Q2:  AIR HUMIDITY AT 1ST LEVEL
C                SFCPRS: SURFACE PRESSURE
C                SMC: VOLUMETRIC SOIL MOISTURE 
C                ZSOIL: SOIL DEPTH
C                NSOIL: NUMBER OF SOIL LAYERS
C                SMCWLT: WILTING POINT
C                SMCREF: REFERENCE SOIL MOISTURE (FIELD CAPACITY)
C
C	OUTPUT:  PC: PLANT COEFFICIENT
C                RC: CANOPY RESISTANCE
C ----------------------------------------------------------------------
C **********************************************************************

      INTEGER   NSOLD, NROOT, NSOIL, K
      PARAMETER (NSOLD = 4)

      REAL SIGMA, RD, CP, SLV
      REAL SOLAR, CH, SFCTMP, Q2, SFCPRS 
      REAL SMC(NSOIL), ZSOIL(NSOIL), PART(NSOLD) 
      REAL SMCWLT, SMCREF, RCMIN, RC, PC, Q2SAT, DQSDT2
      REAL TOPT, RSMAX, RGL, HS, XLAI, RCS, RCT, RCQ, RCSOIL, FF
      REAL P, QS, GX, TAIR4, ST1, SLVCP, RR, DELTA

      PARAMETER (SIGMA=5.67E-8, RD=287.0, CP=1004.5, SLV=2.5E6)

      RCS = 0.0
      RCT = 0.0
      RCQ = 0.0
      RCSOIL = 0.0
      RC = 0.0

C ----------------------------------------------------------------------
C CONTRIBUTION DUE TO INCOMING SOLAR RADIATION.
C ----------------------------------------------------------------------

CC/98/01/05/       FF = 0.55*2.0*SOLAR/RGL
      FF = 0.55*2.0*SOLAR/(RGL*XLAI)
      RCS = (FF + RCMIN/RSMAX) / (1.0 + FF)
      RCS = AMAX1(RCS,0.0001)

C ----------------------------------------------------------------------
C CONTRIBUTION DUE TO THE TEMPERATURE AT THE FIRST MODEL LEVEL.
C ----------------------------------------------------------------------

      RCT = 1.0 - 0.0016*((TOPT-SFCTMP)**2.0)
      RCT = AMAX1(RCT,0.0001)

C ----------------------------------------------------------------------
C CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT THE FIRST MODEL LEVEL.
C ----------------------------------------------------------------------

      P = SFCPRS
      QS = Q2SAT
C RCQ EXPRESSION FROM SSIB 
      RCQ = 1.0/(1.0+HS*(QS-Q2))
      RCQ = AMAX1(RCQ,0.01)

C ----------------------------------------------------------------------
C CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
C DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
C ----------------------------------------------------------------------

      GX = (SMC(1) - SMCWLT) / (SMCREF - SMCWLT)
      IF (GX .GT. 1.) GX = 1.
      IF (GX .LT. 0.) GX = 0.

C****   USING SOIL DEPTH AS WEIGHTING FACTOR
      PART(1) = (ZSOIL(1)/ZSOIL(NROOT)) * GX

C**** USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
CC      PART(1) = RTDIS(1) * GX
      
      DO 10 K = 2, NROOT
        GX = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT)
        IF (GX .GT. 1.) GX = 1.
        IF (GX .LT. 0.) GX = 0.
C****   USING SOIL DEPTH AS WEIGHTING FACTOR        
        PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT)) * GX

C**** USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
CC         PART(K) = RTDIS(K) * GX 
	       
10    CONTINUE

      DO 20 K = 1, NROOT
        RCSOIL = RCSOIL+PART(K)

20    CONTINUE
        RCSOIL = AMAX1(RCSOIL,0.0001)

C ----------------------------------------------------------------------
C         DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS.
C         CONVERT CANOPY RESISTANCE (RC) TO PLANT COEFFICIENT (PC).
C ----------------------------------------------------------------------

CC/98/01/05/          RC = RCMIN/(RCS*RCT*RCQ*RCSOIL)
          RC = RCMIN/(XLAI*RCS*RCT*RCQ*RCSOIL)

          TAIR4 = SFCTMP**4
          ST1 = (4.*SIGMA*RD)/CP
          SLVCP = SLV/CP
          RR = ST1*TAIR4/(SFCPRS*CH) + 1.0
          DELTA = SLVCP*DQSDT2

	  PC = (RR+DELTA)/(RR*(1.+RC*CH)+DELTA)

          RETURN
          END
      FUNCTION CSNOW ( DSNOW )
      
      IMPLICIT NONE

      REAL  C
      REAL  DSNOW
      REAL  CSNOW

      REAL UNIT
      PARAMETER ( UNIT=0.11631 ) 
                                         
C   ****  SIMULATION OF TERMAL SNOW CONDUCTIVITY                   
C   ****  SIMULATION UNITS OF CSNOW IS CAL/(CM*HR* C) 
C   ****  AND IT WILL BE RETURND IN W/(M* C)
C   ****  BASIC VERSION IS DYACHKOVA EQUATION                                

C *****   DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4

      C=0.328*10**(2.25*DSNOW)
      CSNOW=UNIT*C

       
C *****    DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
C       CSNOW=0.0293*(1.+100.*DSNOW**2)
      
C     *****   E. ANDERSEN FROM FLERCHINGER
C     CSNOW=0.021+2.51*DSNOW**2        
      
      RETURN                                                      
      END  
      FUNCTION DEVAP ( ETP1, SMC, ZSOIL, SHDFAC, SMCMAX, B,
     &                 DKSAT, DWSAT, SMCDRY, SMCREF, SMCWLT)

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    NAME:  DIRECT EVAPORATION (DEVAP) FUNCTION  VERSION: N/A
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL B
      REAL DEVAP
      REAL DKSAT
      REAL DWSAT
      REAL ETP1
      REAL FX
      REAL SHDFAC
      REAL SMC
      REAL SMCDRY
      REAL SMCMAX
C      REAL WCND
C      REAL WDF
      REAL ZSOIL
      REAL SMCREF
      REAL SMCWLT

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     GET SOIL H2O DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY GIVEN SMC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C FX      CALL WDFCND ( WDF, WCND, SMC, SMCMAX, B, DKSAT, DWSAT )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FX REPRESENTS SOIL MOISTURE FLUX / ATMOSPHERIC DEMAND
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C FX     FX = ( -2. * WDF * ( SMC - SMCDRY ) / ZSOIL - WCND ) / ETP1
C*** TEST OF BETA METHODE, F. CHEN 7/11/96 ***
C Change 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! formulation to reduce 
C excessive evaporation under little green vegetation:

C      FX = (SMC - SMCWLT) / (SMCREF - SMCWLT)
C      FX = (SMC - SMCWLT) / (SMCMAX - SMCWLT)
      FX = (SMC - SMCDRY) / (SMCMAX - SMCDRY)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FX > 1 REPRESENTS DEMAND CONTROL
C     FX < 1 REPRESENTS FLUX CONTROL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FX = MAX ( MIN ( FX, 1. ) ,0. )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DEVAP = FX * ( 1.0 - SHDFAC ) * ETP1
      RETURN
      END
      FUNCTION DQS (T) 

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  PURPOSE:  TO CALCULATE VALUES OF VAPOR PRESSURE (E)
CC            AND P * DQS/DT (P TIMES CHG IN SAT MXG RATIO WITH RESPECT
CC            TO THE CHG IN TEMP) IN SUBSTITUTION TO THE LOOK-UP TABLES.
CC
CC            SUBSTITUTES LOOK-UP TABLES ASSOCIATED WITH THE DATA
CC            BLOCK  /CHMXR/ .
CC
CC            FORMULAS AND CONSTANTS FROM ROGERS AND YAU, 1989.
CC
CC                         ADDED BY PABLO J. GRUNMANN, 6/30/97.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      REAL  DESDT
      REAL  DQS
CK    REAL  ESD
      REAL  LW
      REAL  T
      REAL  ES
C
CK    REAL CP
CK    REAL CV
CK    REAL CVV
      REAL CPV
      REAL RV
      REAL CW
      REAL EPS
      REAL ESO
      REAL TO
C
CK    PARAMETER  (CP = 1005.)
CK    PARAMETER  (CV = 718.)
CK    PARAMETER (CVV = 1410.)
      PARAMETER (CPV = 1870.)
      PARAMETER  (RV = 461.5)
      PARAMETER  (CW = 4187.)
      PARAMETER (EPS = 0.622)
      PARAMETER (ESO = 611.2)
      PARAMETER  (TO = 273.15)
C 
C     ABOUT THE PARAMETERS:
C      
C     EPS ---------- WATER - DRY AIR MOLECULAR MASS RATIO, EPSILON
C      
C   VALUES FOR SPECIFIC HEAT CAPACITY AND INDIVIDUAL GAS CONSTANTS 
C   IN [JOULES/(KG*KELVIN)] UNITS.
C
C     DRY AIR: 
C             CP, CV
C     WATER VAPOR:
C                 CVV = 1410. 
C                 CPV = 1870.
C                 RV  =  461.5
C     LIQUID WATER:
C                  CW = 4187.
C
C     ESO = ES(T=273.15 K) = SAT. VAPOR PRESSURE (IN PASCAL) AT T=TO
C      TO = 273.15
C      
C     SAT. MIXING  RATIO: QS ~= EPS*ES/P
C     CLAUSIUS-CLAPEYRON: DES/DT = L*ES/(RV*T**2)
C     @QS/@T =  (EPS/P)*DES/DT
C    
          LW = 2.501000E6 - ( CW - CPV ) * ( T - TO )
          ES = ESO*EXP (LW*(1/TO - 1/T)/RV)  
          DESDT = LW*ES/(RV*T*T)
C
C      FOR INSERTION IN DQSDT FUNCTION: 
C      DQSDT = DQS/P , WHERE DQS = EPS*DESDT  
C
          DQS = EPS*DESDT
C     
          RETURN
          END
      FUNCTION DQSDT ( SFCTMP, SFCPRS )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO RETRIEVE THE APPROPRIATE VALUE OF DQSDT (THE CHANGE
CC    =======   OF THE SATURATION MIXING RATIO WITH RESPECT TO THE 
CC              CHANGE IN TEMPERATURE) FROM: 
CC
CC               FORMULAS INTRODUCED IN NEW FUNCTION DQS 
CC                                  (MODIFIED BY PABLO GRUNMANN, 7/9/97).
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL  SFCTMP
      REAL  SFCPRS
      REAL  DQS
      REAL  DQSDT

      IF ((SFCTMP .GE. 173.0) .AND. (SFCTMP  .LE.  373.0)) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       IF THE INPUT SFC AIR TEMP IS BTWN 173 K AND 373 K, USE
C       FUNCTION DQS TO DETERMINE THE SLOPE OF SAT.MIX RATIO FUNCTION
C                                 -ADAPTED TO USE NEW DQS, 7/9/97.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        DQSDT = DQS (SFCTMP) / SFCPRS

      ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       OTHERWISE, SET DQSDT EQUAL TO ZERO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        DQSDT = 0.0

      END IF

      RETURN
      END
      FUNCTION E (T) 

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  PURPOSE:  TO CALCULATE VALUES OF SAT. VAPOR PRESSURE (E)
CC            SUBSTITUTES LOOK-UP TABLES ASSOCIATED WITH THE DATA
CC            BLOCK  /VAPPRS/ .
CC            FORMULAS AND CONSTANTS FROM ROGERS AND YAU, 1989.
CC
CC                         ADDED BY PABLO J. GRUNMANN, 7/9/97.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      REAL  LW
      REAL  T
      REAL  E
C
CK    REAL  EPS
CK    REAL  CP
CK    REAL  CV
CK    REAL  CVV
      REAL  CPV
      REAL  RV
      REAL  CW
      REAL  ESO
      REAL  TO
C
CK    PARAMETER (EPS = 0.622)
CK    PARAMETER  (CP = 1005.)
CK    PARAMETER  (CV = 718.)
CK    PARAMETER (CVV = 1410.)
      PARAMETER (CPV = 1870.)
      PARAMETER  (RV = 461.5)
      PARAMETER  (CW = 4187.)
      PARAMETER (ESO = 611.2)
      PARAMETER  (TO = 273.15)      
C 
C     ABOUT THE PARAMETERS:
C      
C     EPS ---------- WATER - DRY AIR MOLECULAR MASS RATIO, EPSILON
C      
C   VALUES FOR SPECIFIC HEAT CAPACITY AND INDIVIDUAL GAS CONSTANTS 
C   IN [JOULES/(KG*KELVIN)] UNITS.
C
C     DRY AIR: 
C             CP, CV
C     WATER VAPOR:
C                 CVV = 1410. 
C                 CPV = 1870.
C                 RV  =  461.5
C     LIQUID WATER:
C                  CW = 4187.
C
C     ESO = ES(TO) = SAT. VAPOR PRESSURE (IN PASCAL) AT T=TO
C      TO = 273.15
C_______________________________________________________________________
C
C     CLAUSIUS-CLAPEYRON: DES/DT = L*ES/(RV*T**2)
C     
C    
          LW = 2.501000D6 - ( CW - CPV ) * ( T - TO )
          E = ESO*EXP (LW*(1/TO - 1/T)/RV)  
C
     
          RETURN
          END 
      FUNCTION FRH2O(T,SMC,SH2O,SMCMAX,B,PSIS)
 
      IMPLICIT NONE
     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      
CC  PURPOSE:  CALCULATE FREE WATER CONTENT IF TEMPERATURE 
CC  =======   IS BELOW 273.16K (T0). 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL FRH2O, HLICE, GS, CK, DICE, DH2O, T0, ERROR, T
      REAL SMC, SH2O, SMCMAX, B, PSIS, BX, TR, SWL
      REAL SMCM, SMCC, DF, DENOM, SWLK, DSWL

      INTEGER N

      PARAMETER (HLICE=3.335 E5)      
      PARAMETER (GS = 9.81)
      PARAMETER (CK=8.)
      PARAMETER (DICE=920.0)
      PARAMETER (DH2O=1000.0)
      PARAMETER (T0=273.16)
      PARAMETER (ERROR=0.005)
       
C  ***             LIMITS ON PARAMETER B: B < 5.5              ****
C  ***   SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT  ****
C  ***   IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES    ****
C******************************************************************
C
       BX = B
       IF ( B .GT. 5.5 ) BX = 5.5
C------------------------------------------------------------------
        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  ASSUMED THAT FREE WATER CONTENT EQUALS ITS VALUE AT TEMPERATURE 
CC  273.11K IN RANGE OF TEMPERATURES 273.11K - 273.16K TO REDUCE
CC  ITERRATINS PROBLEM IN THIS RANGE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      

      TR=T
      IF(T .GT. 273.11 .AND. T .LE. T0 ) TR=273.11

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  IF TEMPERATURE ABOVE 273.16K SH2O = SMC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       
      IF(TR .GT. T0) THEN
        FRH2O=SMC
        GOTO 77
      ENDIF  

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  CONVERT SOIL MOISTURE VARIABLES INTO %
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
      SWL=(SMC-SH2O)*100.
      SMCM=SMCMAX*100.
      SMCC=SMC*100.                       

      N=0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  START OF ITERRATIONS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

   10  DF = ( -PSIS*GS/HLICE ) * ( ( 1.+CK*SWL*0.01 )**2. ) *
     +      ( SMCM/(SMCC-SWL ) )**BX
        
       DENOM = 2. * CK / ( 1.+CK*SWL*0.01 ) + BX / ( SMCC - SWL )
       SWLK = SWL - (1.-(TR-T0)/(TR*DF))/DENOM
 
       DSWL=ABS(SWLK-SWL)
       SWL=SWLK
       IF(SWL.GE.SMCC) THEN
C        WRITE(*,*) 'WARN2 ',N,SWL,SMCC
       SWL=SMCC*0.8
        SMCC=SMCC*0.97
       ENDIF
       
       N=N+1 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  EXIT IF MORE THAN 900 ITERRATIONS. FREE WATER CONTENT
CC  WILL NOT BE CHANGED
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
      IF(N .GT. 900 ) THEN
        IF(TR . LT. 271.16) THEN
          SWL = ( SMC-0.02 ) * 100.
        ELSE
          SWL=0.
        ENDIF    
        GOTO 20
      ENDIF
        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  END OF ITERRATIONS IF LAST ESTIMATE DIFFERENCE IS LESS THAN 
CC  AN ERROR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
         
      IF ( DSWL .GT. ERROR ) GOTO 10
      
   20 IF(SWL .LT. 0.)  THEN
C        WRITE(*,*) ' SWL < 0. ',SWL
        SWL=0.
      ENDIF  

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  RECALCULATE ICE CONTENT IN % INTO VOLUMETRIC FREE WATER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c this eq was using MAX1(), which yields an integer result
c changed to MAX() - Pablo Grunmann 12-07-98. - !!!!!!!!!!!!
      FRH2O = MAX ( (SMC - SWL*0.01 ) ,  0.02 )

C  ***        FLERCHINGER EQUATION             *******      
CC      FK=(((-HLICE/(GS*PSIS))*((TR-T0)/TR))**(-1/B))*SMCMAX
CC      FRH2O = AMIN1 ( FK, SMC )
CC      IF(FRH2O .LT. 0.) FRH2O = 0.02
C-----------------------------------------------------

   77 CONTINUE
      RETURN
      END      
      SUBROUTINE HRT ( RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,
     +                 SMCWLT, PSISAT, SH2O, DT, B, F1, DF1, QUARTZ )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY
CC    =======   TERM OF THE SOIL THERMAL DIFFUSION EQUATION.  ALSO TO
CC              COMPUTE ( PREPARE ) THE MATRIX COEFFICIENTS FOR THE
CC              TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOLD
      PARAMETER ( NSOLD = 4 )
      INTEGER NSOIL
      INTEGER K
      REAL AI    ( NSOLD )
      REAL BI    ( NSOLD )
      REAL CAIR
      REAL CH2O
      REAL CI    ( NSOLD )
      REAL CSOIL
      REAL DDZ
      REAL DDZ2
      REAL DENOM
      REAL DFMAX1
      REAL DF1
      REAL DFMAXN
      REAL DTSDZ
      REAL DTSDZ2
      REAL F1
      REAL HCPCT
c      REAL MXSMC
c      REAL MXICE
      REAL QUARTZ
      REAL RHSTS ( NSOIL )
      REAL S
      REAL SMC   ( NSOIL )

      REAL SH2O  ( NSOIL )
      REAL SMCMAX
      REAL SMCWLT
      
      REAL STC   ( NSOIL )
      REAL TBOT
      REAL YY
      REAL ZBOT
      REAL ZSOIL ( NSOIL )
      REAL ZZ1

      REAL CICE, T0, TSURF, PSISAT, DT, B, SICE, TBK, TSNSR, TBK1
      REAL TBND, SNKSRC

      INTEGER I

CMIC$ TASKCOMMON ABCI
      COMMON /ABCI/ AI, BI, CI
      PARAMETER(CAIR=1004.0,CH2O=4.2 E6,CSOIL=1.26E6,ZBOT=-3.0)

C   THIS COMMON STATEMENT IS NEED ONLY TO PRINT ESTIMATED SOIL 
C   SURFACE TEMPERATURE, AND CAN BE OMITTED 
      COMMON /PRN_TSOIL/ TSURF
            
C ----------    FROZEN GROUND VERSION    ----------------------------
C    CICE IS ICE HEAT CAPACITY
C
      PARAMETER ( CICE = 2.106 E6, T0 = 273.16 )       
C -------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     RETRIEVE THE THERMAL DIFFUSIVITY FOR THE WETTER OF THE TOP SOIL
C     LAYER OR THE NEXT LOWER SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C ----------    FROZEN GROUND VERSION    ----------------------------
C   SOIL THERMAL CONDUCTIVITY ADJUSTMENT IF IT IS FROZEN
C
      SICE = 0.
      DO I = 1, NSOIL
        SICE = SICE + SMC(I) - SH2O(I)
      END DO

C
C
C
C
C      IF( SMC(1) .GE. SMC(2)) THEN
C	  MXSMC = SMC(1)
C	  MXICE = SMC(1) - SH2O(1)
C      ELSE
C	  MXSMC = SMC(2)
C	  MXICE = SMC(2) - SH2O(2)
c      ENDIF
      CALL TDFCND ( DFMAX1, SMC(1), B, F1,QUARTZ,SMCMAX,SH2O(1))
C__________________________________________________________________
C This IF-statement is no longer used if Peters-Lidard thermal
C conductivity is used for all conditions (frozen,unfrozen):
C      IF ( MXICE .GT. 0.0 ) CALL FRZCND ( DFMAX1, MXICE )
C______________________________________________P.Grunmann_01/99____

C ------------------------------------------------------------
C      MXSMC = MAX ( SMC(1), SMC(2) )
C      CALL TDFCND ( DFMAX1, MXSMC, B, F1 )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      HCPCT = SH2O(1)*CH2O + (1.0-SMCMAX)*CSOIL + (SMCMAX-SMC(1))*CAIR
C VK correction: SMC to SH2O, first term on R.H. side.-P.Grunmann 1/99

C ----------    FROZEN GROUND VERSION    ----------------------------
C   ICE HEAT CAPACUTY IS INCLUDED
C      
     +        + ( SMC(1) - SH2O(1) )*CICE
C -------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
      AI(1) = 0.0
      CI(1) =  ( DFMAX1 * DDZ ) / ( ZSOIL(1) * HCPCT )
      BI(1) = -CI(1) + DF1 / ( 0.5 * ZSOIL(1) * ZSOIL(1)*HCPCT*ZZ1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL
C     LAYERS.  RECALC/ADJUST THE SOIL HEAT FLUX.  USE THE GRADIENT
C     AND FLUX TO CALC RHSTS FOR THE TOP SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C ------------    FROZEN GROUND VERSION    --------------------------
C   SOIL SURFACE TEMPERATURE ADJUSTMENT BY ZZ1 ( THE SAME AS T1 )
C
      TSURF = ( YY + ( ZZ1 - 1 ) * STC(1) ) / ZZ1
C -------------------------------------------------------------------
      
      DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) )
      S = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 )
      RHSTS(1) = ( DFMAX1 * DTSDZ - S ) / ( ZSOIL(1) * HCPCT )

C -----------    FROZEN GROUND VERSION     --------------------------
C   CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
C   AND SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR.
C
      IF ( SICE .NE. 0. .OR. TSURF .LT. T0 ) THEN 
       TBK = TBND ( STC(1), STC(2), ZSOIL, ZBOT, 1, NSOIL )
       TSNSR = SNKSRC ( TSURF, STC(1),TBK, SMC(1), SH2O(1), HCPCT, 
     +             ZSOIL,SMCMAX, PSISAT, B, F1, DT, 1, QUARTZ)
       RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
      ENDIF 
C -------------------------------------------------------------------
            
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     INITIALIZE DDZ2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DDZ2 = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO K = 2, NSOIL

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC THIS SOIL LAYER'S HEAT CAPACITY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        HCPCT = SH2O(K)*CH2O +(1.0-SMCMAX)*CSOIL +(SMCMAX-SMC(K))*CAIR
C VK correction: SMC to SH2O, first term on R.H. side.-P.Grunmann 1/99

C ----------    FROZEN GROUND VERSION    ----------------------------
C   ICE HEAT CAPACUTY IS INCLUDED
C      
     +        + ( SMC(K) - SH2O(K) )*CICE
C -------------------------------------------------------------------   

        IF ( K .NE. NSOIL ) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         RETRIEVE THE THERMAL DIFFUSIVITY FOR THE WETTER OF THE
C         CURRENT OR THE NEXT LOWER SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          MXSMC = MAX ( SMC(K), SMC(K+1) )

C ------------    FROZEN GROUND VERSION    --------------------------
C   SOIL THERMAL DIFFUSIVITY ADJUSTMENT IF SOIL IS FROZEN
C
C      IF ( SMC(K) .GE. SMC(K+1) ) THEN
C	  MXSMC = SMC(K)
C	  MXICE = SMC(K) - SH2O(K)
C      ELSE
C	  MXSMC = SMC(K+1)
C	  MXICE = SMC(K+1) - SH2O(K+1)
C      ENDIF
CC	CALL TDFCND ( DFMAXN, MXSMC, B, F1,QUARTZ,SMCMAX,SH2O(K))
      CALL TDFCND ( DFMAXN, SMC(K), B, F1,QUARTZ,SMCMAX,SH2O(K))
CC__________________________________________________________________
C This IF-statement is no longer used if Peters-Lidard thermal
C conductivity is used for all conditions (frozen,unfrozen):
C      IF ( MXICE .GT. 0.0 ) CALL FRZCND ( DFMAXN, MXICE )
C -------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
          DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
          CI(K) = -DFMAXN * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)

C -----------    FROZEN GROUND VERSION     --------------------------
C   CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
C
          TBK1 = TBND ( STC(K), STC(K+1), ZSOIL, ZBOT, K, NSOIL )
C -------------------------------------------------------------------
          
        ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         RETRIEVE THE THERMAL DIFFUSIVITY FOR THE LOWEST SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C -----------    FROZEN GROUND VERSION     --------------------------
C   SOIL THERMAL DIFFUSIVITY ADJUSTMENT
C   AND  CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
C
C          MXICE = SMC(K) - SH2O(K)
          CALL TDFCND ( DFMAXN, SMC(K), B, F1,QUARTZ,SMCMAX,SH2O(K))
C__________________________________________________________________
C This IF-statement is no longer used if Peters-Lidard thermal
C conductivity is used for all conditions (frozen,unfrozen):
C          IF ( MXICE .GT. 0.0 ) CALL FRZCND ( DFMAXN, MXICE )
          
          TBK1 = TBND ( STC(K), TBOT, ZSOIL, ZBOT, K, NSOIL )
C -------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DTSDZ2 = (STC(K)-TBOT) / (.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         SET MATRIX COEF, CI TO ZERO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          CI(K) = 0.
        END IF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
        RHSTS(K) = ( DFMAXN * DTSDZ2 - DFMAX1 * DTSDZ ) / DENOM

C -----------    FROZEN GROUND VERSION     --------------------------
C   SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR.
C
      IF ( SICE .NE. 0. .OR. TBK .LT. 0. ) THEN
       TSNSR = SNKSRC ( TBK, STC(K),TBK1, SMC(K), SH2O(K), HCPCT, 
     +             ZSOIL,SMCMAX, PSISAT, B, F1, DT, K, QUARTZ)
       RHSTS(K) = RHSTS(K) - TSNSR / DENOM
      ENDIF 
C -------------------------------------------------------------------
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        AI(K) = - DFMAX1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
        BI(K) = -(AI(K) + CI(K))

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C RESET VALUES OF DFMAX1, DTSDZ, AND DDZ FOR LOOP TO NEXT SOIL LYR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C -----------    FROZEN GROUND VERSION     --------------------------
C    REMEMBER UPPER LAYER TEMPERATURE
C
        TBK = TBK1
C -------------------------------------------------------------------
        
        DFMAX1 = DFMAXN
        DTSDZ = DTSDZ2
        DDZ = DDZ2
C   40   CONTINUE
      END DO
      RETURN
      END
      SUBROUTINE HRTICE (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,
     &                 TBOT, B, F1, DF1 )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY
CC    =======   TERM OF THE SOIL THERMAL DIFFUSION EQUATION.  ALSO TO
CC              COMPUTE ( PREPARE ) THE MATRIX COEFFICIENTS FOR THE
CC              TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOLD
      PARAMETER ( NSOLD = 4 )
      INTEGER NSOIL
      REAL AI    ( NSOLD )
      REAL BI    ( NSOLD )
      REAL CAIR
      REAL CH2O
      REAL CI    ( NSOLD )
      REAL CSOIL
      REAL DDZ
      REAL DDZ2
      REAL DENOM
      REAL DFMAX1
      REAL DF1
      REAL DFMAXN
      REAL DTSDZ
      REAL DTSDZ2
      REAL B
      REAL F1
      REAL HCPCT
      INTEGER K
CC      REAL MXSMC
      REAL RHSTS ( NSOIL )
      REAL S
      REAL SMC   ( NSOIL )
      REAL SMCMAX
      REAL STC   ( NSOIL )
      REAL TBOT
      REAL YY
      REAL ZBOT
      REAL ZSOIL ( NSOIL )
      REAL ZZ1
CMIC$ TASKCOMMON ABCI
      COMMON /ABCI/ AI, BI, CI
      PARAMETER(CAIR=1004.0,CH2O=4.2 E6,CSOIL=1.26E6)

      ZBOT=ZSOIL(NSOIL)
      DFMAX1=2.2
      HCPCT=1880.0*917.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
      AI(1) = 0.0
      CI(1) =  ( DFMAX1 * DDZ ) / ( ZSOIL(1) * HCPCT )
      BI(1) = -CI(1) + DF1/( 0.5 * ZSOIL(1) * ZSOIL(1) * HCPCT * ZZ1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL
C     LAYERS.  RECALC/ADJUST THE SOIL HEAT FLUX.  USE THE GRADIENT
C     AND FLUX TO CALC RHSTS FOR THE TOP SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) )
      S = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 )
      RHSTS(1) = ( DFMAX1 * DTSDZ - S ) / ( ZSOIL(1) * HCPCT )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     INITIALIZE DDZ2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DDZ2 = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 40 K = 2, NSOIL

        IF ( K .NE. NSOIL ) THEN
          DFMAXN=2.2

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
          DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
          CI(K) = -DFMAXN * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)

        ELSE
          DFMAXN=2.2

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DTSDZ2 = (STC(K)-TBOT)/(.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         SET MATRIX COEF, CI TO ZERO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          CI(K) = 0.
        END IF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
        RHSTS(K) = ( DFMAXN * DTSDZ2 - DFMAX1 * DTSDZ ) / DENOM

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        AI(K) = - DFMAX1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
        BI(K) = -(AI(K) + CI(K))

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       RESET VALUES OF DFMAX1, DTSDZ, AND DDZ FOR LOOP TO NEXT SOIL LYR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        DFMAX1 = DFMAXN
        DTSDZ = DTSDZ2
        DDZ = DDZ2
   40   CONTINUE
      RETURN
      END
      SUBROUTINE HSTEP ( STCOUT, STCIN, RHSTS, DT, NSOIL )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOLD
      PARAMETER ( NSOLD=4 )
      INTEGER K
      INTEGER NSOIL

      REAL AI    ( NSOLD )
      REAL BI    ( NSOLD )
      REAL CI    ( NSOLD )
      REAL DT
      REAL RHSTS ( NSOIL )

C -----------    FROZEN GROUND VERSION     -------------------------
C      
      REAL STCOUT   ( NSOIL )
      REAL STCIN  ( NSOIL )
C ------------------------------------------------------------------      
      
CMIC$ TASKCOMMON ABCI
      COMMON /ABCI/ AI, BI, CI

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 10 K = 1 , NSOIL
        RHSTS(K) = RHSTS(K) * DT
        AI(K) = AI(K) * DT
        BI(K) = 1. + BI(K) * DT
        CI(K) = CI(K) * DT
 10   CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SOLVE THE TRI-DIAGONAL MATRIX EQUATION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CALL ROSR12 ( CI,AI,BI,CI,RHSTS,RHSTS,NSOIL )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 20 K = 1 , NSOIL
      
C -----------    FROZEN GROUND MODIFICATION    ----------------------
C   NEW TEMPERATURES IN STCOUT FROM OLD STCIN
C       
        STCOUT(K) = STCIN(K) + CI(K)
C -------------------------------------------------------------------
C        STC(K) = STC(K) + CI(K)
        
  20  CONTINUE
      RETURN
      END
      SUBROUTINE JULDATE(JULD,IMONTH,IDAY,JULM)
      
      IMPLICIT NONE

CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC    NAME:  JULIAN DAYS TO MM/DD DATE
CC
CC
CC    PURPOSE:  TO CONVERT FROM JULIAN DAY NUMBER OF THE YEAR TO
CC              CALENDAR DATE  (MONTH, DAY).
CC                                      PABLO J. GRUNMANN, 05/98
CC
CC    VARIABLES:
CC    =========
CC
CC    LABEL	    ............DESCRIPTION...............
CC
CC    IDMONTH(IM)   NUMBER OF DAYS IN MONTH NUMBER "IM", FOR IM=1,12
CC    JULD	    JULIAN DAY
CC    ILDM,ILDMP    JULIAN DAY OF THE LAST DAY OF A PARTICULAR 
CC                  MONTH AND ITS CONSECUTIVE (ILDMP).
CC    IMONTH,IDAY   MONTH, DAY - CALENDAR DATE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INTEGER IM
      INTEGER IDMONTH(12)
      INTEGER JULM(13)
      INTEGER JULD
      INTEGER IDAY
      INTEGER ILDM
      INTEGER ILDMP
      INTEGER IMONTH     
C.................................................................        
C  DEFINE LAST DAY OF FEBRUARY AND MODIFY DATA IDMONTH 
C IF NECESSARY:
C  	       FEB 28
C  	       FEB 29 (LEAP YEAR)
C..................JAN,FEB,MAR,APR,MAY,JUN,JUL,AUG.SEP,OCT,NOV,DEC
      DATA IDMONTH/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
C      DATA JULM/0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 
C     & 334, 365/
C.................................................................        
      ILDMP = 0
      IMONTH = 0
      
C      Print*,JULM
      JULM(1)= 0
      
      DO IM = 1,12
        JULM(IM+1) = JULM(IM) + IDMONTH(IM)
	
        ILDM = JULM(IM) 
        ILDMP = JULM(IM+1)
	IF ((JULD .LE. ILDMP) .AND. (JULD .GT. ILDM)) THEN
	  IMONTH = IM
	  IDAY = JULD - JULM(IM)
	ENDIF
      END DO 
	
      RETURN
      END
      SUBROUTINE MONTH_D(JDAY,XMON,XDAY)
      
      IMPLICIT NONE
      
C    INTERPOLATE TO DAY OF YEAR.
C    THE MONTHLY FIELDS ASSUMED VALID AT 15TH OF MONTH.
C    
C    READ FROM A FILE THAT HAS 13 RECORDS, ONE PER MONTH WITH
C    JANUARY OCCURRING BOTH AS RECORD 1 AND AGAIN AS RECORD 13,
C    THE LATTER TO SIMPLIFY TIME INTERPOLATION FOR DAYS
C    BETWEEN DEC 16 AND JAN 15. WE TREAT JAN 1 TO JAN 15
C    AS JULIAN DAYS 366 TO 380 BELOW, I.E WRAP AROUND YEAR.
C
C   BASED ON SUBROUTINE CNSTS.f WHICH
C   INCLUDED REVISION BY F. CHEN 7/96 TO REFLECT A NEW NESDIS 
C   VEGETATION FRACTION PRODUCT (FIVE-YEAR CLIMATOLOGY WITH 
C   0.144 DEGREE RESOLUTION FROM 89.928S, 180W TO 89.928N, 180E)
C   
C   This is PABLO J. GRUNMANN's version for use with the uncoupled, 
C   column mode Land-Surface Model
C 

C      DATA JULM/0,31,59,90,120,151,181,212,243,273,304,334,365/
C.......getting JULM from subrooutine JULDATE to avoid
C.......problems with leap years.

C
C ****  DO TIME INTERPOLATION ****
C            
       INTEGER IDAY
       INTEGER IMON1
       INTEGER IMON2
       INTEGER IMONTH
       INTEGER JDAY
       INTEGER JULD
       
       INTEGER JULM(13)

      REAL DAY1
      REAL DAY2
      REAL RDAY
      REAL WGHT1
      REAL WGHT2
      REAL XDAY
      
      REAL XMON(13)

       JULD = JDAY
       CALL JULDATE(JULD,IMONTH,IDAY,JULM)
       
       
       IF(JULD.LE.15) JULD=JULD+365
       
       IMON2=IMONTH
       IF (IDAY .GT. 15) IMON2 = IMON2 + 1
       IF (IMON2 .EQ.  1) IMON2 = 13
       
       IMON1 = IMON2 - 1
       
C **** ASSUME DATA VALID AT 15TH OF MONTH 
       DAY2 = REAL(JULM(IMON2)+15)
       DAY1 = REAL(JULM(IMON1)+15)
       RDAY = REAL(JULD)
       WGHT1 = (DAY2-RDAY)/(DAY2-DAY1)
       WGHT2 = (RDAY-DAY1)/(DAY2-DAY1)
C       
         XDAY=WGHT1*XMON(IMON1)+WGHT2*XMON(IMON2)
	 
            
      RETURN
      END
      SUBROUTINE NOPAC ( ETP, ETA, PRCP, SMC, SMCMAX, SMCWLT,
     &                   SMCREF,SMCDRY,CMC,CMCMAX, NSOIL, DT, SHDFAC,
     &                   Q1, Q2, T1, SFCTMP, T24, TH2, F, F1, S, STC,
     &                   EPSCA, B, PC, RCH, RR,  CFACTR, 
     &                   SH2O, SLOPE, KDT, FRZFACT, PSISAT,
     &                   ZSOIL, DKSAT, DWSAT, TBOT, RUNOFF1, RUNOFF2,
     &                   RUNOFF3, EDIR1, EC1, ETT1, NROOT, ICE,RTDIS,
     &                   QUARTZ)

      IMPLICIT NONE

C -----------   FROZEN GROUND VERSION     ---------------------------
C   NEW STATES ADDED:  SH2O,  AND FROZEN GROUND CORRECTION FACTOR
     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE
CC    =======   SOIL MOISTURE CONTENT AND SOIL HEAT CONTENT VALUES FOR
CC              THE CASE WHEN NO SNOW PACK IS PRESENT.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOIL

      REAL B
      REAL BETA
      REAL CFACTR
      REAL CMC
      REAL CMCMAX
      REAL CP
      REAL DEW
      REAL DF1
      REAL DKSAT
      REAL DRIP
      REAL DT
      REAL DWSAT
      REAL EC
      REAL EDIR
      REAL EPSCA
      REAL ETA
      REAL ETA1
      REAL ETP
      REAL ETP1
      REAL ETT
      REAL F
      REAL F1
      REAL FLX1
      REAL FLX2
      REAL FLX3
      REAL KDT
      REAL PC
      REAL PRCP
      REAL PRCP1
      REAL Q2
      REAL RCH
      REAL RHO
      REAL RIB
      REAL RR
      REAL RTDIS (NSOIL)
      REAL RUNOFF,RUNOXX3
      REAL S
      REAL SFCTMP
      REAL SHDFAC
      REAL SIGMA
      REAL SMC   ( NSOIL )
      REAL SH2O  ( NSOIL )
      REAL SMCDRY
      REAL SMCMAX
      REAL SMCREF
      REAL SMCWLT
      REAL STC   ( NSOIL )
      REAL T1
      REAL T24
      REAL TBOT
      REAL TH2
      REAL YY
      REAL YYNUM
      REAL ZSOIL ( NSOIL )
      REAL ZZ1

      REAL Q1, SLOPE, FRZFACT, PSISAT, RUNOFF1, RUNOFF2, RUNOFF3
      REAL EDIR1, EC1, ETT1, QUARTZ, SICE1

      INTEGER NROOT, ICE

CMIC$ TASKCOMMON RITE
      COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF,
     &             DEW,RIB,RUNOXX3
      PARAMETER(CP=1004.5, SIGMA=5.67E-8)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     EXECUTABLE CODE BEGINS HERE.....
C     CONVERT ETP FROM KG M-2 S-1 TO MS-1 AND INITIALIZE DEW.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PRCP1 = PRCP * 0.001
      ETP1 = ETP * 0.001
      DEW = 0.0

C	print*,'NOPAC calling smflx'
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o
	
      IF ( ETP .GT. 0.0 ) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CONVERT PRCP FROM  KG M-2 S-1  TO  M S-1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C ------------    FROZEN GROUND VERSION    --------------------------
C (+) NEW STATES ADDED: SH2O,  AND FROZEN GROUND CORRECTION FACTOR
C
           CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL,
     +             SH2O, SLOPE, KDT, FRZFACT,
     &          SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC,
     &          CMCMAX,SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, 
     &          EDIR1, EC1, ETT1, SFCTMP,Q2,NROOT,RTDIS)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CONVERT MODELED EVAPOTRANSPIRATION FM  M S-1  TO  KG M-2 S-1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        ETA = ETA1 * 1000.0
C	print*,'NOPAC, etp>0, AFTER smflx'
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o

      ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW
C       AND REINITIALIZE ETP1 TO ZERO)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        DEW = -ETP1
        ETP1 = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CONVERT PRCP FROM  KG M-2 S-1  TO  M S-1  AND ADD DEW AMT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        PRCP1 = PRCP1 + DEW
C
        
C ------------    FROZEN GROUND VERSION    --------------------------
C (+) NEW STATES ADDED: SH2O,  AND FROZEN GROUND CORRECTION FACTOR
C
        CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL,
     +             SH2O, SLOPE, KDT, FRZFACT,
     &          SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC,
     &          CMCMAX,SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, 
     &          EDIR1, EC1, ETT1, SFCTMP, Q2, NROOT,RTDIS)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CONVERT MODELED EVAPOTRANSPIRATION FM  M S-1  TO  KG M-2 S-1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        ETA = ETA1 * 1000.0
C	print*,'NOPAC, etp<0, AFTER smflx'
CC	print*,'smc=' , smc
C	print*,'sh2o=', sh2o
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     BASED ON ETP AND E VALUES, DETERMINE BETA
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF ( ETP .LE. 0.0 ) THEN
        BETA = 0.0
        IF ( ETP .LT. 0.0 ) THEN
          BETA = 1.0
          ETA = ETP
        ENDIF
      ELSE
        BETA = ETA / ETP
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC MIXING RATIO AT GRND LEVEL (SKIN)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCC....Q1 = Q2 + ETA * CP / RCH

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
C     CALC ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
C     CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CALL TDFCND ( DF1, SMC(1), B, F1,QUARTZ,SMCMAX,SH2O )
      
C --------------    FROZEN GROUND VERSION     ----------------------
C   THERMAL CONDUCTIVITY ADJUSTMENT
C
      SICE1 = SMC(1) - SH2O(1)
C__________________________________________________________________
C This IF-statement is no longer used if Peters-Lidard thermal
C conductivity is used for all conditions (frozen,unfrozen):
C      IF ( SICE1 .GT. 0.0 ) CALL FRZCND ( DF1, SICE1 )            

      YYNUM = F - SIGMA * T24
      YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-BETA*EPSCA) / RR
      ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0

C ------------    FROZEN GROUND VERSION    --------------------------
C (+) NEW STATES ADDED: SH2O & PARAMETER SMCWLT & PSISAT
C
      CALL SHFLX ( S,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,TBOT,
     +             SMCWLT, PSISAT, SH2O,
     &             B,F1,DF1, ICE, 
     &             QUARTZ )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SET FLX1, AND FLX3 TO ZERO SINCE THEY ARE NOT USED.  FLX2
C     WAS SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FLX1 = 0.0
      FLX3 = 0.0
C
      RETURN
      END
      SUBROUTINE PENMAN(SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,F,T24,S,Q2,
     &                  Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,DQSDT2)

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT.
CC    =======   VARIOUS PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND
CC              PASSED BACK TO THE CALLING ROUTINE FOR LATER USE.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      LOGICAL SNOWNG
      LOGICAL FRZGRA

      REAL A
      REAL BETA
      REAL CH
      REAL CP
      REAL CPH2O
      REAL CPICE
      REAL DELTA
      REAL DEW
      REAL DRIP
      REAL EC
      REAL EDIR
      REAL ELCP
      REAL EPSCA
      REAL ETP
      REAL ETT
      REAL F
      REAL FLX1
      REAL FLX2
      REAL FLX3
      REAL FNET
      REAL LSUBC
      REAL LSUBF
      REAL PRCP
      REAL Q2
      REAL Q2SAT
      REAL R
      REAL RAD
      REAL RCH
      REAL RHO
      REAL RIB
      REAL RR
      REAL RUNOFF,RUNOXX3
      REAL S
      REAL SFCPRS
      REAL SFCTMP
      REAL SIGMA
      REAL T24
      REAL T2V
      REAL TH2
      REAL DQSDT2

CMIC$ TASKCOMMON RITE
      COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF,
     &             DEW,RIB,RUNOXX3

      PARAMETER(CP=1004.6,CPH2O=4.218E+3,CPICE=2.106E+3,R=287.04,
     &   ELCP=2.4888E+3,LSUBF=3.335E+5,LSUBC=2.5E+6,SIGMA=5.67E-8)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     EXECUTABLE CODE BEGINS HERE...
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FLX2 = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DELTA = ELCP * DQSDT2
      T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
      RR = T24 * 6.48E-8 / ( SFCPRS * CH ) + 1.0
      RHO = SFCPRS / ( R * T2V )
      RCH = RHO * CP * CH

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
C     EFFECTS CAUSED BY FALLING PRECIPITATION.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF ( .NOT. SNOWNG ) THEN
        IF ( PRCP .GT. 0.0 ) RR = RR + CPH2O * PRCP / RCH
      ELSE
        RR = RR + CPICE * PRCP / RCH
      ENDIF

      FNET = F - SIGMA * T24 - S

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO
C     ICE ON IMPACT IN THE CALCULATION OF FLX2 AND FNET.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF ( FRZGRA ) THEN
        FLX2 = -LSUBF * PRCP
        FNET = FNET - FLX2
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FINISH PENMAN EQUATION CALCULATIONS.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      RAD = FNET / RCH + TH2 - SFCTMP
      A = ELCP * ( Q2SAT - Q2 )
      EPSCA = ( A * RR + RAD * DELTA ) / ( DELTA + RR )
      ETP = EPSCA * RCH / LSUBC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C    !!!!!!!!!!!! TEST DRAINAGE -> ETP = 0.0 !!!!!!!!!!!!!!!!
C      ETP = 0.0
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      RETURN
      END
       SUBROUTINE PRMSOI(SOLTYP,SBB,SDRYSMC,SF11,SMAXSMC,SREFSMC,
     +            FRZFACT,SSATPSI,SSATDK,SSATDW,SWLTSMC,QTZ)

      IMPLICIT NONE
C
C    SET-UP SOIL PARAMETERS FOR GIVEN SOIL TYPE	
C
C 	INPUT:
C	       SOLTYP: SOIL TYPE (INTEGER INDEX)
C	OUPUT:
C
C	      SOIL PARAMETERS:
C SMAXSMC: MAX SOIL MOISTURE CONTENT (POROSITY)
C SREFSMC: REFERENCE SOIL MOISTURE (ONSET OF SOIL MOISTURE STRESS IN TRANSPIRATION)
C SWLTSMC: WILTING PT SOIL MOISTURE CONTENTS
C SDRYSMC: AIR DRY SOIL MOIST CONTENT LIMITS
C SSATPSI: SATURATED SOIL POTENTIAL
C SSATDK:  SATURATED SOIL HYDRAULIC CONDUCTIVITY 
C SBB:     THE 'B' PARAMETER
C SSATDW:  SATURATED SOIL DIFFUSIVITY
C SF11:    USED TO COMPUTE SOIL DIFFUSIVITY/CONDUCTIVITY
C QUARTZ:  SOIL QUARTZ CONTENT
C *************************************************************************
C
C SOIL TYPES   ZOBLER (1986)	 COSBY ET AL (1984) (quartz cont.(1))
C      1        COARSE  	  LOAMY SAND	    (0.82)
C      2        MEDIUM  	  SILTY CLAY LOAM   (0.10)
C      3        FINE		  LIGHT CLAY	    (0.25)
C      4        COARSE-MEDIUM	  SANDY LOAM	    (0.60)
C      5        COARSE-FINE	  SANDY CLAY	    (0.52)
C      6        MEDIUM-FINE	  CLAY LOAM	    (0.35)
C      7        COARSE-MED-FINE   SANDY CLAY LOAM   (0.60)
C      8        ORGANIC 	  LOAM  	    (0.40)
C      9        LAND ICE	  LOAMY SAND	    (NA using 0.82)
C
C***************************************************************************
C
	INTEGER SOLTYP

        REAL SBB, SDRYSMC, SF11, SMAXSMC, SREFSMC, FRZFACT, SSATPSI
        REAL SSATDK, SSATDW, SWLTSMC, QTZ

        REAL BB(9),DRYSMC(9),F11(9),MAXSMC(9),REFSMC(9),SATPSI(9),
     &       SATDK(9),SATDW(9),WLTSMC(9), QUARTZ(9)

        DATA MAXSMC/0.421, 0.464, 0.468, 0.434, 0.406, 0.465, 0.404,
     &               0.439, 0.421/
        DATA SATPSI/0.04,  0.62, 0.47, 0.14, 0.10, 0.26, 0.14,
     &               0.36,  0.04/
        DATA SATDK /1.41E-5, 0.20E-5, 0.10E-5, 0.52E-5, 0.72E-5,
     &              0.25E-5, 0.45E-5, 0.34E-5, 1.41E-5/
        DATA BB    /4.26,  8.72, 11.55, 4.74, 10.73, 8.17, 6.77,
     &              5.25,  4.26/
        DATA REFSMC/0.283, 0.387, 0.412, 0.312, 0.338, 0.382, 0.315,
     &              0.329, 0.283/
        DATA WLTSMC/0.029, 0.119, 0.139, 0.047, 0.020, 0.103, 0.069,
     &              0.066, 0.029/
c        DATA DRYSMC/0.029, 0.119, 0.139, 0.047, 0.020, 0.103, 0.069,
c     &              0.066, 0.029/
C        DATA DRYSMC/0.07,  0.14,  0.22,  0.08,  0.18,  0.16,  0.12,
C     &               0.10,  0.07/
	DATA DRYSMC/0.07,  0.24,  0.22,  0.08,  0.18,  0.16,  0.12,
     &  	     0.10,  0.07/
        DATA SATDW  /5.71E-6, 2.33E-5, 1.16E-5, 7.95E-6, 1.90E-5,
     &               1.14E-5, 1.06E-5, 1.46E-5, 5.71E-6/
        DATA F11    /-0.999, -1.116, -2.137, -0.572, -3.201, -1.302,
     &               -1.519, -0.329, -0.999/
        DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52, 
     &              0.35, 0.60, 0.40, 0.82/

C NOTE: SATDW = BB*SATDK*(SATPSI/MAXSMC)
C       F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0
C       REFSMC1=MAXSMC*(5.79E-9/SATDK)**(1/(2*BB+3)) 5.79E-9 M/S= 0.5 MM/DAY
C       REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1)
C       WLTSMC1=MAXSMC*(200./SATPSI)**(-1./BB)    (WETZEL AND CHANG, 1987)
C       WLTSMC=WLTSMC1-0.5*WLTSMC1

        SAVE

C    SOIL PARAMETERS

        SBB     = BB(SOLTYP)
	SDRYSMC = DRYSMC(SOLTYP)
	SF11    = F11(SOLTYP)
	SMAXSMC = MAXSMC(SOLTYP)
	SREFSMC = REFSMC(SOLTYP)
        SSATPSI = SATPSI(SOLTYP)
	SSATDK  = SATDK(SOLTYP)
	SSATDW  = SATDW(SOLTYP)
	SWLTSMC = WLTSMC(SOLTYP)
	QTZ     = QUARTZ(SOLTYP)
	
C -------------      FROZEN GROUND VERSION     ---------------------
C    ADJUSTMENT FACTOR OF FROZEN GROUND PARAMETER USING SOIL TYPE 3
C    AS REFERENCE TYPE
C
         FRZFACT = ( SMAXSMC / SREFSMC ) * ( REFSMC(3) / MAXSMC(3) )
C ------------------------------------------------------------------         	
C
	RETURN
	END
      SUBROUTINE PRMVEG(VEGTYP,SHDFAC,RCMIN,RGL,HS,SNUP,
     O     RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL)

      IMPLICIT NONE

C   NEW PARAMETER: SNOW COVER THRESHOLD, SNUP     
C
C   SET-UP VEGETATION PARAMETERS FOR A GIVEN VEGETAION TYPE    
C
C      INPUT:
C             VEGTYP: VEGETATION TYPE (INTEGER INDEX)
C      OUPUT:
C   VEGETATION PARAMETERS:
C     SHDFAC: VEGETATION GREENNESS FRACTION 
C      RCMIN: MIMIMUM STOMATAL RESISTANCE
C   	 RGL: PARAMETER USED IN SOLAR RAD TERM OF CANOPY RESISTANCE FUNCTION
C         HS: PARAMETER USED IN VAPOR PRESSURE DEFICIT TERM OF CANOPY RESISTANCE FUNCTION
C       SNUP: THRESHOLD SNOW DEPTH (IN WATER EQUIVALENT M) THAT IMPLIES 100% SNOW COVER
C
C **********************************************************************
C
C     SSIB VEGETATION TYPES (DORMAN AND SELLERS, 1989; JAM)
C
C    1:   BROADLEAF-EVERGREEN TREES  (TROPICAL FOREST)
C    2:   BROADLEAF-DECIDUOUS TREES
C    3:   BROADLEAF AND NEEDLELEAF TREES (MIXED FOREST)
C    4:   NEEDLELEAF-EVERGREEN TREES
C    5:   NEEDLELEAF-DECIDUOUS TREES (LARCH)
C    6:   BROADLEAF TREES WITH GROUNDCOVER (SAVANNA)
C    7:   GROUNDCOVER ONLY (PERENNIAL)
C    8:   BROADLEAF SHRUBS WITH PERENNIAL GROUNDCOVER
C    9:   BROADLEAF SHRUBS WITH BARE SOIL
C   10:   DWARF TREES AND SHRUBS WITH GROUNDCOVER (TUNDRA)
C   11:   BARE SOIL
C   12:   CULTIVATIONS (THE SAME PARAMETERS AS FOR TYPE 7)
C   13:   GLACIAL
C
C***********************************************************************

	INTEGER VEGTYP
	INTEGER NSOIL
	INTEGER NROOT
	INTEGER I
	
      REAL RTDIS(NSOIL)
      REAL SLDPTH(NSOIL)
      REAL ZSOIL(NSOIL)	

        REAL SHDFAC, RCMIN, RGL, HS, SNUP

	REAL RSMTBL(13),RGLTBL(13),HSTBL(13)
        REAL SNUPX (13)
	
        DATA RSMTBL/150.0,100.0, 125.0, 150.0, 100.0, 70.0,  40.0, 
     *              300.0,400.0, 150.0, 999.0, 40.0,  999.0/ 
        DATA RGLTBL/30.0, 30.0,  30.0,  30.0,  30.0,  65.0,  100.0,
     *              100.0,100.0, 100.0, 999.0, 100.0, 999.0/
        DATA HSTBL/41.69, 54.53, 51.93, 47.35, 47.35, 54.53, 36.35,
     *             42.00, 42.00, 42.00, 999.0, 36.35, 999.0/

C THE FOLLOWING IS THRESHOLD WATER-EQUIV SNOWDEPTH (IN M) 
C IMPLYING 100% SNOW COVERAGE    
        DATA SNUPX/0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.04,
     *             0.04, 0.04, 0.04, 0.025, 0.04, 0.025/    
        
        SNUP = SNUPX ( VEGTYP )

        RCMIN = RSMTBL(VEGTYP)
        RGL = RGLTBL(VEGTYP)
        HS = HSTBL(VEGTYP)
        IF(VEGTYP .EQ. 11) SHDFAC=0.0
C
C CALCULATE ROOT DISTRIBUTION.
C PRESENT VERSION ASSUMES UNIFORM DISTRIBUTION BASED ON SOIL LAYERS.
C	
      DO I=1,NROOT
       RTDIS(I) = -SLDPTH(I)/ZSOIL(NROOT)
      END DO

	RETURN
	END
	
      SUBROUTINE PRTBND(INDI,NOUT,jday,
     & IIDAY,    IITIME,     SFCTMP,  Q2,
     & SFCPRS,  Rg,
     * LW_in, Par_in,  Par_out,  rnet,    
     & S_OBS,  PRCP,    wet,    T1_OBS,
     * T_2cm,   T_4cm,    T_8cm,   T_16cm, 
     & T_32cm,  T_64cm, sm_5cm,
     * sm_20cm, sm_60cm,  w_dir,   u_bar,  eddyuw, 
     * uprim2,  vprim2,   wprim2,  H_OBS,  ETA_OBS,
     * IMONTH,  IDAY,     IREC)

      IMPLICIT NONE

CCC							           CCC
CC							           CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC    NAME:  PRINT BONDVILLE, IL (BND_) DATA:
CC
CC    PURPOSE:  TO WRITE OUT OBSERVATIONAL DATA THAT CAN BE
CC    =======	USED FOR VALIDATION.
CC    
CC    VARIABLES:
CC    =========
C  OBSERVED    SFLX STATE VARIABLES 
C
C   (T1_OBS)                T1: SKIN (GRND SFC) TEMPERATURE (K)
C   (sm_XXcm)    SMC(3 layers): SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
C                               FOR LAYER CENTERED AT XX cm DEPTH.
C   (STC_OBS)    STC(3 layers): SOIL TEMP (K)
C
C  OBSERVED   SFLX OUTPUT 
C  
C (ETA_OBS)  ETA: ACTUAL LATENT HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM SURFACE)
C  (H_OBS)     H: SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM SURFACE)
C  (S_OBS)     S: SOIL HEAT FLUX (W M-2: NEGATIVE, IF DOWNWARD FROM SURFACE)
C             Q1: EFFECTIVE MIXING RATIO AT GRND SFC ( KG KG-1) 
C
CC   ORIGINAL BND_ DATA VARIABLE NAMES
C
C     jday	     Julian Day
C     time	     LST, half hour ending (using IITIME)
C     Ta	     air temperature (C), at 3 m
C     RH	     relative humidity (%) at 3 m
C                    (converted to specific humidity Q2)
C     Pres	     surface pressure in mb
C     Rg	     incoming short wave radiation (W/m2)
C     Par_in	     incoming visible radiation (0.4-0.7 um) in uE/m2/s
C     Par_out	     outgoing or reflected visible light
C     Rnet	     net radiation (W/m2)

C     GHF	     (as S_OBS)  soil or ground heat flux (W/m2)
C       			    
C     RAIN	     total rain for half hour (inches)
C     wet	     wetness sensor (volts), high values indicate wetness.

C     IRT	     (as T1_OBS) surface or skin temp (CONVERTED TO K)

C     T_2cm	     soil temp at 2 cm (C)
C     T_4cm	     soil tmep at 4 cm
C     T_8cm	     soil temp at 8 cm
C     T_16cm	     soil temp at 16 cm
C     T_32cm	     soil temp at 32 cm
C     T_64cm	     soil temp at 64 cm
c
C     STC5           linear interpolations, using the soil
C     STC20          temperatures above, to match the model
C     STC60          depths (5,20 and 60 cm).
C
C     sm_5cm	      soil volumetric water content at 5 cm zone
C     sm_20cm        soil volumetric water content at 20 cm zone
C     sm_60cm	     soil volumetric water content at 60 cm zone
C                    (as SMC_OBS(3 layers))
C
C     w_dir	     wind direction
C     u_bar	     average wind vector speed (m/s), at 6m
C     eddyuw	     (u'w') kinematic shear stress (m2/s2)
C     uprim2	     (u'2) streamwise velocity variance (m2/s2)
C     vprim2	     (v'2) crosswind velocity variance (m2/s2)
C     wprim2	     (w'2) vertical velocity variance  (m2/s2)

C     H 	     (as H_OBS) sensible heat flux (W/m2)
C		     
C     LE	     (as ETA_OBS) latent heat flux (W/m2)


CC The eddy covariance sensors are located at 6 m AGL
CC The bulk density of the soil is 1.4 gm/cm3
CC The site is currently in corn stubble (like it would look
CC  after combining)
CC
CC The units uE/m2/s refer to micro Einsteins per square meter per
CC  second. A uE is 6.02 x 10 (17)  photons.
CC______________________________________________________________________    


      REAL  T0
      PARAMETER (T0=273.16)

      INTEGER INDI, NOUT, NASCII, IREC
      INTEGER IIDAY, IITIME
      INTEGER IMONTH, IDAY, jday
      
      REAL  Q2
      REAL  SFCPRS
      REAL  SDATE
      REAL  DDTIME
      REAL  Rg
      REAL  LW_in
      REAL  Par_in
      REAL  Par_out
      REAL  Rnet
      REAL  S_OBS
      REAL  PRCP
      REAL  wet
      REAL  T1_OBS

      REAL  T_2cm
      REAL  T_4cm
      REAL  T_8cm
      REAL  T_16cm
      REAL  T_32cm
      REAL  T_64cm

      REAL  STC5
      REAL  STC20
      REAL  STC60
      
      REAL  sm_5cm
      REAL  sm_20cm
      REAL  sm_60cm
      REAL  w_dir
      REAL  u_bar
      REAL  eddyuw
      REAL  uprim2
      REAL  vprim2
      REAL  wprim2	
      REAL  H_OBS
      REAL  ETA_OBS
      REAL  RMONTH
      REAL  RIDAY
      REAL  SFCTMP

C ____________________________________________________________
C
C  For now, print soil temperature observations as they are,
C  T_2cm,  T_4cm,  T_8cm,  T_16cm, T_32cm and T_64cm.
C
C  In the future, may use a cubic spline or
C  a simple linear scheme to interpolate
C  to STC_OBS( 5cm ), STC_OBS( 20cm ) and STC_OBS( 60cm ).
C ____________________________________________________________
C 
C  DATA CONVERSIONS TO MATCH MODEL'S VARIABLES

C Soil Temperature in K

       T_2cm = T_2cm + T0
       T_4cm = T_4cm + T0
       T_8cm = T_8cm + T0
      T_16cm = T_16cm + T0
      T_32cm = T_32cm + T0
      T_64cm = T_64cm + T0
      
C SOIL TEMPERATURE AT SOIL MOISTURE MEASUREMENT
C   LEVELS (LINEAR INTERP)
      STC5 = T_4cm + (T_8cm - T_4cm)/4.
      STC20 = T_16cm + 4*(T_32cm - T_16cm)/16.
      STC60 = T_64cm + 4*(T_32cm - T_64cm)/32.


C ____________________________________________________________
C
C  ***** WRITE UNFORMATTED FOR GRADS OUTPUT ********
C

      IF (NOUT .GT. 0) THEN
c convert to real for binary output
      SDATE  = FLOAT(IIDAY)
      DDTIME = FLOAT(IITIME)
      
      RMONTH = FLOAT(IMONTH)
      RIDAY  = FLOAT(IDAY  )
      
      WRITE(NOUT,REC=IREC) SDATE
      IREC = IREC + 1    
      WRITE(NOUT,REC=IREC) DDTIME
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SFCTMP
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) Q2
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SFCPRS
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) Rg
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) LW_in
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) Par_in
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) Par_out
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) rnet
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) S_OBS
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) PRCP
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) wet
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T1_OBS
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T_2cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T_4cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T_8cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T_16cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T_32cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T_64cm
      IREC = IREC + 1
C
      WRITE(NOUT,REC=IREC) STC5
      IREC = IREC + 1		
      WRITE(NOUT,REC=IREC) STC20
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) STC60
      IREC = IREC + 1
C
      WRITE(NOUT,REC=IREC) sm_5cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) sm_20cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) sm_60cm
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) w_dir
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) u_bar
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) eddyuw
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) uprim2
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) vprim2
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) wprim2
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) H_OBS
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETA_OBS
      IREC = IREC + 1

      WRITE(NOUT,REC=IREC) RMONTH
      IREC = IREC + 1	   
      WRITE(NOUT,REC=IREC) RIDAY
      IREC = IREC + 1

      ELSE
            NASCII = -NOUT

      WRITE(NASCII,200) jday, IITIME, SFCTMP, Q2,
     & SFCPRS, Rg, LW_in, Par_in, Par_out,
     & rnet, S_OBS, PRCP, wet, T1_OBS,
     & T_2cm, T_4cm, T_8cm, T_16cm, T_32cm, T_64cm,
     & STC5, STC20, STC60, sm_5cm, sm_20cm, sm_60cm,
     & w_dir, u_bar, eddyuw, uprim2, vprim2, wprim2,
     & H_OBS, ETA_OBS, IMONTH, IDAY
      
      END IF

 200  FORMAT(2(I6,1X),32(F15.4,1X),2(I6,1X))
 
      RETURN
      END
      SUBROUTINE PRTDAILY(INDI,NOUT,NSOIL,jday,
     &                 IIDAY,ITIME,ETADAY,ETADAY_O,
     &                 PCPDAY,SMCDIF,SMCDIF_O,RUNOFFDAY,
     &                 DT,IREC)

      IMPLICIT NONE

CCCC                                                          CCCC
C  From MAIN:                                                    C
C                                                                C
C   CALL PRTDAILY(INDI,NDAILY,NSOIL,IIDAY,ITIME,ETADAY,ETADAY_O, C
C  &                 PCPDAY,SMCDIF,SMCDIF_O,RUNOFFDAY,           C
C  &                 DT,IRECD)                                   C
C                                                                C
CCCC                                                          CCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC    NAME:  PRINT DAILY ACCUMULATED HYDROLOGICAL VARIABLES
CC
CC    VARIABLES:
CC    =========
CC
CC       LABEL     I/O     TYPE    ............DESCRIPTION...............
CC
CC       AET       LOC      RS      ACTUAL EVAPOTRANSPIRATIVE ENERGY
CC                                 (J M-2 S-1)
CC       AET2      LOC      RS      ACTUAL EVAPOTRANSPIRATION (MM)
CC       CH        IOC      RS      DRAG COEF FOR HEAT/MOISTURE
CC       CM        IOC      RS      DRAG COEFFICIENT FOR MOMENTUM
CC       CMC       IOC      RS      CANOPY MOISTURE CONTENT (M)
CC       DEW        IB      RS      DEWFALL AMOUNT  (M S-1)
CC       DRIP       IB      RS      EXCESS CANOPY MOISTURE (M)
CC       EC(EC1)    IB      RS      CANOPY EVAPORATION (M S-1)
CC       EDIR(EDIR1)IB      RS      DIRECT SOIL EVAPORATION (M S-1)
CC       ETA       IOC      RA      FINAL ACTUAL EVAPOTRANSP (KG M-2 S-1)
CC       ETAS      IOC      RA      FINAL ACTUAL EVAPOTRANSP (MM)
CC       ETP       IOC      RA      FINAL POTNTL EVAPOTRANSP (KG M-2 S-1)
CC       ETPS      IOC      RA      FINAL POTNTL EVAPOTRANSP (MM)
CC       ETT(ETT1)  IB      RS      ACCUM PLANT TRANSPIRATION (M S-1)
CC       F         IOC      RS      NET FLUX (TOT DOWNWARD RADIATION)
CC       FLX1       IB      RS      1ST FLUX VALUE (W M-2)
CC       FLX2       IB      RS      2ND FLUX VALUE (W M-2)
CC       FLX3       IB      RS      3RD FLUX VALUE (W M-2)
CC       FOG        IC      LS      INTERMEDIATE HRLY FOG FLAG
CC       FUP        IB      RS      UPWARD GRND LW RADIATION (W M-2)
CC       H         LOC      RS      SENSIBLE HEAT FLUX (W M-2)
CC       HEAT      LOC      RS      SENSIBLE HEAT FLUX SUB-PRODUCT
CC       HEMI       IC      IS      CURRENT HEMISPHERE (1=N, 2=S)
CC       I          IC      IS      1/8 MESH I COORDINATE
CC       ICLAMT     IC      IA      INTERMEDIATE HRLY CLD AMOUNT
CC       ICLTYP     IC      IA      INTERMEDIATE HRLY CLD TYPE
CC       J          IC      IS      1/8 MESH J COORDINATE
CC       K         LOC      IS      LOOP INDEX
CC       NSOIL      IC      IS      SOIL LAYER NUMBER
CC       PET       LOC      RS      POTENTIAL EVAPOTRANSPIRATION (MM)
CC       PRCP      IOC      RA      HALF HOURLY PRECIP AMT (KG M-2 S-1)
CC       Q2        IOC      RS      MIXING RATIO AT 1ST MDL LVL ABV SKIN
CC       Q2SAT     IOC      RS      SAT MXNG RATIO AT 1ST MDL LVL ABV SKIN
CC       RES       LOC      RS      ENERGY BALANCE EQN RESIDUAL (W M-2)
CC       RHO        IB      RS      SFC AIR DENSITY (KG M-3)
CC       RIB        IB      RS      BULK RICHARDSON NUMBER
CC       RLDOWN     IC      RS      DOWNWARD LONGWAVE RADIATION (W M-2)
CC       RR        LOC      RS      SENSIBLE HEAT SUB-PRODUCT
CC       RSOLIN     IC      RS      SOLAR RADIATION (W M-2)
CC       RUNOFF1    IB      RS      GRND SFC RUNOFF (M )
CC       RUNOFF2    IB      RS      UNDERGROUND  RUNOFF (M )
CC       RUNOFF3    IB      RS      RUNOFF WITHIN SOIL LAYERS (M )
CC       SMC       IOC      RA      SOIL MOISTURE CONTENT (VOLUMETRIC)
CC       ZSOIL     IOC      RA      SOIL LAYER DEPTH  ( M )
C      IMPLICIT         DOUBLE PRECISION (A-H,O-Z)

      INTEGER JDAY
      INTEGER IIDAY
      INTEGER ITIME
      INTEGER INDI
      INTEGER IREC
      INTEGER NOUT
      INTEGER NASCII
      INTEGER NSOIL

      REAL  CMCS
      REAL  DT
      REAL  EC1S
      REAL  EDIR1S
      REAL  ETT1S
      REAL  RTIME
      REAL  SDATE

      REAL ETADAY, ETADAY_O, PCPDAY
      REAL SMCDIF, SMCDIF_O, RUNOFFDAY
       

      IF (NOUT .GT. 0) THEN
      SDATE=FLOAT(IIDAY)
      RTIME=FLOAT(ITIME)
      WRITE(NOUT,REC=IREC) SDATE
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RTIME
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) PCPDAY
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETADAY
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETADAY_O
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RUNOFFDAY
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) EDIR1S
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETT1S
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) EC1S
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) CMCS
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SMCDIF
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SMCDIF_O
      IREC = IREC + 1
      
      ELSE
            NASCII = -NOUT

      WRITE(NASCII,200)
     & jday,
     & ITIME,
     & PCPDAY,
     & ETADAY,
     & ETADAY_O,
     & RUNOFFDAY,
     & EDIR1S,
     & ETT1S,
     & EC1S,
     & CMCS,
     & SMCDIF,
     & SMCDIF_O      
      
      
      END IF

 200  FORMAT(I6,1X,I6,10(1x,F15.4))

      RETURN
      END
      SUBROUTINE PRTHMF(INDI,NOUT,NSOIL,jday,
     &                 IIDAY,IITIME,CH,CM,
     &   Z,F,S,CMC,
     &   SMCMAX,SFCTMP,T1,Q1,SFCPRS,
     &   SFCSPD,ETA,ETP,PRCP,
     &   STC,DT,RNET,H,AET,RES,IREC)

      IMPLICIT NONE

CCC							       CCC
C From MAIN:							 C
C      CALL PRTHMF(INDI,NOUT3,NSOIL,IIDAY,itime,CH,CM,Z,F,S,CMC, C
C     & 	 SMCMAX,SFCTMP,T1,Q1,SFCPRS,SFCSPD,ETA,ETP,PRCP, C
C     & 	 STC,DT,RNET,H,AET,IREC3)			 C
CCC							       CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC    NAME:  PRINT THERMODYNAMIC VARIABLES
CC
CC    CPC:   N/A
CC    ===
CC
CC    PURPOSE:  TO CALC AND WRITE OUT PARAMETERS TO AID IN SCIENTIFIC
CC    =======   TUNING AND DEBUGGING OF THE OSU SOIL HYDROLOGY MODEL.
CC              THIS ROUTINE IS ONLY INVOKED FOR POINTS SELECTED BY
CC              SPECIFICATION IN THE CONTROL FILE.
CC
CC    METHOD:   1. CALC NEEDED PARTIAL PRODUCTS AND SUMS.
CC    ======    2. WRITE OUT RESULTS TO PRINT$ FILE.
CC
CC    REFERENCES:  FUNCTIONAL DESCRIPTION, SUBSYSTEM SPECIFICATION
CC
CC    VARIABLES:
CC    =========
CC
CC       LABEL     I/O     TYPE    ............DESCRIPTION...............
CC
CC       AET      LOC      RS      ACTUAL EVAPOTRANSPIRATIVE ENERGY
CC                                 (J M-2 S-1)
CC       AET2     LOC      RS      ACTUAL EVAPOTRANSPIRATION (MM)
CC       CH       IOC      RS      DRAG COEF FOR HEAT/MOISTURE
CC       CM       IOC      RS      DRAG COEFFICIENT FOR MOMENTUM
CC       CMC      IOC      RS      CANOPY MOISTURE CONTENT (M)
CC       DEW       IB      RS      DEWFALL AMOUNT  (M S-1)
CC       DRIP      IB      RS      EXCESS CANOPY MOISTURE (M)
CC       EC        IB      RS      CANOPY EVAPORATION (M S-1)
CC       EDIR      IB      RS      DIRECT SOIL EVAPORATION (M S-1)
CC       ETA      IOC      RA      FINAL ACTUAL EVAPOTRANSP (KG M-2 S-1)
CC       ETP      IOC      RA      FINAL POTNTL EVAPOTRANSP (KG M-2 S-1)
CC       ETPS     IOC      RA      FINAL POTNTL EVAPOTRANSP (J M-2 S-1)
CC       ETT       IB      RS      ACCUM PLANT TRANSPIRATION (M S-1)
CC       F        IOC      RS      NET FLUX (TOT DOWNWARD RADIATION)
CC       FLX1      IB      RS      1ST FLUX VALUE (W M-2)
CC       FLX2      IB      RS      2ND FLUX VALUE (W M-2)
CC       FLX3      IB      RS      3RD FLUX VALUE (W M-2)
CC       FUP       IB      RS      UPWARD GRND LW RADIATION (W M-2)
CC       H        LOC      RS      SENSIBLE HEAT FLUX (W M-2)
CC       HEAT     LOC      RS      SENSIBLE HEAT FLUX SUB-PRODUCT
CC       I         IC      IS      1/8 MESH I COORDINATE
CC       ICLAMT    IC      IA      INTERMEDIATE HRLY CLD AMOUNT
CC       ICLTYP    IC      IA      INTERMEDIATE HRLY CLD TYPE
CC       J         IC      IS      1/8 MESH J COORDINATE
CC       K        LOC      IS      LOOP INDEX
CC       NSOIL     IC      IS      SOIL LAYER NUMBER
CC       PET      LOC      RS      POTENTIAL EVAPOTRANSPIRATION (MM)
CC       PRCP     IOC      RA      HALF HOURLY PRECIP AMT (KG M-2 S-1)
CC       Q2       IOC      RS      MIXING RATIO AT 1ST MDL LVL ABV SKIN
CC       Q2SAT    IOC      RS      SAT MXNG RATIO AT 1ST MDL LVL ABV SKIN
CC       RES      LOC      RS      ENERGY BALANCE EQN RESIDUAL (W M-2)
CC       RHO       IB      RS      SFC AIR DENSITY (KG M-3)
CC       RIB       IB      RS      BULK RICHARDSON NUMBER
CC       RLDOWN    IC      RS      DOWNWARD LONGWAVE RADIATION (W M-2)
CC       RR       LOC      RS      SENSIBLE HEAT SUB-PRODUCT
CC       RSOLIN    IC      RS      SOLAR RADIATION (W M-2)
CC       RUNOFF    IB      RS      GRND SFC RUNOFF (M S-1)
CC       S         IC      RS      GRND SFC FLUX (W M-2)
CC       SFCPRS   IOC      RS      SFC PRESSURE (PASCALS)
CC       SFCSPD    IC      RS      SFC WIND SPEED (M S-1)
CC       SFCTMP   IOC      RS      SFC TEMPERATURE (K)
CC       SMC      IOC      RA      SOIL MOISTURE CONTENT (VOLUMETRIC)
CC       SMCMAX   IOC      RS      MAXIMUM SOIL MOISTURE CONTENT LIMIT
CC       SNODEP   IOC      RA      SNOW DEPTH ( M )
CC       STC      IOC      RA      SOIL TEMPERATURE ( 5 CM AND 95 CM )
CC       T1       IOC      RS      SKIN (GRND SFC) TEMPERATURE (K)
CC       T14      LOC      RS      GRND SFC TEMP TO THE 4TH POWER (K+4)
CC       TIME     IOC      IS      1 HRLY TIME (LOOP INDEX)
CC       Z        IOC      RS      HT ABOVE GRND LVL (M)
CC       ZSOIL    IOC      RA      SOIL LAYER DEPTH  ( M )

      INTEGER IIDAY
      INTEGER jday
      INTEGER IITIME
      INTEGER INDI
      INTEGER IREC
      INTEGER NOUT
      INTEGER NASCII
      INTEGER NSOIL
      
      REAL RNET
      REAL  AET
      REAL  AET2
      REAL  BETA
      REAL  CH
      REAL  CM
      REAL  CMC
      REAL  DEW
      REAL  DRIP
      REAL  DT
      REAL  EC
      REAL  EDIR
      REAL  ETA
      REAL  ETP
      REAL  ETPS
      REAL  ETT
      REAL  F
      REAL  FLX1
      REAL  FLX2
      REAL  FLX3
      REAL  FUP
      REAL  H
      REAL  PRCP
      REAL  Q1
      REAL  RES
      REAL  RHO
      REAL  RIB
      REAL  RUNOFF,RUNOXX3
      REAL  S
      REAL  SDATE
      REAL  RTIME
      REAL  SFCPRS
      REAL  SFCSPD
      REAL  SFCTMP
      REAL  SMCMAX
      REAL  STC      ( 10 )
      REAL  T1
      REAL  T14
      REAL  Z
C
      COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF,
     &             DEW,RIB,RUNOXX3



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CONVERT ACTUAL EVAPOTRANSPIRATION VALUE FROM  KG M-2 S-1  TO
C     W M-2  (BY MULTIPLYING BY 2.5E+6) FOR USE IN ENERGY BALANCE.
C     ALSO CHANGE ACTUAL AND POTENTIAL EVAPORTRANSPIRATION VALUES
C     FROM  KG M-2 S-1  TO  MM  (BY MULTIPLYING BY DT).
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      T14 = T1 * T1 * T1 * T1
      AET = ETA
      ETPS = 2.56E+6 * ETP
      AET2 = DT * ETA
      FUP = 5.67E-8 * T14

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  RESIDUAL OF ALL SURFACE ENERGY BALANCE EQN TERMS:
C     RES = F - H - S - AET - FUP - FLX1 - FLX2 - FLX3
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   NET RADIATION
C
      RNET = F - FUP
C
      IF (NOUT .GT. 0) THEN
      SDATE=FLOAT(IIDAY)
      RTIME=FLOAT(IITIME)
      WRITE(NOUT,REC=IREC) SDATE
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RTIME
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) F
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RNET
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) CH
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) H
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) S
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) AET
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RES
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) T1
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) Q1
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETPS
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) PRCP
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) STC(1)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) STC(2)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) STC(3)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) STC(4)
      IREC = IREC + 1
      
      ELSE
            NASCII = -NOUT

      WRITE(NASCII,200)
     &  jday,
     &  IITIME,
     &  F,
     &  RNET,
     &  CH,
     &  H,
     &  S,
     &  AET,
     &  RES,
     &  T1,
     &  Q1,
     &  ETPS,
     &  PRCP,
     &  STC(1),
     &  STC(2),
     &  STC(3),
     &  STC(4)

      END IF

 200  FORMAT(I6,1X,I6,15(1x,F15.4))
 
      RETURN
      END
      SUBROUTINE PRTHYDF(INDI,NOUT,NSOIL,jday,
     & IIDAY,DDTIME,ETA,
     & ETP,PRCP,SH2O,SMC,ALBEDO,
     & SNEQV, SNMAX, SNOWH, SOILM, SOILW,
     & RUNOFF1,
     & RUNOFF2,
     & RUNOFF3,
     & DT,EDIR1,EC1,ETT1,CMC,IREC)

      IMPLICIT NONE

CCCCC							     CCCCC
C   From MAIN:						     	 C
C							     	 C
C          CALL PRTHYDF(INDI,NOUT1,NSOIL,IIDAY,itime,ETA,	 C
C     & 	      ETP,PRCP,SMC,RUNOFF1, RUNOFF2, RUNOFF3,    C
C     & 	      DT,EDIR ,EC ,ETT ,CMC,IREC1)		 C
C   							         C
CCCCC							     CCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C    NAME:  PRINT HYDROLOGICAL VARIABLES
C
C    VARIABLES:
C    =========
C
C   LABEL     I/O     TYPE    ............DESCRIPTION..............
C
C   AET       LOC      RS      ACTUAL EVAPOTRANSPIRATIVE ENERGY
C			       (J M-2 S-1)
C   AET2      LOC      RS      ACTUAL EVAPOTRANSPIRATION (MM)
C   CH        IOC      RS      DRAG COEF FOR HEAT/MOISTURE
C   CM        IOC      RS      DRAG COEFFICIENT FOR MOMENTUM
C   CMC       IOC      RS      CANOPY MOISTURE CONTENT (M)
C   DEW        IB      RS      DEWFALL AMOUNT  (M S-1)
C   DRIP       IB      RS      EXCESS CANOPY MOISTURE (M)
C   EC(EC1)    IB      RS      CANOPY EVAPORATION (M S-1)
C   EDIR(EDIR1)IB      RS      DIRECT SOIL EVAPORATION (M S-1)
C   ETA       IOC      RA      FINAL ACTUAL EVAPOTRANSP (KG M-2 S-1)
C   ETAS      IOC      RA      FINAL ACTUAL EVAPOTRANSP (MM)
C   ETP       IOC      RA      FINAL POTNTL EVAPOTRANSP (KG M-2 S-1)
C   ETPS      IOC      RA      FINAL POTNTL EVAPOTRANSP (MM)
C   ETT(ETT1)  IB      RS      ACCUM PLANT TRANSPIRATION (M S-1)
C   F	      IOC      RS      NET FLUX (TOT DOWNWARD RADIATION)
C   FLX1       IB      RS      1ST FLUX VALUE (W M-2)
C   FLX2       IB      RS      2ND FLUX VALUE (W M-2)
C   FLX3       IB      RS      3RD FLUX VALUE (W M-2)
C   FOG        IC      LS      INTERMEDIATE HRLY FOG FLAG
C   FUP        IB      RS      UPWARD GRND LW RADIATION (W M-2)
C   H	      LOC      RS      SENSIBLE HEAT FLUX (W M-2)
C   HEAT      LOC      RS      SENSIBLE HEAT FLUX SUB-PRODUCT
C   HEMI       IC      IS      CURRENT HEMISPHERE (1=N, 2=S)
C   I	       IC      IS      1/8 MESH I COORDINATE
C   ICLAMT     IC      IA      INTERMEDIATE HRLY CLD AMOUNT
C   ICLTYP     IC      IA      INTERMEDIATE HRLY CLD TYPE
C   J	       IC      IS      1/8 MESH J COORDINATE
C   K	      LOC      IS      LOOP INDEX
C   NSOIL      IC      IS      SOIL LAYER NUMBER
C   PET       LOC      RS      POTENTIAL EVAPOTRANSPIRATION (MM)
C   PRCP      IOC      RA      HALF HOURLY PRECIP AMT (KG M-2 S-1)
C   Q2        IOC      RS      MIXING RATIO AT 1ST MDL LVL ABV SKIN
C   Q2SAT     IOC      RS      SAT MXNG RATIO AT 1ST MDL LVL ABV SKIN
C   RES       LOC      RS      ENERGY BALANCE EQN RESIDUAL (W M-2)
C   RHO        IB      RS      SFC AIR DENSITY (KG M-3)
C   RIB        IB      RS      BULK RICHARDSON NUMBER
C   RLDOWN     IC      RS      DOWNWARD LONGWAVE RADIATION (W M-2)
C   RR        LOC      RS      SENSIBLE HEAT SUB-PRODUCT
C   RSOLIN     IC      RS      SOLAR RADIATION (W M-2)
C   RUNOFF1    IB      RS      GRND SFC RUNOFF (M )
C   RUNOFF2    IB      RS      UNDERGROUND  RUNOFF (M )
C   RUNOFF3    IB      RS      RUNOFF WITHIN SOIL LAYERS (M )
C   SMC       IOC      RA      SOIL MOISTURE CONTENT (VOLUMETRIC)
C SH2O(NSOIL): UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
C  SNEQV: WATER EQUIVALENT SNOW DEPTH (M) (formerly called SNODEP)
C  SNMAX: SNOW MELT (M) (WATER EQUIVALENT)
C  SNOWH: SNOW DEPTH (M)
C  SOILM: TOTAL SOIL COLUMN WATER CONTENT (M)
C  SOILW: AVAILABLE SOIL MOISTURE (UNITLESS FRACTION) 
C   ZSOIL     IOC      RA      SOIL LAYER DEPTH  ( M )

      INTEGER IIDAY
      INTEGER jday
      INTEGER DDTIME
      INTEGER INDI
      INTEGER IREC
      INTEGER NASCII
      INTEGER NOUT
      INTEGER NSOIL

      REAL  ALBEDO
      REAL  CMC
      REAL  CMCS
      REAL  DT
      REAL  EC1
      REAL  EC1S
      REAL  EDIR1
      REAL  EDIR1S
      REAL  ETA
      REAL  ETP
      REAL  ETT1
      REAL  ETT1S
      REAL  PRCP
      REAL  RTIME
      REAL  RUNOFF1
      REAL  RUNOFF11
      REAL  RUNOFF2
      REAL  RUNOFF22
      REAL  RUNOFF3
      REAL  RUNOFF33
      REAL  SDATE
      REAL  SMC(10)
      REAL  SH2O(10)
      REAL  SNEQV
      REAL  SNMAX
      REAL  SNOWH
      REAL  SOILM
      REAL  SOILW
      REAL  SNEQVMM
      REAL  SNMELT
      REAL  SOILMM
      REAL  TRUNOFF

CC  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
      EC1S=EC1*1000.*DT
      EDIR1S=EDIR1*1000.*DT
      ETT1S = ETT1*1000.*DT
      CMCS=CMC*1000.
      RUNOFF11 = RUNOFF1*DT*1000.
      RUNOFF22 = RUNOFF2*DT*1000.
      RUNOFF33 = RUNOFF3*1000.
      TRUNOFF = RUNOFF11 + RUNOFF22 + RUNOFF33
      
      SNEQVMM = SNEQV*1000.
      SNMELT  = SNMAX*1000.
      SOILMM  = SOILM*1000.
      
      IF (NOUT .GT. 0) THEN
      SDATE=FLOAT(IIDAY)
      RTIME=FLOAT(DDTIME)
      WRITE(NOUT,REC=IREC) SDATE
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RTIME
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) PRCP
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETP
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETA
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RUNOFF11
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RUNOFF22
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) RUNOFF33
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) TRUNOFF
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) EDIR1S
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ETT1S
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) EC1S
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) CMCS
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SH2O(1)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SH2O(2)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SH2O(3)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SH2O(4)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SMC(1)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SMC(2)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SMC(3)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SMC(4)
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) ALBEDO
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SNOWH
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SOILW
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SNEQVMM
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SNMELT
      IREC = IREC + 1
      WRITE(NOUT,REC=IREC) SOILMM
      IREC = IREC + 1
      
      ELSE
            NASCII = -NOUT

      WRITE(NASCII,200) jday, DDTIME, PRCP, ETP,
     & ETA,
     & RUNOFF11, RUNOFF22, RUNOFF33, TRUNOFF,
     & EDIR1S, ETT1S, EC1S, CMCS,
     & SH2O(1), SH2O(2), SH2O(3), SH2O(4),
     & SMC(1), SMC(2), SMC(3), SMC(4),
     & ALBEDO,
     & SNOWH, 
     & SOILW,
     & SNEQVMM,
     & SNMELT,
     & SOILMM

     
      END IF

 200  FORMAT(I6,1X,I6,15(1x,F15.4))

      RETURN
      END
      SUBROUTINE QDATAP (T,P,RH,QD,QS,ES) 

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  PURPOSE:  OBTAIN SPECIFIC HUMIDITY (q) FROM RELATIVE HUMIDITY 
CC            AND GIVEN PRESSURE AND TEMPERATURE.
CC            
CC
CC            FORMULAS AND CONSTANTS FROM ROGERS AND YAU, 1989: 'A 
CC            SHORT COURSE IN CLOUD PHYSICS', PERGAMON PRESS, 3rd ED.
CC
CC                                   Pablo J. Grunmann, 3/6/98.
CC                Updated to eliminate subroutine SVP, 6/24/98.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C----------------------------------------
C In:
C	T    Temperature (K)
C	P    Pressure (Pa)
C	RH   Relative humidity (%)
C-----------------------------------------
C Out:
C	QD   Specific humidity (Kg/Kg)
C	QS   Saturation Specific humidity (Kg/Kg)
C	ES   Saturation vapor pressure for water (Pa)
C----------------------------------------
      REAL  T
      REAL  P
      REAL  RH
      REAL  RHF
      REAL  QD
      REAL  QS
      REAL  ES
      REAL  EP
      REAL  EPS
      REAL E

      PARAMETER (eps=0.622 )
C 
C     ABOUT THE PARAMETER:
C      
C     eps ---------- (Water)/(dry air) molecular mass ratio, epsilon
C _____________________________________________________________________
C
C    function E(T) = Sat. vapor pressure (in Pascal) at 
C                    temperature T (uses Clausius-Clapeyron).
          Es = E(T)
C  CONVERT REL. HUMIDITY (%) TO THE FRACTIONAL VALUE
          RHF = RH/100.

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC      CALCULATE SATURATION MIXING RATIO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C      QS = 0.622 * ES /P   was substituted by a more precise
C formula:                              -PABLO J. GRUNMANN, 05/28/98.
        QS = 0.622 * ES /(P - (1.-0.622)*ES)

C
C  CONVERSION FROM REL. HUMIDITY:
C     (Rogers, pg. 17)
C
	EP = (P*Es*RHF)/(P - Es*(1. - RHF))
	QD = eps*EP/(P - (1. - eps)*EP)
C     
          RETURN
          END
        SUBROUTINE READBND(jday,  itime, SFCTMP,   RH, SFCPRS,  Rg,
     * Par_in, Par_out,  rnet, LW_in, GHF, PRCP, wet, SKN_IRT,
     * T_2cm,  T_4cm,  T_8cm,  T_16cm, T_32cm, T_64cm, sm_5cm,
     * sm_20cm, sm_60cm,  w_dir, u_bar,     eddyuw, 
     * uprim2,   vprim2,      wprim2,         H,     LE,
     * IMONTH,IDAY,NREAD)
     
      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC    NAME:  READBND
CC
CC
CC    PURPOSE:  
CC           THIS IS A MODIFICATION OF V.KOREN'S  
CC           SUBROUTINE READPILPS(ICALB,NREAD,NREAD2,IYEAR,IMONTH,IDAY,
CC                                IHOUR,P,T,Q,U,V,LWDN,RAIN,SOLDN,SOUT,
CC                                TSKNEST)
CC           TO READ FORCING DATA FROM BONDVILLE, IL ATDD SITE.
CC
CC                                      PABLO J. GRUNMANN, 05/98
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                  ARGUMENT LIST IN THE CALL READBND
C                 PREPARED BY PABLO GRUNMANN IN 05/98
C
C       SFCPRS: PRESSURE AT 1ST MDL LVL ABV SKIN (PASCALS)
C         PRCP: PRECIP RATE (KG M-2 S-1)
C       SFCTMP: AIR TEMPERATURE  AT 1ST MDL LVL ABV SKIN (K)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
      INTEGER jday,itime,IMONTH,IDAY,NREAD,I
      INTEGER  JULM(13)
      
      REAL  SFCTMP
      REAL  T
      REAL  RH
      REAL  SFCPRS
      REAL  P
      REAL  Rg
      REAL  Par_in
      REAL  Par_out
      REAL  rnet
      REAL  LW_in
      REAL  GHF
      REAL  PRCP
      REAL  RAIN
      REAL  wet
      REAL  SKN_IRT
      REAL  T_2cm
      REAL  T_4cm
      REAL  T_8cm
      REAL  T_16cm
      REAL  T_32cm
      REAL  T_64cm
      REAL  sm_5cm
      REAL  sm_20cm
      REAL  sm_60cm
      REAL  w_dir
      REAL  u_bar
      REAL  eddyuw
      REAL  uprim2
      REAL  vprim2
      REAL  wprim2
      REAL  w_speed
      REAL  CO2
      REAL  H
      REAL  LE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c Save previous step required forcing data
c to make up in case of missing.
      REAL  XOLD(32)
       do i=1,32
        XOLD(i)=0.0
       end do
       XOLD(2)=   SFCTMP - 273.16
       XOLD(3)=   RH
       XOLD(5)=   SFCPRS/1.E2
       XOLD(6)=   Rg
       XOLD(7)=   LW_in
       XOLD(12)=  1800.*PRCP/25.4
       XOLD(14)=  SKN_IRT
       XOLD(25)=  u_bar
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C      READ BND_ DATA
C  ATTENTION!
C  REMEMBER TO REMOVE TITLE FROM BND_ FILE (1ST LINE), THEN READ DATA

        READ(NREAD,*) jday,  itime, w_speed,  w_dir, 
     * T,   RH,  P,  Rg, Par_in,   Par_out, 
     * rnet,   GHF, RAIN,   wet,  SKN_IRT,
     * T_2cm,  T_4cm,  T_8cm,  T_16cm, T_32cm, T_64cm, 
     * u_bar,  eddyuw, uprim2,   vprim2,   wprim2,     
     * H,  LE, CO2, LW_in, 
     * sm_5cm, sm_20cm, sm_60cm
    
C-----------------------------------------------------------------------
C TILDEN MEYERS' BND_ DATA:
C
C    jday       Julian Day
C    itime       LST, half hour ending
C    Ta         air temperature (C), at 3 m
C    RH         relative humidity (%) at 3 m
C    Pres       surface pressure in mb
C    Rg         incoming short wave radiation (W/m2)
C    Par_in     incoming visible radiation (0.4-0.7 um) in uE/m2/s
C    Par_out    outgoing or reflected visible light
C    Rnet       net radiation (W/m2)
C    GHF        soil or ground heat flux (W/m2)
C    RAIN       total rain for half hour (inches)
C    wet       wetness sensor (voltage - higher values indicatE wetness)
C    SKN_IRT    surface or skin temp (C)
C    T_2_cm     soil temp at 2 cm (C)
C    T_4_cm     soil tmep at 4 cm
C    T_8_cm     soil temp at 8 cm
C    T_16_cm    soil temp at 16 cm
C    T_32_cm    soil temp at 32 cm
C    T_64_cm    soil temp at 64 cm
C    sm_5_cm    soil volumetric water content at 5 cm zone
C    sm_20_cm   soil volumetric water content at 20 cm zone
C    sm_60_cm   soil volumetric water content at 60 cm zone
C    w_dir         wind direction
C    u_bar         average wind vector speed (m/s), at 6m
C    u'w'          kinematic shear stress (m2/s2)
C    u'2           streamwise velocity variance (m2/s2)
C    v'2           crosswind velocity variance (m2/s2)
C    w'2           vertical velocity variance  (m2/s2)
C    H             sensible heat flux (W/m2)
C    LE            latent energy flux (W/m2)
C    
C The eddy covariance sensors are located at 6 m AGL
C The bulk density of the soil is 1.4 gm/cm3
C The site is currently in corn stubble (like it would look 
C  after combining)
C
C The units uE/m2/s refer to micro Einsteins per square meter per
C  second. A uE is 6.02 x 10 (17)  photons.
C
C
C-----------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  THIS PART STILL REQUIRES ATTENTION:                       C
C                                                            C
C   IN LACK OF FORCING DATA, CURRENT PROCEDURE HAS BEEN      C
C (AND IT IS) TO USE PREVIOUS TIME-STEP RECORDED VALUES.     C
C                                                            C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       IF (T       .LT. -73.) THEN
         T       = XOLD(2)
       ENDIF
       IF (RH      .LT.    0. ) THEN
         RH      = XOLD(3)
       ENDIF
       IF (P       .LT.  850. ) THEN
         P       = XOLD(5)
       ENDIF
       IF (Rg      .LT. -100. ) THEN
         Rg      = XOLD(6)
       ENDIF
       IF (LW_in      .LT. 0.) THEN
         LW_in      = XOLD(7)
       ENDIF
       IF (RAIN    .LT.    0. ) THEN
         RAIN    = XOLD(12)
       ENDIF
       IF (u_bar   .LT.    0. ) THEN
C MISSING u_bar, USE W_SPEED WITH LINEAR REGRESSION CONVERSION
         u_bar   = 0.85*w_speed
C BUT, IF EVEN W_SPEED IS NOT AVAILABLE, USE PREVIOUS STEP VALUE
	 IF (u_bar   .LT.    0. ) THEN
	  u_bar   = XOLD(25)
	 ENDIF
       ENDIF


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      CONVERT VARIABLES
C
C OBTAIN DAY AND MONTH: JULIAN DATE SUBROUTINE
      CALL JULDATE(JDAY,IMONTH,IDAY,JULM)
C
C  CONVERT RAIN FROM INCHES IN 30MIN TO KG M^-2 S^-1 
C (SAME AS PRECIP.RATE IN MM/SEC)
C
        PRCP = RAIN*25.4/1800.
C
C  AIR TEMPERATURE IN KELVIN
C
        SFCTMP = T + 273.16
C
C  SFC PRESSURE IN PASCAL
C
        SFCPRS =  P*1.E2
C
C ---- DONE CONVERTING VARIABLES ----------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        RETURN
        END 
      SUBROUTINE READCNTL(CNTRFL,NRUN,DT,NSOIL,NROOT,Z,SLDPTH,
     & SOILTP,VEGTYP,SLOPTYP,ALBEDOM,SHDFACM,MLAI,SNOALB,ICE,
     & TBOT,Z0,CZIL,REFKDT,T1,
     & STC,SMC,SH2O,CMC,SNODEP,SNEQV,FILENAME,
     & LATITUDE,
     & LONGITUDE,
     & JDAY,
     & TIME)
     
      IMPLICIT NONE

C   READS A CONTROL FILE WHICH INCLUDES INPUT PARAMETERS FOR OSU MODEL  -

      INTEGER NRUN
      INTEGER NSOIL
      INTEGER NROOT
      INTEGER SOILTP
      INTEGER VEGTYP
      INTEGER SLOPTYP
      INTEGER ICE

      REAL  SLDPTH(4)
      REAL  STC(4)
      REAL  SMC(4)
      REAL  SH2O(4)
      REAL  DT
      REAL  Z
      REAL  ALBEDOM(13)
      REAL  SHDFACM(13)
      REAL  MLAI(13)
      REAL  SNOALB
      REAL  TBOT
      REAL  Z0
      REAL  T1
      REAL  CMC
      REAL  SNODEP
      REAL  SNEQV
      REAL  CZIL
      REAL  REFKDT
      
      REAL  LATITUDE
      REAL  LONGITUDE
      
      INTEGER  JDAY
      INTEGER  TIME
            
      INTEGER NREAD, I, K

      CHARACTER*100 CNTRFL, LINEFEED
      CHARACTER*9 FILENAME

        NREAD = 21
        OPEN (UNIT=NREAD, FILE=CNTRFL, STATUS='OLD',
     &        FORM='FORMATTED')

C JUMP-READ 2 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED)

      DO I=1,2
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO

C Model Configuration:
        READ (NREAD,  * ) LATITUDE
        READ (NREAD,  * ) LONGITUDE
        READ (NREAD, 100) JDAY
        READ (NREAD, 100) TIME
C
        READ (NREAD, 150) NRUN
        READ (NREAD,  * ) DT
        READ (NREAD, 100) NSOIL
        READ (NREAD, 100) NROOT
        READ (NREAD,  *)  Z
        READ (NREAD,  *) (SLDPTH(K), K=1,NSOIL)

C NEW: NROOT, SLDPTH(K), K=1,NSOIL
C REMOVED:    ZSOIL(K), K=1,NSOIL 

C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED)

      DO I=1,3
        READ(NREAD,'(A)') LINEFEED
        WRITE(*,'(A)') LINEFEED
      END DO

C ATMOSPHERIC DATA FILE TO BE USED FOR FORCING:

        READ(NREAD,'(A)') FILENAME
        PRINT*,'  FILENAME = ', FILENAME
C
C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED)
      DO I=1,3
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO
C  ------------------------
C Land surface characteristics:
C  ------------------------
        READ (NREAD, 100) SOILTP
        READ (NREAD, 100) VEGTYP
        READ (NREAD, 100) SLOPTYP
C ...........SKIP COMMENTS (LINE FEED)
      DO I=1,3
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO
C ...........
        READ (NREAD, *) (ALBEDOM(K), K=1,12)
        ALBEDOM(13) = ALBEDOM(1)
C ...........SKIP COMMENTS (LINE FEED)
      DO I=1,3
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO
C ...........
        READ (NREAD, *) (SHDFACM(K), K=1,12)
        SHDFACM(13) = SHDFACM(1)

C ...........SKIP COMMENTS (LINE FEED)
      DO I=1,3
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO
C ...........
        READ (NREAD, *) (MLAI(K), K=1,12)
        MLAI(13) = MLAI(1)
C ...........SKIP COMMENTS (LINE FEED)
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
C ...........
        READ (NREAD, *) SNOALB
        READ (NREAD, 100) ICE

C NEW: ALBEDO, SHDFAC, ICE, 
C REMOVED: LAND, CFACTR
                
C
C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED)
      DO I=1,3
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO

C ------------------ 
C PHYSICAL PARAMETERS:
C ------------------ 
        READ (NREAD, *) TBOT
        READ (NREAD, *) Z0
        READ (NREAD, *) CZIL
        READ (NREAD, *) REFKDT
C
C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED)
      DO I=1,3
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO

C ------------------------ 
C INITIAL STATE VARIABLES:
C ------------------------ 
        READ (NREAD, *) T1
        READ (NREAD, *) (STC(K), K=1,NSOIL)
        READ (NREAD, *) (SMC(K), K=1,NSOIL)
        READ (NREAD, *) (SH2O(K), K=1,NSOIL)
        READ (NREAD, *) CMC
        READ (NREAD, *) SNODEP
        READ (NREAD, *) SNEQV
C READ LAST 2 LINES (CHARACTER) TO SHOW THAT CONTROLFILE WAS
C COMPLETELY READ
      DO I=1,2
        READ(NREAD,'(A)') LINEFEED
C writ        WRITE(*,'(A)') LINEFEED
      END DO
C NEW: SNEQV
C REMOVED: 
C          ESD
C          KDT
C          DKSAT
C          B
C          DWSAT
C          SMCMAX
C          SMCREF
C          PSISAT
C          SMCWLT
C          DMAX
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C   V. KOREN   5/21/97
C      TWO PARAMETERS TO RUN FROZEN GROUND
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                            C
C        READ (NREAD, *)   CVFRZ             C
C        READ (NREAD, *)   FRZK              C
C                                            C
C         READ (NREAD, *) SLOPE              C
C                                            C
C        CVFRZ = 3.0                         C
C        FRZCL = 0.50                        C
C        SLOPE = 0.0005                      C
C                                            C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

 100    FORMAT(I6)
 150    FORMAT(I7)
C 200    FORMAT(F14.8)
C 400    FORMAT(10(2X, F14.8))
        CLOSE(NREAD)
C
        RETURN
        END     
      SUBROUTINE REDPRM(VEGTYP,SOITYP, SLOPTYP, 
     +     CFACTR, CMCMAX, RSMAX, TOPT,REFKDT,KDT,
     O     Z0,SHDFAC,RCMIN,RGL,HS,
     +     ZBOT, FRZX, PSISAT, SLOPE, SNUP,SALP,
     O     B,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF,
     +     SMCDRY,F1,QUARTZ,
     O     RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL)
      
      IMPLICIT NONE
      
C     ---------------     FROZEN GROUND VERSION    ------------------
C     FROZEN GROUND CORRECTION FACTOR IS ADDED
C     
C     ---------------------------------------------------------------
      
C     THIS SUBROUTINE READS THE SOIL AND VEGETATION PARAMETERS
C     REQUIRED FOR THE EXECUSION OF THE OSU LAND-SURFACE SCHEME. 
C     BY F. CHEN 3/15/95
C     
      INTEGER NROOT,VEGTYP,SOITYP, SLOPTYP
      INTEGER NSOIL

      REAL Z0, SHDFAC, RCMIN, RGL, HS, FRZFACT, PSISAT
      REAL SLOPE, SNUP, SALP
      REAL B, DKSAT, DWSAT, SMCMAX, SMCWLT, SMCREF
      REAL SMCDRY, F1, QUARTZ
      REAL REFDK
      REAL REFKDT
      REAL KDT
      REAL FRZX
      REAL FRZK
      REAL RTDIS(NSOIL)
      REAL SLDPTH(NSOIL)
      REAL ZSOIL(NSOIL)

C     ***   NEW CLASS PARAMETER 'SLOPTYP' WAS INCLUDED TO ESTIMATE
C     ***   LINEAR RESERVOIR COEFFICIENT 'SLOPE' TO THE BASEFLOW RUNOFF 
C     ***   OUT OF THE BOTTOM LAYER. LOWEST CLASS (SLOPTYP=0)MEANS 
C     ***   HIGHEST SLOPE PARAMETER= 1
C     ***   DEFINITION OF SLOPTYP FROM 'ZOBLER' SLOPE TYPE
C     ***  	SLOPE CLASS	PERCENT SLOPE
C     1		0-8
C     2		8-30
C     3		> 30
C     4		0-30
C     5		0-8 & > 30
C     6		8-30 & > 30
C     7		0-8, 8-30, > 30
C     9		GLACIAL ICE
C     BLANK           OCEAN/SEE
C     NOTE:  CLASS 9 FROM 'ZOBLER' FILE SHOULD BE REPLACED BY 8  
C     AND 'BLANK'  9
C     
      
      REAL CFACTR 
      REAL CMCMAX 
      REAL RSMAX
      REAL TOPT
      REAL ZBOT
      
      REAL SLOPES (9)
      DATA SLOPES/0.1, 0.6, 1.0, 0.35, 0.55, 0.8, 0.63, 0., 0./
      
C     Specify depth[m] of lower boundary soil temperature
      
      ZBOT = 3.0
     
C
C SPECIFY SNOW DISTRIBUTION SHAPE PARAMETER  
C SALP   - SHAPE PARAMETER OF DISTRIBUTION FUNCTION
C 	   OF SNOW COVER. FROM ANDERSON'S DATA (HYDRO-17)
C 	   BEST FIT IS WHEN SALP = 2.6

         SALP = 2.6

C     SPECIFY THE VEG/SOIL PARAMETERS
C     Set two canopy water parameters
      
      CFACTR = 0.5
      CMCMAX = 0.5E-3
      
C     Set max. stomatal resistance
      
      RSMAX  = 5000.0
      
C     Set optimum transpiration air temperature
      
      TOPT   = 298.0
      
C     Set drainage parameter based on slope index     
      
      SLOPE = SLOPES ( SLOPTYP)
      
C     SET-UP VEGETATION PARAMETERS
 
      CALL  PRMVEG(VEGTYP,SHDFAC,RCMIN,RGL,HS,SNUP,
     O     RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL)

C     SET-UP SOIL PARAMETERS
      
      CALL PRMSOI(SOITYP,B,SMCDRY,F1,SMCMAX,SMCREF, FRZFACT,
     &     PSISAT,DKSAT,DWSAT,SMCWLT,QUARTZ)
      
C  *****  KDT IS DEFINED BY REFERENCE REFKDT AND DKSAT *****
C    REFDK=2.E-6 IS THE SAT. DK. VALUE FOR THE SOIL TYPE 2 
C
          REFDK = 2.0E-6
          KDT   = REFKDT * DKSAT/REFDK
          
C -----------------------------------------------------------------------   
C    FROZEN GROUND PARAMETER, FRZK, DEFINITION
C    FRZK IS ICE CONTENT THRESHOLD ABOVE WHICH FROZEN SOIL IS IMPERMEABLE
C    REFERENCE VALUE OF THIS PARAMETER FOR THE LIGHT CLAY SOIL (TYPE=3)
C        FRZK = 0.15 M
C
         FRZK = 0.15
         
C    TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
C    FRZFACT IS CALCULATED IN REDPAR SUBROUTINE
C
         FRZX = FRZK * FRZFACT       
C ----------------------------------------------------------------------- 

      RETURN
      END
      SUBROUTINE ROSR12 ( P, A, B, C, D, DELTA, NSOIL )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN
CC    =======   BELOW:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER K
      INTEGER KK
      INTEGER NSOIL
      INTEGER NSOLD
      PARAMETER ( NSOLD = 4 )
      
      REAL P(NSOLD)
      REAL A(NSOLD)
      REAL B(NSOLD)
      REAL C(NSOLD)
      REAL D(NSOLD)
      REAL DELTA(NSOLD)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      C(NSOIL) = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SOLVE THE COEFS FOR THE 1ST SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      P(1) = -C(1) / B(1)
      DELTA(1) = D(1) / B(1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 10 K = 2 , NSOIL
        P(K) = -C(K) * ( 1.0 / (B(K) + A (K) * P(K-1)) )
        DELTA(K) = (D(K)-A(K)*DELTA(K-1))*(1.0/(B(K)+A(K)*P(K-1)))
 10   CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SET P TO DELTA FOR LOWEST SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      P(NSOIL) = DELTA(NSOIL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 20 K = 2 , NSOIL
        KK = NSOIL - K + 1
        P(KK) = P(KK) * P(KK+1) + DELTA(KK)
 20   CONTINUE

      RETURN
      END
      SUBROUTINE SFCDIF(ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS)

      IMPLICIT NONE

C     ******************************************************************
C     *                                                                *
C     *                        SURFACE LAYER                           *
C     *                                                                *
C     ******************************************************************
C-----------------------------------------------------------------------

      REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
      REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, SQVISC
      REAL RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU, PSLHS
      REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
      REAL SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH
      REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
      REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
      REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, RLMA
C      REAL ZTFC

      INTEGER ITRMX, ILECH, ITR

                             PARAMETER
     &(WWST=1.2,WWST2=WWST*WWST,G=9.8,VKRM=0.40,EXCM=0.001
     &,BETA=1./270.,BTG=BETA*G,ELFC=VKRM*BTG
     &,WOLD=.15,WNEW=1.-WOLD,ITRMX=05,PIHF=3.14159265/2.
C-----------------------------------------------------------------------
     &,EPSU2=1.E-4,EPSUST=0.07,EPSIT=1.E-4,EPSA=1.E-8
     &,ZTMIN=-5.,ZTMAX=1.,HPBL=1000.0
     &,SQVISC=258.2)

      PARAMETER
     &(RIC=0.183,RRIC=1.0/RIC,FHNEU=0.8,RFC=0.191
     &,RFAC=RIC/(FHNEU*RFC*RFC))

C-----------------------------------------------------------------------
C **** LECH'S SURFACE FUNCTIONS ****
      PSLMU(ZZ)=-0.96*log(1.0-4.5*ZZ)
      PSLMS(ZZ)=ZZ*RRIC-2.076*(1.-1./(ZZ+1.))
      PSLHU(ZZ)=-0.96*log(1.0-4.5*ZZ)
      PSLHS(ZZ)=ZZ*RFAC-2.076*(1.-1./(ZZ+1.))

C-----------------------------------------------------------------------
C **** PAULSON'S SURFACE FUNCTIONS *****
      PSPMU(XX)=-2.*log((XX+1.)*0.5)-log((XX*XX+1.)*0.5)+2.*ATAN(XX)
     &          -PIHF
      PSPMS(YY)=5.*YY
      PSPHU(XX)=-2.*log((XX*XX+1.)*0.5)
      PSPHS(YY)=5.*YY
C
      ILECH=0
C***********************************************************************
C-----------------------------------------------------------------------
C ZTFC: RATIO OF ZOH/ZOM  LESS OR EQUAL THAN 1
CC......ZTFC=0.1
C CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
C Commented out to allow argument-passed value
C      CZIL=0.2
      
      ZILFC=-CZIL*VKRM*SQVISC
C-----------------------------------------------------------------------
      ZU=Z0
CC.......ZT=Z0*ZTFC
C-----------------------------------------------------------------------
      RDZ=1./ZLM
      CXCH=EXCM*RDZ
C-----------------------------------------------------------------------
      DTHV=THLM-THZ0
C
      DU2=AMAX1(SFCSPD*SFCSPD,EPSU2)
C--------------BELJARS CORRECTION OF USTAR------------------------------
      BTGH=BTG*HPBL
      WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.)
      USTAR=AMAX1(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
C--------------ZILITINKEVITCH APPROACH FOR ZT--------------------------------
      ZT=EXP(ZILFC*SQRT(USTAR*Z0))*Z0
C-----------------------------------------------------------------------
        ZSLU=ZLM+ZU
        ZSLT=ZLM+ZT
C
        RLOGU=log(ZSLU/ZU)
        RLOGT=log(ZSLT/ZT)
C
        RLMO=ELFC*AKHS*DTHV/USTAR**3
C-----------------------------------------------------------------------
            DO 100 ITR=1,ITRMX
C--------------1./MONIN-OBUKKHOV LENGTH-SCALE---------------------------
        ZETALT=AMAX1(ZSLT*RLMO,ZTMIN)
        RLMO=ZETALT/ZSLT
        ZETALU=ZSLU*RLMO
C
        ZETAU=ZU*RLMO
        ZETAT=ZT*RLMO
C-----------------------------------------------------------------------
         IF(ILECH.EQ.0) THEN
          IF(RLMO.LT.0.)THEN
            XLU4=1.-16.*ZETALU
            XLT4=1.-16.*ZETALT
            XU4 =1.-16.*ZETAU
            XT4 =1.-16.*ZETAT
C
            XLU=SQRT(SQRT(XLU4))
            XLT=SQRT(SQRT(XLT4))
            XU =SQRT(SQRT(XU4))
            XT =SQRT(SQRT(XT4))
C
            PSMZ=PSPMU(XU)
            SIMM=PSPMU(XLU)-PSMZ+RLOGU
            PSHZ=PSPHU(XT)
            SIMH=PSPHU(XLT)-PSHZ+RLOGT
          ELSE
            ZETALU=AMIN1(ZETALU,ZTMAX)
            ZETALT=AMIN1(ZETALT,ZTMAX)
C
            PSMZ=PSPMS(ZETAU)
            SIMM=PSPMS(ZETALU)-PSMZ+RLOGU
            PSHZ=PSPHS(ZETAT)
            SIMH=PSPHS(ZETALT)-PSHZ+RLOGT
          ENDIF
         ELSE
C-----------------------------------------------------------------------
C ***** LECH'S FUNCTIONS ****
          IF(RLMO.LT.0.)THEN
            PSMZ=PSLMU(ZETAU)
            SIMM=PSLMU(ZETALU)-PSMZ+RLOGU
            PSHZ=PSLHU(ZETAT)
            SIMH=PSLHU(ZETALT)-PSHZ+RLOGT
          ELSE
            ZETALU=AMIN1(ZETALU,ZTMAX)
            ZETALT=AMIN1(ZETALT,ZTMAX)
C
            PSMZ=PSLMS(ZETAU)
            SIMM=PSLMS(ZETALU)-PSMZ+RLOGU
            PSHZ=PSLHS(ZETAT)
            SIMH=PSLHS(ZETALT)-PSHZ+RLOGT
          ENDIF
         ENDIF
C--------------BELJAARS CORRECTION FOR USTAR----------------------------
        USTAR=AMAX1(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
C--------------ZILITINKEVITCH FIX FOR ZT--------------------------------
        ZT=EXP(ZILFC*SQRT(USTAR*Z0))*Z0

        ZSLT=ZLM+ZT
        RLOGT=log(ZSLT/ZT)
C-----------------------------------------------------------------------
        USTARK=USTAR*VKRM
        AKMS=AMAX1(USTARK/SIMM,CXCH)
        AKHS=AMAX1(USTARK/SIMH,CXCH)
C-----------------------------------------------------------------------
        WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.)
        RLMN=ELFC*AKHS*DTHV/USTAR**3
C-----------------------------------------------------------------------
        RLMA=RLMO*WOLD+RLMN*WNEW
C-----------------------------------------------------------------------
C      IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT)    GO TO 110
C-----------------------------------------------------------------------
        RLMO=RLMA
C-----------------------------------------------------------------------
  100   CONTINUE
C-----------------------------------------------------------------------
  110   CONTINUE
C-----------------------------------------------------------------------
        RETURN
        END   
      SUBROUTINE SFLX(ICE,SATURATED,DT,Z,NSOIL, NROOT, SLDPTH,
     I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2,TBOT,CH,CHKFF,
     I VEGTYP,SOITYP,XLAI,REFKDT,SLOPTYP,SHDFAC,SNOALB,
     O ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX,
     2 T1,CMC,SMC,SH2O,SNOWH,STC,SNEQV,ALBEDO,
     O SOILW,SOILM)

      IMPLICIT NONE
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  SUB-DRIVER FOR FAMILY OF PHYSICS SUBROUTINES OF A
CC              SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL 
CC              MOISTURE, SOIL ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, 
CC              SNOWPACK WATER CONTENT,SNOWDEPTH, AND ALL TERMS
CC              OF THE SURFACE ENERGY BALANCE AND SURFACE WATER
CC              BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF 
CC              DOWNWARD RADIATION AND PRECIP)
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
C ------------    FROZEN GROUND VERSION     -------------------------------
C     NEW STATES: SH2O(NSOIL) - UNFROZEN SOIL MOISTURE
C                 SNOWH       - SNOW DEPTH
C
C -------------------------------------------------------------------------

C  SNOWH actual physical snow depth in m
C  SNEQV liquid water equivalent snow depth in m
C  Ratio SNOWH/SNEQV is the snow density.
C
C                  ARGUMENT LIST IN THE CALL SFLX 
C
C 1. CALLING STATEMENT
C 	CALL SFLX(ICE,LSATURATED,DT,Z,NSOIL,NROOT,SLDPTH,
C  I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2,TBOT,CH,CHKFF,
C  I IVEGTYP,ISOILTP,XLAI,REFKDT,
C  I ISLOPTYP, SHDFAC,SNOALB,
C  O ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX,
C  2 T1,CMC,SMC,
C  2 SH2O,SNOWH,
C  2 STC,SNEQV,ALBEDO,
C  O SOILW,SOILM)
C
C 2. INPUT
C                  *** GENERAL PARAMETERS ***
C
C          ICE: SEA-ICE FLAG  (=1: SEA-ICE, =0: LAND)
C    SATURATED: ATMOSPHERIC SATURATION FLAG 
C              (true: SATURATION -> ETP=0.0 IN SFLX
C               false: OTHERWISE)
C           DT: TIMESTEP (SEC)
C            Z: HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
C        NSOIL: NUMBER OF SOIL LAYERS
C        NROOT: NUMBER OF ROOT-ZONE LAYERS
C       SLDPTH: THE THICKNESS OF EACH SOIL LAYER (M) 
C
C                  *** ATMOSPHERIC VARIABLES ***
C
C         LWDN: LW DOWNWARD RADIATION (W M-2; POSITIVE, IF DOWNWARD)
C        SOLDN: SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, IF DOWNWARD)
C       SFCPRS: PRESSURE AT HEIGHT Z ABOVE GROUND(PASCALS)
C         PRCP: PRECIP RATE (KG M-2 S-1)
C       SFCTMP: AIR TEMPERATURE (K) AT HEIGHT Z ABOVE GROUND 
C          TH2: AIR POTENTIAL TEMPERATURE (K) AT GROUND SURFACE 
C           Q2: MIXING RATIO AT HEIGHT Z ABOVE GROUND (KG KG-1)
C        Q2SAT: SAT MIXING RATIO AT HEIGHT Z ABOVE GROUND (KG KG-1)
C       DQSDT2: SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP (KG KG-1 K-1)
C         TBOT: BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN AIR TEMPERATURE) 
C           CH: SFC EXCH COEF FOR HEAT AND MOISTURE (M S-1)
C        CHKFF: SFC EXCH COEF FOR HEAT AND MOISTURE MULTIPLIED BY RHO*CP
C
C                  *** CANOPY/SOIL CHARACTERISTICS ***
C
C       VEGTYP: VEGETATION TYPE (INTEGER INDEX)
C       SOITYP: SOIL TYPE (INTEGER INDEX)
C         XLAI: LEAF-AREA INDEX (UNITLESS, IN THE RANGE 0.-5.)
C       REFKDT: BASELINE VALUE FOR RUNOFF PARAMETER KDT
C      SLOPTYP: CLASS OF BOTTOM LAYER SLOPE (INTEGER INDEX)
C       SHDFAC: VEGETATION GREENNESS FRACTION 
C       SNOALB: ALBEDO UPPER BOUND OVER DEEP SNOW
C
C 3. INPUT AND OUTPUT
C
C                  *** STATE VARIABLES ***
C
C          T1: SKIN (GRND SFC) TEMPERATURE (K)
C         CMC: CANOPY MOISTURE CONTENT (M)
C  SMC(NSOIL): SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
C SH2O(NSOIL): UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
C       SNOWH:  SNOW DEPTH (M)
C  STC(NSOIL): SOIL TEMP (K)
C       SNEQV: WATER EQUIVALENT SNOW DEPTH (M) (formerly called SNODEP)
C      ALBEDO: SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
C    *** NOTE: THE ABOVE STATE VARIABLES ARE ALSO OUTPUT ***
C
C 4. OUTPUT
C          ETP: POTENTIAL EVAPORATION (W M-2)
C          ETA: ACTUAL LATENT HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM SURFACE)
C            H: SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM SURFACE)
C            S: SOIL HEAT FLUX (W M-2: NEGATIVE, IF DOWNWARD FROM SURFACE)
C      RUNOFF1: SURFACE RUNOFF (M S-1)
C      RUNOFF2: SUBSURFACE RUNOFF (M S-1)
C           Q1: EFFECTIVE MIXING RATIO AT GRND SFC ( KG KG-1) 
C        SNMAX: SNOW MELT (M) (WATER EQUIVALENT)
C        SOILM: TOTAL SOIL COLUMN WATER CONTENT (M)
C        SOILW: AVAILABLE SOIL MOISTURE (UNITLESS FRACTION BETWEEN 
C               SOIL SATURATION AND WILTING POINT) 

      INTEGER NSOLD
      PARAMETER (NSOLD=4)
C
      INTEGER K
      INTEGER ICE
      INTEGER NSOIL,VEGTYP,SOITYP,NROOT
      INTEGER SLOPTYP

      LOGICAL SNOWNG
      LOGICAL FRZGRA
      LOGICAL SATURATED

      REAL B
      REAL CFACTR
      REAL CH
      REAL CMC
      REAL CMCMAX
      REAL DF1
      REAL DKSAT

      REAL DT
      REAL DWSAT
      REAL EPSCA
      REAL ESD
      REAL ETA
      REAL ETP
      REAL F
      REAL F1
      REAL KDT
      REAL LWDN
      REAL PC
      REAL PRCP
      REAL Q1
      REAL Q2
      REAL Q2SAT
      REAL RCH
      REAL REFKDT
      REAL RR
      REAL RTDIS (NSOLD)
      REAL S
      REAL SFCPRS
      REAL SFCTMP
      REAL SHDFAC
      REAL SMC(NSOIL),STC(NSOIL),SLDPTH(NSOIL),SH2O(NSOIL)
      REAL SMCDRY
      REAL SMCMAX
      REAL SMCREF
      REAL SMCWLT
      REAL SNEQV
      REAL SNOWH
      REAL T1
      REAL T1MIN
      REAL T24
      REAL T2V
      REAL TBOT
      REAL TH2
      REAL Z
      REAL Z0
      REAL ZSOIL  ( NSOLD )
      REAL ZBOT

      REAL TOPT, TFREEZ, SOLDN, DQSDT2, CHKFF, XLAI, RUNOFF1, RUNOFF2
      REAL SNMAX, SOILW, RGL, SOILM
      REAL BETA, DRIP, EC, EDIR, ETT, FLX1, FLX2, FLX3, RHO, RUNOF, DEW
      REAL RIB, RUNOFF3, RSMAX, KZ, ALBEDO, RCMIN, HS, PSISAT
      REAL SLOPE, SNUP,SALP,SNOALB
      REAL QUARTZ, FRZX, SNDENS, SNCOND 
      REAL CSNOW, RSNOW, SNOFAC
      REAL SN_NEW, PRCP1, SICE, RC, EDIR1, EC1, ETT1, SOILWM
      REAL SOILWW
      REAL H

      PARAMETER (TFREEZ = 273.16)
      
CMIC$ TASKCOMMON RITE
      COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOF,
     &  	   DEW,RIB,RUNOFF3
C
C   INITIALIZATION

      RUNOFF1 = 0.0
      RUNOFF2 = 0.0
      RUNOFF3 = 0.0
      SNMAX = 0.0

      IF(ICE .EQ. 1) THEN

C SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
        DO KZ = 1, NSOIL
          ZSOIL(KZ)=-3.*FLOAT(KZ)/FLOAT(NSOIL)
        END DO

      ELSE

C CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO 
C BOTTOM OF EACH SOIL LAYER
        ZSOIL(1)=-SLDPTH(1)
        DO KZ = 2, NSOIL
          ZSOIL(KZ)=-SLDPTH(KZ)+ZSOIL(KZ-1)
        END DO

      ENDIF
         
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC   NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, 
CC   INCLUDING SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
        CALL REDPRM(VEGTYP,SOITYP, SLOPTYP, 
     +    CFACTR, CMCMAX, RSMAX, TOPT,REFKDT,KDT,
     O    Z0,SHDFAC,RCMIN,RGL,HS,
     +    ZBOT, FRZX, PSISAT, SLOPE, SNUP, SALP,
     O    B,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF,
     +    SMCDRY,F1,QUARTZ,
     O     RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL)
C
C prt      PRINT*,' 1ST DT=', DT
C prt      PRINT*,' ***********   INPUT PARAMETERS   *************'
C prt      PRINT*,'                                               '
C prt	   PRINT*,'VEGTYP= ',VEGTYP
C prt	   PRINT*,'SOITYP= ',SOITYP
C prt	   PRINT*,'SLOPTYP= ',SLOPTYP
C prt	   PRINT*,'CFACTR= ',CFACTR
C prt	   PRINT*,'CMCMAX= ',CMCMAX
C prt	   PRINT*,'RSMAX= ',RSMAX
C prt	   PRINT*,'TOPT= ',TOPT
C prt	   PRINT*,'REFKDT= ',REFKDT
C prt	   PRINT*,'KDT= ',KDT
C prt	   PRINT*,'Z0= ',Z0
C prt	   PRINT*,'SHDFAC= ',SHDFAC
C prt	   PRINT*,'RCMIN= ',RCMIN
C prt	   PRINT*,'RGL= ',RGL
C prt	   PRINT*,'HS= ',HS
C prt	   PRINT*,'ZBOT= ',ZBOT
C prt	   PRINT*,'FRZX= ',FRZX
C prt	   PRINT*,'PSISAT= ',PSISAT
C prt	   PRINT*,'SLOPE= ',SLOPE
C prt	   PRINT*,'SNUP= ',SNUP
C prt	   PRINT*,'SALP= ',SALP
C prt	   PRINT*,'B= ',B
C prt	   PRINT*,'DKSAT= ',DKSAT
C prt	   PRINT*,'DWSAT= ',DWSAT
C prt	   PRINT*,'SMCMAX= ',SMCMAX
C prt	   PRINT*,'SMCWLT= ',SMCWLT
C prt	   PRINT*,'SMCREF= ',SMCREF
C prt	   PRINT*,'SMCDRY= ',SMCDRY
C prt	   PRINT*,'F1= ',F1
C prt	   PRINT*,'QUARTZ= ',QUARTZ
C prt	   PRINT*,'RTDIS= ',RTDIS
C prt	   PRINT*,'SLDPTH= ',SLDPTH
C prt	   PRINT*,'ZSOIL= ',ZSOIL
C prt	   PRINT*,'NROOT= ',NROOT
C prt      PRINT*,'NSOIL= ',NSOIL

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  INITIALIZE MISC VARIABLES.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      S = 0.0
      SNOWNG = .FALSE.
      FRZGRA = .FALSE.

      ESD = SNEQV 
      IF(ICE .EQ. 1) THEN
        ESD = 0.01
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
C     IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
C     IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
C     TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF ( PRCP .GT. 0.0 ) THEN
        IF ( SFCTMP .LE. TFREEZ ) THEN
          SNOWNG = .TRUE.
        ELSE
          IF ( T1 .LE. TFREEZ ) FRZGRA = .TRUE.
        ENDIF
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     IF EITHER PRCP FLAG IS SET, CONVERT THE PRCP RATE FROM
C     KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH (METERS) AND ACCUM IT.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF(ESD .EQ. 0.0) THEN
        SNDENS = 0.0
        SNOWH = 0.0
	SNCOND = 1.0
      ELSE
        SNDENS=ESD/SNOWH
        SNCOND = CSNOW (SNDENS) 
      ENDIF     

      IF ( ( SNOWNG ) .OR. ( FRZGRA ) ) THEN
      
C ----------   FROZEN GROUND MODIFICATION   ----------------------
C          SNOW DENSITY & CONDUCTIVITY ESTIMATION 
C
      SN_NEW = PRCP * DT * 0.001
      CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
C   **  CSNOW IS FUNCTION TO CALCULATE SNOW CONDUCTIVITY
      SNCOND = CSNOW (SNDENS) 
C
C ----------------------------------------------------------------      
                   
        ESD =ESD + PRCP * DT * 0.001
        PRCP1 = 0.0
      ELSE
        PRCP1 = PRCP
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CHECK: IF ESD < 1.E-6 M. THE SNOW COVER WILL BE IGNORED.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF(ESD.LT.1.E-6) THEN
         SNMAX=ESD
         ESD = 0.0
         
         SNOWH = 0.0
         
         SNOWNG = .FALSE.
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC GRND HEAT FLUX FOR SNOW PACK CASE (WILL EFFECTIVELY
C     CALC THE HEAT FLUX BELOW THE SNOW PACK).  REMEMBER, IF
C     T1 > 273.16 K, THERE WILL BE SNOW MELT.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF ( ESD .GT. 0.0 ) THEN
        T1MIN = AMIN1 ( T1, TFREEZ )

C  ----   FROZEN GROUND MODIFICATION   ----------
         S = SNCOND * (T1MIN - STC(1) ) /SNOWH
C ---------------------------------------------------------
           
CC........S = 0.35 * ( T1MIN - STC(1) ) / ( 10.0 * ESD )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  BOUND THE HEAT FLUX MAGNITUDE UNDER SNOW PACK TO 100 W/M-2
C       SET THE MINIMUM FLUX WITHIN SNOW TO -200 W M-2
C       FOR NUMERICAL REASONS RELATED TO VERY SHALLOW SNOW DEPTH
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        S = MIN ( MAX ( S,-100.0 ), 100.0 )

      ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALCULATE GRND HEAT FLUX FOR CASE WITH NO SNOW PACK
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        CALL TDFCND ( DF1, SMC(1), B, F1,QUARTZ,SMCMAX,SH2O )
        
C ---------   FROZEN GROUND VERSION    ---------------------------
C         SOIL THERMAL DIFFUSIVITY CORRECTION
C
        SICE = SMC(1) - SH2O(1)
C__________________________________________________________________
C This IF-statement is no longer used if Peters-Lidard thermal
C conductivity is used for all conditions (frozen,unfrozen):
C        IF (SICE .NE. 0.0) CALL FRZCND ( DF1, SICE )

C ----------------------------------------------------------------        
                            
        S = DF1 * ( STC(1) - T1 ) / ( 0.5 * ZSOIL(1) )
      ENDIF

C UPDATE ALBEDO AND RADIATION FOR CALL TO PENMAN, NOPAC and SNOPAC
C CALCULATE SNOW ALBEDO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C 
      IF ( (SNUP .GT. 0.0) .AND. (SNEQV .LT. SNUP) ) THEN
            RSNOW = SNEQV/SNUP
            SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP)) 
CCCC........SNOFAC = AMIN1(SNODEP/0.02, 1.0)
      ELSE
            SNOFAC = 1.0
      ENDIF
            
          ALBEDO = ALBEDO+(1.0-SHDFAC)*SNOFAC*(SNOALB-ALBEDO)
	  
          F = SOLDN*(1.0-ALBEDO) + LWDN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALC VIRTUAL TEMPS AND POTENTIAL TEMPS FROM AIR TEMP
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      T2V = SFCTMP * ( 1.0 + 0.61 * Q2 )
CC....TH2 = SFCTMP + ( 0.0097545 * Z )....TH2 is passed in the arg.list

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION
C     (AND OTHER PARTIAL PRODUCTS AND SUMS FOR LATER CALCULATIONS).
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       CALL PENMAN ( SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,F,T24,S,Q2,
     &              Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,DQSDT2)
       IF(SATURATED) ETP = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT 
C     INTO PC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF(SHDFAC .GT. 1.E-6) THEN
      
C  ----  FROZEN GROUND EXTENSION: SMC WAS REPLACED BY SH2O  -----
C      
        CALL CANRES(SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL,
     &            SMCWLT,SMCREF,RCMIN,RC,PC,NROOT,Q2SAT,DQSDT2,
     &            TOPT,RSMAX,RGL,HS,XLAI)
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       IF POINT IS TO BE FULLY PROCESSED....
C       BASED ON WHETHER OR NOT THERE IS A SNOW PACK, CALL SNOPAC OR
C       NOPAC TO CALCULATE AN ESTIMATE OF ACTUAL EVAPOTRANSPIRATION,
C       UPDATE SOIL MOISTURE, AND UPDATE SOIL TEMPERATURE.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        IF ( ESD .EQ. 0.0 ) THEN

          CALL NOPAC ( ETP, ETA, PRCP, SMC, SMCMAX, SMCWLT,
     &                 SMCREF,SMCDRY, CMC, CMCMAX, NSOIL, DT, SHDFAC,
     &                 Q1,Q2,T1,SFCTMP,T24,TH2,F,F1,S,STC,
     &                 EPSCA, B, PC, RCH, RR,  CFACTR,
     +                 SH2O, SLOPE, KDT, FRZX, PSISAT,
     &                 ZSOIL, DKSAT, DWSAT, TBOT,RUNOFF1,RUNOFF2,
     &                 RUNOFF3, EDIR1, EC1, ETT1,NROOT,ICE,RTDIS,QUARTZ)

        ELSE

          CALL SNOPAC ( ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT,
     &                SMCREF, SMCDRY, CMC, CMCMAX, NSOIL, DT, Q1,
     &                Q2,T1,SFCTMP,T24,TH2,F,F1,S,STC,EPSCA,SFCPRS,
     &                B, PC, RCH, RR, CFACTR, SALP, ESD,
     +                SNOWH, SH2O, SLOPE, KDT, FRZX, PSISAT, SNUP,
     &                ZSOIL, DWSAT, DKSAT, TBOT, SHDFAC,RUNOFF1,
     &                RUNOFF2,RUNOFF3,EDIR1,EC1,ETT1,NROOT,SNMAX,ICE,
     &                RTDIS,QUARTZ)
        
        ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  SAVE UPDATED WATER EQUIVALENT SNOW DEPTH FOR RETURN.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SNEQV = ESD

C UPDATE ALBEDO USING UPDATED SNOW FOR RETURN TO PARENT 
C ATMOSPHERIC MODEL 

      IF ( (SNUP .GT. 0.0) .AND. (SNEQV .LT. SNUP) ) THEN
            RSNOW = SNEQV/SNUP
            SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP)) 
CCCC........SNOFAC = AMIN1(SNODEP/0.02, 1.0)
      ELSE
            SNOFAC = 1.0
      ENDIF
            
          ALBEDO = ALBEDO+(1.0-SHDFAC)*SNOFAC*(SNOALB-ALBEDO)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC SURFACE HEAT FLUX FOR RETURN TO PARENT MODEL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      H = -CHKFF * ( TH2 - T1 )
C
C  CONVERT ETA FROM KG M-2 S-1 TO W M-2
C
      ETA = ETA*2.5E+6
      ETP = ETP*2.5E+6

C CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
C         S>0: WARM THE SURFACE  (NIGHT TIME)
C         S<0: COOL THE SURFACE  (DAY TIME)

      S=-1.0*S      
C
C  CONVERT RUNOFF3 FROM M TO M S-1
C
      RUNOFF3 = RUNOFF3/DT
      RUNOFF2 = RUNOFF2+RUNOFF3
C
C TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE 
C SOIL MOISTURE AVAILABILITY FRACTION

      SOILM=-1.0*SMC(1)*ZSOIL(1)
      
      DO K = 2, NSOIL
        SOILM=SOILM+SMC(K)*(ZSOIL(K-1)-ZSOIL(K))
      END DO
      SOILWM=-1.0*(SMCMAX-SMCWLT)*ZSOIL(1)
      SOILWW=-1.0*(SMC(1)-SMCWLT)*ZSOIL(1)
      DO K = 2, NROOT
        SOILWM=SOILWM+(SMCMAX-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
        SOILWW=SOILWW+(SMC(K)-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
      END DO
      SOILW=SOILWW/SOILWM
C
      RETURN
      END
      SUBROUTINE SHFLX ( S,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,TBOT,
     +             SMCWLT, PSISAT, SH2O, B,F1,DF1,ICE,QUARTZ )

      IMPLICIT NONE

C ------------    FROZEN GROUND VERSION    --------------------------
C   NEW STATES ADDED: SH2O & PARAMETER SMCWLT & PSISAT
C
C ------------------------------------------------------------------- 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE SOIL HEAT FLUX.  SOIL TEMPERATURE CONTENT
CC    =======   (AN ARRAY CONTAINING SOIL TEMPERATURES) IS ALSO UPDATED.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOLD
      PARAMETER ( NSOLD = 4 )
      INTEGER NSOIL
      REAL B
      REAL DF1
      REAL DT
      REAL F1
      REAL RHSTS ( NSOLD )
      REAL S
      REAL SMC   ( NSOIL )
      REAL SH2O  ( NSOIL )
      REAL SMCMAX
      REAL SMCWLT
      REAL STC   ( NSOIL )
      REAL STCOUT ( NSOLD )
      REAL T1
      REAL TBOT
      REAL YY
      REAL ZSOIL ( NSOIL )
      REAL ZZ1

      REAL PSISAT, QUARTZ

      INTEGER ICE, IFRZ, I

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF(ICE.EQ.1) THEN
       CALL HRTICE(RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,B,
     *             F1,DF1)
       CALL HSTEP ( STC, STC, RHSTS,DT,NSOIL )
      ELSE
      
       CALL HRT(RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,
       
C -------------    FROZEN GROUND VERSION    -----------------------
C
     +          SMCWLT, PSISAT, SH2O, DT,
C -----------------------------------------------------------------                  
     +          B,F1,DF1,QUARTZ)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     HSTEP ROUTINE CALCS/UPDATES SOIL TEMPS BASED ON RHSTS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        CALL HSTEP ( STCOUT,STC,RHSTS,DT,NSOIL )

C ----------------     FROZEN GROUND VERSION    --------------------
C   DOUBLE CALLING HRT & HSTEP TO REDUCE NOISE FLUCTUATION, THE SAME 
C    AS SOIL MOISTURE CALCULATION. DOUBLE CALL WILL BE ONLY IF THERE
C    IS A ICE ON AT LEAST ONE SOIL LAYER
C
        IFRZ = 0
        DO I = 1,NSOIL
          IF( (SMC(I) - SH2O(I) ) .GT. 0.0 ) IFRZ = 1
        END DO
        
        IF ( IFRZ .EQ. 1 ) THEN     
        DO I = 1,NSOIL
         STC(I) = 0.5 * ( STCOUT(I) + STC(I) )
        END DO
        
        CALL HRT(RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,
       
C -------------    FROZEN GROUND VERSION    -----------------------
C
     +          SMCWLT, PSISAT, SH2O, DT,
C -----------------------------------------------------------------                  
     +          B,F1,DF1,QUARTZ)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     HSTEP ROUTINE CALCS/UPDATES SOIL TEMPS BASED ON RHSTS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        CALL HSTEP ( STC,STC,RHSTS,DT,NSOIL )
       ELSE
        DO I = 1,NSOIL
          STC(I) = STCOUT(I)
        END DO
       ENDIF                  
C ------------------------------------------------------------------
        
       ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE GRND (SKIN) TEMPERATURE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE SFC SOIL HEAT FLUX
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      S = DF1 * (STC(1) - T1) / (0.5 * ZSOIL(1))

      RETURN
      END
      SUBROUTINE SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL,
     &     SH2O, SLOPE, KDT, FRZFACT,
     &     SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC,CMCMAX,
     &     SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, EDIR1, EC1, 
     &     ETT1, SFCTMP,Q2,NROOT,RTDIS)

      IMPLICIT NONE

C ------------    FROZEN GROUND VERSION    --------------------------
C   NEW STATES ADDED: SH2O, AND FROZEN GROUD CORRECTION FACTOR, FRZFACT
C   AND PARAMETER SLOPE 
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE SOIL MOISTURE FLUX.  THE SOIL MOISTURE
CC    =======   CONTENT (SMC - A PER UNIT VOLUME MEASUREMENT) IS A
CC              DEPENDENT VARIABLE THAT IS UPDATED WITH PROGNOSTIC EQNS.
CC              THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOLD
      PARAMETER ( NSOLD = 4 )
      INTEGER K
      INTEGER NSOIL
      REAL B
      REAL BETA
      REAL CFACTR
      REAL CMC
      REAL CMCMAX
      REAL DEW
      REAL DKSAT
      REAL DRIP
      REAL DT
      REAL DWSAT
      REAL EC
      REAL EDIR
      REAL ET     ( NSOLD )
      REAL ETA1
      REAL ETP1
      REAL ETT
      REAL EXCESS
      REAL FLX1
      REAL FLX2
      REAL FLX3
      REAL KDT
      REAL PC
      REAL PCPDRP
      REAL PRCP1
      REAL RHO
      REAL RHSCT
      REAL RHSTT  ( NSOLD )
      REAL RIB
      REAL RTDIS (NSOIL)
      REAL RUNOF
      REAL RUNOFF,RUNOXX3
      REAL SHDFAC
      REAL SMC    ( NSOIL )

C ---------------    FROZEN GROUND VERSION     ---------------------
      
      REAL SH2O   ( NSOIL )
      REAL SICE   ( NSOLD )
      REAL SH2OA   ( NSOLD)
      REAL SH2OFG  ( NSOLD ) 
C -------------------------------------------------------------------
           
      REAL SMCDRY
      REAL SMCMAX
      REAL SMCREF
      REAL SMCWLT
      REAL TRHSCT
      REAL ZSOIL  ( NSOIL )

      REAL SLOPE, FRZFACT, RUNOFF1, RUNOFF2, RUNOFF3, EDIR1, EC1
      REAL ETT1, SFCTMP, Q2, DUMMY, CMC2MS, DEVAP

      INTEGER NROOT, I

CMIC$ TASKCOMMON RITE
      COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOF,
     &             DEW,RIB,RUNOXX3

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     EXECUTABLE CODE BEGINS HERE....IF THE POTENTIAL EVAPOTRANS-
C     PIRATION IS GREATER THAN ZERO...
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DUMMY=0.
      EDIR = 0.
      EC = 0.
      ETT = 0.
      DO 10 K = 1, NSOIL
        ET ( K ) = 0.
   10 CONTINUE

      IF ( ETP1 .GT. 0.0 ) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C --------------     FROZEN GROUND VERSION     ---------------------
C   SMC STATES WERE REPLACED BY SH2O STATES
C
        EDIR = DEVAP ( ETP1, SH2O(1), ZSOIL(1), SHDFAC, SMCMAX,
C ------------------------------------------------------------------        
C        EDIR = DEVAP ( ETP1, SMC(1), ZSOIL(1), SHDFAC, SMCMAX,

     &                B, DKSAT, DWSAT, SMCDRY,SMCREF, SMCWLT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT
C       TRANSPIRATION, AND ACCUMULATE IT FOR ALL SOIL LAYERS.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        ETT = 0.

        IF(SHDFAC.GT.0.0) THEN
        
C --------------     FROZEN GROUND VERSION     ---------------------
C   SMC STATES WERE REPLACED BY SH2O STATES
C
          CALL TRANSP ( ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT,
C ------------------------------------------------------------------
C          CALL TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT,
     &         CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)

          DO 20 K = 1 , NSOIL
            ETT = ETT + ET ( K )
   20     CONTINUE
        ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALCULATE CANOPY EVAPORATION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1

C********  EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE
C          WATER ON THE CANOPY. MODIFIED BY F.CHEN ON 10/18/94
C********

        CMC2MS = CMC / DT
        EC = MIN ( CMC2MS, EC )

      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        EDIR1=EDIR
	EC1=EC
	ETT1=ETT

      ETA1 = EDIR + ETT + EC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      RHSCT = SHDFAC * PRCP1 - EC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CONVERT RHSCT (A RATE) TO TRHSCT (AN AMT) AND ADD IT TO EXISTING
C     CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP
C     AND WILL FALL TO THE GRND.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DRIP = 0.
      TRHSCT = DT * RHSCT
      EXCESS =  CMC + TRHSCT
      IF ( EXCESS .GT. CMCMAX ) DRIP = EXCESS - CMCMAX

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT
C     GOES INTO THE SOIL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PCPDRP = (1. - SHDFAC) * PRCP1 + DRIP / DT

CC       PRINT*,' **************** SMLX ******************'
CC       PRINT*,' PCPDRP=', PCPDRP, ' EDIR=', EDIR,' ET=', ET,
CC     *      'SMC(1)=', SMC(1), 'SMC(2)=', SMC(2), ' PRCP1=', PRCP1

C ---------------     FROZEN GROUND VERSION     --------------------
C    STORE ISE CONTENT AT EACH LAYER BEFORE CALLING SRT & SSTEP
C
          DO I = 1,NSOIL
            SICE(I) = SMC(I) - SH2O(I)
          END DO
C ------------------------------------------------------------------
            
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
C     TENDENCY EQUATIONS. ITERATE TWICE THRU THE CALLS. THIS ELIMINATES
C     2-DELTA-T OSCILLATIONS IN THE SMC VALUES EVEN WITH TIMESTEPS OF
C     UP TO 1 HOUR. THE 1ST SET OF CALLS YIELDS 1ST GUESS SMC VALUES,
C     THE 2ND SET OF CALLS YIELDS FINAL SMC VALUES.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C	print*,'SMFLX, BEFORE srt, sstep'
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o
      IF ( PCPDRP .GT. 0.0 ) THEN

C ---------------    FROZEN GROUND VERSION       ---------------------
C    SMC STATES REPLACED BY SH2O STATES IN SRT SUBR.
C    SH2O & SICE STATES INCLUDED IN SSTEP SUBR.
C    FROZEN GROUND CORRECTION FACTOR, FRZFACT, ADDED
C    ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
C
        CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL,
     &             DWSAT,DKSAT,SMCMAX, B, RUNOFF1, 
     +             RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT, SICE)
     
C	print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER srt - 1'
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o
	
        CALL SSTEP ( SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX,
     &               CMCMAX, RUNOFF3, ZSOIL, SMC, SICE )

C	print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER sstep - 1'
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o
	
        DO 30 K = 1, NSOIL
          SH2OA(K) = ( SH2O(K) + SH2OFG(K) ) * 0.5
   30   CONTINUE
	
        CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL,
     &             DWSAT,DKSAT,SMCMAX, B, RUNOFF1,
     +             RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT, SICE)

C	print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER srt - 2'
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o
	
        CALL SSTEP ( SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX,
     &               CMCMAX, RUNOFF3, ZSOIL,SMC,SICE)
     
C	print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER sstep - 2'
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o

      ELSE

        CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL,
     &             DWSAT,DKSAT,SMCMAX, B, RUNOFF1,
     +             RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT, SICE)
     

C	print*,'SMFLX, AFTER srt '
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o
	
        CALL SSTEP ( SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX,
     &               CMCMAX, RUNOFF3, ZSOIL,SMC,SICE)
     
C	print*,'SMFLX, AFTER sstep '
C	print*,'smc=' , smc
C	print*,'sh2o=', sh2o

C -------------------------------------------------------------------- 

C  *****    BELOW IS OLD CODE    *************************************
C        CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMC,NSOIL,PCPDRP,ZSOIL,
C     &       DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT)        
C        CALL SSTEP ( SMCFG,SMC,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX,
C     &               CMCMAX, RUNOFF3, ZSOIL)
C
C        DO 30 K = 1, NSOIL
C          SMCA(K) = ( SMC(K) + SMCFG(K) ) * 0.5
C   30   CONTINUE
C
C        CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMCA,NSOIL,PCPDRP,ZSOIL,
C     &      DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT)
C        CALL SSTEP ( SMC,SMC,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX,
C     &               CMCMAX, RUNOFF3, ZSOIL)
C      ELSE
C
C        CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMC,NSOIL,PCPDRP,ZSOIL,
C     &        DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT)
C        CALL SSTEP ( SMC,SMC,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX,
C     &               CMCMAX, RUNOFF3, ZSOIL)
C ********************************************************************

      ENDIF

      RUNOF = RUNOFF
      RETURN
      END
      FUNCTION SNKSRC ( TUP,TM,TDN, SMC, SH2O, HCPCT, ZSOIL,
     +                    SMCMAX, PSISAT, B, F1, DT, K, QUARTZ) 
     
      IMPLICIT NONE
     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION
CC    =======   EQUATION. (SH2O) IS AVAILABLE LIQUED WATER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER  K
      
      REAL  B
      REAL  DF
      REAL  DFH2O
      REAL  DFICE
      REAL  DH2O
      REAL  DT
      REAL  DZ
      REAL  DZH
      REAL  F1
      REAL  FREE
      REAL  FRH2O
      REAL  HLICE
      REAL  PSISAT
      REAL  QDN
      REAL  QTOT
      REAL  QUARTZ
      REAL  QUP
      REAL  SH2O
      REAL  SICE
      REAL  SMC
      REAL  SMCMAX
      REAL  SNKSRC
      REAL  T0
      REAL  TAVG
      REAL  TDN
      REAL  TM
      REAL  TUP
      REAL  TZ
      REAL  X0
      REAL  XDN
      REAL  XH2O
      REAL  XUP
      REAL  HCPCT
      REAL  ZSOIL (4)
      
      PARAMETER (HLICE=3.3350 E5)
      PARAMETER (DH2O =1.0000 E3)
      PARAMETER (  T0 =2.7316 E2)
      
C  **************************************************************************
C  ***  OPTION TO RUN NON-FROZEN VERSION:  UNCOMMENT NEXT TWO STATEMENTS ****
C  *******************        WARNING !!!!!        **************************  
C  ***  TO RUN FROZEN GROUND VERSION, NEXT TWO STATEMENTS SHOULD BE      ****
C  ***  COMMENTED                                                        **** 
C
CC      SNKSRC = 0.0
CC      GOTO 77
C  ************************************************************************** 

      SICE=SMC-SH2O      
      IF(K.EQ.1) THEN
        DZ=-ZSOIL(1)
      ELSE
        DZ=ZSOIL(K-1)-ZSOIL(K)
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  CALCULATE DIFFUSIVITIES OF THAWED AND FROZEN GROUND
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  WARNING !!!!!, by Pablo Grunmann, 10/98
C ATTENTION TO THE CALLING ARGUMENT:
C 
C    1. USE SMC instead of SH2O to use VKOREN's Frozen DF (call 
C        and IF statement commented out with C VK)
C    2. USE SH2O in calling argument list for Petrs-Lidard (call 
C        commented out with C P-L)
C                 !!!                      !!!!
C P-L      CALL TDFCND ( DFH2O, SMC, B, F1,QUARTZ,SMCMAX,SH2O) 
C VK       CALL TDFCND ( DFH2O, SMC, B, F1,QUARTZ,SMCMAX, SMC) 
      CALL TDFCND ( DFH2O, SMC, B, F1,QUARTZ,SMCMAX,SH2O)
      DFICE = DFH2O
C Frozen conductivity for original VK approach
C VK      IF ( SICE .GT. 0.0 ) CALL FRZCND ( DFICE, SICE )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  HEAT FLUX FROM THE TOP BOUNDARY OF LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      IF(TUP .LE. T0) THEN
        DF=DFICE
      ELSE
        DF=DFH2O
      ENDIF
      QUP=-DF*(T0-TUP)/(0.5*DZ)
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  HEAT FLUX FROM THE BOTTOM BONDARY OF LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF(TDN .LE. T0) THEN
        DF=DFICE
      ELSE
        DF=DFH2O
      ENDIF
      QDN=-DF*(T0-TDN)/(0.5*DZ)
      QTOT=QUP+QDN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALCULATE POTENTIAL REDUCTION OF LIQUED WATER CONTENT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      XH2O=QTOT*DT/(DH2O*HLICE*DZ) + SH2O
            
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     ESTIMATE UNFROZEN WATER AT TEMPERATURE TAVG,
C     AND CHECK IF CALCULATED WATER CONTENT IS REASONABLE 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
        
C  ****   NEW CALCULATION OF AVERAGE TEMPERATURE (TAVG)   **********
C  ****   IN FREEZING/THAWING LAYER USING UP, DOWN, AND MIDDLE   ***
C  ****   LAYER TEMPERATURES (TUP, TDN, TM)               **********
   
      DZH=DZ*0.5
      IF(TUP .LT. 0.) THEN

        IF(TM .LT. 0.) THEN

          IF(TDN .LT. 0.) THEN
C1  ******   TUP, TM, TDN < 0.  CASE  ************          
            TAVG=(TUP+2.0*TM+TDN)/4.
             
C       C111  ****   END TDN_111  THEN

          ELSE
C2  ******   TUP & TM < 0.,  TDN > 0. CASE   *****
            X0=(T0-TM)*DZH/(TDN-TM)
            TAVG=0.5*(TUP*DZH+TM*(DZH+X0)+T0*(2.*DZH-X0))/DZ
                       
C       C112   ****  END TDN_111  ELSE   ****
          ENDIF      

C   C11  ***********  END TM_1 THEN   ***********
        ELSE
        
          IF(TDN .LT. 0.) THEN
C3  *******  TUP < 0.  TM > 0.  TDN < 0. CASE   ****
            XUP=(T0-TUP)*DZH/(TM-TUP)
            XDN=DZH-(T0-TZ)*DZH/(TDN-TM)
            TAVG=0.5*(TUP*XUP+T0*(2.*DZ-XUP-XDN)+TDN*XDN)/DZ
                      
C       C121   ****   END TDN_121  THEN  *****

          ELSE
C4   ******  TUP < 0  TM > 0  TDN > 0  CASE   ******
            XUP=(T0-TUP)*DZH/(TM-TUP)
            TAVG=0.5*(TUP*XUP+T0*(2.*DZ-XUP))/DZ
                      

C       C122   ***   END TDN_121  ELSE   ***
          ENDIF   
        
C   C12  ***********  END TM_1 ELSE   ********** 
        ENDIF

C1    ********************  END TUP THEN   ***********    
      ELSE

        IF(TM .LT. 0.) THEN

          IF(TDN .LT. 0.) THEN
C5   *****  TUP > 0  TM < 0 TDN < 0 CASE   *********
            XUP=DZH-(T0-TUP)*DZH/(TM-TUP)
            TAVG=0.5*(T0*(DZ-XUP)+TM*(DZH+XUP)+TDN*DZH)/DZ
                      
C       C211   ****   END TDN_211  THEN  *****

          ELSE
C6   *****  TUP > 0  TM < 0  TDN > 0 CASE  *********
            XUP=DZH-(T0-TUP)*DZH/(TM-TUP)
            XDN=(T0-TM)*DZH/(TDN-TM)
            TAVG=0.5*(T0*(2.*DZ-XUP+XDN)+TM*(XUP+XDN))/DZ
                                   

C       C212   ***   END TDN_211  ELSE   ***
          ENDIF   

C   C21   ************   END TM_2  THEN   ************
        ELSE

          IF(TDN .LT. 0.) THEN
C7   *****  TUP > 0  TM > 0  TDN < 0 CASE   ********
            XDN=DZH-(T0-TM)*DZH/(TDN-TM)
            TAVG=(T0*(DZ-XDN)+0.5*(T0+TDN)*XDN)/DZ
                  
C       C221   ****   END TDN_221  THEN  *****

          ELSE
C8   *****  TUP > 0  TM > 0  TDN > 0 CASE   ********
            TAVG=(TUP+2.*TM+TDN)/4.
                      
C       C222   ***   END TDN_221  ELSE   ***
          ENDIF           

C   C22   ************   END TM_2  ELSE   ************
        ENDIF

C   *********************  END TUP ELSE  **************
      ENDIF                      
C 777   
      
100   FREE=FRH2O(TAVG, SMC, SH2O, SMCMAX, B, PSISAT )

      IF ( XH2O .LT. SH2O .AND. XH2O .LT. FREE) THEN 
         IF ( FREE .GT. SH2O ) THEN
              XH2O = SH2O
          ELSE
              XH2O = FREE
          ENDIF
      ENDIF
              
      IF ( XH2O .GT. SH2O .AND. XH2O .GT. FREE )  THEN
         IF ( FREE .LT. SH2O ) THEN
              XH2O = SH2O
          ELSE
              XH2O = FREE
          ENDIF
      ENDIF 
             
      IF(XH2O .LT. 0. ) XH2O=0.
10    IF(XH2O .GT. SMC) XH2O=SMC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALCULATE SINK/SOURCE TERM AND REPLACE PREVIOUS WATER CONTENT 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     
      SNKSRC=-DH2O*HLICE*DZ*(XH2O-SH2O)/DT

        IF(XH2O .LT. 0.) THEN
          XH2O=SH2O
        ENDIF

      SH2O=XH2O
      
77    RETURN
      END
      SUBROUTINE SNOPAC ( ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT,
     &     SMCREF, SMCDRY, CMC, CMCMAX, NSOIL, DT, Q1,
     &     Q2,T1,SFCTMP,T24,TH2,F,F1,S,STC,EPSCA,SFCPRS,
     &     B, PC, RCH, RR, CFACTR, SALP, ESD,
     +     SNOWH, SH2O, SLOPE, KDT, FRZFACT, PSISAT,SNUP,
     &     ZSOIL, DWSAT, DKSAT, TBOT, SHDFAC, RUNOFF1,
     &     RUNOFF2,RUNOFF3,EDIR1,EC1,ETT1,NROOT,SNMAX,ICE,
     &     RTDIS,QUARTZ)

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE
CC    =======   SOIL MOISTURE CONTENT AND SOIL HEAT CONTENT VALUES FOR
CC              THE CASE WHEN A SNOW PACK IS PRESENT.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOIL

      LOGICAL SNOWNG

      REAL B
      REAL BETA
      REAL CFACTR
      REAL CMC
      REAL CMCMAX
      REAL CP
      REAL CPH2O
      REAL CPICE
      REAL DENOM
      REAL DEW
      REAL DF1
      REAL DKSAT
      REAL DRIP
      REAL DT
      REAL DWSAT
      REAL EC
      REAL EDIR
      REAL EPSCA
      REAL ESD
      REAL SNOWH
      REAL ETA
      REAL ETA1
      REAL ETP
      REAL ETP1
      REAL ETP2
      REAL ETT
      REAL EX
      REAL F
      REAL FLX1
      REAL FLX2
      REAL FLX3
      REAL F1
      REAL KDT
      REAL LSUBF
      REAL LSUBC
      REAL LSUBS
      REAL PC
      REAL PRCP
      REAL PRCP1
      REAL Q1
      REAL Q2
      REAL RCH
      REAL RHO
      REAL RIB
      REAL RR
      REAL RTDIS (NSOIL)
      REAL RUNOFF
      REAL S
      REAL S1
      REAL SFCTMP
      REAL SHDFAC
      REAL SIGMA
      REAL SMC     ( NSOIL )
      REAL SH2O    ( NSOIL )
      REAL SMCDRY
      REAL SMCMAX
      REAL SMCREF
      REAL SMCWLT
      REAL SNMAX
      REAL STC     ( NSOIL )
      REAL T1
      REAL T11
      REAL T12
      REAL T12A
      REAL T12B
      REAL T24
      REAL TBOT
      REAL TH2
      REAL YY
      REAL ZSOIL( NSOIL )
      REAL ZZ1

      REAL TFREEZ, SALP, SFCPRS, SLOPE, FRZFACT, PSISAT, SNUP
      REAL RUNOFF1, RUNOFF2, RUNOFF3,RUNOXX3
      REAL EDIR1, EC1, ETT1, QUARTZ
      REAL SNDENS, SNCOND, RSNOW, SNCOVER, QSAT, ETP3, SEH, T14
      REAL SICE1, CSNOW

      INTEGER NROOT, ICE


CMIC$ TASKCOMMON RITE

      COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF,
     &             DEW,RIB,RUNOXX3
      PARAMETER(CP=1004.5,CPH2O=4.218E+3,CPICE=2.106E+3,
     &     LSUBF=3.335E+5,LSUBC=2.5E+6,LSUBS=2.83E+6,SIGMA=5.67E-8)
     
C ------------    FROZEN GROUND VERSION     ------------------------
C   
C  
      PARAMETER ( TFREEZ = 273.16)
C-------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     EXECUTABLE CODE BEGINS HERE...
C     CONVERT ETP FM  KG M-2 S-1  TO  M S-1  AND THEN TO AN AMOUNT
C     (M) AND CALL IT AN EFFECTIVE SNOWPACK REDUCTION AMT, ETP2 (M).
C     THIS IS THE AMT THE SNOWPACK WOULD BE REDUCED DUE TO EVAPORATION
C     FROM THE SNOW SFC DURING THE TIMESTEP.  EVAPORATION WILL
C     PROCEED AT THE POTNTL RATE UNLESS THE SNOW DEPTH IS LESS THAN
C     THE EXPECTED SNOWPACK REDUCTION.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PRCP1=PRCP1*0.001
      ETP2 = ETP * 0.001 * DT
      BETA = 1.0
      IF(ICE.NE.1) THEN
        IF ( ESD .LT. ETP2 ) THEN
          BETA = ESD / ETP2
        ENDIF
      ENDIF
      DEW = 0.0
      IF ( ETP .LT. 0.0 ) THEN
        DEW = -ETP * 0.001
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC AN 'EFFECTIVE SNOW-GRND SFC TEMP' BASED ON HEAT FLUXES
C     BTWN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      T12A = ((F - SIGMA * T24) / RCH+TH2-SFCTMP-BETA*EPSCA) / RR
      
C ---------------    FROZEN GROUND VERSION    -----------------------
C  VARIABLE SNOW DENSITY & CONDUCTIVITY

      SNDENS = ESD / SNOWH
      SNCOND = CSNOW ( SNDENS )
      T12B = SNCOND * STC(1) / ( SNOWH * RR * RCH )
      DENOM = 1.0 + SNCOND / ( SNOWH * RR * RCH )
C -------------------------------------------------------------------      
           
C      T12B = 0.35 * STC(1) / ( 10.0 * ESD * RR * RCH )
C      DENOM = 1.0 + 0.35 / ( 10.0 * ESD * RR * RCH )

      T12 = (SFCTMP + T12A + T12B ) / DENOM

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     RECALCULATE THE HEAT FLUX INTO OR OUT OF THE SOIL, 
C     AND SET THE EFFECTIVE POTENTIAL EVAPOTRANSP TO ZERO.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CC.....S = 0.35 * ( T12 - STC(1) ) / ( 10.0 * ESD )
      
C ---------------    FROZEN GROUND VERSION    -----------------------
C  VARIABLE SNOW DENSITY & CONDUCTIVITY
C
      S = SNCOND * ( T12 - STC(1) ) / ( SNOWH )
C -------------------------------------------------------------------
      
      S = MIN ( MAX ( S,-100.0 ), 100.0 )
      IF ( T12 .LE. TFREEZ ) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BLO FREEZING,
C       NO SNOW MELT WILL OCCUR.  SET THE SKIN TEMP TO THIS
C       EFFECTIVE TEMP AND SET THE EFFECTIVE PRECIP TO ZERO.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        ESD = AMAX1 ( 0.0 , ESD - ETP2 )
        T1 = T12
        FLX3 = 0.0
        EX = 0.0
        SNMAX = 0.0

      ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW
C       MELT WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNMAX.
C       REVISE THE EFFECTIVE SNOW DEPTH.  REVISE THE SKIN TEMP BECAUSE
C       IT WOULD HAVE CHGD DUE TO THE LATENT HEAT RELEASED BY THE
C       MELTING. CALC THE LATENT HEAT RELEASED, FLX3. SET THE EFFECTIVE
C       PRECIP, PRCP1 TO THE SNOW MELT RATE, EX FOR USE IN SMFLX.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C THIS IS THE OLD VERSION OF SNOW MELT
C        EX = 0.001 * RCH * RR * ( T12 - 273.16 ) / LSUBF
C        SNMAX = EX * DT
C        EX = AMIN1 ( ESD / DT, EX )
C        ESD = AMAX1 ( ESD - SNMAX, 0.0 )
C        T1 = AMAX1 ( 273.16, (T12-(LSUBF * EX * 1000.0 / (RCH*RR))))
C        FLX3 = LSUBF * EX * 1000.0
C        PRCP1 = PRCP1 + EX
C BELOW IS THE CORRECTION OF EVAPORATION CONSIDERING SNOW MELTING
C                    F. CHEN (96/3/18)

C        T1 = 273.16

C ------------    FROZEN GROUND VERSION    --------------------------
C   ADJUSTMENT T1 TO ACCOUNT FOR SNOW PATCHES
C
        IF ( SNUP .GT. 0.0 .AND. ESD .LT. SNUP ) THEN
          RSNOW = ESD / SNUP
          SNCOVER = 1.-( EXP(-SALP*RSNOW)-RSNOW*EXP(-SALP))
C          SNCOVER = AMIN1 ( ESD / SNUP, 1.)
        ELSE
          SNCOVER = 1.
        ENDIF
            
        T1 = TFREEZ * SNCOVER + T12 * ( 1.0 - SNCOVER )
C -------------------------------------------------------------------
        
        QSAT = (0.622*6.11E2)/(SFCPRS-0.378*6.11E2)
        ETP = RCH*(QSAT-Q2)/CP
        ETP2 = ETP*0.001*DT
        BETA = 1.0
        IF ( ESD .LE. ETP2 ) THEN
          BETA = ESD / ETP2
          ESD = 0.0
          SNMAX = 0.0
          EX = 0.0
        ELSE
          ESD = MAX ( 0.0 , ESD - ETP2 )
          ETP3 = ETP*LSUBC
          
C ------------    FROZEN GROUND VERSION    --------------------------
C
          S = SNCOND * ( T1 - STC(1) ) / ( SNOWH )
C -------------------------------------------------------------------

C          S = 0.35 * ( T1 - STC(1) ) / ( 10.0 * ESD )

          S = MIN ( MAX ( S,-100.0 ), 100.0 )
          SEH = RCH*(T1-TH2)
          T14 = T1*T1
          T14 = T14*T14
          FLX3 = F - SIGMA*T14 - S - SEH - ETP3
          IF(FLX3.LE.0.0) FLX3=0.0
          EX = FLX3*0.001/LSUBF
          
C ------------    FROZEN GROUND VERSION    --------------------------
C   SNOWMELT REDUCTION DEPENDING ON SNOW COVER
C   IF SNOW COVER LESS THAN 5% NO SNOWMELT REDUCTION
C
          IF ( SNCOVER .GT. 0.05) EX = EX * SNCOVER
C ------------------------------------------------------------------- 
         
          SNMAX = EX * DT
        ENDIF
        
        IF(SNMAX.LT.ESD) THEN
          ESD = ESD - SNMAX
        ELSE
          EX = ESD/DT
          SNMAX = ESD
          ESD = 0.0
          FLX3 = EX*1000.0*LSUBF
C         CALL TDFCND ( DF1, SMC(1), B, F1 )
C         YYNUM = F - SIGMA * T24
C         YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-(FLX3+ETP3)) / RR
C          ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0
C          T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1
C         S = DF1 * ( STC(1) - T1 ) / ( 0.5 * ZSOIL(1) )
        ENDIF
        PRCP1 = PRCP1 + EX

      ENDIF
         
C SET THE EFFECTIVE POTENTIAL EVAPOTRANSP TO ZERO
      ETP1 = 0.0
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SMFLX RETURNS SOIL MOISTURE VALUES AND PRELIMINARY VALUES OF
C     EVAPOTRANSPIRATION.  IN THIS, THE SNOW PACK CASE, THE PRELIM-
C     VALUES (ETA1) ARE NOT USED IN SUBSEQUENT CALCULATION OF E.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF (ICE .NE. 1) THEN
         
C ------------    FROZEN GROUND VERSION    --------------------------
C (+) NEW STATES ADDED: SH2O,  AND FROZEN GROUND CORRECTION FACTOR
C
         CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL,
     +             SH2O, SLOPE, KDT, FRZFACT,
     &             SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,
     &             SMCREF,SHDFAC,CMCMAX,SMCDRY,CFACTR,RUNOFF1,RUNOFF2,
     &             RUNOFF3, EDIR1, EC1, ETT1,SFCTMP,Q2,NROOT,RTDIS)

      ENDIF

      ETA = BETA * ETP

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC/UPDATE THE GRND SFC (SKIN) MIXING RATIO.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C      Q1 = Q2 + ETA * CP / RCH

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     RETRIEVE THE SOIL TEMP DIFFUSIVITY/CONDUCTIVITY.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CALL TDFCND ( DF1, SMC(1), B, F1,QUARTZ,SMCMAX,SH2O )
      
C ------------    FROZEN GROUND VERSION    --------------------------
C   SOIL CONDUCTIVITY ADJUSTMENT
C      
      SICE1 = SMC(1) - SH2O(1)
C__________________________________________________________________
C This IF-statement is no longer used if Peters-Lidard thermal
C conductivity is used for all conditions (frozen,unfrozen):
C      IF ( SICE1 .GT. 0.0 ) CALL FRZCND ( DF1, SICE1 ) 
C -------------------------------------------------------------------
           
      IF(ICE.EQ.1) DF1=2.2

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     THE 'ADJUSTED TOP SOIL LYR TEMP' ( YY ) AND THE 'ADJUSTED SOIL
C     HEAT FLUX' ( ZZ1 ) ARE SET TO THE TOP SOIL LYR TEMP, AND 1,
C     RESPECTIVELY.  THESE ARE CLOSE-ENOUGH APPROXIMATIONS BECAUSE
C     THE SFC HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE
C     THE FLUX AT THE SNOW TOP SURFACE.  T11 IS A DUMMY ARGUEMENT
C     SINCE WE WILL NOT USE ITS VALUE AS REVISED BY SHFLX.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      ZZ1 = 1.0
      YY = STC(1)-0.5*S*ZSOIL(1)*ZZ1/DF1
      T11 = T1

C ------------    FROZEN GROUND VERSION    --------------------------
C   SNOW DEPTH & DENSITY ADJUSTMENT BASED ON SNOW COMPACTION
C      
      CALL SNOWPACK ( ESD, DT, SNOWH, SNDENS, T1, YY )
C -------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE THE SOIL HEAT FLUX
C     ( S1 ) AND THE SKIN TEMP ( T11 ) OUTPUT FROM THIS CALL ARE NOT
C     USED IN ANY SUBSEQUENT CALCULATIONS.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


C ------------    FROZEN GROUND VERSION    --------------------------
C (+) NEW STATES ADDED: SH2O & PARAMETER SMCWLT & PSISAT
C
      CALL SHFLX ( S1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL,TBOT,
     +             SMCWLT, PSISAT, SH2O,
     &             B,F1,DF1,ICE, 
     &             QUARTZ )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     IF PRECIP IS FALLING, CALC HEAT FLUX FROM SNOW SFC TO NEWLY
C     ACCUMULATING PRECIP.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FLX1 = 0.0
      IF ( SNOWNG ) THEN
        FLX1 = CPICE * PRCP * ( T1 - SFCTMP )
      ELSE
        IF ( PRCP .GT. 0.0 ) FLX1 = CPH2O * PRCP * ( T1 - SFCTMP )
      ENDIF
      RETURN
      END
      SUBROUTINE SNOWPACK ( W,DTS,HC,DS,TSNOW,TSOIL )

      IMPLICIT NONE

C **************************************************************
C **        SUBROUTINE TO CALCULATE SNOW COMPACTION          ***
C **  EQUATION OF INCREASING OF SNOW DENSITY WAS OBTAINED    ***
C     AS AN APPROXIMATE SOLUTION OF E. ANDERSON DIFFERENTIAL ***
C     EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19,         ***
C                 BY   VICTOR KOREN   03/25/95               ***
C **************************************************************

C **************************************************************
C  W      IS A WATER EQUIVALENT OF SNOW, IN M                ***
C  DTS    IS A TIME STEP, IN SEC                             ***
C  HC     IS A SNOW DEPTH, IN M                              ***
C  DS     IS A SNOW DENSITY, IN G/CM3                        ***
C  TSNOW  IS A SNOW SURFACE TEMPERATURE, K                   ***
C  TSOIL  IS A SOIL SURFACE TEMPERATURE, K                   ***
C      SUBROUTINE WILL RETURN NEW VALUES OF H AND DS         ***
C **************************************************************

      REAL C1, C2, HC, W, DTS, DS, TSNOW, TSOIL, H, WX
      REAL DT, TSNOWX, TSOILX, TAVG, B, DSX, DW

      INTEGER IPOL, J
      PARAMETER (C1=0.01, C2=21.0)
      REAL PEXP

C **  CONVERSION INTO SIMULATION UNITS   ************************* 

      H=HC*100.
      WX=W*100.
      DT=DTS/3600.
      TSNOWX=TSNOW-273.16
      TSOILX=TSOIL-273.16

C **  CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK              ***

      TAVG=0.5*(TSNOWX+TSOILX)                                    

C **  CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
C              DS=DS0*(EXP(B*W)-1.)/(B*W)
C              B=DT*C1*EXP(0.08*TAVG-C2*DS0)
C **  C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR)) 
C **  C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G

      IF(WX .GT. 1.E-2) THEN
         B=DT*C1*EXP(0.08*TAVG-C2*DS)
C
C  9/03/98   This line was changed to increase an accuracy of exponent
C        DSX=DS*((DEXP(B*WX)-1.)/(B*WX))
C--------------------------------------------------------------------
C  The above line was causing numerical difficulties because of the
C  denominator, that can go to zero (the function
C  has a well defined limit, equal to 1. when B*WX 
C  goes to zero, though)
C
C  Suggestion to substitute the line above by 
C                                       Pablo Grunmann 03-09-98
C--------------------------------------------------------------------
C Remember to declare PEXP

C Precision 
C   desired:  ipol greater than 9 only makes a difference on double
C             precision.
C        ipol=9, for rel.error <~ 1.6 e-6 % (8 significant digits)
C        ipol=8, for rel.error <~ 1.8 e-5 % (7 significant digits)
C        ipol=7, for rel.error <~ 1.8 e-4 % ...

        ipol = 9
	PEXP = 0.
        do j = ipol,1,-1
       	  PEXP = (1. + PEXP)*B*WX/real(j+1) 
        end do 
        PEXP = PEXP + 1.
C	
C  This result imitates the original formula
C  
        DSX=DS*(PEXP)
C--------------------------------------------------------------------
C  End of substitution
C--------------------------------------------------------------------
         IF(DSX .GT. 0.40) DSX=0.40
         IF(DSX .LT. 0.05) DSX=DS
         DS=DSX
         H=WX/DS

C **  CORRECTION OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER
C **  DURING SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED
C **  IN SNOW PER DAY DURING SNOWMELT TILL SNOW DENSITY 0.40

         IF((TSNOWX .GE. 0.) .AND. (H .NE. 0.)) THEN
            DW=0.13*DT/24.
            DS=DS*(1.-DW)+DW
            IF(DS .GT. 0.40) DS=0.40
            H=WX/DS
         ENDIF
                   
      ELSE
         IF(WX .LE. 0.) THEN
            WX=0.
            H=0.
            DS=0.
         ELSE
            H=WX/DS
         ENDIF
 
      ENDIF

C **  REPLACE PREVIOUS LIQUID WATER AND CHANGE DEPTH TO M

C    7 
      HC=H*0.01

      RETURN
      END
      SUBROUTINE SNOW_NEW ( T,P,HC,DS )
      
      IMPLICIT NONE
      
C     CALCULATING SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL
C     T - AIR TEMPERATURE, K
C     P - NEW SNOWFALL, M 
C     HC - SNOW DEPTH, M
C     DS - SNOW DENSITY 
C     NEW VALUES OF SNOW DEPTH & DENSITY WILL BE RETURNED   
      
      REAL HC
      REAL T 
      REAL P
      REAL DS
      REAL H
      REAL PX
      REAL TX
      REAL DS0
      REAL HNEW
      
C     **  CONVERSION INTO SIMULATION UNITS   ************************* 
      
      H=HC*100.
      PX=P*100.
      TX=T-273.16
      
C     ----------------------------------------------------------------------
C     ***  CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE      **
C     ***  EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
C     ***  AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE, 
C     ***  VEMADOLEN, SWEDEN, 1980, 172-177PP.
C-----------------------------------------------------------------------
      
      IF(TX .LE. -15.) THEN
         DS0=0.05
      ELSE                                                      
         DS0=0.05+0.0017*(TX+15.)**1.5
      ENDIF
      
C     **   ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL      
      
      HNEW=PX/DS0
      IF ((H+HNEW) .GT. 0.0) THEN
         DS=(H*DS+HNEW*DS0)/(H+HNEW)
         H=H+HNEW
      ELSE
         DS=0.0
         H=0.0
      END IF
      
      HC=H*0.01
      
      RETURN
      END
      SUBROUTINE SRT (RHSTT,RUNOFF,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,
     &                 ZSOIL,DWSAT,DKSAT,SMCMAX,B, RUNOFF1, 
     +                 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE)      

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY
CC    =======   TERM OF THE SOIL WATER DIFFUSION EQUATION.  ALSO TO
CC              COMPUTE ( PREPARE ) THE MATRIX COEFFICIENTS FOR THE
CC              TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOLD
      PARAMETER ( NSOLD = 4 )
      INTEGER K
      INTEGER NSOIL
      INTEGER CVFRZ
      REAL AI     ( NSOLD )
      REAL B
      REAL BI     ( NSOLD )
      REAL CI     ( NSOLD )
      REAL DMAX   ( NSOLD )
      REAL DDZ
      REAL DDZ2
      REAL DENOM
      REAL DENOM2
      REAL DKSAT
      REAL DSMDZ
      REAL DSMDZ2
      REAL DWSAT
      REAL EDIR
      REAL ET     ( NSOIL )
      REAL INFMAX
      REAL KDT
      REAL MXSMC
      REAL MXSMC2
      REAL NUMER
      REAL PCPDRP
      REAL PDDUM
      REAL RHSTT  ( NSOIL )
      REAL RUNOFF
      
      REAL SH2O   ( NSOIL )
      REAL SH2OA  ( NSOIL )
      REAL SICE   ( NSOIL )
      
      REAL SMCMAX
      REAL WCND
      REAL WCND2
      REAL WDF
      REAL WDF2
      REAL ZSOIL  ( NSOIL )

      REAL RUNOFF1, RUNOFF2, DT, SMCWLT, SLOPE, FRZX, DT1
      REAL SMCAV, DICE, DD, VAL, DDT, PX, FCR, ACRT, SUM
      REAL SSTT, SLOPX

      INTEGER IOHINF, KS, IALP1, J, JJ

CMIC$ TASKCOMMON ABCI
      COMMON /ABCI/ AI, BI, CI

C -----------     FROZEN GROUND VERSION    -------------------------
C   REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
C   AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
C   CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. 
C   BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT
C   CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM.
C   THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6})  
C
        PARAMETER ( CVFRZ = 3 )
C ------------------------------------------------------------------
     
C	print*,'in SRT, Declaration -----------------------'
C	print*,'NSOIL=' , NSOIL
C	print*,'NSOLD=' , NSOLD
        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C ************  INCLUDE THE INFILTRATION FORMULE FROM SCHAAKE AND KOREN 
C           MODEL
C
CC    MODIFIED BY Q DUAN
CC    REPLACE THE MAXIMUM INFILTRATION EQUATION IN OSU MODEL
CC    BY SCHAAKE/KOREN EXPRESSION
CC
CC    INITIALIZE PARAMETERS AND VARIABLES
CC
      
      IOHINF=1

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PDDUM = PCPDRP
      RUNOFF1 = 0.0
      IF ( PCPDRP .NE. 0.0 ) THEN

CC++  MODIFIED BY Q. DUAN, 5/16/94

C        IF (IOHINF .EQ. 1) THEN
  
          DT1 = DT/86400.
          SMCAV = SMCMAX - SMCWLT
          DMAX(1)=-ZSOIL(1)*SMCAV

C -----------     FROZEN GROUND VERSION    ------------------------
C
          DICE = -ZSOIL(1) * SICE(1)
C-------------------------------------------------------------------
          
          DMAX(1)=DMAX(1)*(1.0 - (SH2OA(1)+SICE(1)-SMCWLT)/SMCAV)
          DD=DMAX(1)
      DO KS=2,NSOIL
          
C -----------     FROZEN GROUND VERSION    ------------------------
C
           DICE = DICE + ( ZSOIL(KS-1) - ZSOIL(KS) ) * SICE(KS)
C------------------------------------------------------------------- 
         
           DMAX(KS)=(ZSOIL(KS-1)-ZSOIL(KS))*SMCAV
           DMAX(KS)=DMAX(KS)*(1.0 - (SH2OA(KS)+SICE(KS)-SMCWLT)/SMCAV)
           DD=DD+DMAX(KS)
      END DO
CC          VAL = (1.-EXP(-KDT*SQRT(DT1)))
C REMOVE THE SQRT BY F. CHEN, 07/08/96
          VAL = (1.-EXP(-KDT*DT1))
          DDT = DD*VAL
          PX = PCPDRP*DT  
          IF(PX.LT.0.0) PX = 0.0
          INFMAX = (PX*(DDT/(PX+DDT)))/DT
          
C -----------     FROZEN GROUND VERSION    --------------------------
C    REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
C
         FCR = 1. 
         IF ( DICE .GT. 1.E-2) THEN 
           ACRT = CVFRZ * FRZX / DICE 
           SUM = 1.
           IALP1 = CVFRZ - 1 
           DO J = 1,IALP1
              K = 1
              DO JJ = J+1, IALP1
                K = K * JJ
              END DO   
              SUM = SUM + (ACRT ** ( CVFRZ-J)) / FLOAT (K) 
           END DO 
           FCR = 1. - EXP(-ACRT) * SUM 
         END IF 
         INFMAX = INFMAX * FCR
C -------------------------------------------------------------------

C ************    CORRECTION OF INFILTRATION LIMITATION    **********
C     IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE 
C     VALUE OF HYDROLIC CONDUCTIVITY
C
C         MXSMC = AMAX1 ( SH2OA(1), SH2OA(2) ) 
        MXSMC = SH2OA(1)

C	print*,'SRT, BEFORE WDFCND - 1 ------------------------------'
C	print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
 
         CALL WDFCND ( WDF,WCND,MXSMC,SMCMAX,B,DKSAT,DWSAT )
            INFMAX = AMAX1(INFMAX, WCND)
            INFMAX= AMIN1(INFMAX,PX)

C	print*,'SRT, AFTER WDFCND - 1 ------------------------------'
C	print*,'WDF,WCND=' , WDF,WCND
C	print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
 
C  ******************************************************************
C          INFMAX = MIN(INFMAX, DKSAT)
 
C        ELSE

C          INFMAX = -DWSAT*(SMCMAX-SMCA(1))/(0.5*ZSOIL(1))+DKSAT
C          IF (SMCA(1).GE.SMCMAX) INFMAX = DKSAT
C
C        END IF
C
C*************** RUNOFF1: SURFACE RUNOFF *******
C
          IF ( PCPDRP .GT. INFMAX ) THEN
            RUNOFF1 = PCPDRP - INFMAX
            PDDUM = INFMAX
          END IF

      END IF
C
C*********   BELOW IS THE ORIGINAL OSU MODEL
C
C      PDDUM = PCPDRP
C      PRINT*,' ************** IN SRT *******************'
C      PRINT*,' PCPDRP=', PCPDRP
C      RUNOFF1 = 0.0
C      IF ( PCPDRP .NE. 0.0 ) THEN
C        INFMAX = - DWSAT * (SMCMAX - SMCA(1))/(0.5 * ZSOIL(1)) + DKSAT
C        IF ( SMCA(1) .GE. SMCMAX ) INFMAX = DKSAT
C        IF ( PCPDRP .GT. INFMAX ) THEN
C          RUNOFF1 = PCPDRP - INFMAX
C
C*************** RUNOFF1: SURFACE RUNOFF *******
C
C          PDDUM = INFMAX
C        END IF
C      END IF


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     DETERMINE SOIL MOISTURE DIFFUSIVITY (WDF) & CONDUCTIVITY (WCND)
C     COEFS FOR TOP SOIL LYR USING THE HIGHER OF SMC FOR THE TOP
C     SOIL LAYER OR SMC FOR THE NEXT LOWER SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C_________________________________________________________________
C Comments on August 98 changes: 
C
C This approach disregarded the actual soil moisture in 
C the upper layer if it was less that that of the next, for 
C hydraulic conductivity and diffussivity purposes (supposedly
C to solve problems with the numerical discretization on a very 
C coarse vertical resolution setting - two layer model).
C In our experience, This would actually allow the upper 
C layer to go very dry at a fast rate while the next layer
C keeps  unchanged (except for diffusivity effect)
C
C                                      Pablo J. Grunmann
C 
C      MXSMC = MAX( SH2OA(1), SH2OA(2) )
C___________________________________________________________
C Test a less drastic version of MXSMC: !!!!!!!!!!!!!!!!!!!!
C      MXSMC = MAX( SH2OA(1), (SH2OA(2)+SH2OA(1))/2.0 )
C___________________________________________________________
C_________________________________________________________________
C New approach: dif/conductivity parameters on upper layer
C depend on upper layer only
        MXSMC =  SH2OA(1)

C	print*,'SRT, BEFORE WDFCND - 2'
C	print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
 

      CALL WDFCND ( WDF,WCND,MXSMC,SMCMAX,B,DKSAT,DWSAT )


C	print*,'SRT, AFTER WDFCND - 2'
C	print*,'WDF,WCND=' , WDF,WCND
C	print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DDZ = 1. / ( -.5 * ZSOIL(2) )
      AI(1) = 0.0
      BI(1) = WDF * DDZ / ( -ZSOIL(1) )
      CI(1) = -BI(1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL
C     MOISTURE GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DSMDZ = ( SH2O(1) - SH2O(2) ) / ( -.5 * ZSOIL(2) )
      RHSTT(1) = (WDF * DSMDZ + WCND - PDDUM + EDIR + ET(1))/ZSOIL(1)
      SSTT = WDF * DSMDZ + WCND + EDIR + ET(1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     INITIALIZE DDZ2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DDZ2 = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 10 K = 2 , NSOIL
        DENOM2 = ( ZSOIL(K-1) - ZSOIL(K) )
        IF ( K .NE. NSOIL ) THEN

        SLOPX = 1.
        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       RETRIEVE THE SOIL WATER DIFFUSIVITY FOR THE WETTER OF THE
C       CURRENT OR THE NEXT LOWER SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C          MXSMC2 = MAX ( SMCA(K), SMCA(K+1) )

CC           MXSMC2 = MAX ( SH2OA(K), SH2OA(K+1) )
 
C ************ PILPS TEST *************
          MXSMC2 =  SH2OA(K)

C	print*,'SRT, BEFORE WDFCND - 3'
C	print*,'MXSMC2,SMCMAX=' , MXSMC2,SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
C	print*,'K=' , K
 

          CALL WDFCND ( WDF2,WCND2,MXSMC2,SMCMAX,B,DKSAT,DWSAT )


C	print*,'SRT, AFTER WDFCND - 3'
C	print*,'WDF2,WCND2=' , WDF2,WCND2
C	print*,'MXSMC2,SMCMAX=' , MXSMC2,SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DENOM = ( ZSOIL(K-1) - ZSOIL(K+1) )
          DSMDZ2 = ( SH2O(K) - SH2O(K+1) ) / ( DENOM * 0.5 )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DDZ2 = 2.0 / DENOM
          CI(K) = -WDF2 * DDZ2 / DENOM2
        ELSE

C   SLOPE OF BOTTOM LAYER IS INTRODUCED     ************
C
          SLOPX = SLOPE
C--------------------------------------------------------
          
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC
C         CONDUCTIVITY FOR THIS LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


C	print*,'SRT, BEFORE WDFCND - 4'
C	print*,'SH2OA(NSOIL),SMCMAX=' , SH2OA(NSOIL),SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
C	print*,'K=' , K
 
          CALL WDFCND( WDF2,WCND2,SH2OA(NSOIL),SMCMAX,B,DKSAT,DWSAT )

C	print*,'SRT, AFTER WDFCND - 4'
C	print*,'WDF2,WCND2=' , WDF2,WCND2
C	print*,'SH2OA(NSOIL),SMCMAX=' , SH2OA(NSOIL),SMCMAX
C	print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DSMDZ2 = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         SET MATRIX COEF CI TO ZERO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          CI(K) = 0.0
        END IF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2 - (WDF * DSMDZ) 
     +          - WCND + ET(K)
        RHSTT(K) = NUMER / (-DENOM2)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        AI(K) = -WDF * DDZ / DENOM2
        BI(K) = -( AI(K) + CI(K) )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        IF(K.EQ.NSOIL) THEN
C*************** RUNOFF2: GROUND WATER RUNOFF ***********
           RUNOFF2 = SLOPX * WCND2
        ENDIF

        IF ( K .NE. NSOIL ) THEN
          WDF = WDF2
          WCND = WCND2
          DSMDZ = DSMDZ2
          DDZ = DDZ2
        END IF
 10   CONTINUE

C	print*,'SRT, final Runoff'
C	print*,'RUNOFF1=' , RUNOFF1
C	print*,'RUNOFF2=' , RUNOFF2
 
      RETURN
      END
      SUBROUTINE SSTEP ( SH2OOUT, SH2OIN, CMC, RHSTT, RHSCT, DT,
     &              NSOIL, SMCMAX, CMCMAX, RUNOFF3, ZSOIL,SMC,SICE )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE/UPDATE THE SOIL MOISTURE CONTENT VALUES
CC    =======   AND THE CANOPY MOISTURE CONTENT VALUES.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOLD, I
      PARAMETER ( NSOLD = 4 )
      INTEGER K
      INTEGER NSOIL
      REAL AI     ( NSOLD )
      REAL BI     ( NSOLD )
      REAL CI     ( NSOLD )
      REAL CMC
      REAL CMCMAX
      REAL DT
      REAL RHSCT
      REAL RHSTT  ( NSOIL )

C --------     FROZEN GROUND CORRECTIONS     ------------------      
      REAL SH2OIN  ( NSOIL )
      REAL SH2OOUT ( NSOIL )
      REAL SICE    ( NSOIL )
      REAL SMC     ( NSOIL )
C --------------------------------------------------------------
      
      REAL SMCMAX
      REAL ZSOIL(NSOIL)

      REAL RUNOFF3, RUNOFS, WPLUS, DDZ, STOT, KK11, WFREE, DPLUS


cmic$ TASKCOMMON ABCI
      COMMON /ABCI/ AI, BI, CI

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
C     TRI-DIAGONAL MATRIX ROUTINE.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 10 K = 1 , NSOIL
        RHSTT(K) = RHSTT(K) * DT
        AI(K) = AI(K) * DT
        BI(K) = 1. + BI(K) * DT
        CI(K) = CI(K) * DT
 10   CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CALL ROSR12 ( CI, AI, BI, CI, RHSTT, RHSTT, NSOIL )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
C     NEW VALUE.  MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C   ****************** RUNOFF3: Runoff within soil layers *******

C ***********   WATER BALANCE CALCULATIONS WERE WRONG    ********
C   IT WAS CORRECTED DURING DURING INCLUDING FROZEN GROUND   ****
C   
      RUNOFS = 0.0
      WPLUS = 0.0
      RUNOFF3 = 0.
      DDZ = - ZSOIL(1)
      
      DO 20 K = 1 , NSOIL
        IF ( K .NE. 1 ) DDZ = ZSOIL(K - 1) - ZSOIL(K)
        SH2OOUT(K) = SH2OIN(K) + CI(K) + WPLUS / DDZ
	
C	print*,'IN sstep'
C	print*,'SH2OOUT=', SH2OOUT
	
        STOT = SH2OOUT(K) + SICE(K)
        IF ( STOT .GT. SMCMAX ) THEN
           IF ( K .EQ. 1 ) THEN
              DDZ = -ZSOIL(1)
           ELSE
              KK11 = K - 1
              DDZ = -ZSOIL(K) + ZSOIL(KK11)
           END IF
           WPLUS = ( STOT - SMCMAX ) * DDZ
         ELSE
           WPLUS = 0.
         END IF
         SMC(K) = AMAX1 ( AMIN1( STOT, SMCMAX ), 0.02 )

         SH2OOUT(K) = AMAX1 ( (SMC(K) - SICE(K)), 0.0 )

   20 CONTINUE

C  ***  V. KOREN   9/01/98    ******
C     WATER BALANCE CHECKING UPWARD
C
      IF(WPLUS .GT. 0.) THEN
       DO I=NSOIL-1,1,-1
	IF(I .EQ. 1) THEN
	 DDZ=-ZSOIL(1)
        ELSE
	 DDZ=-ZSOIL(I)+ZSOIL(I-1)
        ENDIF
	WFREE=(SMCMAX-SH2OOUT(I)-SICE(I))*DDZ
	DPLUS=WFREE-WPLUS
	IF(DPLUS .GE. 0.) THEN
	 SH2OOUT(I)=SH2OOUT(I)+WPLUS/DDZ
	 SMC(I)=SH2OOUT(I)+SICE(I)
	 WPLUS=0.
	   
        ELSE
	 SH2OOUT(I)=SH2OOUT(I)+WFREE/DDZ
	 SMC(I)=SH2OOUT(I)+SICE(I)
	 WPLUS=-DPLUS
        ENDIF
       END DO
30     RUNOFF3=WPLUS
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CONVERT RHSCT TO AN 'AMOUNT' VALUE AND ADD TO PREVIOUS
C     CMC VALUE TO GET NEW CMC VALUE.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CMC = CMC + DT * RHSCT
C      CMC = MIN ( MAX ( CMC, 0. ), CMCMAX )
C
C code changed 8-26-98 by Manikin (on original source)
C
      IF (CMC .LT. 1.E-20) CMC=0.0
      CMC = MIN(CMC,CMCMAX)


      RETURN
      END
      FUNCTION TBND (TU, TB, ZSOIL, ZBOT, K, NSOIL)            

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC   PURPOSE:   CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER
CC   =======    BY INTERPOLATION OF THE MIDDLE LAYER TEMPERATURES
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL TBND, T0, TU, TB, ZBOT, ZUP, ZB

      INTEGER NSOIL, K
      
      PARAMETER (T0=273.16)
      REAL ZSOIL (4)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC   USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      IF(K .EQ. 1) THEN
        ZUP=0.
      ELSE
        ZUP=ZSOIL(K-1)
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC   USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
CC   TEMPERATURE INTO THE LAST LAYER BOUNDARY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      IF(K .EQ. NSOIL) THEN
        ZB=2.*ZBOT-ZSOIL(K)
      ELSE
        ZB=ZSOIL(K+1)
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC   LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      TBND=TU+(TB-TU)*(ZUP-ZSOIL(K))/(ZUP-ZB)
      
      RETURN
      END  
      SUBROUTINE TDFCND ( DF, SMC, B, F1, Q,  SMCMAX, SH2O)

      IMPLICIT NONE

C               INCLUDED NEW VARS HERE: ,QUARTZ,SMCMAX,SH2O
C______________________________________________________________
C   NEW CALLING ARGUMENTS in SFLX, HRT   10/98
C  (remove CALL TDFCND (DF, SMC, BB(ISOIL), F11(ISOIL))
C
C	   CALL TDFCND (DF, SMC, BB(ISOIL), F11(ISOIL),
C     & 		 QUARTZ(ISOIL), SMCMAX(ISOIL), SH2O)
C--------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF
CC    =======   THE SOIL FOR A GIVEN POINT AND TIME.
CC
CC    VERSION:  PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
CC    =======
CC                          Added by    PABLO J. GRUNMANN, 08/1998
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       REAL B
       REAL DF
       REAL F1
       REAL GAMMD
       REAL THKDRY
       REAL AKE
       REAL THKICE
       REAL THKO
       REAL THKQTZ
       REAL THKSAT
       REAL THKS
       REAL THKW
       REAL Q
       REAL SATRATIO
       REAL SH2O
       REAL SMC
       REAL SMCMAX
       REAL XU
       REAL XUNFROZ


C IN PRMSOI:
C        DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52, 
C     &              0.35, 0.60, 0.40, 0.82/

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
C     OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
C
C                                  PABLO GRUNMANN, 08/17/98
C REFS.:
C      FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK 
C              AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
C      JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
C              UNIVERSITY OF TRONDHEIM,
C      PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL 
C              CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
C              AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
C              VOL. 55, PP. 1209-1224.
C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C THKW  WATER THERMAL CONDUCTIVITY
C   THKQTZ  THERMAL CONDUCTIVITY FOR QUARTZ
C THKO  THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
C THKS  THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ + OTHER)
C   THKICE  ICE THERMAL CONDUCTIVITY
C   SMCMAX  POROSITY (= SMCMAX)
C        Q  QUARTZ CONTENT(SOIL TYPE)
C

C  NEEDS PARAMETERS
C POROSITY(SOIL TYPE):
C      POROS = SMCMAX
C SATURATION RATIO:
      SATRATIO = SMC/SMCMAX
C      print *, 'SATRATIO=',SATRATIO
C     PARAMETERS  W/(M.K)
      THKICE = 2.2
      THKW = 0.57
      THKO = 2.0
C      IF (Q .LE. 0.2) THKO = 3.0
      THKQTZ = 7.7
C  SOLIDS' CONDUCTIVITY      
      THKS = (THKQTZ**Q)*(THKO**(1.- Q))
C      print *, 'THKS = ',THKS
C  UNFROZEN FRACTION (FROM 1.0, I.E., 100%LIQUID, TO 0.0 (100% FROZEN))
      XUNFROZ=(SH2O + 1.E-9)/(SMC + 1.E-9)
C      print *, '   '
C      print *, 'XUNFROZ = ',XUNFROZ
C      print *, '    '
C  UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
      XU=XUNFROZ*SMCMAX 
C  SATURATED THERMAL CONDUCTIVITY
      THKSAT = THKS**(1.-SMCMAX)*THKICE**(SMCMAX-XU)*THKW**(XU)
C      print *, 'THKSAT = ',THKSAT
C  DRY DENSITY IN KG/M3
      GAMMD = (1. - SMCMAX)*2700.
C      print *, 'GAMMD = ',GAMMD
C  DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
      THKDRY = (0.135*GAMMD + 64.7)/(2700. - 0.947*GAMMD)
C      print *, 'THKDRY = ',THKDRY
C RANGE OF VALIDITY FOR THE KERSTEN NUMBER
      IF ( SATRATIO .GT. 0.1 ) THEN

C    KERSTEN NUMBER (FINE FORMULA, AT LEAST 5% OF PARTICLES<(2.E-6)M)
   	IF ( (XUNFROZ + 0.0005) .LT. SMC ) THEN
C    FROZEN
   	   AKE = SATRATIO
   	ELSE
C    UNFROZEN
   	   AKE = LOG10(SATRATIO) + 1.0
   	ENDIF

      ELSE
        
C USE K = KDRY
	AKE = 0.0
C        print *, 'AKE (ELSE) = ',AKE
      ENDIF
C  THERMAL CONDUCTIVITY
       DF = AKE*(THKSAT - THKDRY) + THKDRY

      RETURN
      END
      SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT,
     &      CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE TRANSPIRATION FROM THE VEGTYP FOR THIS PT.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER I, K
      INTEGER NSOIL
      INTEGER NROOT

      REAL CFACTR
      REAL CMC
      REAL CMCMAX
      REAL ET ( NSOIL )
      REAL ETP1
      REAL ETP1A
      REAL GX (7)
C      REAL PART (4)
      REAL PC
      REAL RTDIS (NSOIL)
      REAL SHDFAC
      REAL SMC    ( NSOIL )
      REAL SMCREF
      REAL SMCWLT
      REAL ZSOIL  ( NSOIL )

      REAL SFCTMP, Q2, SGX, DENOM, RTX

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       INITIALIZE  PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 20 K = 1, NSOIL
         ET(K) = 0.
 20   CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC AN 'ADJUSTED' POTNTL TRANSPIRATION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      ETP1A = SHDFAC * PC * ETP1 * (1.0 - (CMC /CMCMAX) ** CFACTR)
                  
      SGX = 0.0
      DO I = 1, NROOT
         GX(I) = ( SMC(I) - SMCWLT ) / ( SMCREF - SMCWLT )
         GX(I) = MAX ( MIN ( GX(I), 1. ), 0. )
         SGX = SGX + GX (I)
      END DO
      SGX = SGX / NROOT
      
      DENOM = 0.
      DO I = 1,NROOT
         RTX = RTDIS(I) + GX(I) - SGX
         GX(I) = GX(I) * MAX ( RTX, 0. )
         DENOM = DENOM + GX(I)
      END DO   
      IF ( DENOM .LE. 0.0) DENOM = 1.
      
      DO I = 1, NROOT
         ET(I) = ETP1A * GX(I) / DENOM
      END DO  
      
C     ************    BELOW IS OLD CODE    ************
C     ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
C     C        ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
      
C***  USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
C     C        ET(1) = RTDIS(1) * ETP1A
C     C       ET(1) =  ETP1A*PART(1)
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
C     BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
C     ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
C     C        DO 10 K = 2, NROOT
      
C     GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
C     GX = MAX ( MIN ( GX, 1. ), 0. )
C     TEST CANOPY RESISTANCE
C     GX = 1.0
      
C     ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
C     C         ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
      
C***  USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
C     C         ET(K) = RTDIS(K) * ETP1A
C     C         ET(K) = ETP1A*PART(K)
C     C  10    CONTINUE
      
      RETURN
      END
      SUBROUTINE WDFCND ( WDF,WCND,SMC,SMCMAX,B,DKSAT,DWSAT )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE SOIL WATER DIFFUSIVITY AND SOIL
CC    =======   HYDRAULIC CONDUCTIVITY.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL B
      REAL DKSAT
      REAL DWSAT
      REAL EXPON
      REAL FACTR2
      REAL SMC
      REAL SMCMAX
      REAL WCND
      REAL WDF

C	print*,'------------ in WDFCND -------------------------------'
C	print*,'BEFORE WDFCND'
C	print*,'B=',B
C	print*,'DKSAT=',DKSAT
C	print*,'DWSAT=',DWSAT
C	print*,'EXPON=',EXPON
C	print*,'FACTR2=',FACTR2
C	print*,'SMC=',SMC
C	print*,'SMCMAX=',SMCMAX
C	print*,'WCND=',WCND
C	print*,'WDF=',WDF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SMC = SMC
      SMCMAX = SMCMAX
      FACTR2 = SMC / SMCMAX

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      EXPON = B + 2.0
      WDF = DWSAT * FACTR2 ** EXPON
C  Drainage Tests !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C      WDF = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      EXPON = ( 2.0 * B ) + 3.0
      WCND = DKSAT * FACTR2 ** EXPON

C	print*,' WDFCND Results --------------------------------'
C	print*,'B=',B
C	print*,'DKSAT=',DKSAT
C	print*,'DWSAT=',DWSAT
C	print*,'EXPON=',EXPON
C	print*,'FACTR2=',FACTR2
C	print*,'SMC=',SMC
C	print*,'SMCMAX=',SMCMAX
C	print*,'WCND=',WCND
C	print*,'WDF=',WDF
C      PRINT*,' SMC         WDF           WCND             B'
C      PRINT*,SMC,WDF,WCND,B

      RETURN
      END
