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

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

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


Revision 1.5 - (hide annotations) (download)
Fri May 21 13:41:02 2004 UTC (20 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint53c_post, checkpoint53d_pre
Changes since 1.4: +6 -4 lines
snow accumulation not larger than net precipitation

1 jmc 1.5 C $Header: /u/gcmpack/MITgcm/pkg/land/land_stepfwd.F,v 1.4 2004/05/15 20:43:27 jmc Exp $
2 jmc 1.1 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 jmc 1.2 C grd_HeatCp :: Heat capacity of the ground [J/m3/K]
47 jmc 1.3 C enthalpGrdW :: enthalpy of ground water [J/m3]
48 jmc 1.1 C fieldCapac :: field capacity (of water) [m]
49 jmc 1.2 C mWater :: water content of the ground [kg/m3]
50 jmc 1.3 C groundWnp1 :: hold temporary future soil moisture []
51     C grdWexcess :: ground water in excess [m/s]
52 jmc 1.1 C fractRunOff :: fraction of water in excess which leaves as runoff
53 jmc 1.2 C flxkup :: downward flux of water, upper interface (k-1,k)
54     C flxdwn :: downward flux of water, lower interface (k,k+1)
55 jmc 1.3 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 jmc 1.2 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 jmc 1.3 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 jmc 1.2 C ageFac :: snow aging factor [1]
72 jmc 1.3 _RL grd_HeatCp, enthalpGrdW
73     _RL fieldCapac, mWater
74     _RL groundWnp1, grdWexcess, fractRunOff
75 jmc 1.1 _RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76     _RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 jmc 1.3 _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 jmc 1.1 INTEGER i,j,k,kp1
82    
83 jmc 1.2 IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN
84 jmc 1.1 C-- Step forward ground temperature:
85    
86     DO k=1,land_nLev
87     kp1 = MIN(k+1,land_nLev)
88    
89     IF (k.EQ.1) THEN
90     DO j=1,sNy
91     DO i=1,sNx
92     flxkup(i,j) = land_HeatFlx(i,j,bi,bj)
93     ENDDO
94     ENDDO
95     ELSE
96     DO j=1,sNy
97     DO i=1,sNx
98     flxkup(i,j) = flxkdw(i,j)
99     ENDDO
100     ENDDO
101     ENDIF
102    
103     DO j=1,sNy
104     DO i=1,sNx
105     IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
106     C- Thermal conductivity flux, lower interface (k,k+1):
107     flxkdw(i,j) = land_grdLambda*
108     & ( land_groundT(i,j,k,bi,bj)
109     & -land_groundT(i,j,kp1,bi,bj) )
110     & *land_rec_dzC(kp1)
111    
112 jmc 1.2 C- Step forward ground enthalpy, level k :
113     land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
114     & + land_deltaT * (flxkup(i,j)-flxkdw(i,j))/land_dzF(k)
115 jmc 1.1
116     ENDIF
117     ENDDO
118     ENDDO
119    
120     ENDDO
121     C-- step forward ground temperature: end
122     ENDIF
123    
124 jmc 1.2 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
125    
126 jmc 1.4 IF ( land_calc_grW .OR. land_calc_snow ) THEN
127 jmc 1.3 C-- Initialize run-off arrays.
128     DO j=1,sNy
129     DO i=1,sNx
130     land_runOff(i,j,bi,bj) = 0. _d 0
131     land_enRnOf(i,j,bi,bj) = 0. _d 0
132     ENDDO
133     ENDDO
134 jmc 1.4 ENDIF
135    
136     #ifdef LAND_OLD_VERSION
137     IF ( .TRUE. ) THEN
138     #else
139     IF ( land_calc_grW ) THEN
140     #endif
141 jmc 1.2 C-- need (later on) ground temp. to be consistent with updated enthalpy:
142     DO k=1,land_nLev
143     DO j=1,sNy
144     DO i=1,sNx
145     IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
146     mWater = land_rhoLiqW*land_waterCap
147     & *land_groundW(i,j,k,bi,bj)
148     grd_HeatCp = land_heatCs + land_CpWater*mWater
149     temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
150     & / grd_HeatCp
151     temp_af = land_enthalp(i,j,k,bi,bj) / grd_HeatCp
152     land_groundT(i,j,k,bi,bj) =
153     & MIN( temp_bf, MAX(temp_af, 0. _d 0) )
154     ENDIF
155     ENDDO
156     ENDDO
157     ENDDO
158     ENDIF
159    
160     IF ( land_calc_snow ) THEN
161     C-- Step forward Snow thickness (also account for rain temperature)
162     ageFac = 1. _d 0 - land_deltaT/timeSnowAge
163     DO j=1,sNy
164     DO i=1,sNx
165     IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
166     mPmE = land_Pr_m_Ev(i,j,bi,bj)
167     enWfx = land_EnWFlux(i,j,bi,bj)
168     enGr1 = land_enthalp(i,j,1,bi,bj)*land_dzF(1)
169     C- snow aging:
170     land_snowAge(i,j,bi,bj) =
171     & ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac )
172     IF ( enWfx.LT.0. ) THEN
173 jmc 1.5 C- snow precip in excess ( > Evap of snow) or snow prec & Evap of Liq.Water:
174 jmc 1.2 C => start to melt (until ground at freezing point) and then accumulate
175     snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )
176 jmc 1.5 C- snow accumulation cannot be larger that net precip
177     snowPrec = MAX( 0. _d 0 ,
178     & MIN( snowPrec*recip_Lfreez, mPmE ) )
179 jmc 1.2 mPmE = mPmE - snowPrec
180 jmc 1.3 flxEngU(i,j) = enWfx + land_Lfreez*snowPrec
181 jmc 1.2 hNewSnow = land_deltaT * snowPrec / land_rhoSnow
182     C- refresh snow age:
183     land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj)
184     & *EXP( -hNewSnow/hNewSnowAge )
185 jmc 1.3 C- update snow thickness:
186     c land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow
187     C glacier & ice-sheet missing: excess of snow put directly into run-off
188     dhSnowMx = MAX( 0. _d 0,
189     & land_hMaxSnow - land_hSnow(i,j,bi,bj) )
190     dhSnow = MIN( hNewSnow, dhSnowMx )
191     land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow
192     mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT
193     land_runOff(i,j,bi,bj) = mIceDt/land_rhoLiqW
194     land_enRnOf(i,j,bi,bj) = -mIceDt*land_Lfreez
195 jmc 1.2 ELSE
196 jmc 1.5 C- rain precip (whatever Evap is) or Evap of snow exceeds snow precip:
197 jmc 1.2 C => snow melts or sublimates
198     c snowMelt = MIN( enWfx*recip_Lfreez ,
199     c & land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT )
200     mSnow = land_hSnow(i,j,bi,bj)*land_rhoSnow
201     dMsn = enWfx*recip_Lfreez*land_deltaT
202     IF ( dMsn .GE. mSnow ) THEN
203     dMsn = mSnow
204     land_hSnow(i,j,bi,bj) = 0. _d 0
205 jmc 1.3 flxEngU(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT
206 jmc 1.2 ELSE
207 jmc 1.3 flxEngU(i,j) = 0. _d 0
208 jmc 1.2 land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj)
209     & - dMsn / land_rhoSnow
210     ENDIF
211     c IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge
212     mPmE = mPmE + dMsn/land_deltaT
213     ENDIF
214     flxkup(i,j) = mPmE/land_rhoLiqW
215     c land_Pr_m_Ev(i,j,bi,bj) = mPmE
216     IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 )
217     & land_snowAge(i,j,bi,bj) = 0. _d 0
218     C- avoid negative (but very small, < 1.e-34) hSnow that occurs because
219     C of truncation error. Might need to rewrite this part.
220     c IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN
221     c land_hSnow(i,j,bi,bj) = 0. _d 0
222     c land_snowAge(i,j,bi,bj) = 0. _d 0
223     c ENDIF
224     ENDIF
225     ENDDO
226     ENDDO
227     ELSE
228     DO j=1,sNy
229     DO i=1,sNx
230     flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW
231 jmc 1.3 flxEngU(i,j) = 0. _d 0
232 jmc 1.2 ENDDO
233     ENDDO
234     ENDIF
235    
236     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
237    
238 jmc 1.1 IF (land_calc_grW) THEN
239     C-- Step forward ground Water:
240    
241     DO k=1,land_nLev
242     IF (k.EQ.land_nLev) THEN
243     kp1 = k
244     fractRunOff = 1. _d 0
245     ELSE
246     kp1 = k+1
247     fractRunOff = land_fractRunOff
248     ENDIF
249     fieldCapac = land_waterCap*land_dzF(k)
250    
251     DO j=1,sNy
252     DO i=1,sNx
253     IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
254 jmc 1.3
255     #ifdef LAND_OLD_VERSION
256     IF ( .TRUE. ) THEN
257     IF ( k.EQ.land_nLev ) THEN
258     #else
259     IF ( land_groundT(i,j,k,bi,bj).LT.0. _d 0 ) THEN
260     C- Frozen level: only account for upper level fluxes
261     IF ( flxkup(i,j) .LT. 0. _d 0 ) THEN
262     C- Step forward soil moisture (& enthapy), level k :
263     land_groundW(i,j,k,bi,bj) = land_groundW(i,j,k,bi,bj)
264     & + land_deltaT * flxkup(i,j) / fieldCapac
265     IF ( land_calc_snow )
266     & land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
267     & + land_deltaT * flxEngU(i,j) / land_dzF(k)
268     ELSE
269     C- Frozen level: incoming water flux goes directly into run-off
270     land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
271     & + flxkup(i,j)
272     land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
273     & + flxEngU(i,j)
274     ENDIF
275     C- prepare fluxes for next level:
276     flxkup(i,j) = 0. _d 0
277     flxEngU(i,j) = 0. _d 0
278    
279     ELSE
280    
281     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
282 jmc 1.2 C- Diffusion flux of water, lower interface (k,k+1):
283 jmc 1.3 IF ( k.EQ.land_nLev .OR.
284     & land_groundT(i,j,kp1,bi,bj).LT.0. _d 0 ) THEN
285     #endif /* LAND_OLD_VERSION */
286     C- no Diffusion of water if one level is frozen :
287     flxkdw(i,j) = 0. _d 0
288     flxEngL = 0. _d 0
289     ELSE
290     flxkdw(i,j) = fieldCapac*
291     & ( land_groundW(i,j,k,bi,bj)
292     & -land_groundW(i,j,kp1,bi,bj) )
293     & / land_wTauDiff
294     C- energy flux associated with water flux: take upwind Temp
295     IF ( flxkdw(i,j).GE.0. ) THEN
296     flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
297     & *land_groundT(i,j,k,bi,bj)
298     ELSE
299     flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater
300     & *land_groundT(i,j,kp1,bi,bj)
301     ENDIF
302     ENDIF
303 jmc 1.1
304     C- Step forward soil moisture, level k :
305 jmc 1.3 groundWnp1 = land_groundW(i,j,k,bi,bj)
306 jmc 1.2 & + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac
307 jmc 1.3
308     C- Water in excess will leave as run-off or go to level below
309     land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1)
310     grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )
311     & *fieldCapac/land_deltaT
312 jmc 1.1
313     C- Run off: fraction 1-fractRunOff enters level below
314 jmc 1.3 land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
315     & + fractRunOff*grdWexcess
316     C- prepare fluxes for next level:
317     flxkup(i,j) = flxkdw(i,j)
318     & + (1. _d 0-fractRunOff)*grdWexcess
319    
320     IF ( land_calc_snow ) THEN
321     enthalpGrdW = land_rhoLiqW*land_CpWater
322     & *land_groundT(i,j,k,bi,bj)
323     C-- Account for water fluxes in energy budget: update ground Enthalpy
324     land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
325     & + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW
326     & )*land_deltaT/land_dzF(k)
327    
328     land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
329     & + fractRunOff*grdWexcess*enthalpGrdW
330     C- prepare fluxes for next level:
331     flxEngU(i,j) = flxEngL
332     & + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW
333     ENDIF
334 jmc 1.2 ENDIF
335 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
336 jmc 1.2
337 jmc 1.1 ENDIF
338     ENDDO
339     ENDDO
340    
341     ENDDO
342     C-- step forward ground Water: end
343 jmc 1.2 ENDIF
344    
345     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
346    
347 jmc 1.3 IF ( land_calc_grT ) THEN
348     C-- Compute ground temperature from enthalpy (if not already done):
349 jmc 1.2
350     DO k=1,land_nLev
351     DO j=1,sNy
352     DO i=1,sNx
353     C- Ground Heat capacity, layer k:
354     mWater = land_rhoLiqW*land_waterCap
355     & *land_groundW(i,j,k,bi,bj)
356     grd_HeatCp = land_heatCs + land_CpWater*mWater
357     C temperature below freezing:
358     temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
359     & / grd_HeatCp
360     C temperature above freezing:
361     temp_af = land_enthalp(i,j,k,bi,bj) / grd_HeatCp
362     #ifdef LAND_OLD_VERSION
363     land_enthalp(i,j,k,bi,bj) =
364     & grd_HeatCp*land_groundT(i,j,k,bi,bj)
365     #else
366     land_groundT(i,j,k,bi,bj) =
367     & MIN( temp_bf, MAX(temp_af, 0. _d 0) )
368     #endif
369     ENDDO
370     ENDDO
371     ENDDO
372    
373     IF ( land_impl_grT ) THEN
374     DO j=1,sNy
375     DO i=1,sNx
376     IF ( land_hSnow(i,j,bi,bj).GT.0. _d 0 ) THEN
377     land_skinT(i,j,bi,bj) = MIN(land_skinT(i,j,bi,bj), 0. _d 0)
378     ELSE
379     land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
380     ENDIF
381     ENDDO
382     ENDDO
383     ELSE
384     DO j=1,sNy
385     DO i=1,sNx
386     land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
387     ENDDO
388     ENDDO
389     ENDIF
390    
391     C-- Compute ground temperature: end
392 jmc 1.1 ENDIF
393    
394     #endif /* ALLOW_LAND */
395    
396     RETURN
397     END

  ViewVC Help
Powered by ViewVC 1.1.22