! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/soil_temp.F
      SUBROUTINE SOIL_TEMP (DT, G, DRAINAGE, W_FLUX, ETRAN_L, SUBSOIL_L)

      IMPLICIT NONE

!-----------------------------------------------------------------------
! Description:
!     Updates deep soil temperatures, frozen and unfrozen
!     frozen soil water content.
!-----------------------------------------------------------------------

      include "size.h"
      include "mtlm.h"

      INTEGER ITER, ITERMAX             ! Number of iterations and maximum iteration number.
      PARAMETER (ITERMAX=25)

      REAL ERRMAX                       ! Maximum allowable flux error (W/m2).

      PARAMETER (ERRMAX=0.00000005)

      INTEGER I, J, L, N                ! Loop counters.

      LOGICAL NOT_DONE                  ! logical flag to keep iteratating

!     input variables
      REAL DT                           ! Timestep (s)
      REAL G(POINTS)                    ! Net downward surface heat flux (W/m2).
      REAL DRAINAGE(POINTS)
      REAL W_FLUX(POINTS,0:NGND)
      REAL ETRAN_L(POINTS,NGND)
      REAL SUBSOIL_L(POINTS,NGND)

!     local variables
      REAL DHGND(NGND)                  ! Current heat increment to the layer (J/m2)
      REAL DTHU                         ! Rate of change of unfrozen soil moisture concentration with temp (m3 H2O/m3 soil/K)
      REAL DTGND(NGND)
      REAL HCAP_A(NGND)
      REAL D_FLUX(0:NGND)               ! The fluxes of heat between layers (W/m2)
      REAL MF(NGND)                     ! Frozen moisture content of each soil layer (kg/m2)
      REAL MF0(NGND)                    ! Previous value of MF (kg/m2)
      REAL MSAT(NGND)                   ! The saturation moisture content of each layer (kg/m2)
      REAL MU(NGND)                     ! Unfrozen moisture content of each soil layer (kg/m2)
      REAL MU0(NGND)                    ! Previous value of MU (kg/m2).
      REAL TMAX(NGND)                   ! Temperature above which all water is unfrozen (K)
      REAL TGND0(NGND)                  ! Previous value of TGND (K).
      REAL A_FLUX(NGND)                 ! Heat increment to layer due to advection of water (J/m2)

      REAL DZ_EFF(NGND)                 ! Effective thickenss of layer (soil + excess ice) (m)
      REAL HCON(NGND)                   ! Thermal conductivity  (W/mK)
      REAL MAXMELT                      ! Maximum amount of layer ice that can be melted (structural excess) (kg/m2)
      REAL RDZ_GND(NGND)                ! reciprical of DZ_GND
      REAL STH(NGND)                    ! SU + SF
      REAL RSTH(NGND)                   ! reciprical of STH
      REAL ERRMAXDT                     ! ERRMAX converted from flux error to total heat error

      REAL TG(NGND)                     ! one dimensional form of TGND
      REAL HC_B(NGND)                   ! one dimensional form of HCON_B
      REAL W_FL(0:NGND)                 ! one dimensional form of W_FLUX

!     work variables
      REAL DF1_DT1(NGND), DF1_DT2(NGND), DF2_DT1(NGND), DF2_DT2(NGND)
      REAL GA(NGND), GB(NGND), GC(NGND), GD(NGND), GAMCON, XMAX, XMIN

!~Slab excess ice variables~
! ADDROCK              ! Heat to add to rock layer if needed (J/m2)
! DHSB                 ! Heat increment to slab ice (J/m2)
! DTSLAB               ! Temperature increment to slab ice (K)
! DZ_SLAB              ! Thickness of slab (m)
! HCICE                ! Thermal conductivity of ice (W/mK)
! D_FLUXSLAB           ! Heat flux at top / bottom of slab (W/m2)
! MAXHEAT              ! Amount of heat required to fully melt slab (J/m2)
! MSLAB0               ! Initial mass of slab (kg/m2)

! HCSAT                ! The thermal conductivity of the saturated soil at current ratio of ice to liquid water (W/m/K).
! HCAIR                ! Thermal conductivity of air (W/m/K).
! HCICE                ! Thermal conductivity of ice (W/m/K).
! HCWAT                ! Thermal conductivity of liquid water (W/m/K)
! HCON                 ! Dry soil thermal conductivity (W/m/K).
! HCONS                ! The thermal conductivity between adjacent layers including effects of water and ice (W/m/K)
! SU                   ! Fractional saturation of unfrozen water at layer boundaries
! SF                   ! Fractional saturation of frozen water at layer boundaries
! THICE                ! The concentration of ice at saturation  for the current mass fraction of liquid water (m3 H2O/m3 soil)
! THWAT                ! The concentration of liquid water at  saturation for the current mass fraction of liquid water  (m3 H2O/m3 soil).
! V_SAT                ! Volumetric soil moisture concentration at saturation (m3/m3 soil)

      REAL HCAIR, HCICE, HCWAT, THWAT, THICE, HCSAT
      PARAMETER (HCAIR=0.025, HCICE=2.24, HCWAT=0.56)

      REAL HCAP_WmI

