Back to home page

MITgcm
 
 

    

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

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
89992793c5 Jean*0001 #include "LAND_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: LAND_IMPL_TEMP
                0005 C     !INTERFACE:
                0006       SUBROUTINE LAND_IMPL_TEMP(
                0007      I                land_frc, 
                0008      I                dTskin, sFlx,
                0009      O                dTsurf,
                0010      I                bi, bj, myTime, myIter, myThid)
                0011 
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | S/R LAND_IMPL_TEMP
                0015 C     | o solve ground temp. and surface temp. implicitly 
                0016 C     *==========================================================*
                0017 C     | o account for snow layer (with no heat capacity) 
                0018 C     |   and ground water freezing/melting
                0019 C     | o surf. heat flux is linearly dependent on surf. temp.
                0020 C     *==========================================================*
                0021 C     \ev
                0022 
                0023 C     !USES:
                0024       IMPLICIT NONE
                0025 
                0026 C     == Global variables ===
                0027 C-- size for MITgcm & Land package :
                0028 #include "LAND_SIZE.h" 
                0029 
                0030 #include "EEPARAMS.h"
                0031 #include "LAND_PARAMS.h"
                0032 #include "LAND_VARS.h"
                0033 
                0034 C     !INPUT/OUTPUT PARAMETERS:
                0035 C     == Routine arguments ==
                0036 C     land_frc :: land fraction [0-1]
                0037 C     dTskin   :: temp. correction for daily-cycle heating [K]
                0038 C     sFlx     :: surf. heat flux (+=down) function of surface temp. Ts:
                0039 C                 0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
                0040 C     dTsurf   :: surf. temp adjusment: Ts^n+1 - Ts^n
                0041 C     bi,bj    :: Tile index
                0042 C     myTime   :: Current time of simulation ( s )
                0043 C     myIter   :: Current iteration number in simulation
                0044 C     myThid   :: Number of this instance of the routine
                0045       _RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0046       _RL dTskin(sNx,sNy), sFlx(sNx,sNy,0:2) 
                0047       _RL dTsurf(sNx,sNy)
                0048       INTEGER bi, bj, myIter, myThid
                0049       _RL myTime
                0050 CEOP
                0051 
                0052 #ifdef ALLOW_LAND
                0053 
                0054 C     == Local variables ==
                0055 C-  local variables used in solving the ground temp. implicitly :
                0056 C     aLoc         :: ground Conductivity * delT / delZ_12    [J/K]
                0057 C     bLoc         :: minus surf. flux derivative ./. Ts      [W/m2/K]
                0058 C     cLoc         :: temporary value for level.1 heat capacity [J/m2/K]
                0059 C     eLoc         :: temporary value for level.1 ground enthalpy [J/m2]
                0060 C     fLoc         :: temporary value for surface heat flux [W/m2]
                0061 C     alpha        :: snow thicknes / snow conductivity [m2.K/W]
                0062 C     beta         :: local coeff = 1/(1+alpha*bLoc)    [1]
                0063 C     tSurf        :: surf.  temperature   [oC]
                0064 C     tg           :: ground temperature   [oC]    (2 levels)
                0065 C     eg           :: ground enthalpy      [J/m2]  (2 levels)
                0066 C     cg           :: ground heat capacity [J/m2/K](2 levels)
                0067 C     mW           :: ground water mass    [kg/m2] (2 levels)
                0068 C     temp_af      :: ground temperature if above freezing
                0069 C     temp_bf      :: ground temperature if below freezing
                0070 C     mSnow        :: mass of snow         [kg/m2]
                0071 C     dMsn         :: mass of melting snow [kg/m2]
                0072 C     delT         :: time step            [s]
                0073 C     mSnEpsil     :: small snow mass      [kg/m2]
                0074 C     i,j,k        :: loop counters
                0075 C     msgBuf       :: Informational/error meesage buffer
4c9728b121 Jean*0076 C     tmpFlag      :: temp. flag, =.T. until found final groung temp
89992793c5 Jean*0077       _RL aLoc, bLoc, cLoc, eLoc, fLoc, alpha, beta
                0078       _RL tSurf, tg(land_nLev), eg(land_nLev)
                0079       _RL cg(land_nLev), mW(land_nLev)
                0080       _RL temp_af, temp_bf
                0081       _RL mSnow, dMsn, delT
                0082       _RL mSnEpsil
                0083       _RL tmp1, tmp2
                0084       INTEGER i,j,k
4c9728b121 Jean*0085       LOGICAL tmpFlag
                0086 
                0087 #ifdef LAND_DEBUG
