! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/mtlm_state.F

       subroutine MTLM_STATE (DT, G, ETRAN_L, MELTHEAT)

!-----------------------------------------------------------------------
! Routine to update land surface prognostic variables
! (soil moisture, soil temperature, lying snow mass).
! Also diagnoses runoff and snowmelt

!**********************************************************************
! this file is based on code that may have had the following copyright:
! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE.

! Permission has been granted by the authors to the public to copy
! and use this software without charge, provided that this Notice and
! any statement of authorship are reproduced on all copies. Neither the
! Crown nor the U.K. Meteorological Office makes any warranty, express
! or implied, or assumes any liability or responsibility for the use of
! this software.
!**********************************************************************
!-----------------------------------------------------------------------

      implicit none
      include "size.h"
      include "mtlm.h"
      include "mtlm_data.h"

!~ Loop/spatial variables ~
! I, L, N, X =  Loop counters
! Z          =  Depth counter
! NPERMA     =  Soil layer containing permafrost

      integer I, L, N, NPERMA, X
      real Z

!~ Constants ~
! DT        = Timestep (s).
! FILL      = Fill value for permafrost array
! SECSDAY   = Seconds per day
! SLOPEARRAY= Value of L_slope.nc to use for wetlands diagnosis
! WETVAL    = Value of SU that wet areas must exceed to be wetlands

      real DT, WETVAL
      real FILL,  SECSDAY
      integer SLOPEARRAY
      PARAMETER (FILL =  -9.9692099683868690e+36, SLOPEARRAY = 4)
      PARAMETER (SECSDAY = 86400, WETVAL = 0.65)

!~ Other variables ~
! DZ_EFF     = Effective layer thickness (excess ice) (m)
! ETRAN_L    = Evapotranspiration from each layer (kg/m2/s).
! G          = Ground heat flux (W/m2).
! HCAP_A     = Actual heat capacity (dry soil + liquid moisture)
! DRAINAGE   = Flow of water through ground water (kg/m2/s)
! INFIL      = Maximum rate of infiltration into soil
! MELTABLE   = Amount of lying_snow it is possbile to melt (kg/m2)
! MELTHEAT   = Heat available to melt snow (W/m2)
! W_FLUX     = Water flux (kg/m2/s) into layer N

      real ETRAN_L(POINTS, NGND), G(POINTS), INFIL(POINTS)
      real DRAINAGE(POINTS)
      real W_FLUX(POINTS, 0:NGND), HCAP_A
      real MELTABLE, MELTHEAT(POINTS), DZ_EFF(POINTS, NGND)
      real SUBSOIL_L(POINTS, NGND), ADDMELT
      real DELTA, MAXFROZEN(POINTS)

      DZ_EFF(:,:) = 0
! These 4 lines of code went missing when nsoil was integrated into model by
! Eby in autumn 2013. Need to be here for diagnostics and Kovin diffusion to
! work properly
      PERMA(:) = FILL
      PFTHICK(:) = 0.
      ACTIVELT(:) = 0.
      TALIK(:) = 0.

      do I=1,LAND_PTS

        L = LAND_INDEX(I)

        DZ_EFF(L,:) = DZ_GND(:)

!-------------------------------------------------------------------------
! Diagnose snow melt and update the lying snow
! Order of priority for snow fluxes:
! i)   Sublimation
! ii)  Surface melt
! iii) Subsurface melt
!-------------------------------------------------------------------------

!        i) Sublimation + snowfall fluxes
        SNOWMELT(L) = 0.0
        LYING_SNOW(L) = LYING_SNOW(L) + (SNOW(L) - SUBSNOW(L))*DT
        if (LYING_SNOW(L) .lt. 0.) then                  ! If negative snow, change extra sublimation to evaporation
          SUBSNOW(L) = SUBSNOW(L) + LYING_SNOW(L)/DT
          ETRAN_L(L,1) = ETRAN_L(L,1) - LYING_SNOW(L)/DT
          G(L) = G(L) + LHF*LYING_SNOW(L) / DT
          LYING_SNOW(L) = 0.
          MELTABLE = 0.
        else                                             ! Still snow left to melt
