#include "LAND_OPTIONS.h"

CBOP
C     !ROUTINE: LAND_STEPFWD
C     !INTERFACE:
      SUBROUTINE LAND_STEPFWD(
     I                land_frc, bi, bj, myTime, myIter, myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R LAND_STEPFWD
C     | o Land model main S/R: step forward land variables
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     == Global variables ===
C-- size for MITgcm & Land package :
#include "LAND_SIZE.h"

#include "EEPARAMS.h"
#include "LAND_PARAMS.h"
#include "LAND_VARS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     land_frc :: land fraction [0-1]
C     bi,bj    :: Tile index
C     myTime   :: Current time of simulation ( s )
C     myIter   :: Current iteration number in simulation
C     myThid   :: Number of this instance of the routine
      _RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER bi, bj, myIter, myThid
      _RL myTime
CEOP

#ifdef ALLOW_LAND
C     == Local variables ==
C     i,j,k        :: loop counters
C     kp1          :: k+1
C     grd_HeatCp   :: Heat capacity of the ground [J/m3/K]
C     enthalpGrdW  :: enthalpy of ground water [J/m3]
C     fieldCapac   :: field capacity (of water) [m]
C     mWater       :: water content of the ground [kg/m3]
C     groundWnp1   :: hold temporary future soil moisture []
C     grdWexcess   :: ground water in excess [m/s]
C     fractRunOff  :: fraction of water in excess which leaves as runoff
C     flxkup       :: downward flux of water, upper interface (k-1,k)
C     flxdwn       :: downward flux of water, lower interface (k,k+1)
C     flxEngU      :: downward energy flux associated with water flux (W/m2)
C                     upper interface (k-1,k)
C     flxEngL      :: downward energy flux associated with water flux (W/m2)
C                     lower interface (k,k+1)
C     temp_af      :: ground temperature if above freezing
C     temp_bf      :: ground temperature if below freezing
C     mPmE         :: hold temporary (liquid) Precip minus Evap [kg/m2/s]
C     enWfx        :: hold temporary energy flux of Precip [W/m2]
C     enGr1        :: ground enthalpy of level 1  [J/m2]
C     mSnow        :: mass of snow         [kg/m2]
C     dMsn         :: mass of melting snow [kg/m2]
C     snowPrec     :: snow precipitation [kg/m2/s]
C     hNewSnow     :: fresh snow accumulation [m]
C     dhSnowMx     :: potential snow increase [m]
C     dhSnow       :: effective snow increase [m]
C     mIceDt       :: ground-ice growth rate (<- excess of snow) [kg/m2/s]
C     ageFac       :: snow aging factor [1]
      _RL grd_HeatCp, enthalpGrdW
      _RL fieldCapac, mWater
      _RL groundWnp1, grdWexcess, fractRunOff
      _RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL flxEngU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL flxEngL, temp_af, temp_bf, mPmE, enWfx, enGr1
      _RL mSnow, dMsn, snowPrec
      _RL hNewSnow, dhSnowMx, dhSnow, mIceDt, ageFac
      INTEGER i,j,k,kp1

#ifdef LAND_DEBUG
      LOGICAL dBug
      INTEGER iprt,jprt,lprt
      DATA iprt, jprt , lprt / 19 , 20 , 6 /
 1010 FORMAT(A,I3,1P4E11.3)
#endif

      IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN
C--   Step forward ground temperature:

      DO k=1,land_nLev
       kp1 = MIN(k+1,land_nLev)

       IF (k.EQ.1) THEN
        DO j=1,sNy
         DO i=1,sNx
           flxkup(i,j) = land_HeatFlx(i,j,bi,bj)
         ENDDO
        ENDDO
       ELSE
        DO j=1,sNy
         DO i=1,sNx
           flxkup(i,j) = flxkdw(i,j)
         ENDDO
        ENDDO
       ENDIF

       DO j=1,sNy
        DO i=1,sNx
         IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
C-     Thermal conductivity flux, lower interface (k,k+1):
          flxkdw(i,j) = land_grdLambda*
     &             ( land_groundT(i,j,k,bi,bj)
     &              -land_groundT(i,j,kp1,bi,bj) )
     &            *land_rec_dzC(kp1)

C-     Step forward ground enthalpy, level k :
          land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
     &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j))/land_dzF(k)

         ENDIF
        ENDDO
       ENDDO

      ENDDO