89992793c5 Jean*0088       CHARACTER*(MAX_LEN_MBUF) msgBuf
4c9728b121 Jean*0089       LOGICAL dBug, debugFlag
                0090       INTEGER iprt,jprt,lprt
                0091       DATA iprt, jprt , lprt / 19 , 20 , 6 /
                0092       DATA debugFlag / .FALSE. /
                0093  1010 FORMAT(A,I3,1P4E11.3)
                0094 #endif
89992793c5 Jean*0095 
                0096       DATA    mSnEpsil / 1. _d -6 /
                0097 
                0098 C-------------------------------------------------------------------------
                0099 C  solve implicitly the coupled 3 eq. (time level n+1 omitted) :
                0100 C    1a : if hs=0 : Flx = Flx^o + d.F/dT*(Ts - Ts^n) & Ts=Tg1
                0101 C    1b : if hs>0 : Flx = (Ts-Tg1)*Ks/hs =< Flx^o + d.F/dT*(Ts - Ts^n)
                0102 C         & difference used to melt the snow, maintaining Ts=0 
                0103 C    2  : Eg1 - Eg1^n  = delT*Flx - (lambda*delT/delZ_12)*(Tg1-Tg2)
                0104 C    3  : Eg2 - Eg2^n  = (lambda*delT/delZ_12)*(Tg1-Tg2)
                0105 C    were  lambda = ground Conductivity ; Ks = snow Conductivity
                0106 C          k=1,2: Eg_k = Cg_k * Tg_k - Lfreez * mIce_k
                0107 C
                0108 C  using local variables:
                0109 C   a = lambda*delT/delZ_12 ; b = - d.F/dT ;  f = Flx^o + b*Ts^n
                0110 C                         alpha = hs/Ks ;  beta = 1/(1+alpha*b)
                0111 C  3.eq system becomes: 
                0112 C   o if Ts*hs =< 0
                0113 C     1a,b:  Ts = ( Tg1 + alpha*F)*beta
                0114 C      2  : Eg1 + a*(Tg1-Tg2) + (b*delT)*beta*Tg1 = Eg1^n + delT*f*beta
                0115 C      3  : Eg2 + a*(Tg2-Tg1) = Eg2^n
                0116 C   o else: 
                0117 C     1.b : Ts=0 , f = Flx(ts=0) ; snowMelt = (f + Tg1/alpha)/Lfreez
                0118 C      2  : Eg1 + a*(Tg1-Tg2) + (delT/alpha)*Tg1 = Eg1^n
                0119 C      3  : Eg2 + a*(Tg2-Tg1) = Eg2^n
                0120 C-------------------------------------------------------------------------
                0121 
                0122 C---  Solve implicitely for ground temp. & surface temp
                0123 
                0124       delT = land_deltaT
                0125       aLoc = land_grdLambda*land_deltaT*land_rec_dzC(2)
                0126       DO j=1,sNy
                0127        DO i=1,sNx
                0128         IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
                0129 
                0130 C--   initialise local variables
                0131           tmpFlag = .TRUE.
                0132           tSurf = land_skinT(i,j,bi,bj)
                0133           mSnow = land_rhoSnow*land_hSnow(i,j,bi,bj)
                0134           bLoc  = -sFlx(i,j,2)
                0135           fLoc  = sFlx(i,j,1)+bLoc*tSurf
                0136           alpha = land_hSnow(i,j,bi,bj)/diffKsnow
                0137           beta  = 1. _d 0/(1. _d 0+alpha*bLoc)
                0138           DO k=1,land_nLev
                0139             eg(k) = land_dzF(k)*land_enthalp(i,j,k,bi,bj)
                0140             mW(k) = land_dzF(k)*land_groundW(i,j,k,bi,bj)
                0141      &                         *land_waterCap*land_rhoLiqW
4c9728b121 Jean*0142             mW(k) = MAX( mW(k), 0. _d 0 )
89992793c5 Jean*0143             cg(k) = land_dzF(k)*land_heatCs + mW(k)*land_CpWater
                0144             tg(k) = land_groundT(i,j,k,bi,bj)
                0145           ENDDO
4c9728b121 Jean*0146 #ifdef LAND_DEBUG
                0147           dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
                0148           IF (dBug) write(6,1010)
                0149      &         'LAND_IMPL_TEMP: 0 , ts,tg1&2,mSw=',0,tSurf,tg,mSnow
