#include "LAND_OPTIONS.h"

CBOP
C     !ROUTINE: LAND_IMPL_TEMP
C     !INTERFACE:
      SUBROUTINE LAND_IMPL_TEMP(
     I                land_frc, 
     I                dTskin, sFlx,
     O                dTsurf,
     I                bi, bj, myTime, myIter, myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R LAND_IMPL_TEMP
C     | o solve ground temp. and surface temp. implicitly 
C     *==========================================================*
C     | o account for snow layer (with no heat capacity) 
C     |   and ground water freezing/melting
C     | o surf. heat flux is linearly dependent on surf. temp.
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     dTskin   :: temp. correction for daily-cycle heating [K]
C     sFlx     :: surf. heat flux (+=down) function of surface temp. Ts:
C                 0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
C     dTsurf   :: surf. temp adjusment: Ts^n+1 - Ts^n
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)
      _RL dTskin(sNx,sNy), sFlx(sNx,sNy,0:2) 
      _RL dTsurf(sNx,sNy)
      INTEGER bi, bj, myIter, myThid
      _RL myTime
CEOP

#ifdef ALLOW_LAND

C     == Local variables ==
C-  local variables used in solving the ground temp. implicitly :
C     aLoc         :: ground Conductivity * delT / delZ_12    [J/K]
C     bLoc         :: minus surf. flux derivative ./. Ts      [W/m2/K]
C     cLoc         :: temporary value for level.1 heat capacity [J/m2/K]
C     eLoc         :: temporary value for level.1 ground enthalpy [J/m2]
C     fLoc         :: temporary value for surface heat flux [W/m2]
C     alpha        :: snow thicknes / snow conductivity [m2.K/W]
C     beta         :: local coeff = 1/(1+alpha*bLoc)    [1]
C     tSurf        :: surf.  temperature   [oC]
C     tg           :: ground temperature   [oC]    (2 levels)
C     eg           :: ground enthalpy      [J/m2]  (2 levels)
C     cg           :: ground heat capacity [J/m2/K](2 levels)
C     mW           :: ground water mass    [kg/m2] (2 levels)
C     temp_af      :: ground temperature if above freezing
C     temp_bf      :: ground temperature if below freezing
C     mSnow        :: mass of snow         [kg/m2]
C     dMsn         :: mass of melting snow [kg/m2]
C     delT         :: time step            [s]
C     mSnEpsil     :: small snow mass      [kg/m2]
C     i,j,k        :: loop counters
C     msgBuf       :: Informational/error meesage buffer
C     tmpFlag      :: temp. flag, =.T. until found final groung temp
      _RL aLoc, bLoc, cLoc, eLoc, fLoc, alpha, beta
      _RL tSurf, tg(land_nLev), eg(land_nLev)
      _RL cg(land_nLev), mW(land_nLev)
      _RL temp_af, temp_bf
      _RL mSnow, dMsn, delT
      _RL mSnEpsil
      _RL tmp1, tmp2
      INTEGER i,j,k
      LOGICAL tmpFlag

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

      DATA    mSnEpsil / 1. _d -6 /

C-------------------------------------------------------------------------
C  solve implicitly the coupled 3 eq. (time level n+1 omitted) :
C    1a : if hs=0 : Flx = Flx^o + d.F/dT*(Ts - Ts^n) & Ts=Tg1
C    1b : if hs>0 : Flx = (Ts-Tg1)*Ks/hs =< Flx^o + d.F/dT*(Ts - Ts^n)
C         & difference used to melt the snow, maintaining Ts=0 
C    2  : Eg1 - Eg1^n  = delT*Flx - (lambda*delT/delZ_12)*(Tg1-Tg2)
C    3  : Eg2 - Eg2^n  = (lambda*delT/delZ_12)*(Tg1-Tg2)
C    were  lambda = ground Conductivity ; Ks = snow Conductivity
C          k=1,2: Eg_k = Cg_k * Tg_k - Lfreez * mIce_k
C
C  using local variables:
C   a = lambda*delT/delZ_12 ; b = - d.F/dT ;  f = Flx^o + b*Ts^n
C                         alpha = hs/Ks ;  beta = 1/(1+alpha*b)
C  3.eq system becomes: 
C   o if Ts*hs =< 0
C     1a,b:  Ts = ( Tg1 + alpha*F)*beta
C      2  : Eg1 + a*(Tg1-Tg2) + (b*delT)*beta*Tg1 = Eg1^n + delT*f*beta
C      3  : Eg2 + a*(Tg2-Tg1) = Eg2^n
C   o else: 
C     1.b : Ts=0 , f = Flx(ts=0) ; snowMelt = (f + Tg1/alpha)/Lfreez
C      2  : Eg1 + a*(Tg1-Tg2) + (delT/alpha)*Tg1 = Eg1^n
C      3  : Eg2 + a*(Tg2-Tg1) = Eg2^n
C-------------------------------------------------------------------------

