Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:41:42 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
e0a2f8aec4 Jean*0001 #include "LAND_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: LAND_STEPFWD
                0005 C     !INTERFACE:
                0006       SUBROUTINE LAND_STEPFWD(
                0007      I                land_frc, bi, bj, myTime, myIter, myThid)
                0008 
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | S/R LAND_STEPFWD
                0012 C     | o Land model main S/R: step forward land variables
                0013 C     *==========================================================*
                0014 C     \ev
                0015 
                0016 C     !USES:
                0017       IMPLICIT NONE
                0018 
                0019 C     == Global variables ===
                0020 C-- size for MITgcm & Land package :
20cd0df114 Jean*0021 #include "LAND_SIZE.h"
e0a2f8aec4 Jean*0022 
                0023 #include "EEPARAMS.h"
                0024 #include "LAND_PARAMS.h"
                0025 #include "LAND_VARS.h"
                0026 
                0027 C     !INPUT/OUTPUT PARAMETERS:
                0028 C     == Routine arguments ==
                0029 C     land_frc :: land fraction [0-1]
                0030 C     bi,bj    :: Tile index
                0031 C     myTime   :: Current time of simulation ( s )
                0032 C     myIter   :: Current iteration number in simulation
                0033 C     myThid   :: Number of this instance of the routine
                0034       _RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0035       INTEGER bi, bj, myIter, myThid
                0036       _RL myTime
                0037 CEOP
                0038 
                0039 #ifdef ALLOW_LAND
                0040 C     == Local variables ==
                0041 C     i,j,k        :: loop counters
                0042 C     kp1          :: k+1
89992793c5 Jean*0043 C     grd_HeatCp   :: Heat capacity of the ground [J/m3/K]
b30f09edfc Jean*0044 C     enthalpGrdW  :: enthalpy of ground water [J/m3]
e0a2f8aec4 Jean*0045 C     fieldCapac   :: field capacity (of water) [m]
89992793c5 Jean*0046 C     mWater       :: water content of the ground [kg/m3]
b30f09edfc Jean*0047 C     groundWnp1   :: hold temporary future soil moisture []
e0a2f8aec4 Jean*0048 C     grdWexcess   :: ground water in excess [m/s]
b30f09edfc Jean*0049 C     fractRunOff  :: fraction of water in excess which leaves as runoff
89992793c5 Jean*0050 C     flxkup       :: downward flux of water, upper interface (k-1,k)
                0051 C     flxdwn       :: downward flux of water, lower interface (k,k+1)
b30f09edfc Jean*0052 C     flxEngU      :: downward energy flux associated with water flux (W/m2)
                0053 C                     upper interface (k-1,k)
                0054 C     flxEngL      :: downward energy flux associated with water flux (W/m2)
                0055 C                     lower interface (k,k+1)
89992793c5 Jean*0056 C     temp_af      :: ground temperature if above freezing
                0057 C     temp_bf      :: ground temperature if below freezing
                0058 C     mPmE         :: hold temporary (liquid) Precip minus Evap [kg/m2/s]
                0059 C     enWfx        :: hold temporary energy flux of Precip [W/m2]
                0060 C     enGr1        :: ground enthalpy of level 1  [J/m2]
                0061 C     mSnow        :: mass of snow         [kg/m2]
                0062 C     dMsn         :: mass of melting snow [kg/m2]
                0063 C     snowPrec     :: snow precipitation [kg/m2/s]
                0064 C     hNewSnow     :: fresh snow accumulation [m]
b30f09edfc Jean*0065 C     dhSnowMx     :: potential snow increase [m]
                0066 C     dhSnow       :: effective snow increase [m]
                0067 C     mIceDt       :: ground-ice growth rate (<- excess of snow) [kg/m2/s]
89992793c5 Jean*0068 C     ageFac       :: snow aging factor [1]
20cd0df114 Jean*0069       _RL grd_HeatCp, enthalpGrdW
b30f09edfc Jean*0070       _RL fieldCapac, mWater
                0071       _RL groundWnp1, grdWexcess, fractRunOff
e0a2f8aec4 Jean*0072       _RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073       _RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b30f09edfc Jean*0074       _RL flxEngU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075       _RL flxEngL, temp_af, temp_bf, mPmE, enWfx, enGr1
20cd0df114 Jean*0076       _RL mSnow, dMsn, snowPrec
b30f09edfc Jean*0077       _RL hNewSnow, dhSnowMx, dhSnow, mIceDt, ageFac
e0a2f8aec4 Jean*0078       INTEGER i,j,k,kp1
                0079 
