cable_sli_solve.F90 Source File


This file depends on

sourcefile~~cable_sli_solve.f90~~EfferentGraph sourcefile~cable_sli_solve.f90 cable_sli_solve.F90 sourcefile~cable_define_types.f90 cable_define_types.F90 sourcefile~cable_sli_solve.f90->sourcefile~cable_define_types.f90 sourcefile~cable_iovars.f90 cable_iovars.F90 sourcefile~cable_sli_solve.f90->sourcefile~cable_iovars.f90 sourcefile~cable_sli_numbers.f90 cable_sli_numbers.F90 sourcefile~cable_sli_solve.f90->sourcefile~cable_sli_numbers.f90 sourcefile~cable_sli_utils.f90 cable_sli_utils.F90 sourcefile~cable_sli_solve.f90->sourcefile~cable_sli_utils.f90 sourcefile~cable_climate_type_mod.f90 cable_climate_type_mod.F90 sourcefile~cable_define_types.f90->sourcefile~cable_climate_type_mod.f90 sourcefile~cable_iovars.f90->sourcefile~cable_define_types.f90 sourcefile~cable_sli_numbers.f90->sourcefile~cable_define_types.f90 sourcefile~cable_sli_utils.f90->sourcefile~cable_define_types.f90 sourcefile~cable_sli_utils.f90->sourcefile~cable_sli_numbers.f90 sourcefile~cable_common.f90 cable_common.F90 sourcefile~cable_climate_type_mod.f90->sourcefile~cable_common.f90 sourcefile~grid_constants_cbl.f90 grid_constants_cbl.F90 sourcefile~cable_climate_type_mod.f90->sourcefile~grid_constants_cbl.f90 sourcefile~cable_runtime_opts_mod.f90 cable_runtime_opts_mod.F90 sourcefile~cable_common.f90->sourcefile~cable_runtime_opts_mod.f90

Files dependent on this one

sourcefile~~cable_sli_solve.f90~~AfferentGraph sourcefile~cable_sli_solve.f90 cable_sli_solve.F90 sourcefile~cable_sli_main.f90 cable_sli_main.F90 sourcefile~cable_sli_main.f90->sourcefile~cable_sli_solve.f90 sourcefile~cbl_model_driver_offline.f90 cbl_model_driver_offline.F90 sourcefile~cbl_model_driver_offline.f90->sourcefile~cable_sli_main.f90 sourcefile~cable_mpimaster.f90 cable_mpimaster.F90 sourcefile~cable_mpimaster.f90->sourcefile~cbl_model_driver_offline.f90 sourcefile~cable_mpiworker.f90 cable_mpiworker.F90 sourcefile~cable_mpiworker.f90->sourcefile~cbl_model_driver_offline.f90 sourcefile~cable_serial.f90 cable_serial.F90 sourcefile~cable_serial.f90->sourcefile~cbl_model_driver_offline.f90 sourcefile~cable_offline_driver.f90 cable_offline_driver.F90 sourcefile~cable_offline_driver.f90->sourcefile~cable_serial.f90

Source Code

! Original code by P.J. Ross 2001-2007  See: Ross, P.J. (2003)
!     Modeling soil water and solute transport - fast, simplified numerical solutions. Agron. J. 95:1352-1361.
!
! Modified by V. Haverd 2008: matrix expanded by factor of 2 to allow solution of coupled heat/moisture fluxes
! Explicit calculation of heat/moisture fluxes at surface
! Switchable option for litter
! Isotope subroutine (V. Haverd and M. Cuntz, September 2008)
!
! Frozen Soil (V. Haverd, M. Cuntz Sep 2010 - Jan 2011)
! Pond lumped with top soil layer (V. Haverd Jan 2011)
! Snow included in layers -2:0 (V. Haverd Jan 2011)
! Convergence test for dthetaldT (V.Haverd Feb 2011)
! Include heat advection by liquid water flux (V.Haverd Feb 2011)

! Code cleaned and isotope routine updated (V. Haverd, M. Cuntz Aug-Sep 2012)
! lots of changes, especially relating to ponding and freezing together
! Outstanding problems with surface runoff and advection (set h0max large and turn off advection for now)
!
! Sep 17 2012 (V.Haverd)
! remove separate solution for maxpond (h0 > h0max)
! Instead, convert excess pond to runoff at t=tfin, and correct energy stores for loss of pond
! Advection and surface runoff now functioning
!
! Dec 31 2012 (V.Haverd)
! bug fixes to improve energy conservation associated with change in ice,
! especially correction at time of pond disappearance (removal of negative pond)
!
! Jan 5 2013 (V.Haverd)
! Remove iteration over frozen soil within "(while iok==0)" loop
! For frozen soil, at the time of updating variables, evaluate new T and ice content,
! consistent with J0 + deltaJ, where deltaJ is the change in energy obtained from the matrix solution (=LHS_h*dt)
!
! Jan 7-8 2013 (V.Haverd)
! Calls to soil_snow albedo and density
! Improve iteration convergence at time of updating new T and ice content in frozen soil
! Bug fix relating to moisture conservation in snow
!
! Jan 14-15 2013 (V. Haverd)
! Modified definition of Sliq for frozen soil (in hyofS): now defined relative to (theta_sat - theta_r - theta_ice)
! Allow for snow-pack initalisation in the absence of pond
!
! Feb 7 2013 (V.Haverd)
! Revised formulation of dphidT for frozen soil (see hyofS): fixes a few remaining -ve T spikes in frozen soil
!
! Feb 9 2013 (V.Haverd)
! Moved snow layer adjustments to subroutine
!
! Feb 24 2013 (V. Haverd)
! Snow structure expanded to accommodate 3 layers
!
! March 30 (V. Haverd)
!
! Jan 7 2014 (V. Haverd)
! Adjust snow surface bcs such that snow surface temperature can't be more than 0
!
! Jan 8 2014 (V Haverd)
! enable 2nd snow layer so that if snowpack exceeds 3 cm , the excess goes into the layer below.
!
! March 2014 (V. Haverd)
! Extract SEB calcs to subroutine in utils, and trial Force-Restore method
!
! October 2016 (S Wales)
! Restructured code of solve into several subroutines

MODULE sli_solve

  USE cable_def_types_mod, ONLY: r_2, i_d
  USE sli_numbers,         ONLY: &
       experiment, &
       zero, one, two, half, thousand, e3, e5, &  ! numbers
       Tzero, rlambda, lambdaf, lambdas, Dva, rhocp, rhow, gf, hmin, & ! parameters
       csice , cswat, &
       snmin, nsnow_max, fsnowliq_max, &
       params, vars_aquifer, vars_met, vars, vars_snow, solve_type, & ! types
       dSfac, h0min, Smax, h0max, dSmax, dSmaxr, dtmax, dSmax, dSmaxr, & ! numerical limits
       dtmax, dtmin, dTsoilmax, dTLmax, nsteps_ice_max, nsteps_max,tol_dthetaldT, &
       hbot, botbc, &
       Mw, Rgas, & ! boundary condition
       freezefac
  USE sli_utils,           ONLY: &
       x, Sofh, hyofh, hyofS, litter_props, massman_sparse, tri, &
       aquifer_props, Tfrz, thetalmax, Tthetalmax, dthetalmaxdTh, &
       getfluxes_vp, getheatfluxes, flux, sol, &
       csat, slope_csat, potential_evap, tri, setsol, zerovars, &
       esat_ice, slope_esat_ice, Tfrozen, rtbis_Tfrozen, GTFrozen, &
       JSoilLayer, esat, forcerestore, SEB
  USE cable_IO_vars_module, ONLY: logn

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: solve ! solution routine

  INTEGER(i_d), DIMENSION(:), ALLOCATABLE :: nless, n_noconverge ! global counters

  ! Definitions of public entities and private parameters (see above for default
  ! values):
  ! botbc    - bottom boundary condn for water; "constant head", "free drainage",
  !    "seepage", or "zero flux". Constant head means that matric head h
  !    is specified. Free drainage means zero gradient of matric head,
  !    i.e. unit hydraulic gradient. Seepage means zero flux when the
  !    matric head is below zero and an upper limit of zero for the head.
  ! h0max    - max pond depth allowed before runoff.
  ! hbot    - matric head at bottom of profile when botbc set to "constant head".
  ! dSmax    - max change in S (the "effective saturation") of any unsaturated
  !    layer to aim for each time step; controls time step size.
  ! dSmaxr   - maximum negative relative change in S each time step. This
  !    parameter helps avoid very small or negative S.
  ! dtmax    - max time step allowed.
  ! dsmmax   - max solute change per time step (see dSmax); user should set this
  !    according to solute units used. Units for different solutes can be
  !    scaled by the user (e.g. to an expected max of around 1.0).
  ! dSfac    - a change in S of up to dSfac*dSmax is accepted.
  ! dpmaxr   - relative change in matric flux potential (MFP) phi that is
  !    accepted for convergence when finding head h at soil interfaces.
  ! h0min    - min (negative) value for surface pond when it empties.
  ! Smax    - max value for layer saturation to allow some overshoot.
  ! solve    - sub to call to solve RE
  !