C---  Solve implicitely for ground temp. & surface temp

      delT = land_deltaT
      aLoc = land_grdLambda*land_deltaT*land_rec_dzC(2)
      DO j=1,sNy
       DO i=1,sNx
        IF ( land_frc(i,j,bi,bj).GT.0. ) THEN

C--   initialise local variables
          tmpFlag = .TRUE.
          tSurf = land_skinT(i,j,bi,bj)
          mSnow = land_rhoSnow*land_hSnow(i,j,bi,bj)
          bLoc  = -sFlx(i,j,2)
          fLoc  = sFlx(i,j,1)+bLoc*tSurf
          alpha = land_hSnow(i,j,bi,bj)/diffKsnow
          beta  = 1. _d 0/(1. _d 0+alpha*bLoc)
          DO k=1,land_nLev
            eg(k) = land_dzF(k)*land_enthalp(i,j,k,bi,bj)
            mW(k) = land_dzF(k)*land_groundW(i,j,k,bi,bj)
     &                         *land_waterCap*land_rhoLiqW
            mW(k) = MAX( mW(k), 0. _d 0 )
            cg(k) = land_dzF(k)*land_heatCs + mW(k)*land_CpWater
            tg(k) = land_groundT(i,j,k,bi,bj)
          ENDDO
#ifdef LAND_DEBUG
          dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
          IF (dBug) write(6,1010)
     &         'LAND_IMPL_TEMP: 0 , ts,tg1&2,mSw=',0,tSurf,tg,mSnow
          IF (dBug) write(6,1010)
     &         'LAND_IMPL_TEMP: 0 , sFlx=', 0,(sFlx(i,j,k),k=0,2)
#endif

C---   Solve for temp as if no freezing/melting was occuring :
          tg(1) = ( cg(1)*tg(1) + fLoc*delT*beta
     &                          + cg(2)*tg(2)*aLoc/(cg(2)+aLoc)
     &            )
     &          / ( cg(1) + aLoc + bLoc*delT*beta
     &                           - aLoc*aLoc/(cg(2)+aLoc) 
     &            )
          tg(2) = ( cg(2)*tg(2) + aLoc*tg(1) ) / (cg(2)+aLoc)
          tSurf = ( tg(1) + alpha*fLoc ) * beta
          
#ifdef LAND_DEBUG
          IF (dBug) write(6,1010)
     &         'LAND_IMPL_TEMP: 1 , ts,tg1&2,mW1=',1,tSurf,tg,mW(1)
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C---   If melting/freezing (top of snow layer, ground water level 1 or 2)
C      set corresponding temp to freezing point and update enthalpy 
C--------------

          IF ( tg(2)*land_groundT(i,j,2,bi,bj) .LE. 0. _d 0 
     &         .AND. tmpFlag .AND. tSurf*mSnow .LE. 0. _d 0 ) THEN
C--    freezing/melting in level 2: set Tg2 to freezing point
           tmp1 = tg(1)
           tmp2 = tg(2)
           tg(2) = 0.
           eLoc = eg(1) + fLoc*delT*beta
           cLoc = cg(1) + aLoc + bLoc*delT*beta
           temp_bf = (eLoc+land_Lfreez*mW(1))/cLoc
           temp_af =  eLoc/cLoc
           tg(1) = MIN( temp_bf, MAX(temp_af, 0. _d 0) )
           tSurf = ( tg(1) + alpha*fLoc ) * beta
           IF ( tSurf*mSnow .LE. 0. _d 0 ) THEN
             tmpFlag = .FALSE.
             eg(1) = eLoc - (aLoc + bLoc*delT*beta)*tg(1)
             eg(2) = eg(2) + aLoc*tg(1) 