4c9728b121 Jean*0080 #ifdef LAND_DEBUG
                0081       LOGICAL dBug
                0082       INTEGER iprt,jprt,lprt
                0083       DATA iprt, jprt , lprt / 19 , 20 , 6 /
                0084  1010 FORMAT(A,I3,1P4E11.3)
                0085 #endif
                0086 
89992793c5 Jean*0087       IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN
e0a2f8aec4 Jean*0088 C--   Step forward ground temperature:
                0089 
                0090       DO k=1,land_nLev
                0091        kp1 = MIN(k+1,land_nLev)
                0092 
                0093        IF (k.EQ.1) THEN
                0094         DO j=1,sNy
                0095          DO i=1,sNx
                0096            flxkup(i,j) = land_HeatFlx(i,j,bi,bj)
                0097          ENDDO
                0098         ENDDO
                0099        ELSE
                0100         DO j=1,sNy
                0101          DO i=1,sNx
                0102            flxkup(i,j) = flxkdw(i,j)
                0103          ENDDO
                0104         ENDDO
                0105        ENDIF
                0106 
                0107        DO j=1,sNy
                0108         DO i=1,sNx
                0109          IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
                0110 C-     Thermal conductivity flux, lower interface (k,k+1):
                0111           flxkdw(i,j) = land_grdLambda*
                0112      &             ( land_groundT(i,j,k,bi,bj)
                0113      &              -land_groundT(i,j,kp1,bi,bj) )
20cd0df114 Jean*0114      &            *land_rec_dzC(kp1)
                0115 
89992793c5 Jean*0116 C-     Step forward ground enthalpy, level k :
                0117           land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
                0118      &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j))/land_dzF(k)
e0a2f8aec4 Jean*0119 
                0120          ENDIF
                0121         ENDDO
                0122        ENDDO
                0123 
                0124       ENDDO
                0125 C--   step forward ground temperature: end
                0126       ENDIF
                0127 
89992793c5 Jean*0128 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0129 
1665818901 Jean*0130       IF ( land_calc_grW .OR. land_calc_snow ) THEN
b30f09edfc Jean*0131 C--   Initialize run-off arrays.
                0132         DO j=1,sNy
                0133          DO i=1,sNx
                0134            land_runOff(i,j,bi,bj) = 0. _d 0
                0135            land_enRnOf(i,j,bi,bj) = 0. _d 0
                0136          ENDDO
                0137         ENDDO
1665818901 Jean*0138       ENDIF
                0139 
                0140 #ifdef LAND_OLD_VERSION
                0141       IF ( .TRUE. ) THEN
                0142 #else
                0143       IF ( land_calc_grW ) THEN
                0144 #endif
89992793c5 Jean*0145 C--   need (later on) ground temp. to be consistent with updated enthalpy:
                0146         DO k=1,land_nLev
                0147          DO j=1,sNy
                0148           DO i=1,sNx
                0149            IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
                0150             mWater = land_rhoLiqW*land_waterCap
                0151      &              *land_groundW(i,j,k,bi,bj)
4c9728b121 Jean*0152             mWater = MAX( mWater, 0. _d 0 )
89992793c5 Jean*0153             grd_HeatCp = land_heatCs + land_CpWater*mWater
                0154             temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
                0155      &                                           / grd_HeatCp
                0156             temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
20cd0df114 Jean*0157             land_groundT(i,j,k,bi,bj) =
89992793c5 Jean*0158      &              MIN( temp_bf, MAX(temp_af, 0. _d 0) )
4c9728b121 Jean*0159 #ifdef LAND_DEBUG
                0160             dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
                0161             IF (dBug) write(6,1010)
                0162      &        'LAND_STEPFWD: k,temp,af,bf=',
                0163      &       k,land_groundT(i,j,k,bi,bj),temp_af,temp_bf
                0164 #endif
89992793c5 Jean*0165            ENDIF
                0166           ENDDO
                0167          ENDDO
                0168         ENDDO
                0169       ENDIF
                0170 
                0171       IF ( land_calc_snow ) THEN
                0172 C--   Step forward Snow thickness (also account for rain temperature)
                0173         ageFac = 1. _d 0 - land_deltaT/timeSnowAge
                0174         DO j=1,sNy
                0175          DO i=1,sNx
                0176           IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
                0177            mPmE  = land_Pr_m_Ev(i,j,bi,bj)
                0178            enWfx = land_EnWFlux(i,j,bi,bj)
                0179            enGr1 = land_enthalp(i,j,1,bi,bj)*land_dzF(1)
