44 |
C i,j,k :: loop counters |
C i,j,k :: loop counters |
45 |
C kp1 :: k+1 |
C kp1 :: k+1 |
46 |
C grd_HeatCp :: Heat capacity of the ground [J/m3/K] |
C grd_HeatCp :: Heat capacity of the ground [J/m3/K] |
47 |
|
C enthalpGrdW :: enthalpy of ground water [J/m3] |
48 |
C fieldCapac :: field capacity (of water) [m] |
C fieldCapac :: field capacity (of water) [m] |
49 |
C mWater :: water content of the ground [kg/m3] |
C mWater :: water content of the ground [kg/m3] |
50 |
C fractRunOff :: fraction of water in excess which leaves as runoff |
C groundWnp1 :: hold temporary future soil moisture [] |
51 |
C grdWexcess :: ground water in excess [m/s] |
C grdWexcess :: ground water in excess [m/s] |
52 |
C groundWnp1 :: hold temporary future soil moisture |
C fractRunOff :: fraction of water in excess which leaves as runoff |
|
C enthalpGrdW :: enthalpy of ground water [J/m3] |
|
53 |
C flxkup :: downward flux of water, upper interface (k-1,k) |
C flxkup :: downward flux of water, upper interface (k-1,k) |
54 |
C flxdwn :: downward flux of water, lower interface (k,k+1) |
C flxdwn :: downward flux of water, lower interface (k,k+1) |
55 |
C flxEng :: downward energy flux associated with water flux |
C flxEngU :: downward energy flux associated with water flux (W/m2) |
56 |
|
C upper interface (k-1,k) |
57 |
|
C flxEngL :: downward energy flux associated with water flux (W/m2) |
58 |
|
C lower interface (k,k+1) |
59 |
C temp_af :: ground temperature if above freezing |
C temp_af :: ground temperature if above freezing |
60 |
C temp_bf :: ground temperature if below freezing |
C temp_bf :: ground temperature if below freezing |
61 |
C mPmE :: hold temporary (liquid) Precip minus Evap [kg/m2/s] |
C mPmE :: hold temporary (liquid) Precip minus Evap [kg/m2/s] |
65 |
C dMsn :: mass of melting snow [kg/m2] |
C dMsn :: mass of melting snow [kg/m2] |
66 |
C snowPrec :: snow precipitation [kg/m2/s] |
C snowPrec :: snow precipitation [kg/m2/s] |
67 |
C hNewSnow :: fresh snow accumulation [m] |
C hNewSnow :: fresh snow accumulation [m] |
68 |
|
C dhSnowMx :: potential snow increase [m] |
69 |
|
C dhSnow :: effective snow increase [m] |
70 |
|
C mIceDt :: ground-ice growth rate (<- excess of snow) [kg/m2/s] |
71 |
C ageFac :: snow aging factor [1] |
C ageFac :: snow aging factor [1] |
72 |
_RL grd_HeatCp, fieldCapac, mWater |
_RL grd_HeatCp, enthalpGrdW |
73 |
_RL fractRunOff, grdWexcess, groundWnp1, enthalpGrdW |
_RL fieldCapac, mWater |
74 |
|
_RL groundWnp1, grdWexcess, fractRunOff |
75 |
_RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
76 |
_RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
77 |
_RL flxEng(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL flxEngU(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
78 |
_RL temp_af, temp_bf, mPmE, enWfx, enGr1 |
_RL flxEngL, temp_af, temp_bf, mPmE, enWfx, enGr1 |
79 |
_RL mSnow, dMsn, snowPrec, hNewSnow, ageFac |
_RL mSnow, dMsn, snowPrec |
80 |
|
_RL hNewSnow, dhSnowMx, dhSnow, mIceDt, ageFac |
81 |
INTEGER i,j,k,kp1 |
INTEGER i,j,k,kp1 |
82 |
|
|
83 |
IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN |
IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN |
126 |
#ifdef LAND_OLD_VERSION |
#ifdef LAND_OLD_VERSION |
127 |
IF ( .TRUE. ) THEN |
IF ( .TRUE. ) THEN |
128 |
#else |
#else |
129 |
IF ( land_calc_snow ) THEN |
IF ( land_calc_grW ) THEN |
130 |
#endif |
#endif |
131 |
|
C-- Initialize run-off arrays. |
132 |
|
DO j=1,sNy |
133 |
|
DO i=1,sNx |
134 |
|
land_runOff(i,j,bi,bj) = 0. _d 0 |
135 |
|
land_enRnOf(i,j,bi,bj) = 0. _d 0 |
136 |
|
ENDDO |
137 |
|
ENDDO |
138 |
C-- need (later on) ground temp. to be consistent with updated enthalpy: |
C-- need (later on) ground temp. to be consistent with updated enthalpy: |
139 |
DO k=1,land_nLev |
DO k=1,land_nLev |
140 |
DO j=1,sNy |
DO j=1,sNy |
172 |
snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 ) |
snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 ) |
173 |
snowPrec = MAX( snowPrec*recip_Lfreez , 0. _d 0 ) |
snowPrec = MAX( snowPrec*recip_Lfreez , 0. _d 0 ) |
174 |
mPmE = mPmE - snowPrec |
mPmE = mPmE - snowPrec |
175 |
flxEng(i,j) = enWfx + land_Lfreez*snowPrec |
flxEngU(i,j) = enWfx + land_Lfreez*snowPrec |
176 |
hNewSnow = land_deltaT * snowPrec / land_rhoSnow |
hNewSnow = land_deltaT * snowPrec / land_rhoSnow |
|
land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow |
|
177 |
C- refresh snow age: |
C- refresh snow age: |
178 |
land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj) |
land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj) |
179 |
& *EXP( -hNewSnow/hNewSnowAge ) |
& *EXP( -hNewSnow/hNewSnowAge ) |
180 |
|
C- update snow thickness: |
181 |
|
c land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow |
182 |
|
C glacier & ice-sheet missing: excess of snow put directly into run-off |
183 |
|
dhSnowMx = MAX( 0. _d 0, |
184 |
|
& land_hMaxSnow - land_hSnow(i,j,bi,bj) ) |
185 |
|
dhSnow = MIN( hNewSnow, dhSnowMx ) |
186 |
|
land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow |
187 |
|
mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT |
188 |
|
land_runOff(i,j,bi,bj) = mIceDt/land_rhoLiqW |
189 |
|
land_enRnOf(i,j,bi,bj) = -mIceDt*land_Lfreez |
190 |
ELSE |
ELSE |
191 |
C- rain precip (whatever Evap is) or Evap exceeds snow precip : |
C- rain precip (whatever Evap is) or Evap exceeds snow precip : |
192 |
C => snow melts or sublimates |
C => snow melts or sublimates |
197 |
IF ( dMsn .GE. mSnow ) THEN |
IF ( dMsn .GE. mSnow ) THEN |
198 |
dMsn = mSnow |
dMsn = mSnow |
199 |
land_hSnow(i,j,bi,bj) = 0. _d 0 |
land_hSnow(i,j,bi,bj) = 0. _d 0 |
200 |
flxEng(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT |
flxEngU(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT |
201 |
ELSE |
ELSE |
202 |
flxEng(i,j) = 0. _d 0 |
flxEngU(i,j) = 0. _d 0 |
203 |
land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) |
land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) |
204 |
& - dMsn / land_rhoSnow |
& - dMsn / land_rhoSnow |
205 |
ENDIF |
ENDIF |
223 |
DO j=1,sNy |
DO j=1,sNy |
224 |
DO i=1,sNx |
DO i=1,sNx |
225 |
flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW |
flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW |
226 |
flxEng(i,j) = 0. _d 0 |
flxEngU(i,j) = 0. _d 0 |
227 |
ENDDO |
ENDDO |
228 |
ENDDO |
ENDDO |
229 |
ENDIF |
ENDIF |
243 |
ENDIF |
ENDIF |
244 |
fieldCapac = land_waterCap*land_dzF(k) |
fieldCapac = land_waterCap*land_dzF(k) |
245 |
|
|
|
IF (k.EQ.1) THEN |
|
|
DO j=1,sNy |
|
|
DO i=1,sNx |
|
|
land_runOff(i,j,bi,bj) = 0. _d 0 |
|
|
land_enRnOf(i,j,bi,bj) = 0. _d 0 |
|
|
ENDDO |
|
|
ENDDO |
|
|
ELSE |
|
|
DO j=1,sNy |
|
|
DO i=1,sNx |
|
|
flxkup(i,j) = flxkdw(i,j) |
|
|
ENDDO |
|
|
ENDDO |
|
|
ENDIF |
|
|
|
|
246 |
DO j=1,sNy |
DO j=1,sNy |
247 |
DO i=1,sNx |
DO i=1,sNx |
248 |
IF ( land_frc(i,j,bi,bj).GT.0. ) THEN |
IF ( land_frc(i,j,bi,bj).GT.0. ) THEN |
249 |
|
|
250 |
|
#ifdef LAND_OLD_VERSION |
251 |
|
IF ( .TRUE. ) THEN |
252 |
|
IF ( k.EQ.land_nLev ) THEN |
253 |
|
#else |
254 |
|
IF ( land_groundT(i,j,k,bi,bj).LT.0. _d 0 ) THEN |
255 |
|
C- Frozen level: only account for upper level fluxes |
256 |
|
IF ( flxkup(i,j) .LT. 0. _d 0 ) THEN |
257 |
|
C- Step forward soil moisture (& enthapy), level k : |
258 |
|
land_groundW(i,j,k,bi,bj) = land_groundW(i,j,k,bi,bj) |
259 |
|
& + land_deltaT * flxkup(i,j) / fieldCapac |
260 |
|
IF ( land_calc_snow ) |
261 |
|
& land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj) |
262 |
|
& + land_deltaT * flxEngU(i,j) / land_dzF(k) |
263 |
|
ELSE |
264 |
|
C- Frozen level: incoming water flux goes directly into run-off |
265 |
|
land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj) |
266 |
|
& + flxkup(i,j) |
267 |
|
land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj) |
268 |
|
& + flxEngU(i,j) |
269 |
|
ENDIF |
270 |
|
C- prepare fluxes for next level: |
271 |
|
flxkup(i,j) = 0. _d 0 |
272 |
|
flxEngU(i,j) = 0. _d 0 |
273 |
|
|
274 |
|
ELSE |
275 |
|
|
276 |
|
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
277 |
C- Diffusion flux of water, lower interface (k,k+1): |
C- Diffusion flux of water, lower interface (k,k+1): |
278 |
flxkdw(i,j) = fieldCapac* |
IF ( k.EQ.land_nLev .OR. |
279 |
& ( land_groundW(i,j,k,bi,bj) |
& land_groundT(i,j,kp1,bi,bj).LT.0. _d 0 ) THEN |
280 |
& -land_groundW(i,j,kp1,bi,bj) ) |
#endif /* LAND_OLD_VERSION */ |
281 |
& / land_wTauDiff |
C- no Diffusion of water if one level is frozen : |
282 |
|
flxkdw(i,j) = 0. _d 0 |
283 |
|
flxEngL = 0. _d 0 |
284 |
|
ELSE |
285 |
|
flxkdw(i,j) = fieldCapac* |
286 |
|
& ( land_groundW(i,j,k,bi,bj) |
287 |
|
& -land_groundW(i,j,kp1,bi,bj) ) |
288 |
|
& / land_wTauDiff |
289 |
|
C- energy flux associated with water flux: take upwind Temp |
290 |
|
IF ( flxkdw(i,j).GE.0. ) THEN |
291 |
|
flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater |
292 |
|
& *land_groundT(i,j,k,bi,bj) |
293 |
|
ELSE |
294 |
|
flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater |
295 |
|
& *land_groundT(i,j,kp1,bi,bj) |
296 |
|
ENDIF |
297 |
|
ENDIF |
298 |
|
|
299 |
C- Step forward soil moisture, level k : |
C- Step forward soil moisture, level k : |
300 |
groundWnp1 = land_groundW(i,j,k,bi,bj) |
groundWnp1 = land_groundW(i,j,k,bi,bj) |
301 |
& + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac |
& + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac |
302 |
land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1) |
|
303 |
|
C- Water in excess will leave as run-off or go to level below |
304 |
|
land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1) |
305 |
|
grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) ) |
306 |
|
& *fieldCapac/land_deltaT |
307 |
|
|
308 |
C- Run off: fraction 1-fractRunOff enters level below |
C- Run off: fraction 1-fractRunOff enters level below |
309 |
grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) ) |
land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj) |
310 |
& *fieldCapac/land_deltaT |
& + fractRunOff*grdWexcess |
311 |
land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj) |
C- prepare fluxes for next level: |
312 |
& + fractRunOff*grdWexcess |
flxkup(i,j) = flxkdw(i,j) |
313 |
|
& + (1. _d 0-fractRunOff)*grdWexcess |
314 |
IF ( land_calc_snow ) THEN |
|
315 |
C- account for water fluxes in energy budget: |
IF ( land_calc_snow ) THEN |
316 |
enthalpGrdW = land_enthalp(i,j,k,bi,bj) |
enthalpGrdW = land_rhoLiqW*land_CpWater |
317 |
& - land_heatCs*land_groundT(i,j,k,bi,bj) |
& *land_groundT(i,j,k,bi,bj) |
318 |
land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj) |
C-- Account for water fluxes in energy budget: update ground Enthalpy |
319 |
& + fractRunOff*grdWexcess*enthalpGrdW |
land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj) |
320 |
land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj) |
& + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW |
321 |
& + ( flxEng(i,j) - (flxkdw(i,j)+grdWexcess)*enthalpGrdW |
& )*land_deltaT/land_dzF(k) |
322 |
& )*land_deltaT/land_dzF(k) |
|
323 |
ELSE |
land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj) |
324 |
enthalpGrdW = 0. _d 0 |
& + fractRunOff*grdWexcess*enthalpGrdW |
|
ENDIF |
|
325 |
C- prepare fluxes for next level: |
C- prepare fluxes for next level: |
326 |
flxkdw(i,j) = flxkdw(i,j) |
flxEngU(i,j) = flxEngL |
327 |
& + (1. _d 0-fractRunOff)*grdWexcess |
& + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW |
328 |
flxEng(i,j) = flxkdw(i,j)*enthalpGrdW |
ENDIF |
329 |
|
ENDIF |
330 |
|
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
331 |
|
|
332 |
ENDIF |
ENDIF |
333 |
ENDDO |
ENDDO |
339 |
|
|
340 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
341 |
|
|
342 |
IF (land_calc_grT ) THEN |
IF ( land_calc_grT ) THEN |
343 |
C-- Compute ground temperature from enthalpy : |
C-- Compute ground temperature from enthalpy (if not already done): |
344 |
|
|
345 |
DO k=1,land_nLev |
DO k=1,land_nLev |
346 |
DO j=1,sNy |
DO j=1,sNy |