!-----------------------------------------------------------------------------
! Initializations - perform calculations that need to occur only 1x per run.
!-----------------------------------------------------------------------------

      XMAX = 1.0E4
      XMIN = -1.0E4
      HCAP_WmI = HCAP_W - HCAP_I
      HCAP_A(:) = HCAP_ROCK
      do N=1,NGND
        RDZ_GND(N) = 1./DZ_GND(N)
      enddo
      ERRMAXDT = ERRMAX*DT

!-----------------------------------------------------------------------------
! Loop over all land points
!-----------------------------------------------------------------------------
      do I=1,LAND_PTS
        L = LAND_INDEX(I)

!       set local soil column arrays from global arrays
        do N=1,NGND
          TG(N) = TGND(L,N)
          HC_B(N) = HCON_B(L,N)
          STH(N) = SU(L,N) + SF(L,N)
          RSTH(N) = 1./(STH(N) + 1e-30)
        enddo
        do N=0,NGND
          W_FL(N) = W_FLUX(L,N)
        enddo

        do N=1,NSOIL(L)
!         diagnose the frozen and unfrozen water.
          MSAT(N) = RHO_W*DZ_GND(N)*V_SAT(L,N)
          MU(N) = MSAT(N)*SU(L,N)
          MF(N) = MSAT(N)*SF(L,N)

          DZ_EFF(N) = DZ_GND(N)

        enddo

        do N=NSOIL(L)+1,NGND
          DZ_EFF(N) = DZ_GND(N)
        enddo

!-----------------------------------------------------------------------
! Calculate soil heat capacity and thermal conductivity
!-----------------------------------------------------------------------
!       surface conductivity and heat flow calculated in penmon
        do N=1,NSOIL(L)

          TMAX(N) = ZERODEGC - PHIS(L,N)/DPSIDT
     &              *(MSAT(N)/M(L,N) + 1.e-30)**(B(L,N))
          THWAT = V_SAT(L,N)*SU(L,N)*RSTH(N)
          THICE = V_SAT(L,N)*SF(L,N)*RSTH(N)
!          HCSAT = (HCWAT**THWAT)*(HCICE**THICE)/(HCAIR**V_SAT(L,N))
!          HCON(N) = HCON_D(L,N)*(((HCSAT) - 1.)*STH(N) + 1.)
!# if defined O_mtlm_exstruct
!          HCON(N) = (HCON(N)*DZ_GND(N) + HCICE*DZ_STRUCT(N))/DZ_EFF(N)
!# endif

          HCSAT = HCON_D(L,N)*(HCWAT**THWAT)*(HCICE**THICE)
     &            /(HCAIR**V_SAT(L,N))
          HCON(N) = ((HCSAT - HCON_D(L,N))*STH(N) + HCON_D(L,N))

        enddo

!       calculate the thermal conductivity at the layer boundaries
        do N=1,NSOIL(L)-1
          HC_B(N) = (DZ_EFF(N+1)*HCON(N) + DZ_EFF(N)*HCON(N+1))
     &              /(DZ_EFF(N+1) + DZ_EFF(N))
        enddo

        SNOW_HT(L) = LYING_SNOW(L)/RHO_S
        if (SNOW_HT(L) .gt. 0.) then
          if (SNOW_HT(L) .gt. (0.5*DZ_EFF(1))) then
            HC_B(1) = HC_B(1)*(DZ_EFF(1) + DZ_EFF(2))
     &                /(2.*DZ_EFF(1) + DZ_EFF(2)
     &              + (HCON(1)/HCON_SNOW)*(2.*SNOW_HT(L) - DZ_EFF(1)))
          else
            HC_B(1) = HC_B(1)*(DZ_EFF(1) + DZ_EFF(2))
     &                /(DZ_EFF(1) + DZ_EFF(2) + 2.*SNOW_HT(L))
          endif
        endif

! calculate thermal conductivity at boundaries

        HC_B(NSOIL(L)) = (DZ_GND(NSOIL(L)+1)*HCON(NSOIL(L))
     &                 + DZ_EFF(NSOIL(L))*HCON_ROCK)
     &                   /(DZ_GND(NSOIL(L)+1) + DZ_EFF(NSOIL(L)))

!--------------------------------------------------------------------
! Now calculate heat fluxes by diffusion
!--------------------------------------------------------------------
!       calculate soil heat flux

        do N=1,NSOIL(L)
          D_FLUX(N) = -HC_B(N)*2.*(TG(N+1) - TG(N))
     &                /(DZ_EFF(N+1)+DZ_EFF(N))
        enddo

!       calculate rock heat fluxes
        do N=NSOIL(L)+1,NGND-1
          D_FLUX(N) = -HCON_ROCK*2.*(TG(N+1) - TG(N))
     &                /(DZ_GND(N+1) + DZ_GND(N))
        enddo

