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
0004
0005
0006 SUBROUTINE LAND_STEPFWD(
0007 I land_frc, bi, bj, myTime, myIter, myThid)
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 IMPLICIT NONE
0018
0019
0020
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
0028
0029
0030
0031
0032
0033
0034 _RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0035 INTEGER bi, bj, myIter, myThid
0036 _RL myTime
0037
0038
0039 #ifdef ALLOW_LAND
0040
0041
0042
89992793c5 Jean*0043
b30f09edfc Jean*0044
e0a2f8aec4 Jean*0045
89992793c5 Jean*0046
b30f09edfc Jean*0047
e0a2f8aec4 Jean*0048
b30f09edfc Jean*0049
89992793c5 Jean*0050
0051
b30f09edfc Jean*0052
0053
0054
0055
89992793c5 Jean*0056
0057
0058
0059
0060
0061
0062
0063
0064
b30f09edfc Jean*0065
0066
0067
89992793c5 Jean*0068
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
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
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
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
0126 ENDIF
0127
89992793c5 Jean*0128
0129
1665818901 Jean*0130 IF ( land_calc_grW .OR. land_calc_snow ) THEN
b30f09edfc Jean*0131
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
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
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
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
20cd0df114 Jean*0191
89992793c5 Jean*0192 snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )
ecb76f8691 Jean*0193
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
0200 land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj)
0201 & *EXP( -hNewSnow/hNewSnowAge )
b30f09edfc Jean*0202
0203
0204
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
89992793c5 Jean*0219
0220
0221
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
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
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
89992793c5 Jean*0246
0247
0248
0249
0250
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
0264
e0a2f8aec4 Jean*0265 IF (land_calc_grW) THEN
0266
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
0291 IF ( flxkup(i,j) .LT. 0. _d 0 ) THEN
0292
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
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
0306 flxkup(i,j) = 0. _d 0
0307 flxEngU(i,j) = 0. _d 0
0308
0309 ELSE
0310
0311
89992793c5 Jean*0312
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
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
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
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
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
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
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
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
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
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
0387 ENDIF
0388
89992793c5 Jean*0389
0390
b30f09edfc Jean*0391 IF ( land_calc_grT ) THEN
0392
89992793c5 Jean*0393
0394 DO k=1,land_nLev
0395 DO j=1,sNy
0396 DO i=1,sNx
0397
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
0403 temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
0404 & / grd_HeatCp
0405
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
0437 ENDIF
0438
e0a2f8aec4 Jean*0439 #endif /* ALLOW_LAND */
0440
0441 RETURN
0442 END