4c9728b121 Jean*0180 #ifdef LAND_DEBUG
                0181            dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
                0182            IF (dBug) write(6,1010)
                0183      &       'LAND_STEPFWD:mPmE,enWfx,enGr1/dt,hSnow=',0,
                0184      &       mPmE,enWfx,enGr1/land_deltaT,land_hSnow(i,j,bi,bj)
                0185 #endif
89992793c5 Jean*0186 C-    snow aging:
20cd0df114 Jean*0187            land_snowAge(i,j,bi,bj) =
89992793c5 Jean*0188      &         ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac )
                0189            IF ( enWfx.LT.0. ) THEN
ecb76f8691 Jean*0190 C-    snow precip in excess ( > Evap of snow) or snow prec & Evap of Liq.Water:
20cd0df114 Jean*0191 C     => start to melt (until ground at freezing point) and then accumulate
89992793c5 Jean*0192             snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )
ecb76f8691 Jean*0193 C-    snow accumulation cannot be larger that net precip
                0194             snowPrec = MAX( 0. _d 0 ,
                0195      &                      MIN( snowPrec*recip_Lfreez, mPmE ) )
89992793c5 Jean*0196             mPmE = mPmE - snowPrec
b30f09edfc Jean*0197             flxEngU(i,j) = enWfx + land_Lfreez*snowPrec
89992793c5 Jean*0198             hNewSnow = land_deltaT * snowPrec / land_rhoSnow
                0199 C-    refresh snow age:
                0200             land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj)
                0201      &                          *EXP( -hNewSnow/hNewSnowAge )
b30f09edfc Jean*0202 C-    update snow thickness:
                0203 c           land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow
                0204 C     glacier & ice-sheet missing: excess of snow put directly into run-off
20cd0df114 Jean*0205             dhSnowMx = MAX( 0. _d 0,
b30f09edfc Jean*0206      &                      land_hMaxSnow - land_hSnow(i,j,bi,bj) )
                0207             dhSnow = MIN( hNewSnow, dhSnowMx )
                0208             land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow
                0209             mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT
20cd0df114 Jean*0210             land_runOff(i,j,bi,bj) =  mIceDt
b30f09edfc Jean*0211             land_enRnOf(i,j,bi,bj) = -mIceDt*land_Lfreez
4c9728b121 Jean*0212 #ifdef LAND_DEBUG
                0213             IF (dBug) write(6,1010)
                0214      &        'LAND_STEPFWD: 3,snP,mPmE,hNsnw,hSnw=',
                0215      &         3,snowPrec,mPmE,hNewSnow,land_hSnow(i,j,bi,bj)
                0216 #endif
89992793c5 Jean*0217            ELSE
ecb76f8691 Jean*0218 C-    rain precip (whatever Evap is) or Evap of snow exceeds snow precip:
89992793c5 Jean*0219 C     => snow melts or sublimates
                0220 c           snowMelt = MIN( enWfx*recip_Lfreez ,
                0221 c    &                 land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT )
                0222             mSnow = land_hSnow(i,j,bi,bj)*land_rhoSnow
                0223             dMsn = enWfx*recip_Lfreez*land_deltaT
                0224             IF ( dMsn .GE. mSnow ) THEN
                0225               dMsn = mSnow
                0226               land_hSnow(i,j,bi,bj) = 0. _d 0
b30f09edfc Jean*0227               flxEngU(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT
89992793c5 Jean*0228             ELSE
b30f09edfc Jean*0229               flxEngU(i,j) = 0. _d 0
89992793c5 Jean*0230               land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj)
20cd0df114 Jean*0231      &                              - dMsn / land_rhoSnow
89992793c5 Jean*0232             ENDIF
                0233 c           IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge
                0234             mPmE = mPmE + dMsn/land_deltaT
4c9728b121 Jean*0235 #ifdef LAND_DEBUG
                0236             IF (dBug) write(6,1010)
                0237      &        'LAND_STEPFWD: 4,dMsn,mPmE,hSnw,enWfx=',
                0238      &         4,dMsn,mPmE,land_hSnow(i,j,bi,bj),flxEngU(i,j)
                0239 #endif
89992793c5 Jean*0240            ENDIF
                0241            flxkup(i,j) = mPmE/land_rhoLiqW
                0242 c          land_Pr_m_Ev(i,j,bi,bj) = mPmE
                0243            IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 )
                0244      &          land_snowAge(i,j,bi,bj) = 0. _d 0
20cd0df114 Jean*0245 C-    avoid negative (but very small, < 1.e-34) hSnow that occurs because
89992793c5 Jean*0246 C      of truncation error. Might need to rewrite this part.
                0247 c          IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN
                0248 c             land_hSnow(i,j,bi,bj)   = 0. _d 0
                0249 c             land_snowAge(i,j,bi,bj) = 0. _d 0
                0250 c          ENDIF
                0251           ENDIF
                0252          ENDDO
                0253         ENDDO
                0254       ELSE
                0255         DO j=1,sNy
                0256          DO i=1,sNx
                0257            flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW
b30f09edfc Jean*0258            flxEngU(i,j) = 0. _d 0
89992793c5 Jean*0259          ENDDO
                0260         ENDDO
                0261       ENDIF
                0262 
                0263 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0264 
e0a2f8aec4 Jean*0265       IF (land_calc_grW) THEN
                0266 C--   Step forward ground Water:
                0267 
                0268       DO k=1,land_nLev
                0269        IF (k.EQ.land_nLev) THEN
                0270         kp1 = k
                0271         fractRunOff = 1. _d 0
                0272        ELSE
                0273         kp1 = k+1
                0274         fractRunOff = land_fractRunOff
                0275        ENDIF
                0276        fieldCapac = land_waterCap*land_dzF(k)
                0277 
                0278        DO j=1,sNy
                0279         DO i=1,sNx
                0280          IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
4c9728b121 Jean*0281 #ifdef LAND_DEBUG
                0282           dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
                0283 #endif
b30f09edfc Jean*0284 
                0285 #ifdef LAND_OLD_VERSION
                0286           IF ( .TRUE. ) THEN
                0287            IF ( k.EQ.land_nLev ) THEN
                0288 #else
                0289           IF ( land_groundT(i,j,k,bi,bj).LT.0. _d 0 ) THEN
                0290 C-     Frozen level: only account for upper level fluxes
                0291            IF ( flxkup(i,j) .LT. 0. _d 0 ) THEN
                0292 C-     Step forward soil moisture (& enthapy), level k :
                0293             land_groundW(i,j,k,bi,bj) = land_groundW(i,j,k,bi,bj)
                0294      &       + land_deltaT * flxkup(i,j) / fieldCapac
                0295             IF ( land_calc_snow )
                0296      &      land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
                0297      &       + land_deltaT * flxEngU(i,j) / land_dzF(k)
                0298            ELSE
                0299 C-     Frozen level: incoming water flux goes directly into run-off
                0300             land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
20cd0df114 Jean*0301      &                             + flxkup(i,j)*land_rhoLiqW
b30f09edfc Jean*0302             land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
                0303      &                             + flxEngU(i,j)
                0304            ENDIF
                0305 C-     prepare fluxes for next level:
                0306            flxkup(i,j)  = 0. _d 0
                0307            flxEngU(i,j) = 0. _d 0
                0308 
                0309           ELSE
                0310 
                0311 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
89992793c5 Jean*0312 C-     Diffusion flux of water, lower interface (k,k+1):
b30f09edfc Jean*0313            IF ( k.EQ.land_nLev .OR.
                0314      &         land_groundT(i,j,kp1,bi,bj).LT.0. _d 0 ) THEN
                0315 #endif /* LAND_OLD_VERSION */
                0316 C-     no Diffusion of water if one level is frozen :
                0317              flxkdw(i,j) = 0. _d 0
                0318              flxEngL =     0. _d 0
                0319            ELSE
                0320              flxkdw(i,j) = fieldCapac*
                0321      &                   ( land_groundW(i,j,k,bi,bj)
                0322      &                    -land_groundW(i,j,kp1,bi,bj) )
                0323      &                   / land_wTauDiff
                0324 C-     energy flux associated with water flux: take upwind Temp
                0325              IF ( flxkdw(i,j).GE.0. ) THEN
                0326               flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
                0327      &                 *land_groundT(i,j,k,bi,bj)
                0328              ELSE
                0329               flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
                0330      &                 *land_groundT(i,j,kp1,bi,bj)
                0331              ENDIF
                0332            ENDIF
20cd0df114 Jean*0333 
e0a2f8aec4 Jean*0334 C-     Step forward soil moisture, level k :
b30f09edfc Jean*0335            groundWnp1 = land_groundW(i,j,k,bi,bj)
89992793c5 Jean*0336      &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac
b30f09edfc Jean*0337 
4c9728b121 Jean*0338 #ifdef LAND_DEBUG
                0339            IF(dBug)write(6,1010)'LAND_STEPFWD: grdW-1,fx_ku,kd,grdW-1='
                0340      &      ,5,land_groundW(i,j,k,bi,bj)-1.,
                0341      &         flxkup(i,j),flxkdw(i,j),groundWnp1-1.
                0342 #endif
                0343 
b30f09edfc Jean*0344 C-     Water in excess will leave as run-off or go to level below
                0345            land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1)
                0346            grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )
                0347      &                 *fieldCapac/land_deltaT