CONTAINS

  SUBROUTINE get_fluxes_and_derivs( &
       irec, mp, qprec, qprec_snow, n, dx, h0, &
       Tsoil, S, &
       qh, vmet, vlit, vsnow, var, T0, Tsurface, &
       dxL, SL, Tl, &
       plit, par, &
       qali, &
       qvh, &
       qevap, &
       again, getq0,getqn,init, ns, nsat, &
       nsatlast, &
       phip, qpme, hint, phimin, &
       qexd, &
       q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb, qadv, qadvya, qadvyb, &
       qadvTa, qadvTb, vtmp, qliq, qv, qvT, qlya, qlyb, &
       qvya, qvyb, qlTb, qvTa, qvTb, &
       vtop, vbot, &
       lE0, G0, &
       Epot, surface_case, &
       iflux, litter, j, kk, &
       advection, dTqwdTa, dTqwdTb, Tqw, keff, &
       hice &
       )
    IMPLICIT NONE
    INTEGER(i_d)                                           :: irec, mp
    REAL(r_2),      DIMENSION(1:mp)                        :: qprec
    REAL(r_2),      DIMENSION(1:mp)                        :: qprec_snow
    INTEGER(i_d)                                           :: n
    REAL(r_2),      DIMENSION(1:n) :: dx
    REAL(r_2),      DIMENSION(1:mp)                        :: h0
    REAL(r_2),      DIMENSION(1:n) :: Tsoil, S
    REAL(r_2),      DIMENSION(-nsnow_max:n) :: qh
    TYPE(vars_met), DIMENSION(1:mp)                        :: vmet
    TYPE(vars),     DIMENSION(1:mp)                        :: vlit
    TYPE(vars_snow), DIMENSION(1:mp)                       :: vsnow
    TYPE(vars),     DIMENSION(1:n) :: var
    REAL(r_2),      DIMENSION(1:mp)                        :: T0, Tsurface
    REAL(r_2),      DIMENSION(1:mp)                        :: dxL
    REAL(r_2),      DIMENSION(1:mp)                        :: SL, Tl
    TYPE(params),   DIMENSION(1:mp)                        :: plit
    TYPE(params),   DIMENSION(1:n) :: par
    REAL(r_2),      DIMENSION(1:mp), OPTIONAL                        :: qali
    REAL(r_2),      DIMENSION(-nsnow_max:n) :: qvh
    REAL(r_2),    DIMENSION(1:mp)                          :: qevap
    LOGICAL,      DIMENSION(1:mp)                          :: again, getq0,getqn,init
    INTEGER(i_d), DIMENSION(1:mp)                          :: ns, nsat, nsatlast
    REAL(r_2),    DIMENSION(1:mp)                          :: phip
    REAL(r_2),    DIMENSION(1:mp)                          :: qpme
    REAL(r_2),    DIMENSION(1:n) :: hint, phimin, qexd
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qadv, qadvya, qadvyb, qadvTa, qadvTb
    TYPE(vars)                                             :: vtmp
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qliq, qv, qvT, qlya, qlyb, qvya, qvyb, qlTb, qvTa, qvTb
    TYPE(vars),         DIMENSION(1:mp)                    :: vtop, vbot
    REAL(r_2),          DIMENSION(1:mp)                    :: lE0, G0, Epot
    INTEGER(i_d),       DIMENSION(1:mp)                    :: surface_case
    INTEGER(i_d),       DIMENSION(1:mp)                    :: iflux
    LOGICAL                                                :: litter
    INTEGER(i_d)                                           :: j, kk
    INTEGER(i_d)                                           :: advection ! switches
    REAL(r_2)                                              :: dTqwdTa, dTqwdTb, Tqw, keff
    REAL(r_2),          DIMENSION(1:mp)                    :: hice


    !----- get fluxes and derivs
    ! get surface condition
    ! ns==1 if no pond or full pond, i.e. do not solve for change in pond height
    ! ns==0 then change for pond height
    IF ((var(1)%phi <= phip(kk) .AND. h0(kk) <= zero .AND. nsat(kk) < n) .OR. &
         (var(1)%isat==0)) THEN ! no ponding
       ns(kk)    = 1 ! start index for eqns
    ELSE ! ponding
       ns(kk)    = 0
       var(1)%phi = (one+e5)*var(1)%phie+(h0(kk)-hice(kk))*var(1)%Ksat
       TL(kk)    = vmet(kk)%Ta ! initialise pond T
       vtop(kk)%isat    = 1
       vtop(kk)%h       = h0(kk)
       vtop(kk)%phi     = MAX((var(1)%phie -var(1)%he*var(1)%Ksat), &
            (one+e5)*var(1)%phie)+(h0(kk)-hice(kk))*var(1)%Ksat
       vtop(kk)%K       = var(1)%Ksat
       vtop(kk)%rh      = one
       vtop(kk)%lambdav = rlambda
       vtop(kk)%lambdaf = lambdaf
       ! var(1)%phi = var(1)%phie
       ! calculates phi1,eff (pond + top soil layer)   !vh! does this get used???
       CALL flux(par(1), vtop(kk), var(1), half*dx(1), &
            q(0), qya(0), qyb(0), qTa(0), qTb(0))


    ENDIF

    IF (Sl(kk) >= one) THEN
       vlit(kk)%isat = 1
       vlit(kk)%h    = plit(kk)%he
    ENDIF

    ! get bottom boundary condn
    IF (botbc=="seepage") THEN
       vtmp = zerovars()
       vbot = SPREAD(vtmp,1,mp)
       IF (var(n)%h > -half*gf*dx(n)) THEN
          getqn(kk) = .TRUE.
          vbot(kk)%isat    = 1
          vbot(kk)%phi     = (zero-var(n)%he)*var(n)%Ksat+var(n)%phie ! (special for frozen soil)
          vbot(kk)%K       = var(n)%Ksat
          vbot(kk)%rh      = one
          vbot(kk)%lambdav = rlambda
          vbot(kk)%lambdaf = lambdaf
       ELSE
          getqn(kk) = .FALSE.
       ENDIF
    END IF

    ! Fluxes and derivatives at the air/surface interface
    IF (vsnow(kk)%nsnow > 0) THEN
       surface_case(kk) = 2
    ELSE
       surface_case(kk) = 1
    ENDIF
    SELECT CASE (surface_case(kk))
    CASE (1) ! no snow
       CALL SEB(n, par(:), vmet(kk), vsnow(kk), var(:), qprec(kk), qprec_snow(kk), dx(:), &
            h0(kk), Tsoil(:), &
            Tsurface(kk), G0(kk), lE0(kk),Epot(kk),  &
            q(0), qevap(kk), qliq(0), qv(0), &
            qyb(0), qTb(0), qlyb(0), qvyb(0), qlTb(0), qvTb(0), qh(0), &
            qadv(0), qhyb(0), qhTb(0), qadvyb(0), qadvTb(0), irec)
       qya(0)    = zero
       qTa(0)    = zero
       qlya(0)   = zero
       qvya(0)   = zero
       qvTa(0)   = zero
       qhya(0)   = zero
       qhTa(0)   = zero
       qadvya(0) = zero
       qadvTa(0) = zero

    CASE (2) ! snow
       CALL SEB(n, par(:), vmet(kk), vsnow(kk), var(:), qprec(kk), qprec_snow(kk), dx(:), &
            h0(kk), Tsoil(:), &
            Tsurface(kk), G0(kk), lE0(kk), Epot(kk),  &
            q(-vsnow(kk)%nsnow), qevap(kk), qliq(-vsnow(kk)%nsnow), qv(-vsnow(kk)%nsnow), &
            qyb(-vsnow(kk)%nsnow), qTb(-vsnow(kk)%nsnow), qlyb(-vsnow(kk)%nsnow), &
            qvyb(-vsnow(kk)%nsnow), qlTb(-vsnow(kk)%nsnow), qvTb(-vsnow(kk)%nsnow), &
            qh(-vsnow(kk)%nsnow), qadv(-vsnow(kk)%nsnow), qhyb(-vsnow(kk)%nsnow), &
            qhTb(-vsnow(kk)%nsnow), qadvyb(-vsnow(kk)%nsnow), qadvTb(-vsnow(kk)%nsnow), irec)
       qya(-vsnow(kk)%nsnow)    = zero
       qTa(-vsnow(kk)%nsnow)    = zero
       qlya(-vsnow(kk)%nsnow)   = zero
       qvya(-vsnow(kk)%nsnow)   = zero
       qvTa(-vsnow(kk)%nsnow)   = zero
       qhya(-vsnow(kk)%nsnow)   = zero
       qhTa(-vsnow(kk)%nsnow)   = zero
       qadvya(-vsnow(kk)%nsnow) = zero
       qadvTa(-vsnow(kk)%nsnow) = zero

    CASE default
       WRITE(*,*) "solve: illegal surface case."
       STOP 2
    END SELECT ! surface_case
    qpme(kk) = q(0) ! water flux into top of soil column
    ! finished all the surfaces

    ! get moisture fluxes and derivatives (at time t=0, i.e. q0 etc.)

    CALL getfluxes_vp(n, dx(1:n), vtop(kk), vbot(kk), par(1:n), var(1:n), & ! moisture fluxes
         hint(1:n), phimin(1:n), q(0:n), qya(0:n), qyb(0:n), qTa(0:n), qTb(0:n), &
         qliq(0:n), qlya(0:n), qlyb(0:n), qv(0:n), qvT(0:n), qvh(0:n), qvya(0:n), &
         qvyb(0:n), &
         iflux(kk), init(kk), getq0(kk), getqn(kk), Tsoil(1:n), T0(kk), nsat(kk), nsatlast(kk))
    qTa(n) = zero
    qTb(n) = zero
    qvTa(1:n) = qTa(1:n)
    qvTb(1:n) = qTb(1:n)
    qlTb(1:n) = zero


    ! get  fluxes heat and derivatives (at time t=0, i.e. q0 etc.)
    CALL getheatfluxes(n, dx(1:n), dxL(kk), &
         qh(0:n), qhya(0:n), qhyb(0:n), qhTa(0:n), qhTb(0:n), &
         var(1:n), vlit(kk), Tsoil(1:n), TL(kk), litter, &
         q(0:n), qya(0:n), qyb(0:n), qTa(0:n), qTb(0:n), &
         qadv(0:n),qadvya(0:n), qadvyb(0:n), qadvTa(0:n), qadvTb(0:n), &
         advection) ! heat fluxes

    ! get heat and vapour fluxes and derivatives in snow-pack
    IF (vsnow(kk)%nsnow>1) THEN
       ! vapour flux at interfaces between snow layers
       DO j=1, vsnow(kk)%nsnow-1
          keff = 2_r_2*((vsnow(kk)%kE(j)/(thousand*lambdaf))*(vsnow(kk)%kE(j+1)/(thousand*lambdaf))/ &
               ((vsnow(kk)%kE(j)/(thousand*lambdaf))*vsnow(kk)%depth(j+1)+(vsnow(kk)%kE(j+1)/ &
               (thousand*lambdaf))*vsnow(kk)%depth(j)) )
          q(j-vsnow(kk)%nsnow) = keff*(vsnow(kk)%tsn(j)-vsnow(kk)%tsn(j+1))
          qTa(j-vsnow(kk)%nsnow) = MERGE(keff,zero,vsnow(kk)%hliq(j)<=zero)
          qTb(j-vsnow(kk)%nsnow) = MERGE(-keff,zero,vsnow(kk)%hliq(j+1)<=zero)
          qya(j-vsnow(kk)%nsnow) = zero
          qyb(j-vsnow(kk)%nsnow) = zero
          ! conductive heat flux at interface between snow layers
          keff = 2_r_2*(vsnow(kk)%kth(j+1)*vsnow(kk)%kth(j))/ &
               (vsnow(kk)%kth(j+1)*vsnow(kk)%depth(j)+vsnow(kk)%kth(j)*vsnow(kk)%depth(j+1))  ! check this!
          qh(j-vsnow(kk)%nsnow) = keff*(vsnow(kk)%tsn(j)-vsnow(kk)%tsn(j+1))
          qhTa(j-vsnow(kk)%nsnow) = MERGE(zero,keff,vsnow(kk)%hliq(j)>zero)
          qhTb(j-vsnow(kk)%nsnow) = MERGE(zero,-keff,vsnow(kk)%hliq(j+1)>zero)

          ! advective heat flux at interface between snow layers
          Tqw  = MERGE(vsnow(kk)%tsn(j), vsnow(kk)%tsn(j+1), q(j-vsnow(kk)%nsnow)>zero)
          dTqwdTb = MERGE(zero,one, q(j-vsnow(kk)%nsnow)>zero)
          dTqwdTa = MERGE(one,zero, q(j-vsnow(kk)%nsnow)>zero)
          IF (vsnow(kk)%hliq(j)>zero) THEN
             qadv(j-vsnow(kk)%nsnow) = rhow*q(j-vsnow(kk)%nsnow)*cswat*Tqw
             qadvTa(j-vsnow(kk)%nsnow) = zero
          ELSE
             qadv(j-vsnow(kk)%nsnow) = rhow*q(j-vsnow(kk)%nsnow)*cswat*Tqw
             qadvTa(j-vsnow(kk)%nsnow) = rhow*cswat*q(j-vsnow(kk)%nsnow)*dTqwdTa  + &
                  rhow*cswat*Tqw*qTa(j-vsnow(kk)%nsnow)
          ENDIF
          qadvTb(0) = rhow*cswat*q(j-vsnow(kk)%nsnow)*dTqwdTb + rhow*cswat*Tqw*qTb(j-vsnow(kk)%nsnow)

          qh(j-vsnow(kk)%nsnow) = qh(j-vsnow(kk)%nsnow) + qadv(j-vsnow(kk)%nsnow)
          qhTa(j-vsnow(kk)%nsnow) = qhTa(j-vsnow(kk)%nsnow) +  qadvTa(j-vsnow(kk)%nsnow)
          qhTb(j-vsnow(kk)%nsnow) = qhTb(j-vsnow(kk)%nsnow) +  qadvTb(j-vsnow(kk)%nsnow)
       ENDDO
    ENDIF ! end fluxes at snow/snow interfaces

    IF (vsnow(kk)%nsnow.GE.1) THEN
       ! vapour flux at soil/snow interface
       IF (var(1)%isat==1) THEN
          keff = vsnow(kk)%kE(vsnow(kk)%nsnow)/(thousand*lambdaf)/vsnow(kk)%depth(vsnow(kk)%nsnow)/2_r_2
       ENDIF
       IF (var(1)%isat==0) THEN
          keff = 2_r_2*((vsnow(kk)%kE(vsnow(kk)%nsnow)/(thousand*lambdaf))*(var(1)%kE/thousand/var(1)%lambdav))/ &
               ((vsnow(kk)%kE(vsnow(kk)%nsnow)/(thousand*lambdaf))*dx(1)+(var(1)%kE/thousand/var(1)%lambdav)* &
               vsnow(kk)%depth(vsnow(kk)%nsnow))
       ENDIF
       q(0) = keff*(vsnow(kk)%tsn(vsnow(kk)%nsnow)-Tsoil(1))
       qTa(0) = keff
       qTb(0) =-keff

       qya(0) = zero
       qyb(0) = zero
       qv(0)   = q(0)
       qvyb(0) = qyb(0)
       qvTb(0) = qTb(0)
       qvya(0) = qya(0)
       qvTa(0) = qTa(0)
       qliq(0) = zero
       qlyb(0) = zero
       qlTb(0) = zero
       qlya(0) = zero
       qya(0) = zero;

       IF (vsnow(kk)%hliq(vsnow(kk)%nsnow)>zero) THEN
          qhTa(0) = zero
          qTa(0) = zero
       ENDIF

       ! conductive heat flux at snow/soil interface  ! check this!
       keff = 2._r_2*(vsnow(kk)%kth(vsnow(kk)%nsnow)*var(1)%kth)/ &
            (vsnow(kk)%kth(vsnow(kk)%nsnow)*dx(1)+var(1)%kth*vsnow(kk)%depth(vsnow(kk)%nsnow))
       qh(0) = keff*(vsnow(kk)%tsn(vsnow(kk)%nsnow)-Tsoil(1))
       IF (vsnow(kk)%hliq(1)>zero) THEN
          qhTa(0) = zero
       ELSE
          qhTa(0) = keff
       ENDIF
       qhTb(0) = -keff

       ! advective heat flux at snow/soil interface
       Tqw  = MERGE(vsnow(kk)%tsn(vsnow(kk)%nsnow), Tsoil(1), q(0)>zero)
       dTqwdTb = MERGE(zero,one, q(0)>zero)
       dTqwdTa = MERGE(one,zero, q(0)>zero)
       IF (vsnow(kk)%hliq(vsnow(kk)%nsnow)>zero) THEN
          qadv(0) = rhow*q(0)*cswat*Tqw
          qadvTa(0) = zero
       ELSE
          qadv(0) = rhow*q(0)*cswat*Tqw
          qadvTa(0) = rhow*cswat*q(0)*dTqwdTa  +  rhow*cswat*Tqw*qTa(0)
       ENDIF
       qadvTb(0) = rhow*cswat*q(0)*dTqwdTb + rhow*cswat*Tqw*qTb(0)
    ENDIF ! end of heat and vapour fluxes at soil/snow interface

    IF (ns(kk)==0) THEN ! pond included in top soil layer
       ! change qya(1) from dq/dphi (returned by getfluxes) to dq/dh
       qya(1) = var(1)%Ksat*qya(1)
       IF (advection==1) THEN
          qhya(1) = qhya(1) - qadvya(1)
          Tqw  = MERGE(Tsoil(1), Tsoil(2), q(1)>zero)
          qadvya(1) =  rhow*cswat*qya(1)*Tqw  ! apply corrected qya(1) to qadvya(1)
          qhya(1) = qhya(1) + qadvya(1)
       ENDIF
    ENDIF

    ! adjust for bottom boundary condition
    IF (botbc=="zero flux") THEN
       qliq(n) = zero
       qv(n)   = zero
       q(n)    = zero
       qya(n)  = zero
       qlya(n) = zero
       qvya(n) = zero
    ENDIF

    ! specify mositure flux at bottom of soil column (heat flux set to zero)
    IF (botbc /= "constant head") THEN
       SELECT CASE (botbc)
       CASE ("zero flux")
          q(n)   = zero
          qya(n) = zero
       CASE ("free drainage")
          q(n) = gf*var(n)%K
          IF (var(n)%isat == 0) THEN
             qya(n) = gf*var(n)%KS
          ELSE
             qya(n) = zero
          END IF
       CASE ("seepage")
          IF (var(n)%h <= -half*gf*dx(n)) THEN
             q(n)   = zero
             qya(n) = zero
          END IF
       CASE default
          WRITE(*,*) "solve: illegal bottom boundary condition."
          STOP 2
       END SELECT
    END IF
    IF (PRESENT(qali)) THEN
       IF (qali(kk)>zero) THEN
          q(n)   = -qali(kk)
          qya(n) = zero
       END IF
    ENDIF
    IF (experiment==7 .OR. experiment==8) THEN
       q(n)   = q(0)
       qya(n) = zero
    ENDIF
    IF (experiment==8) qh(n) = G0(kk)

    ! adjust lower heat flux for advection
    IF (advection==1) THEN
       qadv(n) = rhow*cswat*(Tsoil(n))*q(n)
       qadvya(n) = rhow*cswat*(Tsoil(n))*qya(n)
       qadvTa(n)= rhow*cswat*q(n)
       qh(n) = qh(n) + qadv(n)
       qhya(n) = qhya(n) + qadvya(n)
       qhTa(n) = qhTa(n) + qadvTa(n)
    ELSE
       qadv(:)   = zero
       qadvya(:) = zero
       qadvyb(:) = zero
       qadvTa(:) = zero
       qadvTb(:) = zero
    ENDIF

    qexd(1:n) = zero ! time derivative for root extraction (assumed fixed at input value)
    again(kk)  = .FALSE. ! flag for recalcn of fluxes (default=false)
    !----- end get fluxes and derivs
  END SUBROUTINE get_fluxes_and_derivs

  SUBROUTINE estimate_timestep( &
       irec, tfin, mp, n, dx, h0, &
       qh, qhTa, qadvTa, qhTb, qadvTb, &
       nsteps, vlit, vsnow, var, &
       dxL, &
       plit, par, &
       deltaTa, &
       qL, qhL, &
       again, ns, nsat, &
       nsatlast, nsteps0, dmax, dt, &
       qpme, t, &
       q, qadv, &
       tmp2d1, tmp2d2, &
       tmp1d1, tmp1d3, &
       iflux, litter, kk, &
       advection, &
       iqex &
       )
    IMPLICIT NONE
    REAL(r_2)                                              :: tfin
    INTEGER(i_d)                                           :: mp, irec
    INTEGER(i_d)                                           :: n
    REAL(r_2),      DIMENSION(1:n) :: dx
    REAL(r_2),      DIMENSION(1:mp)                        :: h0
    REAL(r_2),      DIMENSION(-nsnow_max:n) :: qh
    INTEGER(i_d),   DIMENSION(1:mp)                        :: nsteps
    TYPE(vars),     DIMENSION(1:mp)                        :: vlit
    TYPE(vars_snow), DIMENSION(1:mp)                       :: vsnow
    TYPE(vars),     DIMENSION(1:n) :: var
    REAL(r_2),      DIMENSION(1:mp)                        :: dxL
    TYPE(params),   DIMENSION(1:mp)                        :: plit
    TYPE(params),   DIMENSION(1:n) :: par
    REAL(r_2),      DIMENSION(1:mp)                        :: deltaTa
    REAL(r_2),    DIMENSION(1:mp)                          :: qL, qhL
    LOGICAL,      DIMENSION(1:mp)                          :: again
    INTEGER(i_d), DIMENSION(1:mp)                          :: ns, nsat, nsatlast, nsteps0
    REAL(r_2),    DIMENSION(1:mp)                          :: dmax, dt
    REAL(r_2),    DIMENSION(1:mp)                          :: qpme, t
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: q, qhTa, qhTb
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qadv, qadvTa, qadvTb
    REAL(r_2),    DIMENSION(0:n) :: tmp2d1, tmp2d2,  deltaTmax
    REAL(r_2),          DIMENSION(1:mp)                    :: tmp1d1, tmp1d3
    INTEGER(i_d),       DIMENSION(1:mp)                    :: iflux
    LOGICAL                                                :: litter
    INTEGER(i_d)                                           :: kk
    INTEGER(i_d)                                           :: advection ! switches
    REAL(r_2),          DIMENSION(1:n) :: iqex


    !----- first estimate of time step dt before the calculation
    !      gets revised after the calculation
    dmax(kk)  = zero
    tmp2d1(:) = zero ! 1st order estim of rate of change of moisture storage
    tmp2d2(:) = zero ! 1st order estim of rate of change of temperature
    ! estimate rate of change of moisture storage [m/s]
    WHERE (var(1:n)%isat==0.AND.var(1:n)%iice==0.) tmp2d1(1:n) = &
         ABS(q(1:n)-q(0:n-1)-iqex(1:n))/(par(1:n)%thre*dx(1:n))
    WHERE (var(1:n)%iice==1)  tmp2d1(1:n) =  tmp2d1(1:n)/2._r_2
    ! estimate rate of change of temperature [K/s]
    tmp2d2(1:n) = ABS(qh(1:n)-qh(0:n-1))/(var(1:n)%csoileff*dx(1:n))
    IF (advection==1) THEN
       ! first order estimate of rate of temperature change
       tmp2d2(1:n) = ABS((qh(1:n)-qadv(1:n))-(qh(0:n-1)-qadv(0:n-1)))/(var(1:n)%csoileff*dx(1:n))
       deltaTmax(1:n) = tmp2d2(1:n)*(tfin-t(kk)) ! maximum  T change
    ENDIF
    IF (litter .AND. ns(kk)==1 ) THEN ! litter , no pond
       IF (vlit(kk)%isat==0) THEN ! estimate rate of change of moisture storage [m/s]
          WRITE(*,*) 'Should not be here - QL 01 ', qL(kk)
          tmp2d1(0) = ABS(q(0) - qL(kk))/(plit(kk)%thre*dxL(kk))
       ENDIF
       WRITE(*,*) 'Should not be here - QHL 01 ', qhL(kk)
       tmp2d2(0)  = ABS(qh(0) - qhL(kk))/(vlit(kk)%csoileff*dxL(kk)) ! estimate rate of change of heat storage [K/s]
       tmp1d3(kk) = (dTLmax-ABS(deltaTa(kk))) / tmp2d2(0)
    ELSE
       tmp2d1(0)  = zero
       tmp2d2(0)  = zero
       tmp1d3(kk) = dtmax
    ENDIF

    dmax(kk)   = MAXVAL(tmp2d1(1:n),1) ! max derivative |dS/dt|
    tmp1d1(kk) = MAXVAL(tmp2d2(1:n),1) ! max derivative |dTsoil/dt|
    IF (dmax(kk) > zero) THEN
       dt(kk) = MIN(dSmax/dmax(kk), dTLmax/tmp1d1(kk), tmp1d3(kk)) ! constrained either by moisture or temp
    ELSE ! steady state flow
       IF (qpme(kk)>=q(n)) THEN ! if saturated soil columnn and more precip then drainige -> finish
          dt(kk) = tfin-t(kk) ! step to finish
       ELSE ! otherwise adjust dt because change of pond height
          dt(kk) = -(h0(kk)-half*h0min)/(qpme(kk)-q(n))
       END IF
       dt(kk) = MIN(dt(kk), dTLmax/tmp1d1(kk), tmp1d3(kk)) ! constrained by  temp
    END IF

    ! check that time step is short enough to prevent melting a top snow layer
    IF (vsnow(kk)%nsnow.GT.0 ) THEN
       ! energy required to melt ice in snow-pack
       tmp1d1(kk) = -rhow*(vsnow(kk)%hsnow(1)-vsnow(kk)%hliq(1))*(csice*vsnow(kk)%tsn(1) - lambdaf)
       IF (qh(-vsnow(kk)%nsnow)*dt(kk).GT.tmp1d1(kk)) THEN
          dt(kk) = 0.9_r_2*tmp1d1(kk)/qh(-vsnow(kk)%nsnow)
       ENDIF
    ENDIF

    IF (dt(kk)>dtmax) dt(kk) = dtmax ! user's limit

    ! if initial step, improve phi where S>=1
    ! might be that you get better derivatives, especially at sat/non-sat interfaces
    IF (nsteps(kk)==nsteps0(kk) .AND. nsat(kk)>0 .AND. iflux(kk)==1) THEN
       again(kk) = .TRUE.
       dt(kk)    = dtmin
    END IF
    ! if fully saturated but was not fully saturated before, adjust even within each time step iteration
    IF (nsat(kk)==n .AND. nsatlast(kk)<n .AND. iflux(kk)==1) THEN
       again(kk) = .TRUE. ! profile has just become saturated so adjust phi values
       dt(kk)    = dtmin
    END IF
    ! sprint to the end
    IF (t(kk)+1.1_r_2*dt(kk)>tfin) THEN ! step to finish
       dt(kk) = tfin-t(kk)
       t(kk)  = tfin
    ELSE
       t(kk) = t(kk)+dt(kk) ! tentative update
       IF (again(kk)) t(kk) = t(kk)-dt(kk)
    END IF

    !----- end estimate time step dt


  END SUBROUTINE estimate_timestep

  SUBROUTINE iflux_loop( &
       tfin, irec, mp, qprec, qprec_snow, n, dx, h0, S, thetai, &
       Jsensible, Tsoil, evap, infil, drainage, discharge, &
       qh, nsteps, vmet, vlit, vsnow, var, T0, Tsurface, Hcum, lEcum, &
       Gcum, Qadvcum, Jcol_sensible, &
       Jcol_latent_S, Jcol_latent_T, csoil, kth, phi, dxL, zdelta, SL, Tl, &
       plit, par, wex, ciso_snow, cisoice_snow, &
       qali, &
       qvsig, qlsig, qvTsig, qvh, deltaTa, &
       precip, qevap, qL, qhL, qybL, qTbL, qhTbL, qhybL, rexcol, wcol, &
       again, getq0,getqn,init, again_ice, ih0, iok, itmp, ns, nsat, &
       nsatlast, nsteps0, accel, dmax, dt, dwinfil, dwoff, fac, &
       phip, qpme, rsig, rsigdt, sig, t, hint, phimin, &
       qexd, aa, bb, cc, dd, ee, ff, gg, dy, aah, bbh, cch, ddh, eeh, ffh, ggh, &
       de, q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb, qadv, qadvya, qadvyb, &
       qadvTa, qadvTb, vtmp, qsig, qhsig, qadvsig, qliq, qv, qvT, qlya, qlyb, &
       qvya, qvyb, qlTb, qvTa, qvTb, deltaS, dTsoil, tmp2d1, tmp2d2, &
       cv0, &
       cisoliqice_snow, &
       dthetaldT, thetal, isave, nsteps_ice, imelt, vtop, vbot, v_aquifer, &
       dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge, &
       dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S, &
       deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T, qevapsig, &
       qrunoff, tmp1d1, tmp1d2, tmp1d3,  tmp1d4, deltah0, deltaSL, &
       lE0, G0, &
       Epot, Tfreezing, dtdT, LHS, RHS, LHS_h, RHS_h, surface_case, nns, &
       iflux, litter, i, j, k, kk, condition, littercase, isotopologue, &
       advection, c2, theta, dTqwdTa, dTqwdTb, Tqw, keff, cp, &
       cpeff, hice, h0_0, hice_0, h0_tmp, hice_tmp, qmelt, &
       qtransfer, delta_snowcol, delta_snowT, &
       delta_snowliq, thetai_0, J0, tmp1, tmp2, iqex, &
       nfac1, nfac2, nfac3, nfac4, nfac5, nfac6, nfac7, nfac8, nfac9, &
       nfac10, nfac11, nfac12, J0snow, wcol0snow, h_ex, wpi, err &
       )
    IMPLICIT NONE
    REAL(r_2)                                              :: tfin
    INTEGER(i_d)                                           :: irec, mp
    REAL(r_2),      DIMENSION(1:mp)                        :: qprec
    REAL(r_2),      DIMENSION(1:mp)                        :: qprec_snow
    INTEGER(i_d)                                           :: n
    REAL(r_2),      DIMENSION(1:n) :: dx
    REAL(r_2),      DIMENSION(1:mp)                        :: h0
    REAL(r_2),      DIMENSION(1:n) :: S
    REAL(r_2),      DIMENSION(1:n) :: thetai
    REAL(r_2),      DIMENSION(1:n) :: Jsensible
    REAL(r_2),      DIMENSION(1:n) :: Tsoil
    REAL(r_2),      DIMENSION(1:mp)                        :: evap
    REAL(r_2),      DIMENSION(1:mp)                        :: infil
    REAL(r_2),      DIMENSION(1:mp)                        :: drainage, discharge
    REAL(r_2),      DIMENSION(1:mp,-nsnow_max:n)           :: qh
    INTEGER(i_d),   DIMENSION(1:mp)                        :: nsteps
    TYPE(vars_met), DIMENSION(1:mp)                        :: vmet
    TYPE(vars),     DIMENSION(1:mp)                        :: vlit
    TYPE(vars_snow), DIMENSION(1:mp)                       :: vsnow
    TYPE(vars),     DIMENSION(1:n) :: var
    REAL(r_2),      DIMENSION(1:mp)                        :: T0, Tsurface
    REAL(r_2),      DIMENSION(1:mp)                        :: Hcum, lEcum
    REAL(r_2),      DIMENSION(1:mp)                        :: Gcum, Qadvcum
    REAL(r_2),      DIMENSION(1:mp)                        :: Jcol_sensible, Jcol_latent_S, Jcol_latent_T
    REAL(r_2),      DIMENSION(1:n) :: csoil, kth
    REAL(r_2),      DIMENSION(1:n) :: phi
    REAL(r_2),      DIMENSION(1:mp)                        :: dxL
    REAL(r_2),      DIMENSION(1:mp)                        :: zdelta, SL, Tl
    TYPE(params),   DIMENSION(1:mp)                        :: plit
    TYPE(params),   DIMENSION(1:n) :: par
    REAL(r_2),      DIMENSION(1:mp,1:n), OPTIONAL                    :: wex
    REAL(r_2),      DIMENSION(1:mp,1:nsnow_max)            :: ciso_snow
    REAL(r_2),      DIMENSION(1:mp,1:nsnow_max)            :: cisoice_snow
    REAL(r_2),      DIMENSION(1:mp), OPTIONAL                        :: qali
    REAL(r_2),      DIMENSION(1:mp,-nsnow_max:n)           :: qvsig, qlsig, qvTsig, qvh
    REAL(r_2),      DIMENSION(1:mp)                        :: deltaTa
    REAL(r_2),    DIMENSION(1:mp)                          :: precip, qevap
    REAL(r_2),    DIMENSION(1:mp)                          :: qL, qhL, qybL, qTbL, qhTbL, qhybL, rexcol, wcol
    LOGICAL,      DIMENSION(1:mp)                          :: again, getq0,getqn,init
    LOGICAL,      DIMENSION(1:mp,1:n)                      :: again_ice
    INTEGER(i_d), DIMENSION(1:mp)                          :: ih0, iok, itmp, ns, nsat, nsatlast, nsteps0
    REAL(r_2),    DIMENSION(1:mp)                          :: accel, dmax, dt, dwinfil, dwoff, fac, phip
    REAL(r_2),    DIMENSION(1:mp)                          :: qpme, rsig, rsigdt, sig, t
    REAL(r_2),    DIMENSION(1:n) :: hint, phimin, qexd
    REAL(r_2),    DIMENSION(-nsnow_max+1:n) :: aa, bb, cc, dd, ee, ff, gg, dy
    REAL(r_2),    DIMENSION(-nsnow_max+1:n) :: aah, bbh, cch, ddh, eeh, ffh, ggh, de
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qadv, qadvya, qadvyb, qadvTa, qadvTb
    TYPE(vars)                                             :: vtmp
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qsig, qhsig, qadvsig
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qliq, qv, qvT, qlya, qlyb, qvya, qvyb, qlTb, qvTa, qvTb
    REAL(r_2),    DIMENSION(1:n) :: deltaS, dTsoil
    REAL(r_2),    DIMENSION(0:n) :: tmp2d1, tmp2d2
    REAL(r_2),    DIMENSION(1:n) :: cv0
    REAL(r_2),      DIMENSION(1:nsnow_max) :: cisoliqice_snow
    REAL(r_2),    DIMENSION(1:n) :: dthetaldT, thetal
    INTEGER(i_d), DIMENSION(1:n) :: isave, nsteps_ice, imelt
    TYPE(vars),         DIMENSION(1:mp)                    :: vtop, vbot
    TYPE(vars_aquifer), DIMENSION(1:mp)                    :: v_aquifer
    REAL(r_2),          DIMENSION(1:mp)                    :: dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge
    REAL(r_2),          DIMENSION(1:mp)                    :: dJcol_latent_S, dJcol_latent_T, dJcol_sensible
    REAL(r_2),          DIMENSION(1:n) :: deltaJ_latent_S, deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T
    REAL(r_2),          DIMENSION(1:mp)                    :: qevapsig
    REAL(r_2),          DIMENSION(1:mp)                    :: qrunoff
    REAL(r_2),          DIMENSION(1:mp)                    :: tmp1d1, tmp1d2, tmp1d3,  tmp1d4
    REAL(r_2),          DIMENSION(1:mp)                    :: deltah0
    REAL(r_2),          DIMENSION(1:mp)                    :: deltaSL
    REAL(r_2),          DIMENSION(1:mp)                    :: lE0, G0, Epot
    REAL(r_2),          DIMENSION(1:mp)                    :: Tfreezing
    REAL(r_2),          DIMENSION(1:mp)                    :: dtdT
    REAL(r_2),          DIMENSION(-nsnow_max+1:n) :: LHS, RHS, LHS_h, RHS_h
    INTEGER(i_d),       DIMENSION(1:mp)                    :: surface_case
    INTEGER(i_d),       DIMENSION(1:mp)                    :: nns, iflux
    LOGICAL                                                :: litter
    INTEGER(i_d)                                           :: i, j, k, kk, condition
    INTEGER(i_d)                                           :: littercase, isotopologue, advection ! switches
    REAL(r_2)                                              :: c2, theta
    REAL(r_2)                                              :: dTqwdTa, dTqwdTb, Tqw, keff
    REAL(r_2),          DIMENSION(1:mp)                    :: cp, cpeff, hice, h0_0, hice_0, h0_tmp, hice_tmp
    REAL(r_2),          DIMENSION(nsnow_max)          :: qmelt, hsnow
    REAL(r_2),  DIMENSION(1:mp)                            :: qtransfer
    REAL(r_2),          DIMENSION(nsnow_max) :: delta_snowcol, delta_snowT, delta_snowliq
    REAL(r_2),      DIMENSION(1:n) :: thetai_0, J0
    REAL(r_2)                                              :: tmp1, tmp2
    REAL(r_2),          DIMENSION(1:n) :: iqex
    INTEGER(i_d),       DIMENSION(1:mp)                    :: nfac1, nfac2, nfac3, nfac4, nfac5, &
         nfac6, nfac7, nfac8, nfac9, nfac10, nfac11, nfac12
    REAL(r_2),          DIMENSION(1:mp)                    :: J0snow, wcol0snow
    REAL(r_2), DIMENSION(1:n)                              :: h_ex
    REAL(r_2)                                              :: wpi
    ! Error flag if nstep of SLI > nsteps_max: err=0 -> no error; err/=0 -> error
    INTEGER(i_d),        INTENT(INOUT), OPTIONAL           :: err

    hsnow(:) = zero

    DO WHILE (again(kk)) ! sometimes need twice to adjust phi at satn

       nsatlast(kk) = nsat(kk) ! for detecting onset of profile saturation
       nsat(kk)     = SUM(var(:)%isat,1) ! no. of sat layers
       sig(kk)      = half
       IF (nsat(kk) /= 0) sig(kk) = one ! time weighting sigma
       rsig(kk)     = one/sig(kk)

       ! update variables
       IF (iflux(kk)==1) THEN
          ! Calc flux matric potentials (and derivatives) from S
          ! this set var-structure
          isave(:) = var(:)%isat
          ! Debug for mp=1: remove elemental from hyofS and do loop instead of next line
          CALL hyofS(S(1:n), Tsoil(1:n), par(1:n), var(1:n))
          !do i=1, n
          !   call hyofS(S(i), Tsoil(i), par(i), var(i))
          !end do

          ! End debug hyofS
          cp(kk) = REAL(1-var(1)%iice,r_2)*cswat*rhow & ! heat capacity of pond
               + REAL(var(1)%iice,r_2)*rhow* &
               ((one-var(1)%thetai/par(1)%thre)*cswat + (var(1)%thetai/par(1)%thre)*csice)
          cpeff(kk) = cp(kk) + rhow*lambdaf*var(1)%dthetaldT/par(1)%thre
          ! phip(kk) = max(var(1)%phie-var(1)%he*var(1)%Ksat, (one+e5)*var(1)%phie) !at onset of ponding
          phip(kk) = (one+e5)*var(1)%phie !at onset of ponding
          var(:)%isat  = isave(:)
          thetai(:)    = var(:)%thetai ! ice profile
          thetal(:)    = var(:)%thetal ! liq water profile
          thetai_0(:)  = thetai(:) ! initial ice profile
          dthetaldT(:) = var(:)%dthetaldT
          hice(kk)   = h0(kk)*var(1)%thetai/par(1)%thre
          hice_0(kk) = hice(kk)
          h0_0(kk)   = h0(kk)
          ! sensible heat stored in top layer + pond
          Jsensible(1)   = (var(1)%csoil* dx(1)+h0(kk)*cp(kk))*(Tsoil(1))
          ! sensible heat stored in soil column (2:n)
          Jsensible(2:n) = var(2:n)%csoil*(Tsoil(2:n))* dx(2:n)


          ! calculate this here in case snow pack disappears: ciso of transferred water needs to be retained
          IF ((isotopologue /= 0).AND.(vsnow(kk)%hsnow(1).GT.zero)) THEN
             cisoliqice_snow(1:nsnow_max) = (ciso_snow(kk,1:nsnow_max)*vsnow(kk)%hliq(1:nsnow_max) + &
                  cisoice_snow(kk,1:nsnow_max)*(vsnow(kk)%hsnow(1:nsnow_max)-vsnow(kk)%hliq(1:nsnow_max))) / &
                  vsnow(kk)%hsnow(1:nsnow_max)

          ENDIF

          CALL snow_adjust(irec, mp, n, kk, ns, h0, hice, thetai(:), dx(:), vsnow, var, par, S(:), Tsoil(:), &
               Jcol_latent_S, Jcol_latent_T, Jcol_sensible, deltaJ_sensible_S(:), qmelt(:), qtransfer, j0snow)
          thetai(1) = var(1)%thetai  ! this is the value of thetaice prior to matrix call

          CALL hyofS(S(1), Tsoil(1), par(1), var(1))

       ENDIF ! iflux==1

       ! initialise litter vars
       IF (iflux(kk)==1 .AND. (littercase==1 .OR. littercase == 2)) THEN
          ! call litter_props(Sl(kk), Tl(kk), vlit(kk), plit(kk), h0(kk))
          CALL litter_props(Sl(kk), Tsoil(1), vlit(kk), plit(kk), h0(kk))
       ENDIF

       ! ! phi is solution var at satn, so h calc from phi where S>=1 - done in hyofS above for S<1
       ! where (S(:) >= one) & ! (special for frozen soil)
       !      var(:)%h = var(:)%he + (var(:)%phi-var(:)%phie)/var(:)%Ksat


       CALL get_fluxes_and_derivs( &
            irec, mp, qprec, qprec_snow, n, dx(:), h0, &
            Tsoil(:), S(:), &
            qh(kk,:), vmet, vlit, vsnow, var, T0, Tsurface, &
            dxL, SL, Tl, &
            plit, par, &
            qali, &
            qvh(kk,:), &
            qevap, &
            again, getq0,getqn,init, ns, nsat, &
            nsatlast, &
            phip, qpme, hint(:), phimin(:), &
            qexd(:), &
            q(:), qya(:), qyb(:), qTa(:), qTb(:),qhya(:), qhyb(:), qhTa(:), qhTb(:), qadv(:), qadvya(:), qadvyb(:), &
            qadvTa(:), qadvTb(:), vtmp, qliq(:), qv(:), qvT(:), qlya(:), qlyb(:), &
            qvya(:), qvyb(:), qlTb(:), qvTa(:), qvTb(:), &
            vtop, vbot, &
            lE0, G0, &
            Epot, surface_case, &
            iflux, litter, j, kk, &
            advection, dTqwdTa, dTqwdTb, Tqw, keff, &
            hice &
            )


       CALL estimate_timestep( &
            irec, tfin, mp, n, dx(:), h0, &
            qh(kk,:), qhTa(:), qadvTa(:), qhTb(:), qadvTb(:), &
            nsteps, vlit, vsnow, var, &
            dxL, &
            plit, par, &
            deltaTa, &
            qL, qhL, &
            again, ns, nsat, &
            nsatlast, nsteps0, dmax, dt, &
            qpme, t, &
            q(:), qadv(:), &
            tmp2d1(:), tmp2d2(:), &
            tmp1d1, tmp1d3, &
            iflux, litter, kk, &
            advection, &
            iqex(:) &
            )

       !----- get and solve eqns
       rsigdt(kk) = one/(sig(kk)*dt(kk))
       IF (.NOT. again(kk))  THEN
          dwoff(kk)   = MAX(h0(kk)*(one-var(1)%thetai/par(1)%thre)-h0max,zero)
          qrunoff(kk) = dwoff(kk)/dt(kk)
       ELSE
          qrunoff(kk) = zero
       ENDIF

       ! aa, bb, cc and dd hold coeffs and rhs of linear equation eqn set
       CALL get_and_solve_eqn( &
            irec, mp, qprec,  n, dx(:), h0, S(:), thetai(:), &
            Tsoil(:),       &
            qh(kk,:), nsteps,  vlit, vsnow, var,     &
            dxL,    &
            plit, par,        &
            deltaTa, &
            qevap, qL, qhL, qybL, qTbL, qhTbL, qhybL,    &
            again,  again_ice,  iok, itmp, ns,  &
            accel,  dt,   fac,   &
            phip,  rsig, rsigdt, sig, t,      &
            qexd(:), aa(:), bb(:), cc(:), dd(:), ee(:), ff(:), gg(:), dy(:), &
            aah(:), bbh(:), cch(:), ddh(:), eeh(:), ffh(:), ggh(:), &
            de(:), q(:), qya(:), qyb(:), qTa(:), qTb(:),qhya(:), qhyb(:), qhTa(:), qhTb(:), qadv(:), qadvya(:), qadvyb(:), &
            qadvTa(:), qadvTb(:),  qsig(:), qhsig(:), qadvsig(:),      &
            dTsoil(:), tmp2d1(:),   &
            dthetaldT(:),   nsteps_ice, imelt,     &
            deltaJ_latent_T(:), deltaJ_sensible_S(:),   &
            qrunoff, tmp1d1, tmp1d2, tmp1d3,  tmp1d4,     &
            G0, &
            Tfreezing,   LHS(:), RHS(:), LHS_h(:), RHS_h(:),  nns, &
            iflux, litter, i,   kk, condition,   &
            advection,  c2, theta,     cp, &
            cpeff, hice,    h0_tmp,   hsnow(:), &
            delta_snowT(:), &
            delta_snowliq(:),    J0(:),   iqex(:),  &
            nfac1, nfac2, nfac3, nfac4, nfac5, nfac6, nfac7, nfac8, nfac9, &
            nfac10, nfac11, nfac12, err     &
            )
       IF (err /= 0) RETURN

       CALL update_unknowns( &
            irec, mp, qprec, qprec_snow, n, dx(:), h0, S(:),  &
            Tsoil(:), evap,   infil, drainage, discharge, &
            qh(kk,:), nsteps, vmet, vlit, vsnow, var,   Hcum, lEcum, &
            Gcum, Qadvcum,  &
            csoil(:), kth(:), phi(:),  zdelta, SL, Tl, &
            par,  wex,      &
            qvsig(kk,:), qlsig(kk,:), qvTsig(kk,:),   &
            precip, qevap,       rexcol, wcol,  &
            again,   ih0,   ns,  &
            dt, dwinfil, dwoff,    &
            sig,       &
            qexd(:),        dy(:),        &
            de(:), q(:), qya(:), qyb(:), qTa(:), qTb(:),qhya(:), qhyb(:), qhTa(:), qhTb(:),    &
            qsig(:), qhsig(:), qadvsig(:), qliq(:), qv(:), qvT(:),  qlyb(:), &
            qvya(:), qvyb(:), qlTb(:), qvTa(:), qvTb(:),   dTsoil(:),    &
            v_aquifer,  &
            dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge, &
            dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S(:), &
            deltaJ_latent_T(:), deltaJ_sensible_S(:), deltaJ_sensible_T(:), qevapsig, &
            qrunoff, tmp1d1, tmp1d2, tmp1d3,   deltah0,  deltaSL,  &
            LHS_h(:),  surface_case,  &
            litter, i, j,  kk,    &
            advection,    theta,   Tqw,  cp, &
            hsnow(:), &
            delta_snowcol(:), delta_snowT(:), &
            delta_snowliq(:),    J0(:),   iqex(:),  &
            J0snow, wcol0snow   &
            )

       CALL update_s_t( &
            mp,   n, dx(:), h0, S(:), thetai(:), &
            Tsoil(:),    infil,   &
            nsteps,    var,     &
            par,        &
            qlsig(kk,:),    &
            again,   ih0,   ns,  &
            dt,      &
            sig,       &
            dy(:),        &
            qhsig(:),       &
            deltaS(:), dTsoil(:),    &
            dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S(:), &
            deltaJ_latent_T(:), deltaJ_sensible_S(:), deltaJ_sensible_T(:),  &
            qrunoff, tmp1d1, tmp1d2, tmp1d3,  tmp1d4, deltah0,    &
            Tfreezing,  dtdT,   LHS_h(:),    &
            i, j, k, kk,    &
            theta,     cp, &
            hice,     hice_tmp,   &
            J0(:), tmp1, tmp2,   &
            h_ex, wpi &
            )

       IF (.NOT. again(kk)) THEN
          cv0(1:n)   = var(1:n)%cv  ! save cv before updating
          isave(1:n) = var(1:n)%iice ! save isat before updating
          ! update variables (particularly thetaice)
          ! Debug for mp=1: remove elemental from hyofS and do loop instead of next line
          CALL hyofS(S(1:n), Tsoil(1:n), par(1:n), var(1:n))
          !do i=1, n
          !   call hyofS(S(i), Tsoil(i), par(i), var(i))
          !end do
          ! End debug hyofS
          var(1:n)%iice = isave(1:n)
       ENDIF ! if .not.again

       iflux(kk) = iflux(kk) + 1

    END DO ! do while (again(kk)) ! iflux loop
  END SUBROUTINE iflux_loop

  SUBROUTINE update_s_t( &
       mp,   n, dx, h0, S, thetai, &
       Tsoil,    infil,   &
       nsteps,    var,     &
       par,        &
       qlsig,    &
       again,   ih0,   ns,  &
       dt,      &
       sig,       &
       dy,        &
       qhsig,       &
       deltaS, dTsoil,    &
       dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S, &
       deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T,  &
       qrunoff, tmp1d1, tmp1d2, tmp1d3,  tmp1d4, deltah0,    &
       Tfreezing,  dtdT,   LHS_h,    &
       i, j, k, kk,    &
       theta,     cp, &
       hice,     hice_tmp,   &
       J0, tmp1, tmp2,   &
       h_ex, wpi &
       )
    IMPLICIT NONE
    INTEGER(i_d)                                           :: mp
    INTEGER(i_d)                                           :: n
    REAL(r_2),      DIMENSION(1:n) :: dx
    REAL(r_2),      DIMENSION(1:mp)                        :: h0
    REAL(r_2),      DIMENSION(1:n) :: S
    REAL(r_2),      DIMENSION(1:n) :: thetai
    REAL(r_2),      DIMENSION(1:n) :: Tsoil
    REAL(r_2),      DIMENSION(1:mp)                        :: infil
    INTEGER(i_d),   DIMENSION(1:mp)                        :: nsteps
    TYPE(vars),     DIMENSION(1:n) :: var
    TYPE(params),   DIMENSION(1:n) :: par
    REAL(r_2),      DIMENSION(-nsnow_max:n) :: qlsig
    LOGICAL,      DIMENSION(1:mp)                          :: again
    INTEGER(i_d), DIMENSION(1:mp)                          :: ih0, ns
    REAL(r_2),    DIMENSION(1:mp)                          :: dt
    REAL(r_2),    DIMENSION(1:mp)                          :: sig
    REAL(r_2),    DIMENSION(-nsnow_max+1:n) :: dy
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qhsig
    REAL(r_2),    DIMENSION(1:n) :: deltaS, dTsoil
    REAL(r_2),          DIMENSION(1:mp)                    :: dJcol_latent_S, dJcol_latent_T, dJcol_sensible
    REAL(r_2),          DIMENSION(1:n) :: deltaJ_latent_S, deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T
    REAL(r_2),          DIMENSION(1:mp)                    :: qrunoff
    REAL(r_2),          DIMENSION(1:mp)                    :: tmp1d1, tmp1d2, tmp1d3,  tmp1d4
    REAL(r_2),          DIMENSION(1:mp)                    :: deltah0
    REAL(r_2),          DIMENSION(1:mp)                    :: Tfreezing
    REAL(r_2),          DIMENSION(1:mp)                    :: dtdT
    REAL(r_2),          DIMENSION(-nsnow_max+1:n) :: LHS_h
    INTEGER(i_d)                                           :: i, j, k, kk
    REAL(r_2)                                              :: theta
    REAL(r_2),          DIMENSION(1:mp)                    :: cp, hice, hice_tmp
    REAL(r_2),      DIMENSION(1:n) :: J0
    REAL(r_2)                                              :: tmp1, tmp2

    REAL(r_2), DIMENSION(1:n)                              :: h_ex
    REAL(r_2)                                              :: wpi
    ! update variables (S,T) to end of time step
    ! update variables to sig for use in isotope routine
    DO i=1, n

       IF (.NOT.again(kk)) Tsoil(i)  = Tsoil(i) + dTsoil(i)
       IF (var(i)%isat==0) THEN
          IF (.NOT.again(kk)) THEN
             deltaS(i) = dy(i)  !required for isotope subroutine
             S(i)      = S(i)+dy(i)
             IF (S(i)>one .AND. dy(i)>zero) THEN ! onset of saturation of layer
                var(i)%isat = 1
                var(i)%K    = var(i)%Ksat
                var(i)%phi  = var(i)%phie
             ENDIF
          ENDIF
       ELSE
          deltaS(i)  = zero    !required for isotope subroutine
          IF (i==1) THEN
             var(i)%phi = var(i)%phi + dy(i)*REAL(ns(kk),r_2) ! pond included in top soil layer
          ELSE
             var(i)%phi = var(i)%phi + dy(i)
          ENDIF

          ! var(i)%phi = zero from Ross
          ! VH thinks it is o.k. because hyofS is called later again
          ! MC think that we should probably get rid of -he*Ksat in phip etc.
          !   because we merged the pond with the first layer.
          IF (i==1 .AND. ih0(kk)/=0 .AND. var(i)%phi>=var(i)%phie) var(i)%phi = zero ! pond gone
          IF (var(i)%phi<var(i)%phie .AND. var(i)%iice==0) THEN ! desaturation of layer
             var(i)%isat = 0
             var(i)%K    = var(i)%Ksat
             var(i)%phi  = var(i)%phie
             var(i)%KS   = par(i)%KSe
             var(i)%phiS = par(i)%phiSe
          ELSEIF (var(i)%phi<var(i)%phie .AND. var(i)%iice==1) THEN
             var(i)%isat = 0
             var(i)%K    = var(i)%Ksat
             var(i)%phi  = var(i)%phie
             var(i)%KS   = var(i)%KS
             var(i)%phiS = var(i)%phiS
          ENDIF
       END IF

       IF (.NOT. again(kk)) THEN
