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

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

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

revision 1.2 by jmc, Thu Mar 11 14:42:00 2004 UTC revision 1.3 by jmc, Fri May 14 16:14:48 2004 UTC
# Line 44  C     == Local variables == Line 44  C     == Local variables ==
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]
# Line 62  C     mSnow        :: mass of snow Line 65  C     mSnow        :: mass of snow
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
# Line 118  C---+----1----+----2----+----3----+----4 Line 126  C---+----1----+----2----+----3----+----4
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
# Line 157  C     => start to melt (until ground at Line 172  C     => start to melt (until ground at
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
# Line 173  c    &                 land_hSnow(i,j,bi Line 197  c    &                 land_hSnow(i,j,bi
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
# Line 199  c          ENDIF Line 223  c          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
# Line 219  C--   Step forward ground Water: Line 243  C--   Step forward ground Water:
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
# Line 281  C--   step forward ground Water: end Line 339  C--   step forward ground Water: end
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
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22