C--   step forward ground temperature: end
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF ( land_calc_grW .OR. land_calc_snow ) THEN
C--   Initialize run-off arrays.
        DO j=1,sNy
         DO i=1,sNx
           land_runOff(i,j,bi,bj) = 0. _d 0
           land_enRnOf(i,j,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
      ENDIF

#ifdef LAND_OLD_VERSION
      IF ( .TRUE. ) THEN
#else
      IF ( land_calc_grW ) THEN
#endif
C--   need (later on) ground temp. to be consistent with updated enthalpy:
        DO k=1,land_nLev
         DO j=1,sNy
          DO i=1,sNx
           IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
            mWater = land_rhoLiqW*land_waterCap
     &              *land_groundW(i,j,k,bi,bj)
            mWater = MAX( mWater, 0. _d 0 )
            grd_HeatCp = land_heatCs + land_CpWater*mWater
            temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
     &                                           / grd_HeatCp
            temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
            land_groundT(i,j,k,bi,bj) =
     &              MIN( temp_bf, MAX(temp_af, 0. _d 0) )
#ifdef LAND_DEBUG
            dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
            IF (dBug) write(6,1010)
     &        'LAND_STEPFWD: k,temp,af,bf=',
     &       k,land_groundT(i,j,k,bi,bj),temp_af,temp_bf
#endif
           ENDIF
          ENDDO
         ENDDO
        ENDDO
      ENDIF

      IF ( land_calc_snow ) THEN
C--   Step forward Snow thickness (also account for rain temperature)
        ageFac = 1. _d 0 - land_deltaT/timeSnowAge
        DO j=1,sNy
         DO i=1,sNx
          IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
           mPmE  = land_Pr_m_Ev(i,j,bi,bj)
           enWfx = land_EnWFlux(i,j,bi,bj)
           enGr1 = land_enthalp(i,j,1,bi,bj)*land_dzF(1)
#ifdef LAND_DEBUG
           dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
           IF (dBug) write(6,1010)
     &       'LAND_STEPFWD:mPmE,enWfx,enGr1/dt,hSnow=',0,
     &       mPmE,enWfx,enGr1/land_deltaT,land_hSnow(i,j,bi,bj)
#endif
C-    snow aging:
           land_snowAge(i,j,bi,bj) =
     &         ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac )
           IF ( enWfx.LT.0. ) THEN
C-    snow precip in excess ( > Evap of snow) or snow prec & Evap of Liq.Water:
C     => start to melt (until ground at freezing point) and then accumulate
            snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )
C-    snow accumulation cannot be larger that net precip
            snowPrec = MAX( 0. _d 0 ,
     &                      MIN( snowPrec*recip_Lfreez, mPmE ) )
            mPmE = mPmE - snowPrec
            flxEngU(i,j) = enWfx + land_Lfreez*snowPrec
            hNewSnow = land_deltaT * snowPrec / land_rhoSnow
C-    refresh snow age:
            land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj)
     &                          *EXP( -hNewSnow/hNewSnowAge )
C-    update snow thickness:
c           land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow
C     glacier & ice-sheet missing: excess of snow put directly into run-off
            dhSnowMx = MAX( 0. _d 0,
     &                      land_hMaxSnow - land_hSnow(i,j,bi,bj) )
            dhSnow = MIN( hNewSnow, dhSnowMx )
            land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow
            mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT
            land_runOff(i,j,bi,bj) =  mIceDt
            land_enRnOf(i,j,bi,bj) = -mIceDt*land_Lfreez
#ifdef LAND_DEBUG
            IF (dBug) write(6,1010)
     &        'LAND_STEPFWD: 3,snP,mPmE,hNsnw,hSnw=',
     &         3,snowPrec,mPmE,hNewSnow,land_hSnow(i,j,bi,bj)