!         ii) Surface melt
          MELTABLE = 0.
          MELTABLE = min(MELTHEAT(L)*DT/LHF, LYING_SNOW(L))
          LYING_SNOW(L) = LYING_SNOW(L) - MELTABLE
          SNOWMELT(L) = SNOWMELT(L) + MELTABLE / DT
        endif

        G(L) = G(L) + MELTHEAT(L)- (MELTABLE * LHF / DT)  ! Adjust G in case too much melt heat

!       iii) Subsurface melt
        if ((LYING_SNOW(L).gt.0.0).and.(TGND(L,1).gt.ZERODEGC)) then
          HCAP_A = HCAP_D(L,1) + HCAP_W * RHO_W * V_SAT(L,1) * SU(L,1)
          if ((HCAP_A*DZ_GND(1)*(TGND(L,1)-ZERODEGC)/LHF)
     &      .gt. LYING_SNOW(L)) then
            TGND(L,1) = TGND(L,1) - (LYING_SNOW(L)*LHF)
     &                / (HCAP_A*DZ_GND(1))
            SNOWMELT(L) = SNOWMELT(L) + LYING_SNOW(L)/DT
            LYING_SNOW(L) = 0.!
          else
            SNOWMELT(L) = SNOWMELT(L) + HCAP_A*DZ_GND(1)
     &                  * (TGND(L,1)-ZERODEGC)/(LHF*DT)
            LYING_SNOW(L) = LYING_SNOW(L) - HCAP_A*DZ_GND(1)
     &                    * (TGND(L,1) - ZERODEGC) / LHF
            TGND(L,1) = ZERODEGC
          endif
        endif

        SNOW_HT(L) = LYING_SNOW(L)/RHO_S
        SURF_ROFF(L) = 0.
        W_FLUX(L,0) = RAIN(L) + SNOWMELT(L)

      enddo

!-------------------------------------------------------------------
! Call the soil hydrology code to diagnose water fluxes between soil
! layers and through the bottom of the soil column
!-------------------------------------------------------------------

      call SOIL_HYD (ETRAN_L, DT, DRAINAGE, W_FLUX, SUBSOIL_L)

!-------------------------------------------------------------------
! Call the soil temperature code to diagnose changes in soil layer
! temperatures, heat fluxes between layers and phase changes of water
!-------------------------------------------------------------------
      call SOIL_TEMP (DT, G, DRAINAGE, W_FLUX, ETRAN_L, SUBSOIL_L)

      do I = 1, LAND_PTS
        L = LAND_INDEX(I)

! Add drainage to surface runoff since the ocean can't tell the difference

        SURF_ROFF(L) = SURF_ROFF(L) + DRAINAGE(L)

!-----------------------------------------------------------------
! Update permafrost diagnostics
!-----------------------------------------------------------------
! Update the number of days that the ground has been frozen / unfrozen

        do N=1,NGND
          if (TGND(L,N) .lt. 273.15) then
            DFROZEN_C(L,N) = DFROZEN_C(L,N) + DT/SECSDAY
            DFROZEN_A(L,N) = DFROZEN_A(L,N) + DT/SECSDAY
          else
            DFROZEN_C(L,N) = 0.
          endif
          if (TGND(L,N) .gt. 273.15) then
            DUNFROZEN_C(L,N) = DUNFROZEN_C(L,N) + DT/SECSDAY
          else
            DUNFROZEN_C(L,N) = 0.
          endif
        enddo

        Z = 0.
        NPERMA = 0.
        do N=1,NGND
          if (DFROZEN_C(L,N) .gt. 730.) then  ! Have hit permafrost in the current layer
            PERMA(L) = Z
            NPERMA = N
            Z = Z + DZ_EFF(L,N)
            exit
          else  ! Otherwise, advance to next layer
            Z = Z + DZ_EFF(L,N)
          endif
        enddo
        if (NPERMA .lt. NGND) then ! Have found permafrost, but above bottommost layer
          do N=NPERMA+1, NGND
            if (DFROZEN_C(L,N) .lt. 730.) then
              PFTHICK(L) = Z-PERMA(L)
              exit
            else
              Z = Z + DZ_EFF(L,N)
            endif
            if (N .eq. NGND) then
              PFTHICK(L) = Z - PERMA(L)
              exit
            endif
          enddo
        elseif (DFROZEN_C(L,NGND) .gt. 730) then ! Permafrost just in the bottom layer
          PFTHICK(L) = DZ_EFF(L,NGND)
        endif