!       set boundary conditions (upward geothermal is negative)
        D_FLUX(0) = G(L)
        D_FLUX(NGND) = -GEOHEAT(L)

!       calculate advective heat flux and coefficients
        call SOIL_ADVECT (NGND, NSOIL(L), HCAP_W, HCON_ROCK
     &,                   DZ_EFF, HC_B, TG, W_FL, A_FLUX
     &,                   DF1_DT1, DF1_DT2, DF2_DT1, DF2_DT2)

!       heat flux into a layer
        do N=1,NGND
          DHGND(N) = DT*(D_FLUX(N-1) - D_FLUX(N) + A_FLUX(N))
        enddo

!----------------------------------------------------------------------
! solve temperature for each soil column
!----------------------------------------------------------------------
        ITER = 0
        NOT_DONE = .true.
        do while (NOT_DONE .and. ITER .lt. ITERMAX)
          NOT_DONE = .false.
          ITER = ITER + 1

          do N=1,NSOIL(L)
            TGND0(N) = TG(N)
            MF0(N) = MF(N)
            MU0(N) = MU(N)
            if ((TG(N) .eq. TMAX(N) .and. DHGND(N) .ge. 0.0) .or.
     &        (TG(N) .gt. TMAX(N))) then

              HCAP_A(N) = HCAP_D(L,N) + (HCAP_WmI*MU(N)

     &                  + HCAP_I*M(L,N))*RDZ_GND(N)

            else
!             phase changes
              DTHU = DPSIDT*MSAT(N)/(B(L,N)*PHIS(L,N)*RHO_W*DZ_GND(N))
     &               *(-DPSIDT*(TG(N) - ZERODEGC)/PHIS(L,N))
     &               **(-1./B(L,N) - 1.)
               HCAP_A(N) = HCAP_D(L,N) + (HCAP_WmI*MU(N)

     &                  + HCAP_I*M(L,N))*RDZ_GND(N)

     &                  + RHO_W*DTHU*(HCAP_WmI*TG(N) + LHF)

            endif
          enddo

          do N=1,NGND
            TGND0(N) = TG(N)
            GAMCON = DT/(HCAP_A(N)*DZ_EFF(N))
            GA(N) = -GAMCON*DF1_DT1(N)
            GB(N) = 1. - GAMCON*(DF1_DT2(N) - DF2_DT1(N))
            GC(N) = GAMCON*DF2_DT2(N)
            GD(N) = DHGND(N)/(HCAP_A(N)*DZ_EFF(N))
          enddo

          CALL GAUSSIAN (NGND, GA, GB, GC, GD, XMIN, XMAX, DTGND)

!         recalculate frozen and unfrozen moisture
          do N=1,NSOIL(L)
            DHGND(N) = HCAP_A(N)*DZ_EFF(N)*DTGND(N)
            TG(N) = TG(N) + DTGND(N)
            if ((TG(N) - TMAX(N))*(TGND0(N) - TMAX(N)) .LT. 0.) then
              DTGND(N) = DTGND(N) + TMAX(N) - TG(N)
              TG(N) = TMAX(N)
            endif
            if (TG(N) .ge. TMAX(N)) then
              MU(N) = M(L,N)
              MF(N) = 0.
            else
              MU(N) = MSAT(N)*(-DPSIDT*(TG(N) - ZERODEGC)
     &                /PHIS(L,N))**(-1./B(L,N))
              MF(N) = M(L,N) - MU(N)
            endif
!           calculate the error in heat conservation
            DHGND(N) = DHGND(N) - (HCAP_D(L,N)*DZ_GND(N)*DTGND(N)
     &               + HCAP_W*(MU(N)*TG(N) - MU0(N)*TGND0(N))
     &               + HCAP_I*(MF(N)*TG(N) - MF0(N)*TGND0(N)))
     &               + (MF(N) - MF0(N))*LHF

            if (ABS(DHGND(N)) .gt. ERRMAXDT) then
              NOT_DONE = .true.
            else
              DHGND(N) = 0.
            endif
          enddo

          do N = NSOIL(L)+1,NGND
            DHGND(N) = HCAP_A(N)*DZ_EFF(N)*DTGND(N)
            TG(N) = TG(N) + DTGND(N)
            DHGND(N) = DHGND(N) - (HCAP_A(N)*DZ_GND(N)*DTGND(N))
          enddo

        enddo ! End of iteration do while

!-----------------------------------------------------------------------
! Diagnose fractional values of unfrozen an  frozen water.
!-----------------------------------------------------------------------
        do N=1,NSOIL(L)
          SU(L,N) = MU(N)/MSAT(N)
          SF(L,N) = MAX(MF(N)/MSAT(N),0.)
        enddo

        do N=1,NGND
          TGND(L,N) = TG(N)
          HCON_B(L,N) = HC_B(N)
        enddo
        do N=0,NGND
          W_FLUX(L,N) = W_FL(N)
        enddo

      enddo

      return
      end