#endif
           ELSE
C-    rain precip (whatever Evap is) or Evap of snow exceeds snow precip:
C     => snow melts or sublimates
c           snowMelt = MIN( enWfx*recip_Lfreez ,
c    &                 land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT )
            mSnow = land_hSnow(i,j,bi,bj)*land_rhoSnow
            dMsn = enWfx*recip_Lfreez*land_deltaT
            IF ( dMsn .GE. mSnow ) THEN
              dMsn = mSnow
              land_hSnow(i,j,bi,bj) = 0. _d 0
              flxEngU(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT
            ELSE
              flxEngU(i,j) = 0. _d 0
              land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj)
     &                              - dMsn / land_rhoSnow
            ENDIF
c           IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge
            mPmE = mPmE + dMsn/land_deltaT
#ifdef LAND_DEBUG
            IF (dBug) write(6,1010)
     &        'LAND_STEPFWD: 4,dMsn,mPmE,hSnw,enWfx=',
     &         4,dMsn,mPmE,land_hSnow(i,j,bi,bj),flxEngU(i,j)
#endif
           ENDIF
           flxkup(i,j) = mPmE/land_rhoLiqW
c          land_Pr_m_Ev(i,j,bi,bj) = mPmE
           IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 )
     &          land_snowAge(i,j,bi,bj) = 0. _d 0
C-    avoid negative (but very small, < 1.e-34) hSnow that occurs because
C      of truncation error. Might need to rewrite this part.
c          IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN
c             land_hSnow(i,j,bi,bj)   = 0. _d 0
c             land_snowAge(i,j,bi,bj) = 0. _d 0
c          ENDIF
          ENDIF
         ENDDO
        ENDDO
      ELSE
        DO j=1,sNy
         DO i=1,sNx
           flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW
           flxEngU(i,j) = 0. _d 0
         ENDDO
        ENDDO
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF (land_calc_grW) THEN
C--   Step forward ground Water:

      DO k=1,land_nLev
       IF (k.EQ.land_nLev) THEN
        kp1 = k
        fractRunOff = 1. _d 0
       ELSE
        kp1 = k+1
        fractRunOff = land_fractRunOff
       ENDIF
       fieldCapac = land_waterCap*land_dzF(k)

       DO j=1,sNy
        DO i=1,sNx
         IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
#ifdef LAND_DEBUG
          dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
#endif

#ifdef LAND_OLD_VERSION
          IF ( .TRUE. ) THEN
           IF ( k.EQ.land_nLev ) THEN
#else
          IF ( land_groundT(i,j,k,bi,bj).LT.0. _d 0 ) THEN
C-     Frozen level: only account for upper level fluxes
           IF ( flxkup(i,j) .LT. 0. _d 0 ) THEN
C-     Step forward soil moisture (& enthapy), level k :
            land_groundW(i,j,k,bi,bj) = land_groundW(i,j,k,bi,bj)
     &       + land_deltaT * flxkup(i,j) / fieldCapac
            IF ( land_calc_snow )
     &      land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
     &       + land_deltaT * flxEngU(i,j) / land_dzF(k)
           ELSE
C-     Frozen level: incoming water flux goes directly into run-off
            land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
     &                             + flxkup(i,j)*land_rhoLiqW
            land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
     &                             + flxEngU(i,j)
           ENDIF
C-     prepare fluxes for next level:
           flxkup(i,j)  = 0. _d 0
           flxEngU(i,j) = 0. _d 0

          ELSE

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-     Diffusion flux of water, lower interface (k,k+1):
           IF ( k.EQ.land_nLev .OR.
     &         land_groundT(i,j,kp1,bi,bj).LT.0. _d 0 ) THEN
#endif /* LAND_OLD_VERSION */
C-     no Diffusion of water if one level is frozen :
             flxkdw(i,j) = 0. _d 0
             flxEngL =     0. _d 0
           ELSE
             flxkdw(i,j) = fieldCapac*
     &                   ( land_groundW(i,j,k,bi,bj)
     &                    -land_groundW(i,j,kp1,bi,bj) )
     &                   / land_wTauDiff