#ifdef LAND_DEBUG
           ELSEIF ( debugFlag ) THEN
             WRITE(msgBuf,'(A,2I4,2I3,I10)')
     &             'LAND_IMPL_TEMP: i,j,bi,bj,Iter=',
     &             i,j,bi,bj,myIter
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
             WRITE(msgBuf,'(A,1P4E12.4)')
     &             'LAND_IMPL_TEMP: groundT,t2,ts=',
     &             land_groundT(i,j,1,bi,bj),land_groundT(i,j,2,bi,bj),
     &             tmp2,(tmp1+alpha*fLoc)*beta
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
             WRITE(msgBuf,'(A,1P4E12.4)')
     &             'LAND_IMPL_TEMP: Tg,tSurf,mSnw=',
     &              tg,tSurf,mSnow
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
             WRITE(msgBuf,'(A,1P4E14.6)')
     &             'LAND_IMPL_TEMP: eg,mW=', eg,mW
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
#endif
           ENDIF

C-  if tg2_new*tg2_old < 0 : end
          ENDIF

C--------------

          IF ( tg(1)*land_groundT(i,j,1,bi,bj) .LE. 0. _d 0 
     &         .AND. tmpFlag .AND. tSurf*mSnow .LE. 0. _d 0 ) THEN
C--    freezing/melting in level 1: set Tg1 to freezing point
           tmp1 = tg(1)
           tg(1) = 0.
           tg(2) = cg(2)*tg(2) / (cg(2)+aLoc)
           tSurf = alpha*fLoc * beta
           IF ( tSurf*mSnow .LE. 0. _d 0 ) THEN
             tmpFlag = .FALSE.
             eg(2) = eg(2) - aLoc*tg(2) 
             eg(1) = eg(1) + aLoc*tg(2) + fLoc*delT*beta
             IF ( eg(1)*mSnow .GT. 0. _d 0 ) THEN
C-           melt snow from bottom
              dMsn = MIN( mSnow, eg(1)*recip_Lfreez )
              land_Pr_m_Ev(i,j,bi,bj) = dMsn/delT
              land_hSnow(i,j,bi,bj) = (mSnow - dMsn)/land_rhoSnow 
              eg(1) = eg(1) - dMsn*land_Lfreez
#ifdef LAND_DEBUG
              IF (dBug) write(6,1010)
     &         'LAND_IMPL_TEMP: Bot-Melt : dMsn,dEg1,eg1=',
     &         1, dMsn, -dMsn*land_Lfreez, eg(1) 
#endif
             ENDIF
#ifdef LAND_DEBUG
           ELSEIF ( debugFlag ) THEN
             WRITE(msgBuf,'(A,2I4,2I3,I10)')
     &             'LAND_IMPL_TEMP: i,j,bi,bj,Iter=',
     &             i,j,bi,bj,myIter
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
             WRITE(msgBuf,'(A,4F11.6)')
     &             'LAND_IMPL_TEMP: groundT,t1,ts=',
     &             land_groundT(i,j,1,bi,bj),land_groundT(i,j,2,bi,bj),
     &             tmp1,(tmp1+alpha*fLoc)*beta
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
             WRITE(msgBuf,'(A,4F12.7)')
     &             'LAND_IMPL_TEMP: Tg,tSurf,mSnow=',
     &             tg,tSurf,mSnow
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
             WRITE(msgBuf,'(A,1P4E14.6)')
     &             'LAND_IMPL_TEMP: eg,mW=', eg,mW
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
             WRITE(msgBuf,'(A)')  
     &             'LAND_IMPL_TEMP: snow with ts > 0 ! but continue'
             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
#endif
           ENDIF

C-  if tg1_new*tg1_old < 0 : end
          ENDIF

C--------------

          IF ( tmpFlag .AND. tSurf*mSnow .GT. 0. _d 0 ) THEN
C--    snow is melting at the surface: set ts=0 & use fixed heat flux Flx(ts=0)
#ifdef LAND_DEBUG
              IF (dBug) write(6,1010)
     &         'LAND_IMPL_TEMP: Top-Melt : fx0, fx1, fx1-b*Ts =',
     &         1, sFlx(i,j,0), fLoc, fLoc-bLoc*tSurf