1cd785a34e Jean*0150           IF (dBug) write(6,1010)
                0151      &         'LAND_IMPL_TEMP: 0 , sFlx=', 0,(sFlx(i,j,k),k=0,2)
4c9728b121 Jean*0152 #endif
                0153 
89992793c5 Jean*0154 C---   Solve for temp as if no freezing/melting was occuring :
                0155           tg(1) = ( cg(1)*tg(1) + fLoc*delT*beta
                0156      &                          + cg(2)*tg(2)*aLoc/(cg(2)+aLoc)
                0157      &            )
                0158      &          / ( cg(1) + aLoc + bLoc*delT*beta
                0159      &                           - aLoc*aLoc/(cg(2)+aLoc) 
                0160      &            )
                0161           tg(2) = ( cg(2)*tg(2) + aLoc*tg(1) ) / (cg(2)+aLoc)
                0162           tSurf = ( tg(1) + alpha*fLoc ) * beta
                0163           
4c9728b121 Jean*0164 #ifdef LAND_DEBUG
                0165           IF (dBug) write(6,1010)
                0166      &         'LAND_IMPL_TEMP: 1 , ts,tg1&2,mW1=',1,tSurf,tg,mW(1)
                0167 #endif
89992793c5 Jean*0168 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0169 C---   If melting/freezing (top of snow layer, ground water level 1 or 2)
                0170 C      set corresponding temp to freezing point and update enthalpy 
                0171 C--------------
                0172 
                0173           IF ( tg(2)*land_groundT(i,j,2,bi,bj) .LE. 0. _d 0 
                0174      &         .AND. tmpFlag .AND. tSurf*mSnow .LE. 0. _d 0 ) THEN
                0175 C--    freezing/melting in level 2: set Tg2 to freezing point
                0176            tmp1 = tg(1)
                0177            tmp2 = tg(2)
                0178            tg(2) = 0.
                0179            eLoc = eg(1) + fLoc*delT*beta
                0180            cLoc = cg(1) + aLoc + bLoc*delT*beta
                0181            temp_bf = (eLoc+land_Lfreez*mW(1))/cLoc
                0182            temp_af =  eLoc/cLoc
                0183            tg(1) = MIN( temp_bf, MAX(temp_af, 0. _d 0) )
                0184            tSurf = ( tg(1) + alpha*fLoc ) * beta
                0185            IF ( tSurf*mSnow .LE. 0. _d 0 ) THEN
                0186              tmpFlag = .FALSE.
                0187              eg(1) = eLoc - (aLoc + bLoc*delT*beta)*tg(1)
                0188              eg(2) = eg(2) + aLoc*tg(1) 
4c9728b121 Jean*0189 #ifdef LAND_DEBUG
89992793c5 Jean*0190            ELSEIF ( debugFlag ) THEN
                0191              WRITE(msgBuf,'(A,2I4,2I3,I10)')
                0192      &             'LAND_IMPL_TEMP: i,j,bi,bj,Iter=',
                0193      &             i,j,bi,bj,myIter
                0194              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0195      &                    SQUEEZE_RIGHT, myThid )
                0196              WRITE(msgBuf,'(A,1P4E12.4)')
                0197      &             'LAND_IMPL_TEMP: groundT,t2,ts=',
                0198      &             land_groundT(i,j,1,bi,bj),land_groundT(i,j,2,bi,bj),
                0199      &             tmp2,(tmp1+alpha*fLoc)*beta
                0200              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0201      &                    SQUEEZE_RIGHT, myThid )
                0202              WRITE(msgBuf,'(A,1P4E12.4)')
                0203      &             'LAND_IMPL_TEMP: Tg,tSurf,mSnw=',
                0204      &              tg,tSurf,mSnow
                0205              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0206      &                    SQUEEZE_RIGHT, myThid )
                0207              WRITE(msgBuf,'(A,1P4E14.6)')
                0208      &             'LAND_IMPL_TEMP: eg,mW=', eg,mW
                0209              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0210      &                    SQUEEZE_RIGHT, myThid )
