text:       SUBROUTINE HRT1 ( RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,TBOT,
       +                 ZBOT, PSISAT, SH2O, DT, B, SH, SR, SF, 
  C
  C12/02  V.K.  ADDED NEW FORMAL PARAMETER - IVERS
  CC     +                 STYPE, QUARTZ, SMCMR )
       +                 STYPE, QUARTZ, SMCMR, IVERS, CKSOIL)
       
  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.
  CC     TSS - SNOW-BARE SOIL WEGHTED TEMPERATURE DURING SNOWMELT
  CC     SOILT - SOIL SURFACE TEMPERATURE
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
        INTEGER NSOIL, K
        REAL CAIR, CH2O, CSOIL, DDZ, DDZ2, DENOM, DFMAX1, DFMAXN,
       +     DTSDZ, DTSDZ2, HCPCT, MXICE, S, SMCMAX, STYPE, 
       +     TBOT, ZBOT, YY, SH, SR, SF, SOILT, TSS
        REAL SMC   ( NSOIL )
        REAL SH2O  ( NSOIL )
        REAL STC   ( NSOIL )
        REAL ZSOIL ( NSOIL )
        REAL RHSTS ( NSOIL )
        
        PARAMETER ( NSZ_ABC = 10 )
        REAL AI    ( NSZ_ABC )
        REAL BI    ( NSZ_ABC )
        REAL CI    ( NSZ_ABC )
        COMMON /ABC/ AI, BI, CI
        
  c  test only
        common /test1/ iday
              
        PARAMETER( CAIR=1004.0,CH2O=4.2 E6,CSOIL=1.26E6 )
        PARAMETER ( CICE = 2.106 E6, T0 = 273.16 )       
  
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  C     RETRIEVE THE THERMAL DIFFUSIVITY 
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  
  C  ACCOUNT FOR RESIDUE EFFECT ON CONDUCTIVITY
        SICE = 0.
        DO I = 1, NSOIL
          SICE = SICE + SMC(I) - SH2O(I)
        ENDDO
  
  c3/10 surfwater
        smc0=smcmr*0.15
        sh2o0=smc0
        MXICE = SMC0 - SH2O0
        dfmax1 = cnd_johns(stype,sh2o0,mxice,smcmr,0.)
  c3/10 surfwater end
  
  C  CALCULATE SNOW/SOIL FLUX AND SOIL SURFACE TEMPERATURE IF SNOW
        SOILT=YY
        TSS=-99. 
        IF(SH .NE. 0.) THEN
         TSS=YY
         IF(YY .GT. T0) TSS=SF*T0+(1-SF)*YY
         DFSNOW=CSNOW(SR)
  
  C  CALCULATION OF SOIL TEMPERATURE IS BASED ON Qsnow=Qsoil
         ALP=1./(1. + DFMAX1*SH/(-DFSNOW*0.5*ZSOIL(1)))
         SOILT=ALP*TSS+(1-ALP)*STC(1)
        ENDIF 
        S = DFMAX1 * ( STC(1) - SOILT ) / ( 0.5 * ZSOIL(1) )
        S = MIN ( MAX ( S,-100.0 ), 100.0 )
  
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  C     CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  
        HCPCT = SH2O0*CH2O + (1.0-SMCMR)*CSOIL + (SMCMR-SMC0)*CAIR
       +        + ( SMC0 - SH2O0 ) * CICE
  
  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)
  
  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) )
        RHSTS(1) = ( DFMAX1 * DTSDZ - S ) / ( ZSOIL(1) * HCPCT )
  
  C   CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
  C   AND SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR.
        
        TBK = TBND ( STC(1), STC(2), ZSOIL, ZBOT, 1, NSOIL )
        IF ( SICE .NE. 0. .OR. SOILT .LT. T0 ) THEN 
         TSNSR = SNKSRC ( SOILT, STC(1),TBK, SMC(1), SH2O(1), HCPCT, 
  C
  C12/02  V.K.  ADDED NEW FORMAL PARAMETER - IVERS
  CC     +             ZSOIL, SMCMR, PSISAT, B, DT, 1, STYPE, 0.)
  c 3/10 surface water
       +     ZSOIL, SMCMR, PSISAT, B, DT, 1, STYPE, 0.,IVERS,CKSOIL,smc0)
       
         RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
        ENDIF 
        DDZ2 = 0.0
  
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  C     LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  
        DO 40 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
       +        + ( SMC(K) - SH2O(K) ) * CICE
  
          IF ( K .NE. NSOIL ) THEN
  
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  C         RETRIEVE THE THERMAL DIFFUSIVITY
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  
          MXICE = SMC(K) - SH2O(K)
          dfmaxn=cnd_johns(stype,sh2o(k),mxice,smcmax,quartz)
  
  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   CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
            TBK1 = TBND ( STC(K), STC(K+1), ZSOIL, ZBOT, K, NSOIL )
            
          ELSE
  
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  C         RETRIEVE THE THERMAL DIFFUSIVITY FOR THE LOWEST SOIL LAYER
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  
  C   SOIL THERMAL DIFFUSIVITY ADJUSTMENT
  C   AND  CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
  C
            MXICE = SMC(K) - SH2O(K)
            dfmaxn=cnd_johns(stype,sh2o(k),mxice,smcmax,quartz)
            TBK1 = TBND ( STC(K), TBOT, ZSOIL, ZBOT, K, NSOIL )
  
  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   SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR.
        IF ( SICE .NE. 0. .OR. TBK .LT. T0 ) THEN
         TSNSR = SNKSRC ( TBK, STC(K),TBK1, SMC(K), SH2O(K), HCPCT, 
  C
  C12/02  V.K.  ADDED NEW FORMAL PARAMETER - IVERS
  CC     +             ZSOIL,SMCMAX, PSISAT, B, DT, K,STYPE,QUARTZ)
  c 3/10 surface water
       +   ZSOIL,SMCMAX, PSISAT, B, DT, K,STYPE,QUARTZ,IVERS,CKSOIL,smc0)
       
         RHSTS(K) = RHSTS(K) - TSNSR / DENOM
        ENDIF 
        
  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
          TBK = TBK1
          DFMAX1 = DFMAXN
          DTSDZ = DTSDZ2
          DDZ = DDZ2
     40   CONTINUE
  
        RETURN
        END
  	  
  	  
        SUBROUTINE HSTEP ( STCOUT, STCIN, RHSTS, DT, NSOIL )
  
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  CC    PURPOSE:  TO CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        
        INTEGER NSOIL, K
        REAL DT
        REAL RHSTS ( NSOIL )
        REAL STCOUT   ( NSOIL )
        REAL STCIN  ( NSOIL )
  
        PARAMETER ( NSZ_ABC = 10 )
        REAL AI    ( NSZ_ABC )
        REAL BI    ( NSZ_ABC )
        REAL CI    ( NSZ_ABC )
        COMMON /ABC/ 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
          STCOUT(K) = STCIN(K) + CI(K)
    20  CONTINUE
        RETURN
        END
  	  
  	  
        SUBROUTINE ROSR12 ( P, A, B, C, D, DELTA, NSOIL )
  
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  CC    PURPOSE:  TO INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN
  CC    =======   BELOW:
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  
        INTEGER K, KK, NSOIL
        REAL P(NSOIL)
        REAL A(NSOIL)
        REAL B(NSOIL)
        REAL C(NSOIL)
        REAL D(NSOIL)
        REAL DELTA(NSOIL)
  
  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
     
        
  