! find active layer thickness and talik thickness (for taliks above first PF layer)
        if (NPERMA .gt. 1) then   ! Must have found permafrost, but there will be an ALT
          do N=1,NPERMA-1
            if (DFROZEN_C(L,N) .lt. 730. .and.
     &        DUNFROZEN_C(L,N) .lt. 730.) then
              ACTIVELT(L) = ACTIVELT(L) + DZ_EFF(L,N)
            else
              exit
            endif
          enddo
          do N=1,NPERMA-1
            if (DUNFROZEN_C(L,N) .gt. 730.) then
              TALIK(L) = TALIK(L) + DZ_EFF(L,N)
            endif
          enddo
        endif

      enddo

!------------------------------------------------------------------
! Diagnose wetlands
!------------------------------------------------------------------
      YEARCOUNT = YEARCOUNT + DT/SECSDAY

      do I = 1, LAND_PTS
        L = LAND_INDEX(I)
        if (SU(L,1) .ge. WETVAL) then
          DWET(L) = DWET(L) + DT/SECSDAY
        endif
        if (LYING_SNOW(L) .gt. 0.) then
          DSNOW(L) = DSNOW(L) + DT/SECSDAY
        endif
      enddo

      if ((YEARCOUNT .ge. 365.)) then
        WETMAP(:,:) = 0.
        MAXFROZEN = MAXVAL(DFROZEN_A,2)
        FROZENMAP(:,:) = 0.
        do I = 1,LAND_PTS
          L = LAND_INDEX(I)
          DWET(L) = DWET(L) * 365./YEARCOUNT
          if ((DWET(L) .gt. 0).and.(DWET(L).le. 365)) then
          ANN_DWET(L) = DWET(L)
          else
          ANN_DWET(L) = 0.

          endif
          if (DWET(L) .gt. 0.) then
             WETMAP(L,1) = SLOPEMAP(L,SLOPEARRAY)
             if (DWET(L) .lt. 360.) then
               WETMAP(L,3) = SLOPEMAP(L,SLOPEARRAY)
             else
               WETMAP(L,2) = SLOPEMAP(L,SLOPEARRAY)
             endif
          endif

          if (PERMA(L) .gt. 0.) then
            FROZENMAP(L,1) = 1.
          elseif (MAXFROZEN(L) .gt. 15.) then
            FROZENMAP(L,2) = 1.
          elseif (MAXFROZEN(L) .ge. 1.) then
            FROZENMAP(L,3) = 1.
          elseif (MAXFROZEN(L) .gt. 0.) then
            FROZENMAP(L,4) = 1.
          endif
          if (DSNOW(L) .gt. 15.) then
            FROZENMAP(L,5) = 1.
          elseif (DSNOW(L) .ge. 1.) then
            FROZENMAP(L,6) = 1.
          elseif (DSNOW(L). gt. 0.) then
            FROZENMAP(L,7) = 1.
          endif

          DWET(L) = 0.
          DFROZEN_A(L,:) = 0.
          DSNOW(L) = 0.

        enddo

        YEARKOV_SW = 1

        YEARCOUNT = 0.
      endif

!------------------------------------------------------------------
! Update Schaefer Dmin
!-------------------------------------------------------------------
      do I = 1, LAND_PTS
        L = LAND_INDEX(I)
        if(INT(SCHDMIN(L)) .lt. SCHDMAX) then
	  do N = 1,SCHDMAX
	     if(INT(SCHDMIN(L)) .lt. N) then
	       if(TGND(L,N) .gt. 273.15) then
	          SCHDMIN(L)=REAL(N)
	       endif
	     endif
	  enddo
	endif
      enddo

!-------------------------------------------------------------------
! Reboot Schaefer Dmin during spinup b/c memory quantity
!-------------------------------------------------------------------

      return
      end