4c9728b121 Jean*0211 #endif
89992793c5 Jean*0212            ENDIF
                0213 
                0214 C-  if tg2_new*tg2_old < 0 : end
                0215           ENDIF
                0216 
                0217 C--------------
                0218 
                0219           IF ( tg(1)*land_groundT(i,j,1,bi,bj) .LE. 0. _d 0 
                0220      &         .AND. tmpFlag .AND. tSurf*mSnow .LE. 0. _d 0 ) THEN
                0221 C--    freezing/melting in level 1: set Tg1 to freezing point
                0222            tmp1 = tg(1)
                0223            tg(1) = 0.
                0224            tg(2) = cg(2)*tg(2) / (cg(2)+aLoc)
                0225            tSurf = alpha*fLoc * beta
                0226            IF ( tSurf*mSnow .LE. 0. _d 0 ) THEN
                0227              tmpFlag = .FALSE.
                0228              eg(2) = eg(2) - aLoc*tg(2) 
                0229              eg(1) = eg(1) + aLoc*tg(2) + fLoc*delT*beta
                0230              IF ( eg(1)*mSnow .GT. 0. _d 0 ) THEN
                0231 C-           melt snow from bottom
                0232               dMsn = MIN( mSnow, eg(1)*recip_Lfreez )
                0233               land_Pr_m_Ev(i,j,bi,bj) = dMsn/delT
                0234               land_hSnow(i,j,bi,bj) = (mSnow - dMsn)/land_rhoSnow 
                0235               eg(1) = eg(1) - dMsn*land_Lfreez
4c9728b121 Jean*0236 #ifdef LAND_DEBUG
                0237               IF (dBug) write(6,1010)
                0238      &         'LAND_IMPL_TEMP: Bot-Melt : dMsn,dEg1,eg1=',
                0239      &         1, dMsn, -dMsn*land_Lfreez, eg(1) 
                0240 #endif
89992793c5 Jean*0241              ENDIF
4c9728b121 Jean*0242 #ifdef LAND_DEBUG
89992793c5 Jean*0243            ELSEIF ( debugFlag ) THEN
                0244              WRITE(msgBuf,'(A,2I4,2I3,I10)')
                0245      &             'LAND_IMPL_TEMP: i,j,bi,bj,Iter=',
                0246      &             i,j,bi,bj,myIter
                0247              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0248      &                    SQUEEZE_RIGHT, myThid )
                0249              WRITE(msgBuf,'(A,4F11.6)')
                0250      &             'LAND_IMPL_TEMP: groundT,t1,ts=',
                0251      &             land_groundT(i,j,1,bi,bj),land_groundT(i,j,2,bi,bj),
                0252      &             tmp1,(tmp1+alpha*fLoc)*beta
                0253              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0254      &                    SQUEEZE_RIGHT, myThid )
                0255              WRITE(msgBuf,'(A,4F12.7)')
                0256      &             'LAND_IMPL_TEMP: Tg,tSurf,mSnow=',
                0257      &             tg,tSurf,mSnow
                0258              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0259      &                    SQUEEZE_RIGHT, myThid )
                0260              WRITE(msgBuf,'(A,1P4E14.6)')
                0261      &             'LAND_IMPL_TEMP: eg,mW=', eg,mW
                0262              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0263      &                    SQUEEZE_RIGHT, myThid )
                0264              WRITE(msgBuf,'(A)')  
                0265      &             'LAND_IMPL_TEMP: snow with ts > 0 ! but continue'
                0266              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0267      &                    SQUEEZE_RIGHT, myThid )
4c9728b121 Jean*0268 #endif
89992793c5 Jean*0269            ENDIF
                0270 
                0271 C-  if tg1_new*tg1_old < 0 : end
                0272           ENDIF
                0273 
                0274 C--------------
                0275 
                0276           IF ( tmpFlag .AND. tSurf*mSnow .GT. 0. _d 0 ) THEN
                0277 C--    snow is melting at the surface: set ts=0 & use fixed heat flux Flx(ts=0)
4c9728b121 Jean*0278 #ifdef LAND_DEBUG
                0279               IF (dBug) write(6,1010)
                0280      &         'LAND_IMPL_TEMP: Top-Melt : fx0, fx1, fx1-b*Ts =',
                0281      &         1, sFlx(i,j,0), fLoc, fLoc-bLoc*tSurf
                0282 #endif