e0a2f8aec4 Jean*0348 
                0349 C-     Run off: fraction 1-fractRunOff enters level below
b30f09edfc Jean*0350            land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
20cd0df114 Jean*0351      &                        + fractRunOff*grdWexcess*land_rhoLiqW
b30f09edfc Jean*0352 C-     prepare fluxes for next level:
                0353            flxkup(i,j) = flxkdw(i,j)
                0354      &              + (1. _d 0-fractRunOff)*grdWexcess
                0355 
                0356            IF ( land_calc_snow ) THEN
                0357             enthalpGrdW = land_rhoLiqW*land_CpWater
                0358      &                   *land_groundT(i,j,k,bi,bj)
                0359 C--    Account for water fluxes in energy budget: update ground Enthalpy
                0360             land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
20cd0df114 Jean*0361      &         + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW
b30f09edfc Jean*0362      &           )*land_deltaT/land_dzF(k)
                0363 
                0364             land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
                0365      &                        + fractRunOff*grdWexcess*enthalpGrdW
89992793c5 Jean*0366 C-     prepare fluxes for next level:
b30f09edfc Jean*0367             flxEngU(i,j) = flxEngL
                0368      &              + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW
                0369            ENDIF
4c9728b121 Jean*0370 #ifdef LAND_DEBUG
                0371            IF (dBug) write(6,1010) 'LAND_STEPFWD: Temp,FlxE,FlxW=',
                0372      &      7, land_groundT(i,j,k,bi,bj), flxEngU(i,j), flxkup(i,j)
                0373 #endif
b30f09edfc Jean*0374           ENDIF
                0375 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
4c9728b121 Jean*0376 #ifdef LAND_DEBUG
                0377            IF (dBug) write(6,1010) 'LAND_STEPFWD: RO,enRO=',
                0378      &      8, land_runOff(i,j,bi,bj),land_enRnOf(i,j,bi,bj)
                0379 #endif
89992793c5 Jean*0380 
e0a2f8aec4 Jean*0381          ENDIF
                0382         ENDDO
                0383        ENDDO
                0384 
                0385       ENDDO
                0386 C--   step forward ground Water: end
                0387       ENDIF
                0388 
89992793c5 Jean*0389 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0390 
b30f09edfc Jean*0391       IF ( land_calc_grT ) THEN
                0392 C--   Compute ground temperature from enthalpy (if not already done):
89992793c5 Jean*0393 
                0394        DO k=1,land_nLev
                0395         DO j=1,sNy
                0396          DO i=1,sNx
                0397 C-     Ground Heat capacity, layer k:
                0398           mWater = land_rhoLiqW*land_waterCap
                0399      &            *land_groundW(i,j,k,bi,bj)
4c9728b121 Jean*0400           mWater = MAX( mWater, 0. _d 0 )
89992793c5 Jean*0401           grd_HeatCp = land_heatCs + land_CpWater*mWater
                0402 C         temperature below freezing:
                0403           temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
                0404      &                                         / grd_HeatCp
                0405 C         temperature above freezing:
                0406           temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
                0407 #ifdef LAND_OLD_VERSION
20cd0df114 Jean*0408           land_enthalp(i,j,k,bi,bj) =
89992793c5 Jean*0409      &          grd_HeatCp*land_groundT(i,j,k,bi,bj)
                0410 #else
20cd0df114 Jean*0411           land_groundT(i,j,k,bi,bj) =
89992793c5 Jean*0412      &            MIN( temp_bf, MAX(temp_af, 0. _d 0) )
                0413 #endif
                0414          ENDDO
                0415         ENDDO
                0416        ENDDO
                0417 
                0418        IF ( land_impl_grT ) THEN
                0419         DO j=1,sNy
                0420          DO i=1,sNx
                0421           IF ( land_hSnow(i,j,bi,bj).GT.0. _d 0 ) THEN
                0422            land_skinT(i,j,bi,bj) = MIN(land_skinT(i,j,bi,bj), 0. _d 0)
                0423           ELSE
                0424            land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
                0425           ENDIF
                0426          ENDDO
                0427         ENDDO
20cd0df114 Jean*0428        ELSE
89992793c5 Jean*0429         DO j=1,sNy
                0430          DO i=1,sNx
                0431            land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
                0432          ENDDO
                0433         ENDDO
                0434        ENDIF
                0435 
                0436 C--   Compute ground temperature: end
                0437       ENDIF
                0438 
e0a2f8aec4 Jean*0439 #endif /* ALLOW_LAND */
                0440 
                0441       RETURN
                0442       END