MODULE sli_main_mod

CONTAINS

  SUBROUTINE sli_main(ktau, dt, veg, soil, ssnow, met, canopy, air, rad, SEB_only)

    ! Main subroutine for Soil-litter-iso soil model
    ! Vanessa Haverd, CSIRO Marine and Atmospheric Research
    ! and Matthias Cuntz, UFZ - Helmholtz Centre for Environmental Research, 2010
    ! Modified to operate for multiple veg tiles but a single soil column, March 2011
    ! Rewritten for same number of soil columns as veg tiles May 2012
    USE cable_def_types_mod,       ONLY: veg_parameter_type, soil_parameter_type, soil_snow_type, met_type, &
         canopy_type, air_type, radiation_type, ms, mp, r_2, i_d
    USE cable_common_module , ONLY: cable_user
    USE sli_numbers,        ONLY:  zero, half, one, two, four, thousand, & ! numbers
         Tzero, experiment, &                                       ! variables
         vars_met, vars, params, vars_snow, &                                  ! types
         MW, snmin, Rgas, Lambdas, lambdaf, csice, cswat, rhow, nsnow_max, e5, &
         freezefac, topmodel, alpha, fsat_max, botbc
    USE sli_utils,          ONLY: x, dx, par, setpar, setpar_Loetsch, setx, plit, dxL, setlitterpar, esat, &
         esat_ice, slope_esat_ice, thetalmax, Tfrz,  hyofS, SEB
    USE sli_roots,          ONLY: setroots, getrex
    USE sli_solve,          ONLY: solve

    USE  cable_IO_vars_module, ONLY: verbose

    IMPLICIT NONE
    REAL,                      INTENT(IN)    :: dt
    TYPE(veg_parameter_type),  INTENT(INOUT) :: veg     ! all r_1
    TYPE(soil_parameter_type), INTENT(INOUT) :: soil    ! all r_1
    TYPE(soil_snow_type),      INTENT(INOUT) :: ssnow   ! r_1, r_2 desaster
    TYPE(met_type),            INTENT(INOUT) :: met     ! all r_1
    TYPE(canopy_type),         INTENT(INOUT) :: canopy  ! all r_1
    TYPE(air_type),            INTENT(INOUT) :: air     ! all r_1
    TYPE (radiation_type),     INTENT(IN)    :: rad
    INTEGER,                   INTENT(IN)    :: ktau ! integration step number
    INTEGER,                   INTENT(IN)    :: SEB_only ! integration step number

    REAL(r_2), PARAMETER :: emsoil=0.97
    REAL(r_2), PARAMETER :: rhocp=1.1822e3
    REAL(r_2), PARAMETER :: Dva = 2.17e-5
    INTEGER(i_d) :: i, k, kk, setroot
    REAL(r_2)    :: ti, tf
    TYPE(vars_met),  DIMENSION(1:mp)      :: vmet ! Meteorology above soil
    TYPE(vars),      DIMENSION(1:mp)      :: vlit
    TYPE(vars),      DIMENSION(1:mp,1:ms) :: var
    TYPE(vars_snow), DIMENSION(1:mp)      :: vsnow
    INTEGER(i_d),    DIMENSION(1:mp)      :: nsteps
    REAL(r_2),       DIMENSION(1:mp,1:ms) :: Tsoil, S, thetai, Jsensible
    REAL(r_2),       DIMENSION(1:mp)      :: SL, TL, T0
    REAL(r_2),       DIMENSION(1:mp)      :: drn, evap, infil, qprec, qprec_snow,qprec_snow_tmp, qprec_tot, runoff, runoff_sat
    REAL(r_2),       DIMENSION(1:mp)      :: win, wp, wpi, h0, deltah0, h0old,hsnowold, discharge
    REAL(r_2),       DIMENSION(1:mp)      :: ip, ipi ! volumetric ice content of profile (final and initial)
    REAL(r_2),       DIMENSION(1:mp,1:ms) :: wex, csoil, qex, kth, phi, thetal_max, Sliq, Ksat
    REAL(r_2),       DIMENSION(1:mp,1:ms) :: FS
    ! surface temperature (top of top soil layer or top of litter layer)
    REAL(r_2),       DIMENSION(1:mp)      :: Tsurface
    REAL(r_2),       DIMENSION(1:mp)      :: gr, grc
    REAL(r_2),       DIMENSION(1:mp)      :: Etrans
    REAL(r_2),       DIMENSION(1:mp)      :: gamm
    REAL(r_2),       DIMENSION(1:mp)      :: G0, H, lE, Epot
    REAL(r_2),       DIMENSION(1:mp,-nsnow_max:ms) :: qh, qvsig, qlsig, qvTsig, qvh
    REAL(r_2),       DIMENSION(1:mp)      :: deltaTa, lE_old !, SA, SB, wpAi, wpBi, wpA, wpB
    REAL(r_2),       DIMENSION(1:mp)      :: evap_pot, deltaice_cum_T, deltaice_cum_S, zdelta
    REAL(r_2),       DIMENSION(1:mp)      :: fws
    REAL(r_2),       DIMENSION(1:mp)      :: Qadvcum, Jcol_sensible, Jcol_latent_S, Jcol_latent_T
    REAL(r_2),       DIMENSION(1:mp)      :: tmp1d1, deltaEsnow
    REAL(r_2),       DIMENSION(1:mp)      :: hice, phie
    REAL(r_2)                             :: tmp1d1a, tmp1d2, tmp1d3, tmp1d4, &
         tmp1d5, tmp1d6, tmp1d7, tmp1d8, tmp1d9,tmp1d10, tmp1d11, &
         tmp1d12,tmp1d13, tmp1d14, tmp1d15, tmp1d16
    REAL(r_2) :: rbw, rbh, rrc ! resistances for output
    INTEGER(i_d), DIMENSION(1:mp) :: index

    ! Topmodel
    REAL(r_2), DIMENSION(1:mp) :: fsat ! topmodel saturated area
    REAL(r_2), DIMENSION(1:mp) :: qb   ! topmodel baseflow

    ! Model switches
    INTEGER(i_d), PARAMETER :: litter       = 2 ! which litter model
    ! 0: no litter
    ! 1: full litter
    ! 2: litter resistance
    INTEGER(i_d), PARAMETER :: advection    = 1 ! heat advection by water
    INTEGER(i_d), PARAMETER :: isotopologue = 0 ! which isotope
    ! 0: no isotope calculations
    ! 1: HDO
    ! 2: H218O
    ! 3: HDO & H218O
    ! 0: normal run
    INTEGER(i_d), PARAMETER :: septs        = 0 ! coupled or uncoupled energy and water calculation
    ! 0: coupled calc
    ! 1: uncoupled energy (T) and moisture (S)
    INTEGER(i_d), PARAMETER :: condition    = 3 ! condition matrix before solving
    ! 0: no conditioning
    ! 1: condition columns
    ! 2: condition lines
    ! 3: condition first lines then columns
    LOGICAL, SAVE :: first = .TRUE.
    INTEGER(i_d), SAVE  :: counter

    ! Error flag if nstep of SLI > nsteps_max: err=0 -> no error; err/=0 -> error
    INTEGER(i_d), DIMENSION(1:mp) :: err



    REAL(r_2), DIMENSION(ms) :: zmm,dzmm
    REAL(r_2), DIMENSION(mp) :: zaq

    ! initialise cumulative variables
    ! Jcol_sensible = zero
    ! Jcol_latent_S = zero
    ! Jcol_latent_T = zero
    ! deltaice_cum_T = zero
    ! deltaice_cum_S = zero
    ! Jsensible = zero
    drn    = zero
    ! discharge = zero
    infil  = zero
    evap   = zero
    ! evap_pot  = zero
    runoff = zero
    ! qh = zerocanopy%fevc(i) = ecx(i)*(1.0-canopy%fwet(i))
    ! H      = zero
    ! G0      = zero
    ! lE     = zero
    ! csoil = zero
    ! kth = zero
    ! Qadvcum  = zero
    wex    = zero
    ! qlsig        = zero
    ! qvsig        = zero
    ! qvh        = zero
    ! qvtsig        = zero
    thetal_max = zero
    Sliq       = zero
    Ksat       = zero
    phie       = zero
    phi        = zero
    hice       = zero
    err        = 0

    ! output files for testing purposes
    IF (first .AND. verbose) THEN
       OPEN (unit=332,file="vh08.out",status="replace",position="rewind")
       OPEN (unit=334,file="S.out",status="replace",position="rewind")
       OPEN (unit=336,file="Tsoil.out",status="replace",position="rewind")
       OPEN (unit=335,file="SEB.out",status="replace",position="rewind")
       OPEN (unit=337,file="soil_log.out",status="replace",position="rewind")
       OPEN(unit=338, file="thetai.out", status="replace", position="rewind")
       OPEN(unit=340, file="snow.out", status="replace", position="rewind")
       OPEN(unit=345, file="diags.out",status="replace", position="rewind")
       OPEN(unit=369, file="vmet.out", status="replace", position="rewind", recl=20*20)
       OPEN(unit=370, file="qex.out",status="replace", position="rewind")
       OPEN(unit=371, file="q.out",status="replace", position="rewind")

       !open(unit=339, file="latlong.out",status="replace", position="rewind")
       ! write(339,"(20000f8.2)") rad%latitude
       !write(339,"(20000f8.2)") rad%longitude
       counter = 0
    ENDIF

    counter = counter + 1

    ! Save soil / snow surface temperature from last time step:
    ssnow%otss = ssnow%tss

    ! set layer thicknesses
    IF (.NOT. ALLOCATED(x)) CALL setx(mp, ms, soil)

    ! Set root density distribution (leave in for sli offline)
    setroot = 0  ! reset rooting depths
    IF (setroot == 1) THEN
       CALL setroots(x*100.0_r_2, REAL(veg%F10,r_2), REAL(veg%ZR,r_2)*100.0_r_2, FS)
       veg%froot = FS
    ELSE
       FS = REAL(veg%froot,r_2)
    ENDIF
    ! set required soil hydraulic params
    IF (.NOT. ALLOCATED(par)) THEN
       index=(/(i,i=1,mp,1)/)
       IF (experiment == 16) THEN
          CALL setpar_Loetsch(mp, ms, x-half*dx)
       ELSE
          CALL setpar(mp, ms, soil, index)
       ENDIF
    ENDIF

    ! If we want solutes:
!     if (.not. allocated(bd)) allocate(bd(soil%nhorizons(k)))

    ! Litter parameters:
    IF (.NOT. ALLOCATED(plit)) THEN
       index=(/(i,i=1,mp,1)/)
       CALL setlitterpar(mp, veg, index)
    ENDIF

    ! Met data above soil:
    IF (first) THEN
       vmet%rha   = zero
       vmet%rrc   = zero
       vmet%Rn    = zero
       vmet%Rnsw = zero
       vmet%cva   = zero
       vmet%civa  = zero
       vmet%phiva = zero
    ENDIF

    vmet%Ta  = REAL(met%Tvair,r_2) - Tzero
    vmet%Da  = REAL(met%dva,r_2)

    IF (cable_user%or_evap .AND. SEB_only==1) THEN
       vmet%rbh = ssnow%rtsoil +  (one - &
            ssnow%satfrac(:))*ssnow%rtevap_unsat(:) + ssnow%satfrac(:)*ssnow%rtevap_sat(:)
       vmet%rbw = vmet%rbh + canopy%sublayer_dz(:)/(0.27_r_2/1189.8)
    ELSE
       vmet%rbh = ssnow%rtsoil
       vmet%rbw = vmet%rbh
    END IF

    gr       = four * emsoil * (vmet%Ta+Tzero)**3 *5.67e-8_r_2 ! radiation conductance Wm-2K-1
    grc      = one/vmet%rbh   + gr/rhocp
    vmet%rrc = one/grc                            ! resistance to radiative and convective heat transfer

    vmet%rha   = MAX(MIN((esat(vmet%Ta)-vmet%Da)/esat(vmet%Ta),one),0.1_r_2)
    vmet%cva   = vmet%rha * esat(vmet%Ta)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) ! m3 H2O (liq) m-3 (air)
    vmet%phiva = Dva * vmet%cva
    vmet%Rn    = canopy%fns
    ! vmet%Rnsw  = rad%qssabs  ! shortwave radiation absorbed
    vmet%Rnsw = zero ! all radiation absorbed at snow surface
    ! vmet%Rnsw = vmet%Rn ! all radiation absorbed beneath snow surface
    Etrans     = MAX(canopy%fevc/air%rlam/thousand, zero) ! m s-1
    WHERE (canopy%fevc .LT. zero)
       canopy%fevw = canopy%fevw+canopy%fevc
       canopy%fevc = zero
    END WHERE
    h0         = ssnow%h0

    ! zero runoff here, in case error is returned to avoid excessive runoff from previous time-step. (Runoff is multipled by dt in cable_serial.F90)
    ssnow%rnof1 = 0.0
    ssnow%rnof2 = 0.0
    ssnow%runoff = 0.0
    ! Set isotopes to zero
    vmet%civa = zero

    ! Initialisations:
    IF (first) THEN
       DO kk=1, mp
          vsnow(kk)%hsnow      = zero
          vsnow(kk)%depth      = zero
          vsnow(kk)%totdepth   = zero
          vsnow(kk)%wcol       = zero
          vsnow(kk)%dens       = 120._r_2 ! snow density kg m-3
          vsnow(kk)%tsn        = zero
          vsnow(kk)%kH         = 0.16_r_2 ! snow thermal cond (W m-2 K-1)
          vsnow(kk)%Dv         = Dva      ! m2 s-1
          vsnow(kk)%sl         = zero
          vsnow(kk)%kE         = zero
          vsnow(kk)%kth        = vsnow(kk)%kH
          vsnow(kk)%cv         = zero
          vsnow(kk)%hliq       = zero
          vsnow(kk)%melt       = zero
          vsnow(kk)%nsnow      = 0
          vsnow(kk)%nsnow_last = 0
          ssnow%cls(kk)        = one
          vsnow(kk)%J                      = zero
          vsnow(kk)%fsnowliq_max           = 0.03_r_2
          vsnow(kk)%deltaJlatent           = zero
          vsnow(kk)%deltaJsensible         = zero
          vsnow(kk)%Qadv_snow              = zero
          vsnow(kk)%Qadv_rain              = zero
          vsnow(kk)%Qadv_melt              = zero
          vsnow(kk)%Qadv_vap               = zero
          vsnow(kk)%Qcond_net              = zero
          vsnow(kk)%Qadv_transfer          = zero
          vsnow(kk)%Qmelt                  = zero
          vsnow(kk)%Qtransfer              = zero
          vsnow(kk)%FluxDivergence         = zero
          vsnow(kk)%deltaJ                 = zero
          vsnow(kk)%Qvap                   = zero
          vsnow(kk)%MoistureFluxDivergence = zero
          vsnow(kk)%Qprec                  = zero
          vsnow(kk)%Qevap                  = zero
          vsnow(kk)%deltawcol              = zero
          vsnow(kk)%Jsensible              = zero
          vsnow(kk)%Jlatent                = zero
       ENDDO
       first = .FALSE.
    ENDIF


    Tsurface = ssnow%Tsurface
    T0       = ssnow%Tsurface
    deltaTa  = zero
    lE_old   = ssnow%lE
    zdelta   = ssnow%zdelta



    SL    = 0.5_r_2   ! degree of litter saturation
    Tsoil = ssnow%Tsoil
    TL(:) = Tsoil(:,1) ! litter T

    !thetai      = ssnow%thetai
    !var%thetai = thetai
    ssnow%smelt = zero
    ssnow%cls   = one

    !mrd561 UM changes
    thetai      =ssnow%wbice! ssnow%thetai
    var%thetai = thetai
    DO k=1,ms
       DO i=1,mp
          ssnow%S(i,k) = ssnow%wb(i,k)/(par(i,k)%thr+(par(i,k)%the-par(i,k)%thr))
       END DO
    END DO
    S           = ssnow%S                ! degree of soil saturation


    ! ----------------------------------------------------------------
    ! Iinitialise phi where it is (frozen and saturated) and where (pond >zero)

    ! Test steep freezing curve
    WHERE (Tsoil<Tfrz(S, par%he, one/(par%lambc*freezefac)))
       par%lam = par%lambc * freezefac
    ELSEWHERE
       par%lam = par%lambc
    END WHERE
    par%eta = two/par%lam + two + one

    ! where (Tsoil<Tfrz(S,par%he,one/par%lam))
    WHERE (Tsoil<Tfrz(S, par%he, one/(par%lambc*freezefac)))
       thetal_max = thetalmax(Tsoil,S,par%he,one/par%lam,par%thre,par%the)
       Sliq       = (thetal_max - (par%the-par%thre))/par%thre
       Ksat       = par%Ke*EXP(par%eta*LOG(Sliq))
    ELSEWHERE
       Sliq       = S
       Ksat       = par%Ke
    endwhere

    ! where ((Tsoil<Tfrz(S,par%he,one/par%lam)) .and. (S>=1))
    WHERE ((Tsoil<Tfrz(S, par%he, one/(par%lambc*freezefac))) .AND. (S>=1))
       phi = par%phie*EXP(-LOG(Sliq)/par%lam)*EXP(par%eta*LOG(Sliq))
    endwhere

    WHERE (h0>zero)
       hice(:)  = h0(:)*(S(:,1)-Sliq(:,1))
       phie(:)  = par(:,1)%phie*EXP(-LOG(Sliq(:,1))/par(:,1)%lam)*EXP(par(:,1)%eta*LOG(Sliq(:,1)))
       !phi(:,1) = max((phie -par(:,1)%he*Ksat(:,1)), (one+e5)*phie)+(h0(:)-hice(:))*Ksat(:,1)
       phi(:,1) = (one+e5)*phie + (h0(:)-hice(:))*Ksat(:,1)
    endwhere
    var%phi = phi

    ! ----------------------------------------------------------------

    DO kk=1, mp
       vsnow(kk)%nsnow = ssnow%nsnow(kk)
       vsnow(kk)%wcol  = ssnow%snowd(kk)/thousand ! diagnostic SWE (with or without dedicated snow pack)
       IF ( vsnow(kk)%nsnow > 0)  THEN ! dedicated snow pack
          ! define variables associated with snow
          ssnow%isflag(kk)   = 1
          vsnow(kk)%hsnow(:) = ssnow%smass(kk,1:nsnow_max)/thousand ! SWE (state variable)
          vsnow(kk)%depth(:) = ssnow%sdepth(kk,1:nsnow_max)  ! depth of snow pack (m)
          vsnow(kk)%dens(:)  = ssnow%ssdn(kk,1:nsnow_max)
          WHERE (vsnow(kk)%dens <= 200._r_2)
             vsnow(kk)%fsnowliq_max = 0.03_r_2
          ELSEWHERE
             vsnow(kk)%fsnowliq_max = 0.03_r_2 + (0.1_r_2 - 0.03_r_2)*(vsnow(kk)%dens-200._r_2)/vsnow(kk)%dens
          endwhere
          vsnow(kk)%fsnowliq_max = 0.1 !MC! ???
          vsnow(kk)%tsn(:) = REAL(ssnow%tggsn(kk,1:nsnow_max),r_2) - Tzero
          vsnow(kk)%kH(:)  = ssnow%sconds(kk,1:nsnow_max)
          vsnow(kk)%Dv(:)  = Dva*(MAX(REAL(ssnow%tggsn(kk,1:nsnow_max),r_2)/Tzero,0.0_r_2))**1.88_r_2 ! m2 s-1
          vsnow(kk)%sl     = slope_esat_ice(vsnow(kk)%tsn) * Mw/thousand/Rgas/(vsnow(kk)%tsn+Tzero)
          vsnow(kk)%kE     = vsnow(kk)%Dv*vsnow(kk)%sl*thousand*lambdaf
          vsnow(kk)%kth    = vsnow(kk)%kE + vsnow(kk)%kH
          vsnow(kk)%cv     = esat_ice(vsnow(kk)%tsn)*Mw/thousand/Rgas/(vsnow(kk)%tsn+Tzero) ! m3 m-3
          vsnow(kk)%hliq   = ssnow%snowliq(kk,1:nsnow_max)/thousand ! amount of liq snow water
          vsnow(kk)%melt   = zero ! amount of melted snow leaving each snowlayer (mm/dt)
       ELSE
          ssnow%isflag(kk)       = 0
          vsnow(kk)%hsnow        = zero
          vsnow(kk)%depth        = zero
          vsnow(kk)%dens         = 120_r_2 ! snow density kg m-3
          vsnow(kk)%tsn          = zero
          vsnow(kk)%kH           = 0.16_r_2    ! snow thermal cond (W m-2 K-1)
          vsnow(kk)%Dv           = Dva ! *(Tzero/Tzero)**1.88_r_2 ! m2 s-1
          vsnow(kk)%sl           = zero
          vsnow(kk)%kE           = zero
          vsnow(kk)%kth          = vsnow(kk)%kH
          vsnow(kk)%cv           = zero
          vsnow(kk)%hliq         = zero
          vsnow(kk)%melt         = zero
          vsnow(kk)%Jlatent      = zero
          vsnow(kk)%Jsensible    = zero
          vsnow(kk)%J            = zero
          vsnow(kk)%nsnow        = 0
          vsnow(kk)%fsnowliq_max = 0.03_r_2
       ENDIF
       ! heat stored in snowpack
       WHERE (vsnow(kk)%hsnow > zero)
          vsnow(kk)%Jsensible = (vsnow(kk)%hsnow-vsnow(kk)%hliq)*rhow*csice*(vsnow(kk)%Tsn) + &
               vsnow(kk)%hliq*rhow*cswat*(vsnow(kk)%Tsn)
          vsnow(kk)%Jlatent   = (vsnow(kk)%hsnow-vsnow(kk)%hliq)*rhow*(-lambdaf)
       ELSEWHERE
          vsnow(kk)%Jsensible = zero
          vsnow(kk)%Jlatent   = zero
       END WHERE
       vsnow(kk)%totdepth               = SUM(vsnow(kk)%depth)
       vsnow(kk)%J                      = SUM(vsnow(kk)%Jsensible+vsnow(kk)%Jlatent)
       vsnow(kk)%deltaJ                 = zero
       vsnow(kk)%deltawcol              = zero
       vsnow(kk)%deltaJlatent           = zero
       vsnow(kk)%deltaJsensible         = zero
       vsnow(kk)%Qadv_snow              = zero
       vsnow(kk)%Qadv_rain              = zero
       vsnow(kk)%Qadv_melt              = zero
       vsnow(kk)%Qadv_vap               = zero
       vsnow(kk)%Qcond_net              = zero
       vsnow(kk)%Qadv_transfer          = zero
       vsnow(kk)%Qmelt                  = zero
       vsnow(kk)%Qtransfer              = zero
       vsnow(kk)%FluxDivergence         = zero
       vsnow(kk)%Qvap                   = zero
       vsnow(kk)%MoistureFluxDivergence = zero
       vsnow(kk)%Qprec                  = zero
       vsnow(kk)%Qevap                  = zero
       vsnow(kk)%nsnow_last             = vsnow(kk)%nsnow
       ssnow%osnowd(kk) = ssnow%snowd(kk)
    ENDDO

    deltaTa = zero
    lE_old  = ssnow%lE
    gamm    = REAL(veg%gamma,r_2)

    IF (cable_user%test_new_gw) THEN
       WHERE (canopy%through>=met%precip_sn)
          qprec      = ssnow%fwtop/thousand            ! liq precip rate (m s-1)
          qprec_snow = (met%precip_sn)/thousand/dt
       ELSEWHERE
          qprec = MAX(ssnow%fwtop/thousand, zero)
          qprec_snow = zero
       endwhere
    ELSE
       WHERE (canopy%through>=met%precip_sn)
          qprec      = MAX((canopy%through-met%precip_sn)/thousand/dt , zero) ! liq precip rate (m s-1)
          qprec_snow = (met%precip_sn)/thousand/dt
       ELSEWHERE
          qprec = MAX(canopy%through, zero)
          qprec_snow = zero
       endwhere
    END IF

    ! re-calculate qprec_snow and qprec based on total precip and air T (ref Jin et al. Table II, Hyd Proc, 1999
    ! qprec_tot = qprec + qprec_snow
    ! where (vmet%Ta > 2.5_r_2)
    !    qprec_snow = zero
    !    qprec = qprec_tot
    ! elsewhere ((vmet%Ta <= 2.5_r_2) .and. (vmet%Ta > 2.0_r_2))
    !    qprec_snow = 0.6_r_2 * qprec_tot
    !    qprec = qprec_tot - qprec_snow
    ! elsewhere ((vmet%Ta <= 2.0_r_2) .and. (vmet%Ta > zero))
    !    qprec_snow = (1._r_2 - (54.62_r_2 - 0.2_r_2 *(vmet%Ta + Tzero)))*qprec_tot
    !    qprec = qprec_tot - qprec_snow
    ! elsewhere (vmet%Ta <= zero)
    !    qprec = zero
    !    qprec_snow = qprec_tot
    ! endwhere
    qprec_snow_tmp = qprec_snow
    h0old = ssnow%h0 ! pond height

    DO kk=1, mp
       hsnowold(kk) = SUM(vsnow(kk)%hsnow(1:nsnow_max))
    ENDDO

    ! Heat balance variables

    ! Water balance variables:
    ipi = SUM(ssnow%thetai*dx,2)  + h0*ssnow%thetai(:,1)/par(:,1)%thre       ! ice in profile initially
    ! water in profile initially
    wpi = SUM((par%thr + (par%the-par%thr)*S)*dx,2)  + plit%thre*SL*dxL

    nsteps = 0
    ti     = zero
    tf     = dt ! initial and final times

    win    = zero ! water input (total precip)
    evap   = zero
    runoff = zero
    infil  = zero
    drn    = zero

    IF (SEB_only == 0) THEN

       IF (cable_user%fwsoil_switch.NE.'Haverd2013') THEN

          DO kk=1, mp
             CALL getrex(ssnow%S(kk,:), ssnow%rex(kk,:), fws(kk), FS(kk,:), par(kk,:)%the, &
                  par(kk,:)%thw, Etrans(kk), gamm(kk), dx(kk,:), REAL(dt,r_2))
          ENDDO

       ELSE
          fws = canopy%fwsoil

       ENDIF

       DO k=1,ms
          qex(:,k)= ssnow%rex(:,k)
       ENDDO

       IF (cable_user%test_new_gw) THEN

          qex(:,:) = qex(:,:) + ssnow%qhlev(:,1:ms)/thousand
          qex(:,ms) = qex(:,ms) + ssnow%Qrecharge(:)/thousand
          botbc = "zero flux"

       END IF


    ENDIF

    ! topmodel - 0:no; 1: only sat; 2: only base; 3: sat and base
    zdelta     = zero
    fsat       = zero
    runoff_sat = zero
    qb         = zero
    IF (topmodel > 0) THEN
       zdelta = MAX( SUM((par%thr + (par%the-par%thr))*dx,2) - wpi, zero )
       IF (MOD(topmodel, 2) == 1) THEN
          ! saturated fraction
          fsat       = MIN( par(:,ms)%fsatmax * EXP(-zdelta), one )
          runoff_sat = fsat * qprec * thousand * dt
          qprec      = qprec * (one-fsat)  ! precip available for infiltration
       ENDIF
       IF (topmodel > 1) THEN
          botbc = "zero flux"
          ! calculate topmodel base flow
          qb        = alpha * par(:,ms)%Ke*EXP(-(par(:,ms)%zeta*zdelta))
          qex(:,ms) =  qex(:,ms) + qb
       ENDIF
    ENDIF



    IF (SEB_only == 1) THEN

       DO kk=1, mp

          ! call hyofS(S(kk,:), Tsoil(kk,:), par(kk,:), var(kk,:))
          CALL hyofS(S(kk,1), Tsoil(kk,1), par(kk,1), var(kk,1))
          CALL SEB(ms, par(kk,:), vmet(kk), vsnow(kk), var(kk,:), qprec(kk), qprec_snow(kk), dx(kk,:), &
               h0(kk), Tsoil(kk,:), &
               Tsurface(kk), G0(kk), lE(kk),Epot(kk), &
               tmp1d1a, tmp1d2, tmp1d3, tmp1d4, &
               tmp1d5, tmp1d6, tmp1d7, tmp1d8, tmp1d9,tmp1d10, tmp1d11, &
               tmp1d12,tmp1d13, tmp1d14, tmp1d15, tmp1d16, ktau)

       ENDDO
       canopy%ga  = REAL(G0)
       canopy%fes = REAL(lE)
       canopy%fhs = canopy%fns - canopy%ga - REAL(canopy%fes)
       ssnow%tss  = REAL(Tsurface + Tzero)
       ssnow%potev  = REAL(Epot)

    ELSE ! full SLI
       ! save for output, because they get changed with litter in solve
       rbw = vmet(1)%rbw
       rbh = vmet(1)%rbh
       rrc = vmet(1)%rrc
       !write(*,*), 'b4 solve', ktau
       CALL solve( ti, tf, ktau, mp, qprec, qprec_snow, ms, dx, &
            h0, S, thetai, Jsensible, Tsoil, evap, &
            evap_pot, runoff, infil, drn, discharge, qh, &
            nsteps, vmet, vlit, vsnow, var, csoil, kth, phi, T0, Tsurface, &
            H, lE, G0, Qadvcum, Jcol_sensible, &
            Jcol_latent_S, Jcol_latent_T, deltaice_cum_T, &
            deltaice_cum_S, dxL, zdelta, SL, TL, &
            plit, par, qex=qex, &
            wex=wex, qvsig=qvsig, qlsig=qlsig, qvTsig=qvTsig, qvh=qvh, &
            deltaTa=deltaTa, lE_old=lE_old, &
            dolitter=litter, doisotopologue=isotopologue, docondition=condition, &
            doadvection=advection, err=err)

       qprec_snow = qprec_snow_tmp
       H             = (H/(tf-ti))
       lE            = lE/(tf-ti)
       G0            = G0/(tf-ti)
       Jcol_latent_S = Jcol_latent_S/(tf-ti)
       Jcol_latent_T = Jcol_latent_T/(tf-ti)
       Jcol_sensible = Jcol_sensible/(tf-ti)
       Qadvcum       = Qadvcum/(tf-ti)

       DO kk=1, mp
          tmp1d1(kk) = (SUM(vsnow(kk)%Jsensible) + SUM(vsnow(kk)%Jlatent))
          ! heat stored in snowpack
          WHERE (vsnow(kk)%hsnow > zero)
             vsnow(kk)%Jsensible = (vsnow(kk)%hsnow-vsnow(kk)%hliq)*rhow*csice*(vsnow(kk)%Tsn) + &
                  vsnow(kk)%hliq*rhow*cswat*(vsnow(kk)%Tsn)
             vsnow(kk)%Jlatent = (vsnow(kk)%hsnow-vsnow(kk)%hliq)*rhow*(-lambdaf)
          ELSEWHERE
             vsnow(kk)%Jsensible = zero
             vsnow(kk)%Jlatent = zero
          endwhere
          deltaEsnow(kk) = SUM(vsnow(kk)%Jsensible) + SUM(vsnow(kk)%Jlatent) - tmp1d1(kk)
          deltah0(kk) = h0(kk)-h0old(kk)+SUM(vsnow(kk)%hsnow(1:vsnow(kk)%nsnow))-hsnowold(kk)
       ENDDO
       DO k=1, ms
          WHERE (err(:)==0) ssnow%thetai(:,k) = thetai(:,k)
       END DO
       ip  = SUM(ssnow%thetai*dx,2) + h0*ssnow%thetai(:,1)/par(:,1)%thre   ! ice in profile at tf
       ! water at tf
       wp  = SUM((par%thr + (par%the-par%thr)*S)*dx,2) + plit%thre*SL*dxL
       win = win + (qprec+qprec_snow)*(tf-ti)

       IF (verbose) THEN
          k=1
          WRITE(332,"(i8,i8,18e16.6)") ktau, nsteps(k), wp(k)-wpi(k), infil(k)-drn(k), runoff(k), &
               win(k)-(wp(k)-wpi(k)+deltah0(k)+runoff(k)+evap(k)+drn(k))-Etrans(k)*dt, wp(k), &
               evap(k), evap_pot(k), infil(k), &
               drn(k), h0(k), Etrans(k)*dt, discharge(k), fws(k), (ip(k)-ipi(k)), fsat(k), runoff_sat(k), qb(k)

!if (ktau==5) then
!
! write(*,*) win(k), (wp(k)-wpi(k)), deltah0(k),runoff(k), evap(k), drn(k), Etrans(k)*dt, canopy%fwsoil(k), canopy%fevc(k)
!stop
!endif


          WRITE(334,"(100f15.6)") S(k,:), S(k,:)*par(k,:)%thre+par(k,:)%thr
          WRITE(336,"(100f15.6)") Tsoil(k,:)
          WRITE(335,"(100e20.12)") vmet(k)%Ta, Tsurface(k), zero, H(k), lE(k), &
               G0(k),Jcol_sensible(k),Jcol_latent_S(k), Jcol_latent_T(k), &
               vmet(k)%Rn, TL(k), SL(k), deltaice_cum_T(k), &
               deltaice_cum_S(k), zero, Tsurface(k), vmet(k)%rha, &
               Qadvcum(k), SUM((Jsensible(k,:)-ssnow%gammzz(k,:)),1)
          WRITE(338,"(100f18.6)") thetai(k,:)
          WRITE(369,"(20e20.12)") vmet(k)%Ta, vmet(k)%rha, rbw, &
               rbh, rrc, vmet(k)%Rn, &
               vmet(k)%Da, vmet(k)%cva, vmet(k)%civa, &
               vmet(k)%phiva, Etrans(k), qprec(k), qprec_snow(k), rad%latitude(k),  rad%longitude(k)
          WRITE(370,"(20e20.12)")  qex
          WRITE(371,"(20e20.12)")  qvsig+qlsig

       ENDIF

       IF (cable_user%test_new_gw) THEN

          DO i=1,mp
             ssnow%GWwb(i) = ssnow%GWwb(i)  + (ssnow%Qrecharge(i)-ssnow%qhlev(i,ms+1))*dt/soil%GWdz(i)/thousand

             IF (ssnow%GWwb(i) .GT. soil%GWssat_vec(i)) THEN
                ssnow%qhlev(i,ms+1) = ssnow%qhlev(i,ms+1)  + (ssnow%GWwb(i) - soil%GWssat_vec(i))*soil%GWdz(i)*thousand/dt
                ssnow%GWwb(i) = soil%GWssat_vec(i)
                ssnow%qhz(i) = SUM(ssnow%qhlev(i,1:ms+1),dim=1)
             END IF

             IF (ssnow%GWwb(i) .LT. soil%GWwatr(i)) THEN
                ssnow%qhlev(i,ms+1) = ssnow%qhlev(i,ms+1)  - (ssnow%GWwb(i) - soil%GWwatr(i))*soil%GWdz(i)*thousand/dt
                ssnow%GWwb(i) = soil%GWwatr(i)
                ssnow%qhz(i) = SUM(ssnow%qhlev(i,1:ms+1),dim=1)
             END IF
          END DO
       END IF


       ! Update variables for output:
       WHERE (err(1:mp) == 0)
          ssnow%tss      = REAL(Tsurface + Tzero)
          ssnow%wbtot    = REAL(wp*thousand)
          canopy%ga      = REAL(G0)
          canopy%fes     = lE
          canopy%fhs     = canopy%fns - canopy%ga - canopy%fes
          ssnow%rnof1    = REAL(runoff*thousand/dt )
          ssnow%rnof2    = REAL(drn*thousand/dt )
          ssnow%runoff = ssnow%rnof1 + ssnow%rnof2
          ssnow%zdelta   = zdelta
          ssnow%SL       = SL
          ssnow%TL       = TL
          ssnow%delwcol  = (wp-wpi+deltah0)*thousand  ! includes change in snow pack via deltah0
          ssnow%Tsurface = Tsurface
          ssnow%lE       = lE
          ssnow%evap     = evap*thousand
          ssnow%nsteps   = REAL(nsteps)
          ssnow%h0       = h0
       endwhere
       DO k=1, ms
          WHERE (err(:) == 0)
             ssnow%tgg(:,k)    = REAL(Tsoil(:,k) + Tzero)
             ssnow%wb(:,k)     = REAL(S(:,k)*(par(:,k)%thr+(par(:,k)%the-par(:,k)%thr)))
             ssnow%wbice(:,k)  = thetai(:,k)
             ssnow%S(:,k)      = S(:,k)
             ssnow%Tsoil(:,k)  = Tsoil(:,k)
             ssnow%rex(:,k)    = wex(:,k)*thousand
             ssnow%kth(:,k)    = kth(:,k)
             ssnow%gammzz(:,k) = Jsensible(:,k)
          END WHERE
       END DO

       IF (cable_user%fwsoil_switch.NE.'Haverd2013') THEN
          WHERE (err(1:mp) == 0) canopy%fwsoil = REAL(fws)
       ENDIF

       IF (litter==0) THEN
          ssnow%rlitt = zero
       ELSE
          WHERE (err(1:mp) == 0) ssnow%rlitt = dxL/vlit%Dv
       ENDIF

       ! update CABLE snow variables
       DO kk=1, mp
          IF (err(kk) == 0) THEN
             ssnow%snowd(kk)               = REAL(vsnow(kk)%wcol*thousand)     ! amount of snow  (mm liq water eq)
             ! amount of snow in dedicated snow pack (mm liq water eq)
             ssnow%smass(kk,1:nsnow_max)   = REAL(vsnow(kk)%hsnow(:)*thousand)
             ssnow%sdepth(kk,1:nsnow_max)  = REAL(vsnow(kk)%depth(:))          ! depth of snow pack (m)
             ssnow%ssdn(kk,1:nsnow_max)    = REAL(vsnow(kk)%dens(:))           ! density of snow (kg m-3)
             ssnow%tggsn(kk,1:nsnow_max)   = REAL(vsnow(kk)%tsn(:) + Tzero)    ! abs T of snowpack
             ssnow%sconds(kk,1:nsnow_max)  = REAL(vsnow(kk)%kH(:))             ! thermal conductivty of snowpack
             ssnow%snowliq(kk,1:nsnow_max) = REAL(vsnow(kk)%hliq(:)*thousand)  ! amount of liq snow water
             ! amount of melted snow leaving bottom of snow pack (mm/dt)
             ssnow%smelt(kk)               = REAL(vsnow(kk)%Qmelt*thousand/dt)
             ssnow%nsnow(kk)               = vsnow(kk)%nsnow
             IF (SUM(ssnow%sdepth(kk,1:nsnow_max)) > zero) THEN
                ssnow%ssdnn(kk) = ssnow%snowd(kk)/SUM(ssnow%sdepth(kk,1:nsnow_max))
             ENDIF
          ENDIF
       ENDDO

       WHERE (err(1:mp) == 0) ssnow%isflag = 0

       ! snow output
       IF (1 == 0) THEN
          k = 1
          WRITE(340,"(100e16.6)") SUM(vsnow(k)%hsnow(1:vsnow(k)%nsnow)), vsnow(k)%tsn(1),SUM(vsnow(k)%hliq(1:vsnow(k)%nsnow)), &
               qprec_snow(k)*dt, vsnow(k)%Qmelt, qprec(k)*dt, &
               vsnow(k)%Qevap,vsnow(k)%Qvap,ssnow%albsoilsn(k,1), ssnow%albsoilsn(k,2), ssnow%sconds(k,1), &
               vsnow(k)%dens(1),SUM(vsnow(k)%depth(1:vsnow(k)%nsnow)), vsnow(k)%J, &
               vsnow(k)%MoistureFluxDivergence, vsnow(k)%FluxDivergence, vsnow(k)%dens(nsnow_max), vsnow(k)%tsn(nsnow_max), &
               qh(k,0), vmet(k)%rbh
       ENDIF

       WHERE (err(1:mp) == 0)
          canopy%ofes = canopy%fes
          ! Update total latent heat to reflect updated soil component:
          canopy%fe = canopy%fev + REAL(canopy%fes)
          ! Update total sensible heat to reflect updated soil component:
          canopy%fh = REAL(canopy%fhv) + canopy%fhs
       END WHERE

    ENDIF ! SEB only

  END SUBROUTINE sli_main
END MODULE sli_main_mod
