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
0004
0005
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
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 IMPLICIT NONE
0025
0026
0027
0028 #include "LAND_SIZE.h"
0029
0030 #include "EEPARAMS.h"
0031 #include "LAND_PARAMS.h"
0032 #include "LAND_VARS.h"
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
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
0051
0052 #ifdef ALLOW_LAND
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
4c9728b121 Jean*0076
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
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
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
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
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
0169
0170
0171
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
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
0215 ENDIF
0216
0217
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
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
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
0272 ENDIF
0273
0274
0275
0276 IF ( tmpFlag .AND. tSurf*mSnow .GT. 0. _d 0 ) THEN
0277
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
0293
0294 dMsn = mSnow
0295 tg(1) = 0. _d 0
0296 tg(2) = cg(2)*tg(2) / (cg(2)+aLoc)
0297 ELSE
0298
0299
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
0317
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
0328 ELSEIF ( tmpFlag ) THEN
0329
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
0337
0338
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