/[MITgcm]/MITgcm/pkg/land/land_stepfwd.F
ViewVC logotype

Contents of /MITgcm/pkg/land/land_stepfwd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.6 - (show annotations) (download)
Thu Jun 3 16:43:14 2004 UTC (19 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint57t_post, checkpoint57o_post, checkpoint58e_post, checkpoint57v_post, checkpoint53d_post, checkpoint58u_post, checkpoint58w_post, checkpoint54a_pre, checkpoint57m_post, checkpoint55c_post, checkpoint54e_post, checkpoint57s_post, checkpoint54a_post, checkpoint57k_post, checkpoint55d_pre, checkpoint57d_post, checkpoint57g_post, checkpoint57b_post, checkpoint57c_pre, checkpoint58r_post, checkpoint55j_post, checkpoint56b_post, checkpoint57i_post, checkpoint57y_post, checkpoint57e_post, checkpoint55h_post, checkpoint58n_post, checkpoint58x_post, checkpoint57g_pre, checkpoint54b_post, checkpoint55b_post, checkpoint58t_post, checkpoint58h_post, checkpoint54d_post, checkpoint56c_post, checkpoint57y_pre, checkpoint55, checkpoint57f_pre, checkpoint57a_post, checkpoint54, checkpoint58q_post, checkpoint54f_post, checkpoint55g_post, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint55f_post, checkpoint59c, checkpoint59b, checkpoint59h, checkpoint57r_post, checkpoint59, checkpoint58, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, eckpoint57e_pre, checkpoint57h_done, checkpoint58f_post, checkpoint53g_post, checkpoint57x_post, checkpoint57n_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint57f_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint57c_post, checkpoint58y_post, checkpoint55e_post, checkpoint58k_post, checkpoint58v_post, checkpoint53f_post, checkpoint55a_post, checkpoint54c_post, checkpoint58s_post, checkpoint58p_post, checkpoint57j_post, checkpoint58b_post, checkpoint57h_pre, checkpoint58m_post, checkpoint57l_post, checkpoint57h_post, checkpoint56a_post, checkpoint55d_post
Changes since 1.5: +49 -1 lines
add option "LAND_DEBUG" ; fix Pb in implicit Temp. solver when
  Ts arround freezing (was causing grdW < 0)

1 C $Header: /u/gcmpack/MITgcm/pkg/land/land_stepfwd.F,v 1.5 2004/05/21 13:41:02 jmc Exp $
2 C $Name: $
3
4 #include "LAND_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: LAND_STEPFWD
8 C !INTERFACE:
9 SUBROUTINE LAND_STEPFWD(
10 I land_frc, bi, bj, myTime, myIter, myThid)
11
12 C !DESCRIPTION: \bv
13 C *==========================================================*
14 C | S/R LAND_STEPFWD
15 C | o Land model main S/R: step forward land variables
16 C *==========================================================*
17 C \ev
18
19 C !USES:
20 IMPLICIT NONE
21
22 C == Global variables ===
23 C-- size for MITgcm & Land package :
24 #include "LAND_SIZE.h"
25
26 #include "EEPARAMS.h"
27 #include "LAND_PARAMS.h"
28 #include "LAND_VARS.h"
29
30 C !INPUT/OUTPUT PARAMETERS:
31 C == Routine arguments ==
32 C land_frc :: land fraction [0-1]
33 C bi,bj :: Tile index
34 C myTime :: Current time of simulation ( s )
35 C myIter :: Current iteration number in simulation
36 C myThid :: Number of this instance of the routine
37 _RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
38 INTEGER bi, bj, myIter, myThid
39 _RL myTime
40 CEOP
41
42 #ifdef ALLOW_LAND
43 C == Local variables ==
44 C i,j,k :: loop counters
45 C kp1 :: k+1
46 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]
49 C mWater :: water content of the ground [kg/m3]
50 C groundWnp1 :: hold temporary future soil moisture []
51 C grdWexcess :: ground water in excess [m/s]
52 C fractRunOff :: fraction of water in excess which leaves as runoff
53 C flxkup :: downward flux of water, upper interface (k-1,k)
54 C flxdwn :: downward flux of water, lower interface (k,k+1)
55 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
60 C temp_bf :: ground temperature if below freezing
61 C mPmE :: hold temporary (liquid) Precip minus Evap [kg/m2/s]
62 C enWfx :: hold temporary energy flux of Precip [W/m2]
63 C enGr1 :: ground enthalpy of level 1 [J/m2]
64 C mSnow :: mass of snow [kg/m2]
65 C dMsn :: mass of melting snow [kg/m2]
66 C snowPrec :: snow precipitation [kg/m2/s]
67 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]
72 _RL grd_HeatCp, enthalpGrdW
73 _RL fieldCapac, mWater
74 _RL groundWnp1, grdWexcess, fractRunOff
75 _RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 _RL flxEngU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78 _RL flxEngL, temp_af, temp_bf, mPmE, enWfx, enGr1
79 _RL mSnow, dMsn, snowPrec
80 _RL hNewSnow, dhSnowMx, dhSnow, mIceDt, ageFac
81 INTEGER i,j,k,kp1
82
83 #ifdef LAND_DEBUG
84 LOGICAL dBug
85 INTEGER iprt,jprt,lprt
86 DATA iprt, jprt , lprt / 19 , 20 , 6 /
87 1010 FORMAT(A,I3,1P4E11.3)
88 #endif
89
90 IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN
91 C-- Step forward ground temperature:
92
93 DO k=1,land_nLev
94 kp1 = MIN(k+1,land_nLev)
95
96 IF (k.EQ.1) THEN
97 DO j=1,sNy
98 DO i=1,sNx
99 flxkup(i,j) = land_HeatFlx(i,j,bi,bj)
100 ENDDO
101 ENDDO
102 ELSE
103 DO j=1,sNy
104 DO i=1,sNx
105 flxkup(i,j) = flxkdw(i,j)
106 ENDDO
107 ENDDO
108 ENDIF
109
110 DO j=1,sNy
111 DO i=1,sNx
112 IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
113 C- Thermal conductivity flux, lower interface (k,k+1):
114 flxkdw(i,j) = land_grdLambda*
115 & ( land_groundT(i,j,k,bi,bj)
116 & -land_groundT(i,j,kp1,bi,bj) )
117 & *land_rec_dzC(kp1)
118
119 C- Step forward ground enthalpy, level k :
120 land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
121 & + land_deltaT * (flxkup(i,j)-flxkdw(i,j))/land_dzF(k)
122
123 ENDIF
124 ENDDO
125 ENDDO
126
127 ENDDO
128 C-- step forward ground temperature: end
129 ENDIF
130
131 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
132
133 IF ( land_calc_grW .OR. land_calc_snow ) THEN
134 C-- Initialize run-off arrays.
135 DO j=1,sNy
136 DO i=1,sNx
137 land_runOff(i,j,bi,bj) = 0. _d 0
138 land_enRnOf(i,j,bi,bj) = 0. _d 0
139 ENDDO
140 ENDDO
141 ENDIF
142
143 #ifdef LAND_OLD_VERSION
144 IF ( .TRUE. ) THEN
145 #else
146 IF ( land_calc_grW ) THEN
147 #endif
148 C-- need (later on) ground temp. to be consistent with updated enthalpy:
149 DO k=1,land_nLev
150 DO j=1,sNy
151 DO i=1,sNx
152 IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
153 mWater = land_rhoLiqW*land_waterCap
154 & *land_groundW(i,j,k,bi,bj)
155 mWater = MAX( mWater, 0. _d 0 )
156 grd_HeatCp = land_heatCs + land_CpWater*mWater
157 temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
158 & / grd_HeatCp
159 temp_af = land_enthalp(i,j,k,bi,bj) / grd_HeatCp
160 land_groundT(i,j,k,bi,bj) =
161 & MIN( temp_bf, MAX(temp_af, 0. _d 0) )
162 #ifdef LAND_DEBUG
163 dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
164 IF (dBug) write(6,1010)
165 & 'LAND_STEPFWD: k,temp,af,bf=',
166 & k,land_groundT(i,j,k,bi,bj),temp_af,temp_bf
167 #endif
168 ENDIF
169 ENDDO
170 ENDDO
171 ENDDO
172 ENDIF
173
174 IF ( land_calc_snow ) THEN
175 C-- Step forward Snow thickness (also account for rain temperature)
176 ageFac = 1. _d 0 - land_deltaT/timeSnowAge
177 DO j=1,sNy
178 DO i=1,sNx
179 IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
180 mPmE = land_Pr_m_Ev(i,j,bi,bj)
181 enWfx = land_EnWFlux(i,j,bi,bj)
182 enGr1 = land_enthalp(i,j,1,bi,bj)*land_dzF(1)
183 #ifdef LAND_DEBUG
184 dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
185 IF (dBug) write(6,1010)
186 & 'LAND_STEPFWD:mPmE,enWfx,enGr1/dt,hSnow=',0,
187 & mPmE,enWfx,enGr1/land_deltaT,land_hSnow(i,j,bi,bj)
188 #endif
189 C- snow aging:
190 land_snowAge(i,j,bi,bj) =
191 & ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac )
192 IF ( enWfx.LT.0. ) THEN
193 C- snow precip in excess ( > Evap of snow) or snow prec & Evap of Liq.Water:
194 C => start to melt (until ground at freezing point) and then accumulate
195 snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )
196 C- snow accumulation cannot be larger that net precip
197 snowPrec = MAX( 0. _d 0 ,
198 & MIN( snowPrec*recip_Lfreez, mPmE ) )
199 mPmE = mPmE - snowPrec
200 flxEngU(i,j) = enWfx + land_Lfreez*snowPrec
201 hNewSnow = land_deltaT * snowPrec / land_rhoSnow
202 C- refresh snow age:
203 land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj)
204 & *EXP( -hNewSnow/hNewSnowAge )
205 C- update snow thickness:
206 c land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow
207 C glacier & ice-sheet missing: excess of snow put directly into run-off
208 dhSnowMx = MAX( 0. _d 0,
209 & land_hMaxSnow - land_hSnow(i,j,bi,bj) )
210 dhSnow = MIN( hNewSnow, dhSnowMx )
211 land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow
212 mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT
213 land_runOff(i,j,bi,bj) = mIceDt/land_rhoLiqW
214 land_enRnOf(i,j,bi,bj) = -mIceDt*land_Lfreez
215 #ifdef LAND_DEBUG
216 IF (dBug) write(6,1010)
217 & 'LAND_STEPFWD: 3,snP,mPmE,hNsnw,hSnw=',
218 & 3,snowPrec,mPmE,hNewSnow,land_hSnow(i,j,bi,bj)
219 #endif
220 ELSE
221 C- rain precip (whatever Evap is) or Evap of snow exceeds snow precip:
222 C => snow melts or sublimates
223 c snowMelt = MIN( enWfx*recip_Lfreez ,
224 c & land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT )
225 mSnow = land_hSnow(i,j,bi,bj)*land_rhoSnow
226 dMsn = enWfx*recip_Lfreez*land_deltaT
227 IF ( dMsn .GE. mSnow ) THEN
228 dMsn = mSnow
229 land_hSnow(i,j,bi,bj) = 0. _d 0
230 flxEngU(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT
231 ELSE
232 flxEngU(i,j) = 0. _d 0
233 land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj)
234 & - dMsn / land_rhoSnow
235 ENDIF
236 c IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge
237 mPmE = mPmE + dMsn/land_deltaT
238 #ifdef LAND_DEBUG
239 IF (dBug) write(6,1010)
240 & 'LAND_STEPFWD: 4,dMsn,mPmE,hSnw,enWfx=',
241 & 4,dMsn,mPmE,land_hSnow(i,j,bi,bj),flxEngU(i,j)
242 #endif
243 ENDIF
244 flxkup(i,j) = mPmE/land_rhoLiqW
245 c land_Pr_m_Ev(i,j,bi,bj) = mPmE
246 IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 )
247 & land_snowAge(i,j,bi,bj) = 0. _d 0
248 C- avoid negative (but very small, < 1.e-34) hSnow that occurs because
249 C of truncation error. Might need to rewrite this part.
250 c IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN
251 c land_hSnow(i,j,bi,bj) = 0. _d 0
252 c land_snowAge(i,j,bi,bj) = 0. _d 0
253 c ENDIF
254 ENDIF
255 ENDDO
256 ENDDO
257 ELSE
258 DO j=1,sNy
259 DO i=1,sNx
260 flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW
261 flxEngU(i,j) = 0. _d 0
262 ENDDO
263 ENDDO
264 ENDIF
265
266 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
267
268 IF (land_calc_grW) THEN
269 C-- Step forward ground Water:
270
271 DO k=1,land_nLev
272 IF (k.EQ.land_nLev) THEN
273 kp1 = k
274 fractRunOff = 1. _d 0
275 ELSE
276 kp1 = k+1
277 fractRunOff = land_fractRunOff
278 ENDIF
279 fieldCapac = land_waterCap*land_dzF(k)
280
281 DO j=1,sNy
282 DO i=1,sNx
283 IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
284 #ifdef LAND_DEBUG
285 dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
286 #endif
287
288 #ifdef LAND_OLD_VERSION
289 IF ( .TRUE. ) THEN
290 IF ( k.EQ.land_nLev ) THEN
291 #else
292 IF ( land_groundT(i,j,k,bi,bj).LT.0. _d 0 ) THEN
293 C- Frozen level: only account for upper level fluxes
294 IF ( flxkup(i,j) .LT. 0. _d 0 ) THEN
295 C- Step forward soil moisture (& enthapy), level k :
296 land_groundW(i,j,k,bi,bj) = land_groundW(i,j,k,bi,bj)
297 & + land_deltaT * flxkup(i,j) / fieldCapac
298 IF ( land_calc_snow )
299 & land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
300 & + land_deltaT * flxEngU(i,j) / land_dzF(k)
301 ELSE
302 C- Frozen level: incoming water flux goes directly into run-off
303 land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
304 & + flxkup(i,j)
305 land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
306 & + flxEngU(i,j)
307 ENDIF
308 C- prepare fluxes for next level:
309 flxkup(i,j) = 0. _d 0
310 flxEngU(i,j) = 0. _d 0
311
312 ELSE
313
314 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
315 C- Diffusion flux of water, lower interface (k,k+1):
316 IF ( k.EQ.land_nLev .OR.
317 & land_groundT(i,j,kp1,bi,bj).LT.0. _d 0 ) THEN
318 #endif /* LAND_OLD_VERSION */
319 C- no Diffusion of water if one level is frozen :
320 flxkdw(i,j) = 0. _d 0
321 flxEngL = 0. _d 0
322 ELSE
323 flxkdw(i,j) = fieldCapac*
324 & ( land_groundW(i,j,k,bi,bj)
325 & -land_groundW(i,j,kp1,bi,bj) )
326 & / land_wTauDiff
327 C- energy flux associated with water flux: take upwind Temp
328 IF ( flxkdw(i,j).GE.0. ) THEN
329 flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
330 & *land_groundT(i,j,k,bi,bj)
331 ELSE
332 flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
333 & *land_groundT(i,j,kp1,bi,bj)
334 ENDIF
335 ENDIF
336
337 C- Step forward soil moisture, level k :
338 groundWnp1 = land_groundW(i,j,k,bi,bj)
339 & + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac
340
341 #ifdef LAND_DEBUG
342 IF(dBug)write(6,1010)'LAND_STEPFWD: grdW-1,fx_ku,kd,grdW-1='
343 & ,5,land_groundW(i,j,k,bi,bj)-1.,
344 & flxkup(i,j),flxkdw(i,j),groundWnp1-1.
345 #endif
346
347 C- Water in excess will leave as run-off or go to level below
348 land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1)
349 grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )
350 & *fieldCapac/land_deltaT
351
352 C- Run off: fraction 1-fractRunOff enters level below
353 land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
354 & + fractRunOff*grdWexcess
355 C- prepare fluxes for next level:
356 flxkup(i,j) = flxkdw(i,j)
357 & + (1. _d 0-fractRunOff)*grdWexcess
358
359 IF ( land_calc_snow ) THEN
360 enthalpGrdW = land_rhoLiqW*land_CpWater
361 & *land_groundT(i,j,k,bi,bj)
362 C-- Account for water fluxes in energy budget: update ground Enthalpy
363 land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
364 & + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW
365 & )*land_deltaT/land_dzF(k)
366
367 land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
368 & + fractRunOff*grdWexcess*enthalpGrdW
369 C- prepare fluxes for next level:
370 flxEngU(i,j) = flxEngL
371 & + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW
372 ENDIF
373 #ifdef LAND_DEBUG
374 IF (dBug) write(6,1010) 'LAND_STEPFWD: Temp,FlxE,FlxW=',
375 & 7, land_groundT(i,j,k,bi,bj), flxEngU(i,j), flxkup(i,j)
376 #endif
377 ENDIF
378 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
379 #ifdef LAND_DEBUG
380 IF (dBug) write(6,1010) 'LAND_STEPFWD: RO,enRO=',
381 & 8, land_runOff(i,j,bi,bj),land_enRnOf(i,j,bi,bj)
382 #endif
383
384 ENDIF
385 ENDDO
386 ENDDO
387
388 ENDDO
389 C-- step forward ground Water: end
390 ENDIF
391
392 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
393
394 IF ( land_calc_grT ) THEN
395 C-- Compute ground temperature from enthalpy (if not already done):
396
397 DO k=1,land_nLev
398 DO j=1,sNy
399 DO i=1,sNx
400 C- Ground Heat capacity, layer k:
401 mWater = land_rhoLiqW*land_waterCap
402 & *land_groundW(i,j,k,bi,bj)
403 mWater = MAX( mWater, 0. _d 0 )
404 grd_HeatCp = land_heatCs + land_CpWater*mWater
405 C temperature below freezing:
406 temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
407 & / grd_HeatCp
408 C temperature above freezing:
409 temp_af = land_enthalp(i,j,k,bi,bj) / grd_HeatCp
410 #ifdef LAND_OLD_VERSION
411 land_enthalp(i,j,k,bi,bj) =
412 & grd_HeatCp*land_groundT(i,j,k,bi,bj)
413 #else
414 land_groundT(i,j,k,bi,bj) =
415 & MIN( temp_bf, MAX(temp_af, 0. _d 0) )
416 #endif
417 ENDDO
418 ENDDO
419 ENDDO
420
421 IF ( land_impl_grT ) THEN
422 DO j=1,sNy
423 DO i=1,sNx
424 IF ( land_hSnow(i,j,bi,bj).GT.0. _d 0 ) THEN
425 land_skinT(i,j,bi,bj) = MIN(land_skinT(i,j,bi,bj), 0. _d 0)
426 ELSE
427 land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
428 ENDIF
429 ENDDO
430 ENDDO
431 ELSE
432 DO j=1,sNy
433 DO i=1,sNx
434 land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
435 ENDDO
436 ENDDO
437 ENDIF
438
439 C-- Compute ground temperature: end
440 ENDIF
441
442 #endif /* ALLOW_LAND */
443
444 RETURN
445 END

  ViewVC Help
Powered by ViewVC 1.1.22