89992793c5 Jean*0283            tSurf = 0. _d 0
                0284            fLoc = sFlx(i,j,0)
                0285            dTsurf(i,j) = 1000.
                0286            tg(1) = land_groundT(i,j,1,bi,bj)
                0287            tg(2) = land_groundT(i,j,2,bi,bj)
                0288 
                0289            eLoc = cg(1)*tg(1)
                0290      &          + delT*fLoc - land_Lfreez*mSnow + aLoc*tg(2)
                0291            IF ( eLoc .GT. 0. _d 0 .OR. mSnow.LT.mSnEpsil ) THEN
                0292 C-     all snow melt: do not solve diffusion of heat in snow layer
                0293 C      but put surf. heat flux directly to 1rst level and set tg1=0
                0294              dMsn = mSnow
                0295              tg(1) = 0. _d 0
                0296              tg(2) = cg(2)*tg(2) / (cg(2)+aLoc)
                0297            ELSE
                0298 C-     solve diffusion of heat in snow layer ; heat flux difference 
                0299 C      (surf.Flx - diffusion.Flx) is used to melt the snow from top.
                0300             tg(1) = ( cg(1)*tg(1) + cg(2)*tg(2)*aLoc/(cg(2)+aLoc) )
                0301      &       / ( cg(1)+aLoc + delT/alpha - aLoc*aLoc/(cg(2)+aLoc) )
                0302             tg(2) = ( cg(2)*tg(2) + aLoc*tg(1) ) / (cg(2)+aLoc)
                0303             IF ( tg(2)*land_groundT(i,j,2,bi,bj).LE.0. _d 0 ) THEN
                0304               tg(2) = 0.
                0305               tg(1) = cg(1)*tg(1) / ( cg(1)+aLoc + delT/alpha )
                0306             ELSEIF ( tg(1)*land_groundT(i,j,1,bi,bj).LE.0. _d 0 ) THEN
                0307               tg(1) = 0.
                0308               tg(2) = cg(2)*tg(2) / (cg(2)+aLoc)
                0309             ENDIF
                0310             dMsn = ( fLoc+tg(1)/alpha )*delT*recip_Lfreez
4c9728b121 Jean*0311 #ifdef LAND_DEBUG
                0312               IF (dBug) write(6,1010)
                0313      &         'LAND_IMPL_TEMP: Surf-Melt: dMsn,fLoc,tg1/alpha=',
                0314      &          2, dMsn, fLoc,tg(1)/alpha
                0315 #endif
                0316 C-  note: Fx0 < -tg(1)/alpha can happen (due to non-linearity in Fx(Ts)), 
                0317 C-     => do not melt nor accumulate snow but put d.Flx in Eg1
                0318             dMsn = MIN( MAX(dMsn, 0. _d 0), mSnow )
89992793c5 Jean*0319            ENDIF
                0320            tmpFlag = .FALSE.
                0321            eg(2) = eg(2) + aLoc*(tg(1)-tg(2))
                0322            eg(1) = eg(1) - aLoc*(tg(1)-tg(2))
                0323      &           + delT*fLoc - land_Lfreez*dMsn
                0324            land_Pr_m_Ev(i,j,bi,bj) = dMsn/delT
                0325            land_hSnow(i,j,bi,bj) = (mSnow - dMsn)/land_rhoSnow 
                0326 
                0327 C-  if ts*hSnow > 0 , else:
                0328           ELSEIF ( tmpFlag ) THEN
                0329 C--   snow is not melting & no freezing/melting in ground level 1 & 2
                0330            eg(2) = eg(2) + aLoc*(tg(1)-tg(2))
                0331            eg(1) = eg(1) - aLoc*(tg(1)-tg(2))
                0332      &           + delT*(fLoc-bLoc*Tsurf)
                0333            tmpFlag = .FALSE.
                0334           ENDIF
                0335 
                0336 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0337 
                0338 C---  Save new values :
                0339           IF ( dTsurf(i,j) .LE. 999. ) 
                0340      &         dTsurf(i,j) = tSurf - land_skinT(i,j,bi,bj)
                0341           land_skinT(i,j,bi,bj) = tSurf
                0342           DO k=1,land_nLev
                0343             land_enthalp(i,j,k,bi,bj) = eg(k)/land_dzF(k)
                0344             land_groundT(i,j,k,bi,bj) = tg(k)
                0345           ENDDO
                0346 
4c9728b121 Jean*0347 #ifdef LAND_DEBUG
                0348           IF (dBug) write(6,1010) 'LAND_IMPL_TEMP: 9, ts,tg1&2,dTs=',9,
                0349      &         tSurf, tg, dTsurf(i,j)
                0350           IF (dBug) write(6,1010) 'LAND_IMPL_TEMP: 9, Eg1,Eg2,mPmE=',9,
                0351      &       (land_enthalp(i,j,k,bi,bj),k=1,2), land_Pr_m_Ev(i,j,bi,bj)
                0352 #endif
                0353 
89992793c5 Jean*0354         ENDIF
                0355        ENDDO
                0356       ENDDO
                0357 
                0358 #endif /* ALLOW_LAND */
                0359 
                0360       RETURN
                0361       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!