C $Header: /u/gcmpack/MITgcm/pkg/land/land_impl_temp.F,v 1.3 2004/06/21 16:50:21 jmc Exp $
C $Name: $
#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