#endif
           tSurf = 0. _d 0
           fLoc = sFlx(i,j,0)
           dTsurf(i,j) = 1000.
           tg(1) = land_groundT(i,j,1,bi,bj)
           tg(2) = land_groundT(i,j,2,bi,bj)

           eLoc = cg(1)*tg(1)
     &          + delT*fLoc - land_Lfreez*mSnow + aLoc*tg(2)
           IF ( eLoc .GT. 0. _d 0 .OR. mSnow.LT.mSnEpsil ) THEN
C-     all snow melt: do not solve diffusion of heat in snow layer
C      but put surf. heat flux directly to 1rst level and set tg1=0
             dMsn = mSnow
             tg(1) = 0. _d 0
             tg(2) = cg(2)*tg(2) / (cg(2)+aLoc)
           ELSE
C-     solve diffusion of heat in snow layer ; heat flux difference 
C      (surf.Flx - diffusion.Flx) is used to melt the snow from top.
            tg(1) = ( cg(1)*tg(1) + cg(2)*tg(2)*aLoc/(cg(2)+aLoc) )
     &       / ( cg(1)+aLoc + delT/alpha - aLoc*aLoc/(cg(2)+aLoc) )
            tg(2) = ( cg(2)*tg(2) + aLoc*tg(1) ) / (cg(2)+aLoc)
            IF ( tg(2)*land_groundT(i,j,2,bi,bj).LE.0. _d 0 ) THEN
              tg(2) = 0.
              tg(1) = cg(1)*tg(1) / ( cg(1)+aLoc + delT/alpha )
            ELSEIF ( tg(1)*land_groundT(i,j,1,bi,bj).LE.0. _d 0 ) THEN
              tg(1) = 0.
              tg(2) = cg(2)*tg(2) / (cg(2)+aLoc)
            ENDIF
            dMsn = ( fLoc+tg(1)/alpha )*delT*recip_Lfreez
#ifdef LAND_DEBUG
              IF (dBug) write(6,1010)
     &         'LAND_IMPL_TEMP: Surf-Melt: dMsn,fLoc,tg1/alpha=',
     &          2, dMsn, fLoc,tg(1)/alpha
#endif
C-  note: Fx0 < -tg(1)/alpha can happen (due to non-linearity in Fx(Ts)), 
C-     => do not melt nor accumulate snow but put d.Flx in Eg1
            dMsn = MIN( MAX(dMsn, 0. _d 0), mSnow )
           ENDIF
           tmpFlag = .FALSE.
           eg(2) = eg(2) + aLoc*(tg(1)-tg(2))
           eg(1) = eg(1) - aLoc*(tg(1)-tg(2))
     &           + delT*fLoc - land_Lfreez*dMsn
           land_Pr_m_Ev(i,j,bi,bj) = dMsn/delT
           land_hSnow(i,j,bi,bj) = (mSnow - dMsn)/land_rhoSnow 

C-  if ts*hSnow > 0 , else:
          ELSEIF ( tmpFlag ) THEN
C--   snow is not melting & no freezing/melting in ground level 1 & 2
           eg(2) = eg(2) + aLoc*(tg(1)-tg(2))
           eg(1) = eg(1) - aLoc*(tg(1)-tg(2))
     &           + delT*(fLoc-bLoc*Tsurf)
           tmpFlag = .FALSE.
          ENDIF

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

C---  Save new values :
          IF ( dTsurf(i,j) .LE. 999. ) 
     &         dTsurf(i,j) = tSurf - land_skinT(i,j,bi,bj)
          land_skinT(i,j,bi,bj) = tSurf
          DO k=1,land_nLev
            land_enthalp(i,j,k,bi,bj) = eg(k)/land_dzF(k)
            land_groundT(i,j,k,bi,bj) = tg(k)
          ENDDO

#ifdef LAND_DEBUG
          IF (dBug) write(6,1010) 'LAND_IMPL_TEMP: 9, ts,tg1&2,dTs=',9,
     &         tSurf, tg, dTsurf(i,j)
          IF (dBug) write(6,1010) 'LAND_IMPL_TEMP: 9, Eg1,Eg2,mPmE=',9,
     &       (land_enthalp(i,j,k,bi,bj),k=1,2), land_Pr_m_Ev(i,j,bi,bj)
#endif

        ENDIF
       ENDDO
      ENDDO

#endif /* ALLOW_LAND */

      RETURN
      END