C-     energy flux associated with water flux: take upwind Temp
             IF ( flxkdw(i,j).GE.0. ) THEN
              flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
     &                 *land_groundT(i,j,k,bi,bj)
             ELSE
              flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
     &                 *land_groundT(i,j,kp1,bi,bj)
             ENDIF
           ENDIF

C-     Step forward soil moisture, level k :
           groundWnp1 = land_groundW(i,j,k,bi,bj)
     &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac

#ifdef LAND_DEBUG
           IF(dBug)write(6,1010)'LAND_STEPFWD: grdW-1,fx_ku,kd,grdW-1='
     &      ,5,land_groundW(i,j,k,bi,bj)-1.,
     &         flxkup(i,j),flxkdw(i,j),groundWnp1-1.
#endif

C-     Water in excess will leave as run-off or go to level below
           land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1)
           grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )
     &                 *fieldCapac/land_deltaT

C-     Run off: fraction 1-fractRunOff enters level below
           land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
     &                        + fractRunOff*grdWexcess*land_rhoLiqW
C-     prepare fluxes for next level:
           flxkup(i,j) = flxkdw(i,j)
     &              + (1. _d 0-fractRunOff)*grdWexcess

           IF ( land_calc_snow ) THEN
            enthalpGrdW = land_rhoLiqW*land_CpWater
     &                   *land_groundT(i,j,k,bi,bj)
C--    Account for water fluxes in energy budget: update ground Enthalpy
            land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
     &         + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW
     &           )*land_deltaT/land_dzF(k)

            land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
     &                        + fractRunOff*grdWexcess*enthalpGrdW
C-     prepare fluxes for next level:
            flxEngU(i,j) = flxEngL
     &              + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW
           ENDIF
#ifdef LAND_DEBUG
           IF (dBug) write(6,1010) 'LAND_STEPFWD: Temp,FlxE,FlxW=',
     &      7, land_groundT(i,j,k,bi,bj), flxEngU(i,j), flxkup(i,j)
#endif
          ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef LAND_DEBUG
           IF (dBug) write(6,1010) 'LAND_STEPFWD: RO,enRO=',
     &      8, land_runOff(i,j,bi,bj),land_enRnOf(i,j,bi,bj)
#endif

         ENDIF
        ENDDO
       ENDDO

      ENDDO
C--   step forward ground Water: end
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF ( land_calc_grT ) THEN
C--   Compute ground temperature from enthalpy (if not already done):

       DO k=1,land_nLev
        DO j=1,sNy
         DO i=1,sNx
C-     Ground Heat capacity, layer k:
          mWater = land_rhoLiqW*land_waterCap
     &            *land_groundW(i,j,k,bi,bj)
          mWater = MAX( mWater, 0. _d 0 )
          grd_HeatCp = land_heatCs + land_CpWater*mWater
C         temperature below freezing:
          temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
     &                                         / grd_HeatCp
C         temperature above freezing:
          temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
#ifdef LAND_OLD_VERSION
          land_enthalp(i,j,k,bi,bj) =
     &          grd_HeatCp*land_groundT(i,j,k,bi,bj)
#else
          land_groundT(i,j,k,bi,bj) =
     &            MIN( temp_bf, MAX(temp_af, 0. _d 0) )
#endif
         ENDDO
        ENDDO
       ENDDO

       IF ( land_impl_grT ) THEN
        DO j=1,sNy
         DO i=1,sNx
          IF ( land_hSnow(i,j,bi,bj).GT.0. _d 0 ) THEN
           land_skinT(i,j,bi,bj) = MIN(land_skinT(i,j,bi,bj), 0. _d 0)
          ELSE
           land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
          ENDIF
         ENDDO
        ENDDO
       ELSE
        DO j=1,sNy
         DO i=1,sNx
           land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
         ENDDO
        ENDDO
       ENDIF

C--   Compute ground temperature: end
      ENDIF

#endif /* ALLOW_LAND */

      RETURN
      END