!$                   if (h0(kk)<zero .and. var(1)%isat==0 .and. (i==1)) then ! start negative pond correction
!$                      infil(kk)     = infil(kk)+h0(kk)
!$                      ! negative pond converted to loss of soil moisture in top two layers
!$                      S(1) = S(1) + h0(kk)/(dx(1)*par(1)%thre)* dx(1)/(dx(1)+dx(2))
!$                      deltaS(1) = h0(kk)/(dx(1)*par(1)%thre)* dx(1)/(dx(1)+dx(2))
!$                      ! adjust flux between layer 1 and layer 2
!$                      qlsig(1) = qlsig(1) + h0(kk)* dx(2)/(dx(1)+dx(2))/dt(kk)
!$                      deltaS(2) =   h0(kk)/(dx(2)*par(2)%thre)* dx(2)/(dx(1)+dx(2))
!$                      if (var(2)%isat.eq.0) then
!$                         dy(2) = dy(2) + deltaS(2)
!$                      else
!$                         dy(2) = deltaS(2)
!$                         var(2)%isat = 0
!$                      endif
!$
!$                      deltah0(kk) = -(h0(kk)-deltah0(kk)) ! whole pond lost (new change = - original pond height)
!$                      h0(kk) = zero  ! zero pond remaining
!$                   endif

          ! water in profile initially
          wpi = SUM((par(:)%thr + (par(:)%the-par(:)%thr)*S(:))*dx(:))
          IF (h0(kk)<zero .AND. var(1)%isat==0 .AND. (i==1)) THEN ! start negative pond correction
             wpi = SUM((par(:)%thr + (par(:)%the-par(:)%thr)*S(:))*dx(:))

             infil(kk)     = infil(kk)+h0(kk)
             DO j=1,n
                h_ex(j) = -h0(kk)* (par(j)%the-par(j)%thr)*S(j)* dx(j)/wpi
             ENDDO

             DO j=1,n-1
                qlsig(j) = qlsig(j) +SUM(h_ex(j:n))
             ENDDO

             DO j = 1,n
                IF (j<n .AND. ( S(j)*(dx(j)*par(j)%thre) + h_ex(j)).LT.zero) THEN
                   h_ex(j+1) = h_ex(j+1) +(h_ex(j)- S(j)*(dx(j)*par(j)%thre))
                   h_ex(j) = S(j)*(dx(j)*par(j)%thre)
                ENDIF
                deltaS(j)= -MIN(h_ex(j)/(dx(j)*par(j)%thre),S(j))
                IF (var(j)%isat.EQ.0) THEN
                   dy(j) = dy(j) + deltaS(j)
                ELSE
                   dy(j) = deltaS(j)
                   var(j)%isat = 0
                ENDIF
                IF (j==1) S(j)=S(j)+deltaS(j)

             ENDDO
             deltah0(kk) = -(h0(kk)-deltah0(kk)) ! whole pond lost (new change = - original pond height)
             h0(kk) = zero  ! zero pond remaining

          ENDIF


          IF (S(i)<zero .AND. i<n) THEN ! start correction for negative soil moisture in any layer
             wpi = SUM((par(i+1:n)%thr + (par(i+1:n)%the-par(i+1:n)%thr)*S(i+1:n))*dx(i+1:n))

             DO j=i+1,n
                h_ex(j) = -deltaS(i)*dx(i)*par(i)%thre * &
                     (par(j)%thr + (par(j)%the-par(j)%thr)*S(j))*dx(j) /wpi
             ENDDO

             DO j=i+1,n
                qlsig(j) = qlsig(j) +SUM(h_ex(j:n))
             ENDDO

             DO j = i+1,n
                IF (j<n .AND. ( S(j)*(dx(j)*par(j)%thre) - h_ex(j)).LT.zero) THEN
                   h_ex(j+1) = h_ex(j+1) +(h_ex(j)- S(j)*(dx(j)*par(j)%thre))
                   h_ex(j) = S(j)*(dx(j)*par(j)%thre)
                ENDIF
                deltaS(j)= -MIN(h_ex(j)/(dx(j)*par(j)%thre),S(j))
                IF (var(j)%isat.EQ.0) THEN
                   dy(j) = dy(j) + deltaS(j)
                ELSE
                   dy(j) = deltaS(j)
                   var(j)%isat = 0
                ENDIF
                IF (j==i) S(j)=S(j)+deltaS(j)

             ENDDO

             S(i)=S(i)-deltaS(i)
             deltaS(i) = zero

          ENDIF

          ! corrections for onset of freezing and melting
          ! calculate freezing point temperature at new S
          Tfreezing(kk) = Tfrz(S(i), par(i)%he, one/(par(i)%lambc*freezefac))
          ! check for onset of freezing and adjust (increase) temperature to account for latent heat
          ! release by ice formation

          IF (experiment /= 184) THEN
             IF (Tsoil(i) < Tfreezing(kk) .AND. var(i)%iice==0) THEN  ! start correction for onset of freezing
                dtdT(kk)      = dthetalmaxdTh(Tfreezing(kk), S(i), par(i)%he, one/(par(i)%lambc*freezefac), &
                     par(i)%thre,par(i)%the)
                tmp1d3(kk) = 1000._r_2
                tmp1d2(kk) = Tsoil(i) - dTsoil(i)
                k = 1
                IF ((i==1).AND.(h0(kk)- deltah0(kk))>zero) THEN
                   h0(kk) = h0(kk) - deltah0(kk)  ! calculate stored heat using h0 at t-dt

                   DO WHILE ((k<nsteps_ice_max).AND.(tmp1d3(kk)>tol_dthetaldT))
                      tmp1d1(kk)  = LHS_h(i)*dt(kk)-deltaJ_sensible_S(i) + tmp1d2(kk)* &
                           (var(i)%csoil*dx(i)+cp(kk)*h0(kk)) + &
                           Tfreezing(kk)*rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * &
                           dtdT(kk)*(dx(i)+h0(kk)/par(i)%thre)
                      tmp1d1(kk)  = tmp1d1(kk) /( (var(i)%csoil*dx(i)+cp(kk)*h0(kk)) + &
                           rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * &
                           dtdT(kk)*(dx(i)+h0(kk)/par(i)%thre))
                      dtdT(kk) = (thetalmax(tmp1d1(kk), S(i), par(i)%he, one/(par(i)%lambc*freezefac), &
                           par(i)%thre, par(i)%the) &
                           - (S(i)*par(i)%thre + (par(i)%the-par(i)%thre)))/ &
                           (tmp1d1(kk) - Tfreezing(kk))
                      tmp1d3(kk) = ABS((tmp1d1(kk)- Tfreezing(kk))*dtdT(kk))
                      k = k + 1
                   ENDDO  ! end of do while
                   dTsoil(i) = tmp1d1(kk) - tmp1d2(kk)
                   Tsoil(i) = tmp1d1(kk)
                   deltaJ_sensible_T(i) = (var(i)%csoil*dx(i)+ cp(kk)*h0(kk))*dTsoil(i) + &
                        (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*(dx(1)+h0(kk)/par(1)%thre)* &
                        ((cswat-csice)*(tmp1d2(kk)))
                   deltaJ_latent_T(i) = (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*lambdaf* &
                        (dx(1)+h0(kk)/par(1)%thre)
                   h0(kk) = h0(kk) + deltah0(kk)

                   ! change in heat stored in soil column
                   dJcol_latent_S(kk) = SUM(deltaJ_latent_S(1:n))
                   dJcol_latent_T(kk) = SUM(deltaJ_latent_T(1:n))
                   dJcol_sensible(kk) = SUM(deltaJ_sensible_T(1:n)) + SUM(deltaJ_sensible_S(1:n))
                   ! tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + &
                   !      dJcol_sensible(kk) - (qhsig(0)-qhsig(6))*dt(kk) + &
                   !      qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk)
                   tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + &
                        dJcol_sensible(kk) - (qhsig(0)-qhsig(n))*dt(kk) + &
                        qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk)
                   ! if (abs(tmp1d1(kk)).gt.10.) then
                   !    write(*,*) 'E imbalance check 2'
                   ! endif
                ELSE
                   DO WHILE ((k<nsteps_ice_max).AND.(tmp1d3(kk)>tol_dthetaldT))
                      ! zeroth estimate of corrected temperature
                      tmp1d1(kk)  = LHS_h(i)*dt(kk)-deltaJ_sensible_S(i) + tmp1d2(kk)* &
                           (var(i)%csoil*dx(i)) + Tfreezing(kk)*rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * &
                           dtdT(kk)*dx(i)
                      tmp1d1(kk)  = tmp1d1(kk) /( (var(i)%csoil*dx(i)) + &
                           rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * dtdT(kk)*dx(i))
                      dtdT(kk) =  (thetalmax(tmp1d1(kk), S(i), par(i)%he, &
                           one/(par(i)%lambc*freezefac), par(i)%thre, par(i)%the) &
                           - (S(i)*par(i)%thre + (par(i)%the-par(i)%thre)))/ &
                           (tmp1d1(kk) - Tfreezing(kk))
                      tmp1d3(kk) = ABS((tmp1d1(kk)- Tfreezing(kk))*dtdT(kk))
                      k = k + 1
                   ENDDO  ! end of do while
                   dTsoil(i) = tmp1d1(kk)- (Tsoil(i)-dTsoil(i))
                   Tsoil(i) = tmp1d1(kk)
                   deltaJ_sensible_T(i) = (var(i)%csoil*dx(i))*dTsoil(i) + &
                        (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*dx(i)* &
                        ((cswat-csice)*(tmp1d2(kk)))
                   deltaJ_latent_T(i) = (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*lambdaf*dx(i)
                ENDIF
             ENDIF  ! end correction for onset of freezing
          ENDIF

          IF (var(i)%iice==1) THEN
             theta         = S(i)*(par(i)%thre) + (par(i)%the - par(i)%thre)
             tmp1d1(kk) = Tsoil(i) - dTsoil(i)
             ! energy for complete melting
             IF ((i==1).AND.(h0(kk)>zero)) THEN
                tmp1d3(kk) = thetai(i)*dx(i)*rhow*lambdaf + &
                     thetai(i)*(h0(kk)- deltah0(kk))/par(i)%thre*rhow*lambdaf

                tmp1d2(kk) = (Tfreezing(kk))*(rhow*cswat*(theta*dx(1)+h0(kk))+ &
                     par(1)%rho*par(1)%css*dx(1)) - &
                     (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(1)+(h0(kk)- deltah0(kk))* &
                     (one-thetai(1)/par(1)%thre))+par(1)%rho*par(1)%css*dx(1)) - &
                     (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(1)+(h0(kk)- deltah0(kk))* &
                     thetai(1)/par(1)%thre))
             ELSE
                tmp1d3(kk) = thetai(i)*dx(i)*rhow*lambdaf

                tmp1d2(kk) = (Tfreezing(kk))*(rhow*cswat*(theta*dx(i))+par(i)%rho*par(i)%css*dx(i)) - &
                     (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(i))+par(i)%rho*par(i)%css*dx(i)) - &
                     (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(i)))
             ENDIF
          ENDIF
          ! check for onset of thawing and decrease temperature to account for latent heat required to melt ice

          IF ((var(i)%iice==1).AND.((J0(i) + LHS_h(i)*dt(kk)).GE.JSoilLayer(Tfreezing(kk), &
               dx(i), theta,par(i)%css, par(i)%rho, &
               MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, &
               par(i)%he, one/(par(i)%lambc*freezefac)))) THEN
             ! start correction for onset of thawing
             ! Correct for "overthawing". This extra energy is used to heat the soil even more
             ! The correction comes from equating the energy balances before and after correction.
             tmp1d1(kk) = Tsoil(i) - dTsoil(i)

             IF ((i==1) .AND. (h0(kk)>zero)) THEN
                deltaJ_latent_T(i)=      tmp1d3(kk)
                deltaJ_latent_S(i) = zero
                h0(kk) = h0(kk) - deltah0(kk)

                deltaJ_sensible_S(i) = (tmp1d1(kk))*(rhow*cswat*(theta*dx(1)+(h0(kk)+deltah0(kk)))+ &
                     par(1)%rho*par(1)%css*dx(1)) - &
                     (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(1)+(h0(kk))* &
                     (one-thetai(1)/par(1)%thre))+par(1)%rho*par(1)%css*dx(1)) - &
                     (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(1)+(h0(kk))* &
                     thetai(1)/par(1)%thre))

                dTsoil(i) = (LHS_h(i)*dt(kk) - (deltaJ_latent_T(i) + deltaJ_sensible_S(i)))/ &
                     (rhow*cswat*(theta*dx(1)+(h0(kk)+deltah0(kk)))+par(1)%rho*par(1)%css*dx(1))

                deltaJ_sensible_T(i) = dTsoil(i)* &
                     (rhow*cswat*(theta*dx(1)+(h0(kk)+deltah0(kk)))+par(1)%rho*par(1)%css*dx(1))

                h0(kk) = h0(kk) + deltah0(kk)
             ELSE
                deltaJ_latent_T(i) = thetai(i)*dx(i)*rhow*lambdaf
                deltaJ_latent_S(i) = zero

                theta         = S(i)*(par(i)%thre) + (par(i)%the - par(i)%thre)
                var(i)%csoil = cswat*theta*rhow + par(i)%css*par(i)%rho

                deltaJ_sensible_S(i) = (tmp1d1(kk))*(rhow*cswat*(theta*dx(i))+ &
                     par(i)%rho*par(i)%css*dx(i)) - &
                     (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(i))+par(i)%rho*par(i)%css*dx(i)) - &
                     (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(i)))
                dTsoil(i) = (LHS_h(i)*dt(kk) - (deltaJ_latent_T(i) + deltaJ_sensible_S(i)))/ &
                     (dx(i)*var(i)%csoil)
                deltaJ_sensible_T(i) = var(i)%csoil*dTsoil(i)*dx(i)
             ENDIF
             Tsoil(i) = tmp1d1(kk) + dTsoil(i)
             ! correct thetai and T in frozen soil
          ELSEIF ((var(i)%iice==1).AND.((J0(i) + LHS_h(i)*dt(kk))<JSoilLayer(Tfreezing(kk), &
               dx(i), theta,par(i)%css, par(i)%rho, &
               MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, &
               par(i)%he, one/(par(i)%lambc*freezefac)))) THEN

             tmp1d2(kk) = J0(i) + LHS_h(i)*dt(kk) ! total energy in  soil layer
             IF (Tsoil(i) .LT. -40.0_r_2) THEN
                ! assume no lquid below min temp threshhold
                var(i)%thetal = 0.0_r_2
                var(i)%thetai = theta
                IF (i.EQ.1) hice(kk) = h0(kk)
                tmp1d3(kk) = (tmp1d2(kk) + rhow*lambdaf*(theta*dx(i) +  MERGE(h0(kk),zero,i==1)))/ &
                     (dx(i)*par(i)%css*par(i)%rho + rhow*csice*(theta*dx(i) + &
                     MERGE(h0(kk),zero,i==1)))
             ELSE
                !check there is a zero
                tmp1 = GTfrozen(Tsoil(i)-dTsoil(i)-50._r_2, tmp1d2(kk), dx(i), theta,par(i)%css, par(i)%rho, &
                     MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, par(i)%he, one/(par(i)%lambc*freezefac))

                tmp2 = GTFrozen(Tfreezing(kk), tmp1d2(kk), dx(i), theta,par(i)%css, par(i)%rho, &
                     MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, par(i)%he, one/(par(i)%lambc*freezefac))
                ! there is a zero in between
                IF ((tmp1*tmp2) < zero) THEN
                   tmp1d3(kk) = rtbis_Tfrozen(tmp1d2(kk), dx(i), theta,par(i)%css, par(i)%rho, &
                        MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, &
                        par(i)%he, one/(par(i)%lambc*freezefac), &
                        Tsoil(i)-dTsoil(i)-50._r_2, Tfreezing(kk), 0.0001_r_2)
                   tmp1d4(kk) = thetalmax(tmp1d3(kk), S(i), par(i)%he, one/(par(i)%lambc*freezefac), &
                        par(i)%thre, par(i)%the) ! liquid content at solution for Tsoil
                ELSE
                   WRITE(logn,*) "Found no solution for Tfrozen 1. ", kk, i
                   WRITE(logn,*) "Assume soil is totally frozen"
                   var(i)%thetal = 0.0_r_2
                   var(i)%thetai = theta
                   IF (i.EQ.1) hice(kk) = h0(kk)
                   tmp1d3(kk) = (tmp1d2(kk) + rhow*lambdaf*(theta*dx(i) +  MERGE(h0(kk),zero,i==1)))/ &
                        (dx(i)*par(i)%css*par(i)%rho + rhow*csice*(theta*dx(i) + &
                        MERGE(h0(kk),zero,i==1)))
                   WRITE(logn,*) "frozen soil temperature: ", tmp1d3(kk)
                   WRITE(logn,*) nsteps(kk), S(i), Tsoil(i), dTsoil(i), h0(kk), tmp1, tmp2, tmp1d2(kk), theta, &
                        JSoilLayer(Tfreezing(kk), &
                        dx(i), theta,par(i)%css, par(i)%rho, &
                        MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, &
                        par(i)%he, one/(par(i)%lambc*freezefac)), J0(i) + LHS_h(i)*dt(kk), Tfreezing(kk)
                ENDIF
                var(i)%thetal = tmp1d4(kk)
                var(i)%thetai = theta - tmp1d4(kk)
             ENDIF ! Tsoil <-40
             IF (i==1) THEN
                hice_tmp(kk) = hice(kk)
                hice(kk) = h0(kk)*var(1)%thetai/par(1)%thre

                ! correct total energy stored in pond + soil
                deltaJ_latent_S(i) = - rhow*lambdaf*((hice(kk)-hice_tmp(kk)) + &
                     dx(1)*(var(1)%thetai-thetai(1)))
             ELSE
                ! correct total energy stored in  soil
                deltaJ_latent_S(i) = - rhow*lambdaf*( + &
                     dx(i)*(var(i)%thetai-thetai(i)))


             ENDIF

             deltaJ_latent_T(i) = zero

             deltaJ_sensible_T(i) = JSoilLayer(tmp1d3(kk), &
                  dx(i), theta,par(i)%css, par(i)%rho, &
                  MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, &
                  par(i)%he, one/(par(i)%lambc*freezefac)) - &
                  J0(i)- &
                  deltaJ_latent_S(i)

             deltaJ_sensible_S(i) = 0._r_2

             thetai(i) = var(i)%thetai
             Tsoil(i) = tmp1d3(kk)

          ENDIF ! soil remains frozen

       ENDIF ! if .not.again

    END DO ! i=1, n => update variables S,T
  END SUBROUTINE update_s_t


  SUBROUTINE update_unknowns( &
       irec, mp, qprec, qprec_snow, n, dx, h0, S,  &
       Tsoil, evap,   infil, drainage, discharge, &
       qh, nsteps, vmet, vlit, vsnow, var,   Hcum, lEcum, &
       Gcum, Qadvcum,  &
       csoil, kth, phi,  zdelta, SL, Tl, &
       par,  wex,      &
       qvsig, qlsig, qvTsig,   &
       precip, qevap,       rexcol, wcol,  &
       again,   ih0,   ns,  &
       dt, dwinfil, dwoff,    &
       sig,       &
       qexd,        dy,        &
       de, q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb,    &
       qsig, qhsig, qadvsig, qliq, qv, qvT,  qlyb, &
       qvya, qvyb, qlTb, qvTa, qvTb,   dTsoil,    &
       v_aquifer,  &
       dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge, &
       dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S, &
       deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T, qevapsig, &
       qrunoff, tmp1d1, tmp1d2, tmp1d3,   deltah0,  deltaSL,  &
       LHS_h,  surface_case,  &
       litter, i, j,  kk,    &
       advection,    theta,   Tqw,  cp, &
       hsnow, &
       delta_snowcol, delta_snowT, &
       delta_snowliq,    J0,   iqex,  &
       J0snow, wcol0snow   &
       )
    IMPLICIT NONE
    INTEGER(i_d)                                           :: irec, mp
    REAL(r_2),      DIMENSION(1:mp)                        :: qprec
    REAL(r_2),      DIMENSION(1:mp)                        :: qprec_snow
    INTEGER(i_d)                                           :: n
    REAL(r_2),      DIMENSION(1:n) :: dx
    REAL(r_2),      DIMENSION(1:mp)                        :: h0
    REAL(r_2),      DIMENSION(1:n) :: S
    REAL(r_2),      DIMENSION(1:n) :: Tsoil
    REAL(r_2),      DIMENSION(1:mp)                        :: evap
    REAL(r_2),      DIMENSION(1:mp)                        :: infil
    REAL(r_2),      DIMENSION(1:mp)                        :: drainage, discharge
    REAL(r_2),      DIMENSION(-nsnow_max:n) :: qh
    INTEGER(i_d),   DIMENSION(1:mp)                        :: nsteps
    TYPE(vars_met), DIMENSION(1:mp)                        :: vmet
    TYPE(vars),     DIMENSION(1:mp)                        :: vlit
    TYPE(vars_snow), DIMENSION(1:mp)                       :: vsnow
    TYPE(vars),     DIMENSION(1:n) :: var
    REAL(r_2),      DIMENSION(1:mp)                        :: Hcum, lEcum
    REAL(r_2),      DIMENSION(1:mp)                        :: Gcum, Qadvcum
    REAL(r_2),      DIMENSION(1:n) :: csoil, kth
    REAL(r_2),      DIMENSION(1:n) :: phi
    REAL(r_2),      DIMENSION(1:mp)                        :: zdelta, SL, Tl
    TYPE(params),   DIMENSION(1:n) :: par
    REAL(r_2),      DIMENSION(1:mp,1:n), OPTIONAL                    :: wex
    REAL(r_2),      DIMENSION(-nsnow_max:n) :: qvsig, qlsig, qvTsig
    REAL(r_2),    DIMENSION(1:mp)                          :: precip, qevap
    REAL(r_2),    DIMENSION(1:mp)                          :: rexcol, wcol
    LOGICAL,      DIMENSION(1:mp)                          :: again
    INTEGER(i_d), DIMENSION(1:mp)                          :: ih0, ns
    REAL(r_2),    DIMENSION(1:mp)                          :: dt, dwinfil, dwoff
    REAL(r_2),    DIMENSION(1:mp)                          :: sig
    REAL(r_2),    DIMENSION(1:n) :: qexd
    REAL(r_2),    DIMENSION(-nsnow_max+1:n) :: dy
    REAL(r_2),    DIMENSION(-nsnow_max+1:n) :: de
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qsig, qhsig, qadvsig
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qliq, qv, qvT, qlyb, qvya, qvyb, qlTb, qvTa, qvTb
    REAL(r_2),    DIMENSION(1:n) :: dTsoil
    TYPE(vars_aquifer), DIMENSION(1:mp)                    :: v_aquifer
    REAL(r_2),          DIMENSION(1:mp)                    :: dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge
    REAL(r_2),          DIMENSION(1:mp)                    :: dJcol_latent_S, dJcol_latent_T, dJcol_sensible
    REAL(r_2),          DIMENSION(1:n) :: deltaJ_latent_S, deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T
    REAL(r_2),          DIMENSION(1:mp)                    :: qevapsig
    REAL(r_2),          DIMENSION(1:mp)                    :: qrunoff
    REAL(r_2),          DIMENSION(1:mp)                    :: tmp1d1, tmp1d2, tmp1d3
    REAL(r_2),          DIMENSION(1:mp)                    :: deltah0
    REAL(r_2),          DIMENSION(1:mp)                    :: deltaSL
    REAL(r_2),          DIMENSION(-nsnow_max+1:n) :: LHS_h
    INTEGER(i_d),       DIMENSION(1:mp)                    :: surface_case
    LOGICAL                                                :: litter
    INTEGER(i_d)                                           :: i, j, kk
    INTEGER(i_d)                                           :: advection ! switches
    REAL(r_2)                                              :: theta
    REAL(r_2)                                              :: Tqw
    REAL(r_2),          DIMENSION(1:mp)                    :: cp
    REAL(r_2),          DIMENSION(nsnow_max) :: hsnow
    REAL(r_2),          DIMENSION(nsnow_max) :: delta_snowcol, delta_snowT, delta_snowliq
    REAL(r_2),      DIMENSION(1:n) :: J0
    REAL(r_2),          DIMENSION(1:n) :: iqex

    REAL(r_2),          DIMENSION(1:mp)                    :: J0snow, wcol0snow
    !----- update unknowns
    ! i.e. update state variables to end of time step
    !      cumulate some surface fluxes
    ih0(kk) = 0
    IF (.NOT. again(kk)) THEN
       ! get surface fluxes at sigma of the time step
       dwoff(kk)  = zero
       deltah0(kk) = zero
       precip(kk) = precip(kk) + (qprec(kk)+qprec_snow(kk))*dt(kk)
       dwcol(kk) = SUM(par(1:n)%thre*dx(1:n)*dy(1:n)*(ABS(var(1:n)%isat-1)),1) + &
            REAL(1-ns(kk),r_2)*dy(1)

       IF (vsnow(kk)%nsnow==1) dwcol(kk) = dwcol(kk) + dy(0)

       ! absolute heat stored in soil column before update
       DO j=1, n
          theta  =  S(j)*(par(j)%thre) + (par(j)%the - par(j)%thre)
          J0(j) = JSoilLayer(Tsoil(j), &
               dx(j), theta,par(j)%css, par(j)%rho, &
               MERGE(h0(kk),zero,j.EQ.1), par(j)%thre, par(j)%the, &
               par(j)%he, one/(par(j)%lambc*freezefac))
       END DO

       ! change in heat stored in soil column
       DO j=1, n
          IF (j==1) THEN
             ! pond included in top soil layer
             deltaJ_latent_S(j) = -dx(j)*dy(j)*par(j)%thre*REAL(var(j)%iice,r_2)* &
                  REAL(1-var(j)%isat,r_2)*rhow*lambdaf - &
                  REAL(1-ns(kk),r_2)*dy(1)*REAL(var(j)%iice,r_2)*rhow*lambdaf*(var(1)%thetai/par(1)%thre)

             deltaJ_latent_T(j) = var(j)%dthetaldT*dx(j)*dTsoil(j)*rhow*lambdaf* &
                  REAL(var(j)%iice,r_2) + &
                  h0(kk)/par(1)%thre*var(j)%dthetaldT*dTsoil(j)*rhow*lambdaf*REAL(var(j)%iice,r_2)

             deltaJ_sensible_T(j) = var(j)%csoil*dx(j)*dTsoil(j) + &
                  cp(kk)*h0(kk)*dTsoil(j) + &
                  var(1)%iice*var(j)%dthetaldT*(Tsoil(j))*rhow*(cswat-csice)*dx(j)*dTsoil(j) + &
                  var(j)%iice*var(j)%dthetaldT*(Tsoil(j))*rhow*(cswat-csice)* &
                  h0(kk)/par(j)%thre*dTsoil(j)

             IF (advection==1) THEN
                deltaJ_sensible_S(j) = rhow*cswat*(Tsoil(j))*dx(j)*par(j)%thre*dy(j)* &
                     REAL(1-var(j)%isat,r_2)*REAL(1-var(j)%iice,r_2)  + &
                     rhow*csice*(Tsoil(j))*dx(j)*par(j)%thre*dy(j)* &
                     REAL(1-var(j)%isat,r_2)*REAL(var(j)%iice,r_2) + &
                     REAL(1-ns(kk),r_2)*REAL(1-var(1)%iice,r_2)* &
                     dy(1)*rhow*cswat*(Tsoil(1))  + &
                     REAL(1-ns(kk),r_2)*REAL(var(1)%iice,r_2)* &
                     dy(1)*rhow*csice*(Tsoil(1))*(var(1)%thetai/par(1)%thre) + &
                     REAL(1-ns(kk),r_2)*REAL(var(1)%iice,r_2)* &
                     dy(1)*rhow*cswat*(Tsoil(1))*(one-(var(1)%thetai/par(1)%thre))
             ENDIF
          ELSE ! (j>1)
             deltaJ_latent_S(j) = -dx(j)*dy(j)*par(j)%thre*REAL(var(j)%iice,r_2)* &
                  REAL(1-var(j)%isat,r_2)*rhow*lambdaf
             deltaJ_sensible_T(j) = var(j)%csoil*dx(j)*dTsoil(j) + &
                  var(j)%iice*var(j)%dthetaldT*(Tsoil(j))*rhow*(cswat-csice)*dx(j)*dTsoil(j)


             IF (advection==1) THEN
                deltaJ_sensible_S(j) = rhow*cswat*(Tsoil(j))*dx(j)*dy(j)*par(j)%thre* &
                     REAL(1-var(j)%isat,r_2)*REAL(1-var(j)%iice,r_2) &
                     + rhow*csice*(Tsoil(j))*dx(j)*dy(j)*par(j)%thre &
                     *REAL(1-var(j)%isat,r_2)*REAL(var(j)%iice,r_2)
             ENDIF

             deltaJ_latent_T(j) = var(j)%dthetaldT*dx(j)*dTsoil(j)*rhow*lambdaf*REAL(var(j)%iice,r_2)
          ENDIF
       END DO

       ! change in heat stored in soil column
       dJcol_latent_S(kk) = SUM(deltaJ_latent_S(1:n))
       dJcol_latent_T(kk) = SUM(deltaJ_latent_T(1:n))
       dJcol_sensible(kk) = SUM(deltaJ_sensible_T(1:n)) + SUM(deltaJ_sensible_S(1:n))
       ! tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + dJcol_sensible(kk) - (qhsig(0)-qhsig(6))*dt(kk) &
       !      + qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk)
       tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + dJcol_sensible(kk) - (qhsig(0)-qhsig(n))*dt(kk) &
            + qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk)

       IF (ns(kk)<1) THEN ! change in pond height
          h0(kk) = h0(kk) + dy(1)
          deltah0(kk) = dy(1)
          IF (h0(kk)<zero .AND. dy(1)<zero) THEN
             ih0(kk)=1 ! pond gone
          ENDIF
       ENDIF

       ! cumulate evaporation from top of soil column or top of litter/pond or top of snow pack
       SELECT CASE (surface_case(kk))
       CASE(1)  ! no snow
          evap(kk)     = evap(kk) +qevap(kk)*dt(kk) - sig(kk)*(qyb(0)*dy(1)+qTb(0)*dTsoil(1))*dt(kk)
          qevapsig(kk) = qevap(kk) - sig(kk)*(qyb(0)*dy(1)+qTb(0)*dTsoil(1))
          Gcum(kk) = Gcum(kk)+(qhsig(0)-qadvsig(0))*dt(kk)
          Qadvcum(kk) = Qadvcum(kk) + qadvsig(0)*dt(kk) - qadvsig(n)*dt(kk)
          lEcum(kk)    = lEcum(kk) + qevapsig(kk)*thousand*var(1)%lambdav*dt(kk)
          Hcum(kk) = Hcum(kk) + (vmet(kk)%Rn*dt(kk)-(qhsig(0)-qadvsig(0))*dt(kk)- &
               qevapsig(kk)*thousand*var(1)%lambdav*dt(kk))
          dwinfil(kk) = (qprec(kk)+qprec_snow(kk)-qevapsig(kk))*dt(kk)
       CASE(2) ! dedicated snow layer
          evap(kk)     = evap(kk) +qevap(kk)*dt(kk) - sig(kk)*(qyb(-vsnow(kk)%nsnow)*dy(1-vsnow(kk)%nsnow)+ &
               qTb(-vsnow(kk)%nsnow)*de(1-vsnow(kk)%nsnow))*dt(kk)
          qevapsig(kk) = qevap(kk) - sig(kk)*(qyb(-vsnow(kk)%nsnow)*dy(1-vsnow(kk)%nsnow)+ &
               qTb(-vsnow(kk)%nsnow)*de(1-vsnow(kk)%nsnow))
          Gcum(kk) = Gcum(kk)+(qhsig(-vsnow(kk)%nsnow)-qadvsig(-vsnow(kk)%nsnow))*dt(kk)
          Qadvcum(kk) = Qadvcum(kk) + qadvsig(-vsnow(kk)%nsnow)*dt(kk) - qadvsig(n)*dt(kk)
          IF (vsnow(kk)%hliq(1).GT.zero) THEN
             lEcum(kk)    = lEcum(kk) + qevapsig(kk)*thousand*lambdaf*dt(kk)
             Hcum(kk) = Hcum(kk) + (vmet(kk)%Rn*dt(kk)-(qhsig(-vsnow(kk)%nsnow)-qadvsig(-vsnow(kk)%nsnow))*dt(kk)- &
                  qevapsig(kk)*thousand*lambdaf*dt(kk))
          ELSE
             lEcum(kk)    = lEcum(kk) + qevapsig(kk)*thousand*lambdas*dt(kk)
             Hcum(kk) = Hcum(kk) + (vmet(kk)%Rn*dt(kk)-(qhsig(-vsnow(kk)%nsnow)-qadvsig(-vsnow(kk)%nsnow))*dt(kk)- &
                  qevapsig(kk)*thousand*lambdas*dt(kk))
          ENDIF
          dwinfil(kk) = (qprec(kk)+qprec_snow(kk)-qevapsig(kk))*dt(kk)
       END SELECT ! surface_case

       dwdrainage(kk)     = q(n)*dt(kk) +sig(kk)*dt(kk)*(qya(n)*dy(n)+qTa(n)*dTsoil(n))
       IF (botbc=="aquifer" .AND. v_aquifer(kk)%isat==0) THEN
          dwdischarge(kk) = v_aquifer(kk)%discharge*dt(kk)
       ELSE
          dwdischarge(kk) = dwdrainage(kk)
       ENDIF

       drexcol(kk) = SUM(iqex(:),1)*dt(kk)
       ! cumulative fluxes
       infil(kk)     = infil(kk) + dwinfil(kk)
       inlit(kk)     = inlit(kk) + dwinlit(kk)
       wcol(kk)      = wcol(kk) + dwcol(kk)
       drainage(kk)  = drainage(kk) + dwdrainage(kk)
       discharge(kk) = discharge(kk) + qrunoff(kk)*dt(kk)
       rexcol(kk)    = rexcol(kk) + drexcol(kk)
       Qadvcum(kk) = Qadvcum(kk) - SUM(iqex(:)*(Tsoil(:)),1)*dt(kk)*rhow*cswat - &
            qrunoff(kk)*(Tsoil(1))*dt(kk)*rhow*cswat
       ! if (irec.eq.109) then
       !     write(*,"(1i8,16e16.6)") irec, infil(kk)-(wcol(kk)+discharge(kk)+drainage(kk))-rexcol(kk), &
       !         wcol(kk), deltah0(kk), h0(kk), discharge(kk), evap(kk), &
       !         drainage(kk), rexcol(kk)
       !     write(*,"(1i8,16e16.6)") irec, dwinfil(kk)-dwcol(kk)-dwdrainage(kk)-drexcol(kk)-qrunoff(kk)*dt(kk)
       ! endif

       ! evaluate soil fluxes at sigma of time step for use in isotope routine

       SELECT CASE (surface_case(kk))
       CASE(1)  ! no snow
          qsig(0)  = q(0)  + sig(kk)*(qyb(0)*dy(1)  + qya(0)*dy(0)  + qTb(0)*de(1) &
               + qTa(0)*de(0))
          qhsig(0) = qh(0) + sig(kk)*(qhyb(0)*dy(1) + qhya(0)*dy(0) + qhTb(0)*de(1) &
               + qhTa(0)*de(0))
          qvsig(0) = qv(0)+sig(kk)*qvyb(0)*dy(1)+sig(kk)*qvTb(0)*de(1)
          qlsig(0) = qliq(0)+sig(kk)*qlyb(0)*dy(1)+sig(kk)*qlTb(0)*de(1)

          qsig(1:n-1)   = q(1:n-1) + sig(kk)*(qya(1:n-1)*dy(1:n-1) + qyb(1:n-1)*dy(2:n) &
               + qTa(1:n-1)*de(1:n-1) + qTb(1:n-1)*de(2:n))
          qsig(n)       = q(n) + sig(kk)*(qya(n)*dy(n) + qTa(n)*de(n))
          qhsig(1:n-1)  = qh(1:n-1) + sig(kk)*(qhya(1:n-1)*dy(1:n-1) + qhyb(1:n-1)*dy(2:n) &
               + qhTa(1:n-1)*de(1:n-1) + qhTb(1:n-1)*de(2:n))
          qhsig(n)      = qh(n) + sig(kk)*(qhya(n)*dy(n) + qhTa(n)*de(n))
          qvsig(1:n-1)  = qv(1:n-1) + sig(kk)*(qvya(1:n-1)*dy(1:n-1) + qvyb(1:n-1)*dy(2:n) &
               + qvTa(1:n-1)*de(1:n-1) + qvTb(1:n-1)*de(2:n))
          qvsig(n)      = zero

          qvTsig(0)     = qvT(0) + sig(kk)*qvTb(0)*de(1)
          qvTsig(1:n-1) = qvT(1:n-1) + sig(kk)*(qvTa(1:n-1)*de(1:n-1) + qvTb(1:n-1)*de(2:n))
          qvTsig(n)     = zero
          qlsig(1:n-1)  = qsig(1:n-1) - qvsig(1:n-1)
          qlsig(n)      = qsig(n)

       CASE(2) ! dedicated snow layer
          qsig(-vsnow(kk)%nsnow)  = q(-vsnow(kk)%nsnow)  + &
               sig(kk)*(qyb(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow+1)  + &
                                ! qya(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow)  + &
               qTb(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow+1)) ! &
          !  + qTa(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow))
          qhsig(-vsnow(kk)%nsnow) = qh(-vsnow(kk)%nsnow) + &
               sig(kk)*(qhyb(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow+1) + &
                                ! qhya(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow) + &
               qhTb(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow+1)) !&
          !  + qhTa(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow))
          qvsig(-vsnow(kk)%nsnow) = qv(-vsnow(kk)%nsnow) + &
               sig(kk)*qvyb(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow+1) + &
               sig(kk)*qvTb(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow+1)

          qlsig(-vsnow(kk)%nsnow:-1) = zero
          qlsig(0) = qliq(0)+sig(kk)*qlyb(0)*dy(1)+sig(kk)*qlTb(0)*de(1)

          qsig(-vsnow(kk)%nsnow+1:n-1)   = q(-vsnow(kk)%nsnow+1:n-1) + &
               sig(kk)*(qya(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+1:n-1) + &
               qyb(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+2:n) &
               + qTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + &
               qTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n))
          qsig(n)       = q(n) + sig(kk)*(qya(n)*dy(n) + qTa(n)*de(n))
          qhsig(-vsnow(kk)%nsnow+1:n-1)  = qh(-vsnow(kk)%nsnow+1:n-1) + &
               sig(kk)*(qhya(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+1:n-1) + &
               qhyb(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+2:n) &
               + qhTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + &
               qhTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n))
          qhsig(n)      = qh(n) + sig(kk)*(qhya(n)*dy(n) + qhTa(n)*de(n))
          qvsig(-vsnow(kk)%nsnow+1:n-1)  = qv(-vsnow(kk)%nsnow+1:n-1) + &
               sig(kk)*(qvya(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+1:n-1) + &
               qvyb(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+2:n) &
               + qvTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + &
               qvTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n))
          qvsig(n)      = zero

          qvTsig(0)     = qvT(0) + sig(kk)*qvTb(0)*de(1)
          qvTsig(-vsnow(kk)%nsnow+1:n-1) = qvT(-vsnow(kk)%nsnow+1:n-1) + &
               sig(kk)*(qvTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + &
               qvTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n))
          qvTsig(n)     = zero
          qlsig(-vsnow(kk)%nsnow+1:n-1)  = qsig(-vsnow(kk)%nsnow+1:n-1) - qvsig(-vsnow(kk)%nsnow+1:n-1)
          qlsig(n)      = qsig(n)
       END SELECT ! surface_case

       IF (botbc=="aquifer") THEN ! update aquifer props
          IF (v_aquifer(kk)%isat==0) THEN
             IF (v_aquifer(kk)%WA+(dwdrainage(kk)-dwdischarge(kk))*dt(kk) > &
                  (v_aquifer(kk)%zzero-v_aquifer(kk)%zsoil)*v_aquifer(kk)%Sy) THEN
                v_aquifer(kk)%isat = 1  ! aquifer saturated
                S(n) = S(n) + (-(v_aquifer(kk)%WA+(dwdrainage(kk)-dwdischarge(kk))*dt(kk)) &
                     +(v_aquifer(kk)%zzero-v_aquifer(kk)%zsoil)*v_aquifer(kk)%Sy)/(dx(n)*par(n)%thre)
             ENDIF
             v_aquifer(kk)%WA        = MIN(v_aquifer(kk)%WA+(dwdrainage(kk)-dwdischarge(kk))*dt(kk), &
                  (v_aquifer(kk)%zzero-v_aquifer(kk)%zsoil)*v_aquifer(kk)%Sy)
             v_aquifer(kk)%zdelta    = v_aquifer(kk)%zzero - v_aquifer(kk)%Wa/v_aquifer(kk)%Sy
             ! new discharge rate
             v_aquifer(kk)%discharge = v_aquifer(kk)%Rsmax*EXP(-v_aquifer(kk)%f*v_aquifer(kk)%zdelta)
          ELSEIF (v_aquifer(kk)%isat==1) THEN
             ! check for desat of aquifer
             IF (dwdrainage(kk) < v_aquifer(kk)%Rsmax*EXP(-v_aquifer(kk)%f*v_aquifer(kk)%zsoil)) THEN
                v_aquifer(kk)%isat      = 0
                v_aquifer(kk)%zdelta    = v_aquifer(kk)%zsoil
                ! new discharge rate
                v_aquifer(kk)%discharge = v_aquifer(kk)%Rsmax*EXP(-v_aquifer(kk)%f*v_aquifer(kk)%zdelta)
             ENDIF
          ENDIF
          zdelta(kk) = v_aquifer(kk)%zdelta
       ENDIF       ! end update aquifer props

       IF (botbc=="constant head") THEN
          drn(kk) = drn(kk)+(q(n)+sig(kk)*qya(n)*dy(n))*dt(kk)
       ELSE
          drn(kk) = drn(kk)+(q(n)+sig(kk)*qya(n)*dy(n))*dt(kk)
       END IF

       IF (PRESENT(wex)) THEN
          wex(kk,1:n) = wex(kk,1:n)+(iqex(1:n)+SPREAD(sig(kk),1,n)*qexd(1:n)*dy(1:n))*SPREAD(dt(kk),1,n)
       END IF

       IF (litter .AND. ns(kk)==1 ) THEN ! adjust litter moisture content if no ponding
          SL(kk)      = SL(kk) + dy(0)
          deltaSL(kk) = dy(0)
       ELSE
          deltaSL(kk) = zero
       ENDIF

       TL(kk) = TL(kk) + de(0)   ! Tlitter assumed the same as pond T
       IF (SL(kk) >= one) vlit(kk)%isat = 1   ! check for litter saturation

       csoil(:) = var(1:n)%csoil
       kth(:) = var(1:n)%kth
       phi(:) = var(1:n)%phi

       ! update snow pack water content and temperature (or liquid water)
       IF (vsnow(kk)%nsnow>0) THEN
          vsnow(kk)%Qadv_snow = vsnow(kk)%Qadv_snow + rhow*(qprec_snow(kk))* &
               (csice*(MIN(vmet(kk)%Ta,zero))-lambdaf)*dt(kk)
          vsnow(kk)%Qadv_rain = vsnow(kk)%Qadv_rain + rhow*(qprec(kk))*cswat*(MAX(vmet(kk)%Ta,zero))*dt(kk)

          ! update heat flux components
          Tqw  = MERGE(vmet(kk)%Ta, vsnow(kk)%tsn(vsnow(kk)%nsnow), -qevap(kk)>zero)

          !vsnow(kk)%Qadv_vap = vsnow(kk)%Qadv_vap + &
          !   rhow*(-qevapsig(kk))*cswat*Tqw*dt(kk) - &
          !    qadvsig(0)*dt(kk)

          vsnow(kk)%Qadv_vap =vsnow(kk)%Qadv_vap + (qadvsig(-vsnow(kk)%nsnow)-qadvsig(0))*dt(kk) - &
               rhow*(qprec_snow(kk))*(csice*(MIN(vmet(kk)%Ta,zero))-lambdaf)*dt(kk) - &
               rhow*(qprec(kk))*cswat*(vmet(kk)%Ta)*dt(kk)

          vsnow(kk)%Qcond_net = vsnow(kk)%Qcond_net + (qhsig(-vsnow(kk)%nsnow) -  qhsig(0))*dt(kk) - &
               (qadvsig(-vsnow(kk)%nsnow) -  qadvsig(0))*dt(kk)

          DO i=1, vsnow(kk)%nsnow
             IF (vsnow(kk)%hliq(i)>zero) THEN

                IF ((vsnow(kk)%hliq(i) + de(i-vsnow(kk)%nsnow))>(dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i))) THEN
                   ! account here for new hliq exceeding new hsnow
                   tmp1d1(kk) = vsnow(kk)%hliq(i)*cswat*rhow*zero + &
                        (vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow*(zero*csice-lambdaf)
                   tmp1d2(kk) = ((tmp1d1(kk) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk)))/ &
                        ((dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i))*cswat*rhow)
                   vsnow(kk)%tsn(i) = tmp1d2(kk)

                   IF (vsnow(kk)%tsn(i)>zero) THEN
                      WRITE(*,*) 'T>zero ', irec, nsteps(kk), i, vsnow(kk)%tsn(i)
                   ENDIF

                   vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) + &
                        lambdaf*(vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow
                   vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk) - &
                        lambdaf*(vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow
                   vsnow(kk)%hliq(i) = vsnow(kk)%hsnow(i) + dy(i-vsnow(kk)%nsnow)
                ELSEIF ((vsnow(kk)%hliq(i)+de(i-vsnow(kk)%nsnow))<zero) THEN
                   ! account here for new hliq less than zero
                   tmp1d1(kk) = (csice*(vsnow(kk)%tsn(i))-lambdaf)*rhow*(vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i)) + &
                        cswat*(vsnow(kk)%tsn(i))*rhow*vsnow(kk)%hliq(i) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk)
                   delta_snowT(i)  = (tmp1d1(kk)/rhow/(dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i))+lambdaf)/csice
                   vsnow(kk)%tsn(i) = delta_snowT(i)

                   vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) - lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + &
                        lambdaf*vsnow(kk)%hliq(i)*rhow
                   vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk) - &
                        (- lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + lambdaf*vsnow(kk)%hliq(i)*rhow)
                   vsnow(kk)%hliq(i) = zero
                ELSE
                   vsnow(kk)%hliq(i) = vsnow(kk)%hliq(i) + de(i-vsnow(kk)%nsnow)
                   delta_snowliq(i) = de(i-vsnow(kk)%nsnow)
                   delta_snowT(i) = zero
                   vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) - &
                        lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + lambdaf*de(i-vsnow(kk)%nsnow)*rhow
                   vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + &
                        dy(i-vsnow(kk)%nsnow)*rhow*csice*(vsnow(kk)%tsn(i)) + &
                        de(i-vsnow(kk)%nsnow)*rhow*(cswat-csice)*(vsnow(kk)%tsn(i))
                ENDIF
             ELSE !hliq = zero
                IF ((vsnow(kk)%tsn(i) + de(i-vsnow(kk)%nsnow))<zero) THEN
                   vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) - lambdaf*dy(i-vsnow(kk)%nsnow)*rhow
                   vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + &
                        dy(i-vsnow(kk)%nsnow)*rhow*csice*(vsnow(kk)%tsn(i)) + &
                        de(i-vsnow(kk)%nsnow)*csice*rhow*hsnow(i)
                   vsnow(kk)%tsn(i) = vsnow(kk)%tsn(i) + de(i-vsnow(kk)%nsnow)
                   delta_snowT(i) = de(i-vsnow(kk)%nsnow)
                   delta_snowliq(i) = zero
                   vsnow(kk)%hliq(i) = zero
                ELSE ! vsnow(kk)%tsn + de(0))>zero
                   delta_snowT(i) = zero - vsnow(kk)%tsn(i)
                   delta_snowliq(i) = (LHS_h(i-vsnow(kk)%nsnow)*dt(kk) - &
                        delta_snowT(i)*csice*hsnow(i) *rhow - &
                        dy(i-vsnow(kk)%nsnow)*rhow*(csice*(vsnow(kk)%tsn(i))-lambdaf))/ &
                        (rhow*(zero*(cswat-csice)+lambdaf))

                   vsnow(kk)%hliq(i) = vsnow(kk)%hliq(i) + delta_snowliq(i)
                   vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) &
                        - lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + &
                        lambdaf*rhow*delta_snowliq(i)
                   vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i)  &
                        + dy(i-vsnow(kk)%nsnow)*rhow*csice*(vsnow(kk)%tsn(i)) + &
                        delta_snowT(i)*csice*rhow*hsnow(i)
                   vsnow(kk)%tsn(i) = vsnow(kk)%tsn(i) + delta_snowT(i)
                ENDIF ! vsnow(kk)%tsn + de(0))>zero
             ENDIF ! hliq = zero
             vsnow(kk)%hsnow(i) = dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i)
             vsnow(kk)%Jsensible(i) = (vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow*csice*(vsnow(kk)%Tsn(i))
             vsnow(kk)%Jlatent(i) = (vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow*(-lambdaf)

             delta_snowcol(i) = dy(i-vsnow(kk)%nsnow)
          ENDDO ! loop over snow layers
          ! total internal energy and water content of snowpack
          vsnow(kk)%wcol = SUM(vsnow(kk)%hsnow(1:vsnow(kk)%nsnow))
          vsnow(kk)%deltawcol = vsnow(kk)%wcol - wcol0snow(kk)
          vsnow(kk)%Qprec = vsnow(kk)%Qprec + qprec_snow(kk)*dt(kk) + qprec(kk)*dt(kk)
          vsnow(kk)%Qevap = vsnow(kk)%Qevap + qevapsig(kk)*dt(kk)
          vsnow(kk)%Qvap = vsnow(kk)%Qvap + qsig(0)*dt(kk)

          vsnow(kk)%J = SUM(vsnow(kk)%Jsensible(1:vsnow(kk)%nsnow)+vsnow(kk)%Jlatent(1:vsnow(kk)%nsnow))

          vsnow(kk)%deltaJ = vsnow(kk)%J - J0snow(kk)

          IF (vsnow(kk)%nsnow.LT.nsnow_max) THEN
             vsnow(kk)%hsnow(vsnow(kk)%nsnow+1:nsnow_max) = zero
          ENDIF

          ! update densities and depths of snow layers

          ! effect of new snowfall

          ! snowfall density (tmp1d1), LaChapelle 1969
          IF (vmet(kk)%Ta > 2.0_r_2) THEN
             tmp1d1(kk) = 189.0_r_2
          ELSEIF ((vmet(kk)%Ta > -15.0_r_2) .AND. (vmet(kk)%Ta <= 2.0_r_2)) THEN
             tmp1d1(kk) = 50.0_r_2 + 1.7_r_2*(vmet(kk)%Ta+15.0_r_2)**1.5_r_2
          ELSE
             tmp1d1(kk) = 50.0_r_2
          ENDIF

          IF ((vsnow(kk)%hsnow(1)-qprec_snow(kk)*dt(kk)) > 1.e-3_r_2) THEN
             vsnow(kk)%dens(1) = (vsnow(kk)%dens(1)*(vsnow(kk)%hsnow(1)-qprec_snow(kk)*dt(kk)) &
                  + tmp1d1(kk)*qprec_snow(kk)*dt(kk))/vsnow(kk)%hsnow(1)
          ENDIF

          DO i=1, vsnow(kk)%nsnow
             vsnow(kk)%depth(i) = vsnow(kk)%hsnow(i)/(vsnow(kk)%dens(i)/rhow)

             ! adjust depth for compaction
             IF (vsnow(kk)%hliq(i) .GT. zero) THEN
                tmp1d1(kk) = MERGE(one, EXP(-0.06_r_2*(vsnow(kk)%dens(i)-150._r_2)), vsnow(kk)%dens(i)<150._r_2)
                tmp1d2(kk) = 2
             ELSE
                tmp1d1(kk) = MERGE(one, EXP(-0.046_r_2*(vsnow(kk)%dens(i)-150._r_2)), vsnow(kk)%dens(i)<150._r_2)
                tmp1d2(kk) = 1
             ENDIF
             tmp1d3(kk) =  vsnow(kk)%hsnow(i)*rhow/2._r_2    ! overburden kg/m2
             IF (i.GT.1) THEN
                tmp1d3(kk) = tmp1d3(kk) + SUM(vsnow(kk)%hsnow(1:i-1)*rhow)
             ENDIF

             IF (vsnow(kk)%depth(i)*dt(kk)* &
                  (2.8e-6_r_2*tmp1d1(kk)*tmp1d2(kk)*EXP(-vsnow(kk)%tsn(i)/25.0_r_2) + &
                  tmp1d3(kk)/3.6e6_r_2*EXP(-0.08_r_2*vsnow(kk)%tsn(i))*EXP(-0.023_r_2*vsnow(kk)%dens(i))) &
                  .LT.0.5_r_2*vsnow(kk)%depth(i)) THEN
                vsnow(kk)%depth(i) = vsnow(kk)%depth(i) - vsnow(kk)%depth(i)*dt(kk)* &
                     (2.8e-6_r_2*tmp1d1(kk)*tmp1d2(kk)*EXP(-vsnow(kk)%tsn(i)/25.0_r_2) + & ! metamorphism
                     tmp1d3(kk)/3.6e6_r_2*EXP(-0.08_r_2*vsnow(kk)%tsn(i))*EXP(-0.023_r_2*vsnow(kk)%dens(i))) ! overburden
             ENDIF
             vsnow(kk)%dens(i) = vsnow(kk)%hsnow(i)*rhow/vsnow(kk)%depth(i)
          ENDDO

       ENDIF ! end update snowlayers
       vsnow(kk)%MoistureFluxDivergence = vsnow(kk)%Qprec - vsnow(kk)%Qevap - vsnow(kk)%Qvap &
            - vsnow(kk)%Qmelt +  vsnow(kk)%Qtransfer
       vsnow(kk)%FluxDivergence = vsnow(kk)%Qadv_rain + vsnow(kk)%Qadv_snow + vsnow(kk)%Qadv_vap + &
            vsnow(kk)%Qadv_melt + vsnow(kk)%Qcond_net  +  vsnow(kk)%Qadv_transfer

    ENDIF ! .not. again

  END SUBROUTINE update_unknowns

  SUBROUTINE get_and_solve_eqn( &
       irec, mp, qprec,  n, dx, h0, S, thetai, &
       Tsoil,       &
       qh, nsteps,  vlit, vsnow, var,     &
       dxL,    &
       plit, par,        &
       deltaTa, &
       qevap, qL, qhL, qybL, qTbL, qhTbL, qhybL,    &
       again,  again_ice,  iok, itmp, ns,  &
       accel,  dt,   fac,   &
       phip,  rsig, rsigdt, sig, t,      &
       qexd, aa, bb, cc, dd, ee, ff, gg, dy, aah, bbh, cch, ddh, eeh, ffh, ggh, &
       de, q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb, qadv, qadvya, qadvyb, &
       qadvTa, qadvTb,  qsig, qhsig, qadvsig,      &
       dTsoil, tmp2d1,   &
       dthetaldT,   nsteps_ice, imelt,     &
       deltaJ_latent_T, deltaJ_sensible_S,   &
       qrunoff, tmp1d1, tmp1d2, tmp1d3,  tmp1d4,     &
       G0, &
       Tfreezing,   LHS, RHS, LHS_h, RHS_h,  nns, &
       iflux, litter, i,   kk, condition,   &
       advection,  c2, theta,     cp, &
       cpeff, hice,    h0_tmp,   hsnow, &
       delta_snowT, &
       delta_snowliq,    J0,   iqex,  &
       nfac1, nfac2, nfac3, nfac4, nfac5, nfac6, nfac7, nfac8, nfac9, &
       nfac10, nfac11, nfac12, err  &
       )
    IMPLICIT NONE
    INTEGER(i_d)                                           :: irec, mp
    REAL(r_2),      DIMENSION(1:mp)                        :: qprec
    INTEGER(i_d)                                           :: n
    REAL(r_2),      DIMENSION(1:n) :: dx
    REAL(r_2),      DIMENSION(1:mp)                        :: h0
    REAL(r_2),      DIMENSION(1:n) :: S
    REAL(r_2),      DIMENSION(1:n) :: thetai
    REAL(r_2),      DIMENSION(1:n) :: Tsoil
    REAL(r_2),      DIMENSION(-nsnow_max:n) :: qh
    INTEGER(i_d),   DIMENSION(1:mp)                        :: nsteps
    TYPE(vars),     DIMENSION(1:mp)                        :: vlit
    TYPE(vars_snow), DIMENSION(1:mp)                       :: vsnow
    TYPE(vars),     DIMENSION(1:n) :: var
    REAL(r_2),      DIMENSION(1:mp)                        :: dxL
    TYPE(params),   DIMENSION(1:mp)                        :: plit
    TYPE(params),   DIMENSION(1:n) :: par
    REAL(r_2),      DIMENSION(1:mp)                        :: deltaTa
    REAL(r_2),    DIMENSION(1:mp)                          :: qevap
    REAL(r_2),    DIMENSION(1:mp)                          :: qL, qhL, qybL, qTbL, qhTbL, qhybL
    LOGICAL,      DIMENSION(1:mp)                          :: again
    LOGICAL,      DIMENSION(1:mp,1:n)                      :: again_ice
    INTEGER(i_d), DIMENSION(1:mp)                          :: iok, itmp, ns
    REAL(r_2),    DIMENSION(1:mp)                          :: accel, dt, fac, phip
    REAL(r_2),    DIMENSION(1:mp)                          :: rsig, rsigdt, sig, t
    REAL(r_2),    DIMENSION(1:n) :: qexd
    REAL(r_2),    DIMENSION(-nsnow_max+1:n) :: aa, bb, cc, dd, ee, ff, gg, dy
    REAL(r_2),    DIMENSION(-nsnow_max+1:n) :: aah, bbh, cch, ddh, eeh, ffh, ggh, de
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qadv, qadvya, qadvyb, qadvTa, qadvTb
    REAL(r_2),    DIMENSION(-nsnow_max:n) :: qsig, qhsig, qadvsig
    REAL(r_2),    DIMENSION(1:n) :: dTsoil
    REAL(r_2),    DIMENSION(0:n) :: tmp2d1
    REAL(r_2),    DIMENSION(1:n) :: dthetaldT
    INTEGER(i_d), DIMENSION(1:n) :: nsteps_ice, imelt
    REAL(r_2),          DIMENSION(1:n) :: deltaJ_latent_T, deltaJ_sensible_S
    REAL(r_2),          DIMENSION(1:mp)                    :: qrunoff
    REAL(r_2),          DIMENSION(1:mp)                    :: tmp1d1, tmp1d2, tmp1d3,  tmp1d4
    REAL(r_2),          DIMENSION(1:mp)                    :: G0
    REAL(r_2),          DIMENSION(1:mp)                    :: Tfreezing
    REAL(r_2),          DIMENSION(-nsnow_max+1:n) :: LHS, RHS, LHS_h, RHS_h
    INTEGER(i_d),       DIMENSION(1:mp)                    :: nns, iflux
    LOGICAL                                                :: litter
    INTEGER(i_d)                                           :: i, kk, condition
    INTEGER(i_d)                                           :: advection ! switches
    REAL(r_2)                                              :: c2, theta
    REAL(r_2),          DIMENSION(1:mp)                    :: cp, cpeff, hice, h0_tmp
    REAL(r_2),          DIMENSION(nsnow_max) :: hsnow
    REAL(r_2),          DIMENSION(nsnow_max) :: delta_snowT, delta_snowliq
    REAL(r_2),      DIMENSION(1:n) :: J0
    REAL(r_2),          DIMENSION(1:n) :: iqex
    INTEGER(i_d),       DIMENSION(1:mp)                    :: nfac1, nfac2, nfac3, nfac4, nfac5, &
         nfac6, nfac7, nfac8, nfac9, nfac10, nfac11, nfac12
    ! Error flag if nstep of SLI > nsteps_max: err=0 -> no error; err/=0 -> error
    INTEGER(i_d),       INTENT(INOUT), OPTIONAL            :: err

    iok(kk)         = 0 ! flag for time step test
    itmp(kk)        = 0 ! counter to abort if not getting solution
    again_ice(kk,:) = .FALSE.
    nsteps_ice(1:n) = 0
    h0_tmp(kk) = h0(kk)
    DO WHILE (iok(kk)==0) ! keep reducing time step until all ok
       itmp(kk)  = itmp(kk) + 1
       accel(kk) = one - 0.05_r_2*REAL(MIN(10,MAX(0,itmp(kk)-4)),r_2) ! acceleration [0.5,1], start with 1
       IF (itmp(kk) > 1000) THEN
          WRITE(logn,*) "Solve: too many iterations of equation solution"
          WRITE(logn,*) " irec, kk, S"
          WRITE(logn,*) irec, kk, S(:)
          WRITE(logn,*) " irec, kk, Tsoil"
          WRITE(logn,*) irec, kk, Tsoil(:)
          WRITE(logn,*) " irec, kk, qex"
          WRITE(logn,*) irec, kk, iqex(:)
          WRITE(logn,*) " irec, kk, h0, hsnow, hsnowliq"
          WRITE(logn,*) irec, kk, h0(kk), vsnow(kk)%hsnow, vsnow(kk)%hliq
          WRITE(logn,*) nfac1(kk), nfac2(kk), nfac3(kk), nfac4(kk), nfac5(kk), &
               nfac6(kk), nfac7(kk), nfac8(kk), nfac9(kk), nfac10(kk), nfac11(kk), nfac12(kk)
          err = 1
          RETURN
       END IF

       ! prelim estimate of new top snow layer depth for use in energy cons eq'n             !
       IF ((vsnow(kk)%nsnow>0))  THEN
          hsnow(1:vsnow(kk)%nsnow) = vsnow(kk)%hsnow(1:vsnow(kk)%nsnow) + &
               (-q(1-vsnow(kk)%nsnow:0) + q(-vsnow(kk)%nsnow:-1))*dt(kk)
       ENDIF

       aa(1:n)   =  qya(0:n-1)
       ee(0:n-1) = -qyb(0:n-1)
       bb(1:n)   =  qTa(0:n-1)
       ff(0:n-1) = -qTb(0:n-1)
       gg(1:n) = -(q(0:n-1)-q(1:n)-iqex(1:n))*rsig(kk)
       gg(1) = gg(1)+qrunoff(kk)*rsig(kk)

       aah(1:n)   =  qhya(0:n-1)
       eeh(0:n-1) = -qhyb(0:n-1)
       bbh(1:n)   =  qhTa(0:n-1)
       ffh(0:n-1) = -qhTb(0:n-1)
       ggh(1:n) =  -(qh(0:n-1)-qh(1:n))*rsig(kk)

       IF (advection==1) THEN
          ggh(1:n) = ggh(1:n) + iqex(1:n)*rsig(kk)*(Tsoil(1:n))*cswat*rhow
          ggh(1) = ggh(1) + qrunoff(kk)*rsig(kk)*(Tsoil(1))*cswat*rhow
       ENDIF

       IF (litter) THEN ! full litter model: litter in zeroth layer
          ! only use zeroth layer for litter (pond included in layer 1)
          WRITE(*,*) 'Should not be here - QYBL 01 ', qybl(kk)
          cc(0) = -qya(0) - rsigdt(kk)*plit(kk)%thre*dxL(kk) + qybL(kk)
          gg(0)  = -(qprec(kk)-qevap(kk)-q(0))*rsig(kk)
          ggh(0) = -(G0(kk)-qh(0))*rsig(kk)
          dd(0)  = -qTa(0)
       ENDIF
       aa(0)  = qya(0)
       aah(0) = zero
       bbh(0) = zero

       WHERE (var(1:n)%isat==0) ! unsaturated layers
          cc(1:n) = qyb(0:n-1) - qya(1:n) - par(1:n)%thre*dx(1:n)*rsigdt(kk) - qexd(1:n)
       ELSEWHERE ! saturated layers
          cc(1:n) = qyb(0:n-1) - qya(1:n) - qexd(1:n)
       endwhere

       IF (ns(kk)<1) THEN ! pond included in top soil layer, solving for change in pond height
          cc(1) = -qya(1)-rsigdt(kk) -qexd(1)
       ENDIF

       cch(1:n) = qhyb(0:n-1)-qhya(1:n) +   &
            REAL(var(1:n)%iice,r_2)*REAL(1-var(1:n)%isat,r_2)*rhow*lambdaf*par(1:n)%thre*dx(1:n)*rsigdt(kk)

       IF (ns(kk)==0) THEN ! change in pond height (top layer saturated)
          cch(1) = cch(1) +  &
               REAL(var(1)%iice,r_2)*rhow*lambdaf*rsigdt(kk)*(var(1)%thetai/par(1)%thre)
       ENDIF

       ! modification to cch for advection
       IF (advection==1) THEN
          IF (ns(kk)==0) THEN  ! changing pond included in top soil layer
             cch(1) = cch(1) -rhow*cswat*(Tsoil(1))*rsigdt(kk)*REAL(1-var(1)%iice,r_2) &
                  -rhow*csice*(Tsoil(1))*rsigdt(kk)*REAL(var(1)%iice,r_2) &
                  *(var(1)%thetai/par(1)%thre) &
                  -rhow*cswat*(Tsoil(1))*rsigdt(kk)*REAL(var(1)%iice,r_2)* &
                  (one-(var(1)%thetai/par(1)%thre))
          ELSE ! no pond
             cch(1) = cch(1) -rhow*(Tsoil(1))*par(1)%thre*dx(1)*rsigdt(kk)* &
                  REAL(1-var(1)%isat,r_2)*(cswat*REAL(1-var(1)%iice,r_2) &
                  +csice*REAL(var(1)%iice,r_2))
          ENDIF

          cch(2:n) = cch(2:n) -rhow*(Tsoil(2:n))*par(2:n)%thre*dx(2:n)*rsigdt(kk)* &
               REAL(1-var(2:n)%isat,r_2)*(cswat*REAL(1-var(2:n)%iice,r_2) &
               +csice*REAL(var(2:n)%iice,r_2))
       ENDIF

       dd(1:n)  = qTb(0:n-1)-qTa(1:n)

       ddh(1:n) = qhTb(0:n-1)-qhTa(1:n) - &
                                ! Only apply latent heat component of heat capacity to total deltaT if soil remains frozen
            var(1:n)%csoileff*dx(1:n)*rsigdt(kk)- &
            (cswat-csice)*dx(1:n)*var(1:n)%dthetaldt*rhow*(Tsoil(1:n))* &
            rsigdt(kk)*REAL(var(1:n)%iice,r_2)
       ! modification to ddh(1) for pond
       ddh(1) = ddh(1) - cpeff(kk)*h0_tmp(kk)*rsigdt(kk) - &
            (cswat-csice)*h0(kk)/par(1)%thre*var(1)%dthetaldt*rhow*(Tsoil(1))* &
            rsigdt(kk)*REAL(var(1)%iice,r_2)

       IF (advection==1) THEN
          ddh(1:n) = ddh(1:n) - iqex(1:n)*cswat*rhow
          ddh(1) = ddh(1) - qrunoff(kk)*cswat*rhow
       ENDIF

       ! modification of matrix to incorporate single snow layer
       IF (vsnow(kk)%nsnow==1) THEN
          cc(0) = qyb(-1)-qya(0)-rsigdt(kk)
          dd(0) =  qTb(-1)-qTa(0)
          ee(0) =   -qyb(0)
          ff(0) =   -qTb(0)
          gg(0) = -(q(-1)-q(0))*rsig(kk)
          IF (vsnow(kk)%hliq(1)>zero) THEN ! liquid phase present, solve for change in liq content
             ddh(0) =  - rhow*((vsnow(kk)%tsn(1))*(cswat-csice)+lambdaf)*rsigdt(kk)
          ELSE ! solid phase only; solve for change in snow t
             ddh(0) = qhTb(-1) - qhTa(0) - rhow*csice*hsnow(1)*rsigdt(kk)
          ENDIF
          cch(0) = qhyb(-1)-qhya(0) - rhow*(csice*(vsnow(kk)%tsn(1))-lambdaf)*rsigdt(kk)
          eeh(0) = -qhyb(0)
          ffh(0) = -qhTb(0)
          ggh(0) = -(qh(-1)-qh(0))*rsig(kk)
       ENDIF
       ! modification of matrix to incorporate more than one snow layer
       IF (vsnow(kk)%nsnow>1) THEN
          aa(1-vsnow(kk)%nsnow:0)   =  qya(-vsnow(kk)%nsnow:-1)
          bb(1-vsnow(kk)%nsnow:0)   =  qTa(-vsnow(kk)%nsnow:-1)
          cc(1-vsnow(kk)%nsnow:0) = qyb(-vsnow(kk)%nsnow:-1)-qya(1-vsnow(kk)%nsnow:0)-rsigdt(kk)
          dd(1-vsnow(kk)%nsnow:0) =  qTb(-vsnow(kk)%nsnow:-1)-qTa(1-vsnow(kk)%nsnow:0)
          ee(1-vsnow(kk)%nsnow:0) =   -qyb(1-vsnow(kk)%nsnow:0)
          ff(1-vsnow(kk)%nsnow:0) =   -qTb(1-vsnow(kk)%nsnow:0)
          gg(1-vsnow(kk)%nsnow:0) = -(q(-vsnow(kk)%nsnow:-1)-q(1-vsnow(kk)%nsnow:0))*rsig(kk)
          ddh(1-vsnow(kk)%nsnow:0) = MERGE( & ! liquid phase present, solve for change in liq content
               -rhow*((vsnow(kk)%tsn(1:vsnow(kk)%nsnow))*(cswat-csice)+lambdaf)*rsigdt(kk), &
               qhTb(-vsnow(kk)%nsnow:-1) - qhTa(1-vsnow(kk)%nsnow:0) - &
               rhow*csice*hsnow(1:vsnow(kk)%nsnow)*rsigdt(kk), & ! solid phase only; solve for change in snow t
               vsnow(kk)%hliq(1:vsnow(kk)%nsnow) > zero)
          aah(1-vsnow(kk)%nsnow:0)   =  qhya(-vsnow(kk)%nsnow:-1)
          bbh(1-vsnow(kk)%nsnow:0)   =  qhTa(-vsnow(kk)%nsnow:-1)

          cch(1-vsnow(kk)%nsnow:0) = qhyb(-vsnow(kk)%nsnow:-1)-qhya(1-vsnow(kk)%nsnow:0) - &
               rhow*(csice*(vsnow(kk)%tsn(1:vsnow(kk)%nsnow))-lambdaf)*rsigdt(kk)
          eeh(1-vsnow(kk)%nsnow:0) = -qhyb(1-vsnow(kk)%nsnow:0)
          ffh(1-vsnow(kk)%nsnow:0) = -qhTb(1-vsnow(kk)%nsnow:0)
          ggh(1-vsnow(kk)%nsnow:0) = -(qh(-vsnow(kk)%nsnow:-1)-qh(1-vsnow(kk)%nsnow:0))*rsig(kk)
       ENDIF
       IF (litter .AND. ns(kk)==1) THEN ! litter and no pond
          ! watch for deltaTa
          WRITE(*,*) 'Should not be here - QYBL 02 ', qybl(kk)
          cc(0)  = qybL(kk) -qya(0) -rsigdt(kk)*plit(kk)%thre*dxL(kk)
          dd(0)  = qTbL(kk) - qTa(0)
          ee(0)  = -qyb(0)
          ff(0)  = -qTb(0)
          gg(0)  = (q(0) - qL(kk))*rsig(kk) + deltaTa(kk)*(qTbL(kk)-qTa(0))
          gg(1)  = -(q(0)-q(1)-iqex(1))*rsig(kk) + deltaTa(kk)*qTa(0)
          cch(0) = qhybL(kk) - qhya(0)
          ddh(0) = -qhTa(0)-vlit(kk)%csoil*dxL(kk)*rsigdt(kk)+qhTbL(kk)
          eeh(0) = -qhyb(0)
          ffh(0) = -qhTb(0)
          ggh(0) = (qh(0)-qhL(kk))*rsig(kk) + deltaTa(kk)*(qhTbL(kk) - qhTa(0))
          ggh(1) = -(qh(0)-qh(1))*rsig(kk) + deltaTa(kk)*qhTa(0)
       ENDIF

       ! litter and pond !vh! need to check this now that pond is lumped with top soil layer
       IF (litter .AND. ns(kk)==0) THEN
          ddh(0) = -qhTa(0)-cswat*rhow*h0(kk)*rsigdt(kk) -vlit(kk)%csoil*dxL(kk)*rsigdt(kk)
       ENDIF

       nns(kk) = 1  ! pond included in top soil layer
       IF (vsnow(kk)%nsnow>0) THEN
          nns(kk) = 1-vsnow(kk)%nsnow
       ENDIF

       CALL massman_sparse(aa(nns(kk)+1:n), aah(nns(kk)+1:n), bb(nns(kk)+1:n), bbh(nns(kk)+1:n), &
            cc(nns(kk):n), cch(nns(kk):n), dd(nns(kk):n), ddh(nns(kk):n), ee(nns(kk):n-1), &
            eeh(nns(kk):n-1), &
            ff(nns(kk):n-1), ffh(nns(kk):n-1), gg(nns(kk):n), ggh(nns(kk):n), &
            dy(nns(kk):n), de(nns(kk):n), condition=condition, err=err)
       IF (err /= 0) THEN
          WRITE(logn,*) "Sparse matrix solution failed ", irec, kk
          WRITE(logn,*) Tsoil(1), S(1)
          RETURN
       ENDIF

       dTsoil(1:n) = de(1:n)


       WHERE (vsnow(kk)%hliq(:)>zero)
          delta_snowT(:) = zero
       ELSEWHERE
          delta_snowT(:) = de(LBOUND(de,1):0)
       endwhere

       ! evaluate soil fluxes at sigma of time step
       qsig(0)  = q(0)+sig(kk)*qyb(0)*dy(1) + sig(kk)*qya(0)*dy(0) + &
            sig(kk)*qTb(0)*dTsoil(1) + sig(kk)*qTa(0)*de(0)

       qhsig(0) = qh(0) + sig(kk)*qhyb(0)*dy(1) + sig(kk)*qhya(0)*dy(0) + &
            sig(kk)*qhTb(0)*dTsoil(1) + sig(kk)*qhTa(0)*de(0)

       qadvsig(0) = qadv(0) + sig(kk)*qadvyb(0)*dy(1) &
            + sig(kk)*qadvya(0)*dy(0) + sig(kk)*qadvTb(0)*dTsoil(1) + sig(kk)*qadvTa(0)*de(0)

       qsig(1:n-1) = q(1:n-1) + sig(kk)*(qya(1:n-1)*dy(1:n-1)+qyb(1:n-1)*dy(2:n) &
            +qTa(1:n-1)*dTsoil(1:n-1)+qTb(1:n-1)*dTsoil(2:n))
       qsig(n)     = q(n) + sig(kk)*(qya(n)*dy(n)+qTa(n)*dTsoil(n))

       qhsig(1:n-1) = qh(1:n-1) + sig(kk)*(qhya(1:n-1)*dy(1:n-1)+qhyb(1:n-1)*dy(2:n) &
            +qhTa(1:n-1)*dTsoil(1:n-1)+qhTb(1:n-1)*dTsoil(2:n))
       qhsig(n) = qh(n) + sig(kk)*(qhya(n)*dy(n)+qhTa(n)*dTsoil(n))

       qadvsig(1:n-1) = qadv(1:n-1) + sig(kk)*(qadvya(1:n-1)*dy(1:n-1)+qadvyb(1:n-1)*dy(2:n) &
            +qadvTa(1:n-1)*dTsoil(1:n-1)+qadvTb(1:n-1)*dTsoil(2:n))
       qadvsig(n) = qadv(n) + sig(kk)*(qadvya(n)*dy(n)+qadvTa(n)*dTsoil(n))

       LHS_h(1) = (dx(1)*var(1)%csoileff+h0_tmp(kk)*cpeff(kk))*dTsoil(1)/dt(kk) - &
            dy(1)*par(1)%thre*dx(1)*rhow*lambdaf/dt(kk)*var(1)%iice* &
            REAL(1-var(1)%isat,r_2) - &
            dy(1)*(var(1)%thetai/par(1)%thre)*rhow*lambdaf/dt(kk)*REAL(1-ns(kk),r_2)* &
            var(1)%iice  + &
            var(1)%dthetaldT*(Tsoil(1))*rhow*(cswat-csice)*dx(1)*dTsoil(1)/dt(kk)* &
            var(1)%iice + &
            var(1)%dthetaldT*(Tsoil(1))*rhow*(cswat-csice)*h0(kk)/par(1)%thre* &
            dTsoil(1)/dt(kk)*var(1)%iice

       LHS_h(2:n) = (dx(2:n)*var(2:n)%csoileff)*dTsoil(2:n)/dt(kk)  - &
            dy(2:n)*par(2:n)%thre*dx(2:n)*rhow*lambdaf/dt(kk)*var(2:n)%iice* &
            REAL(1-var(2:n)%isat,r_2) + &
            var(2:n)%dthetaldT*(Tsoil(2:n))*rhow*(cswat-csice)*dx(2:n)* &
            dTsoil(2:n)/dt(kk)*var(2:n)%iice

       RHS_h(1:n) = qhsig(0:n-1) -qhsig(1:n)

       IF (advection==1) THEN
          LHS_h(1:n) = LHS_h(1:n) &
               + REAL(1-var(1:n)%iice,r_2) * REAL(1-var(1:n)%isat,r_2) * &
               dx(1:n) * dy(1:n)/dt(kk) * par(1:n)%thre * (Tsoil(1:n)) * rhow * cswat &
               + REAL(var(1:n)%iice,r_2) * REAL(1-var(1:n)%isat,r_2) * &
               dx(1:n) * dy(1:n)/dt(kk) * par(1:n)%thre * rhow * csice * (Tsoil(1:n))

          LHS_h(1) = LHS_h(1) &
               + REAL(1-ns(kk),r_2) * REAL(1-var(1)%iice,r_2) * &
               dy(1)/dt(kk) * rhow * cswat * (Tsoil(1)) &
               + REAL(1-ns(kk),r_2) * REAL(var(1)%iice,r_2)* &
               dy(1)/dt(kk) * rhow * csice * (Tsoil(1)) * (var(1)%thetai/par(1)%thre) &
               + REAL(1-ns(kk),r_2) * REAL(var(1)%iice,r_2)  * &
               dy(1)/dt(kk) * rhow * cswat * (Tsoil(1)) * (one-(var(1)%thetai/par(1)%thre))

          RHS_h(1:n) = RHS_h(1:n) - iqex(1:n)*cswat*rhow*(