/[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.1 by jmc, Thu Jun 12 17:54:22 2003 UTC revision 1.2 by jmc, Thu Mar 11 14:42:00 2004 UTC
# Line 27  C-- size for MITgcm & Land package : Line 27  C-- size for MITgcm & Land package :
27  #include "LAND_PARAMS.h"  #include "LAND_PARAMS.h"
28  #include "LAND_VARS.h"  #include "LAND_VARS.h"
29    
 c #include "PARAMS.h"  
 c #include "GRID.h"  
 c #include "DYNVARS.h"  
   
30  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
31  C     == Routine arguments ==  C     == Routine arguments ==
32  C     land_frc :: land fraction [0-1]  C     land_frc :: land fraction [0-1]
# Line 47  CEOP Line 43  CEOP
43  C     == Local variables ==  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  C     grd_HeatCp   :: Heat capacity of the ground [J/m3/K]
47  C     fieldCapac   :: field capacity (of water) [m]  C     fieldCapac   :: field capacity (of water) [m]
48  C     ground_dTdt  :: ground temperature tendency  C     mWater       :: water content of the ground [kg/m3]
 C     ground_dWdt  :: soil moisture tendency  
 C     flxkup       :: downward flux, upper interface (k-1,k)  
 C     flxdwn       :: downward flux, lower interface (k,k+1)  
49  C     fractRunOff  :: fraction of water in excess which leaves as runoff  C     fractRunOff  :: fraction of water in excess which leaves as runoff
50  C     grdWexcess   :: ground water in excess [m/s]  C     grdWexcess   :: ground water in excess [m/s]
51        _RL grd_HeatCp, fieldCapac, ground_dTdt, ground_dWdt  C     groundWnp1   :: hold temporary future soil moisture
52        _RL fractRunOff, grdWexcess, groundWnp1  C     enthalpGrdW  :: enthalpy of ground water [J/m3]
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     flxEng       :: downward energy flux associated with water flux
56    C     temp_af      :: ground temperature if above freezing
57    C     temp_bf      :: ground temperature if below freezing
58    C     mPmE         :: hold temporary (liquid) Precip minus Evap [kg/m2/s]
59    C     enWfx        :: hold temporary energy flux of Precip [W/m2]
60    C     enGr1        :: ground enthalpy of level 1  [J/m2]
61    C     mSnow        :: mass of snow         [kg/m2]
62    C     dMsn         :: mass of melting snow [kg/m2]
63    C     snowPrec     :: snow precipitation [kg/m2/s]
64    C     hNewSnow     :: fresh snow accumulation [m]
65    C     ageFac       :: snow aging factor [1]
66          _RL grd_HeatCp, fieldCapac, mWater
67          _RL fractRunOff, grdWexcess, groundWnp1, enthalpGrdW
68        _RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69        _RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70          _RL flxEng(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71          _RL temp_af, temp_bf, mPmE, enWfx, enGr1
72          _RL mSnow, dMsn, snowPrec, hNewSnow, ageFac
73        INTEGER i,j,k,kp1        INTEGER i,j,k,kp1
74    
75        IF (land_calc_grT) THEN        IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN
76  C--   Step forward ground temperature:  C--   Step forward ground temperature:
77    
78        DO k=1,land_nLev        DO k=1,land_nLev
# Line 90  C-     Thermal conductivity flux, lower Line 101  C-     Thermal conductivity flux, lower
101       &              -land_groundT(i,j,kp1,bi,bj) )       &              -land_groundT(i,j,kp1,bi,bj) )
102       &            *land_rec_dzC(kp1)       &            *land_rec_dzC(kp1)
103    
104  C-     Ground Heat capacity, layer k:  C-     Step forward ground enthalpy, level k :
105            grd_HeatCp = land_heatCs            land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
106       &               + land_heatCw*land_groundW(i,j,k,bi,bj)       &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j))/land_dzF(k)
      &                            *land_waterCap  
   
 C-     Net temperature tendency  
           ground_dTdt = (flxkup(i,j)-flxkdw(i,j))  
      &                / (grd_HeatCp*land_dzF(k))  
           
 C-     Step forward ground temperature, level k :  
           land_groundT(i,j,k,bi,bj) = land_groundT(i,j,k,bi,bj)  
      &                              + land_deltaT*ground_dTdt  
107    
108           ENDIF           ENDIF
109          ENDDO          ENDDO
# Line 111  C-     Step forward ground temperature, Line 113  C-     Step forward ground temperature,
113  C--   step forward ground temperature: end  C--   step forward ground temperature: end
114        ENDIF        ENDIF
115    
116    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
117    
118    #ifdef LAND_OLD_VERSION
119          IF ( .TRUE. ) THEN
120    #else
121          IF ( land_calc_snow ) THEN
122    #endif
123    C--   need (later on) ground temp. to be consistent with updated enthalpy:
124            DO k=1,land_nLev
125             DO j=1,sNy
126              DO i=1,sNx
127               IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
128                mWater = land_rhoLiqW*land_waterCap
129         &              *land_groundW(i,j,k,bi,bj)
130                grd_HeatCp = land_heatCs + land_CpWater*mWater
131                temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
132         &                                           / grd_HeatCp
133                temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
134                land_groundT(i,j,k,bi,bj) =
135         &              MIN( temp_bf, MAX(temp_af, 0. _d 0) )
136               ENDIF
137              ENDDO
138             ENDDO
139            ENDDO
140          ENDIF
141    
142          IF ( land_calc_snow ) THEN
143    C--   Step forward Snow thickness (also account for rain temperature)
144            ageFac = 1. _d 0 - land_deltaT/timeSnowAge
145            DO j=1,sNy
146             DO i=1,sNx
147              IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
148               mPmE  = land_Pr_m_Ev(i,j,bi,bj)
149               enWfx = land_EnWFlux(i,j,bi,bj)
150               enGr1 = land_enthalp(i,j,1,bi,bj)*land_dzF(1)
151    C-    snow aging:
152               land_snowAge(i,j,bi,bj) =
153         &         ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac )
154               IF ( enWfx.LT.0. ) THEN
155    C-    snow precip in excess (Snow > Evap) :
156    C     => start to melt (until ground at freezing point) and then accumulate
157                snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )
158                snowPrec = MAX( snowPrec*recip_Lfreez , 0. _d 0 )
159                mPmE = mPmE - snowPrec
160                flxEng(i,j) = enWfx + land_Lfreez*snowPrec
161                hNewSnow = land_deltaT * snowPrec / land_rhoSnow
162                land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow
163    C-    refresh snow age:
164                land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj)
165         &                          *EXP( -hNewSnow/hNewSnowAge )
166               ELSE
167    C-    rain precip (whatever Evap is) or Evap exceeds snow precip :
168    C     => snow melts or sublimates
169    c           snowMelt = MIN( enWfx*recip_Lfreez ,
170    c    &                 land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT )
171                mSnow = land_hSnow(i,j,bi,bj)*land_rhoSnow
172                dMsn = enWfx*recip_Lfreez*land_deltaT
173                IF ( dMsn .GE. mSnow ) THEN
174                  dMsn = mSnow
175                  land_hSnow(i,j,bi,bj) = 0. _d 0
176                  flxEng(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT
177                ELSE
178                  flxEng(i,j) = 0. _d 0
179                  land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj)
180         &                              - dMsn / land_rhoSnow
181                ENDIF
182    c           IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge
183                mPmE = mPmE + dMsn/land_deltaT
184               ENDIF
185               flxkup(i,j) = mPmE/land_rhoLiqW
186    c          land_Pr_m_Ev(i,j,bi,bj) = mPmE
187               IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 )
188         &          land_snowAge(i,j,bi,bj) = 0. _d 0
189    C-    avoid negative (but very small, < 1.e-34) hSnow that occurs because
190    C      of truncation error. Might need to rewrite this part.
191    c          IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN
192    c             land_hSnow(i,j,bi,bj)   = 0. _d 0
193    c             land_snowAge(i,j,bi,bj) = 0. _d 0
194    c          ENDIF
195              ENDIF
196             ENDDO
197            ENDDO
198          ELSE
199            DO j=1,sNy
200             DO i=1,sNx
201               flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW
202               flxEng(i,j) = 0. _d 0
203             ENDDO
204            ENDDO
205          ENDIF
206    
207    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
208    
209        IF (land_calc_grW) THEN        IF (land_calc_grW) THEN
210  C--   Step forward ground Water:  C--   Step forward ground Water:
211    
# Line 127  C--   Step forward ground Water: Line 222  C--   Step forward ground Water:
222         IF (k.EQ.1) THEN         IF (k.EQ.1) THEN
223          DO j=1,sNy          DO j=1,sNy
224           DO i=1,sNx           DO i=1,sNx
            flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)  
225             land_runOff(i,j,bi,bj) = 0. _d 0             land_runOff(i,j,bi,bj) = 0. _d 0
226               land_enRnOf(i,j,bi,bj) = 0. _d 0
227           ENDDO           ENDDO
228          ENDDO          ENDDO
229         ELSE         ELSE
# Line 142  C--   Step forward ground Water: Line 237  C--   Step forward ground Water:
237         DO j=1,sNy         DO j=1,sNy
238          DO i=1,sNx          DO i=1,sNx
239           IF ( land_frc(i,j,bi,bj).GT.0. ) THEN           IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
240  C-     Diffusion flux of soil moisture, lower interface (k,k+1):  C-     Diffusion flux of water, lower interface (k,k+1):
241            flxkdw(i,j) = fieldCapac*            flxkdw(i,j) = fieldCapac*
242       &                ( land_groundW(i,j,k,bi,bj)       &                ( land_groundW(i,j,k,bi,bj)
243       &                 -land_groundW(i,j,kp1,bi,bj) )       &                 -land_groundW(i,j,kp1,bi,bj) )
244       &                / land_wTauDiff       &                / land_wTauDiff
245    
 C-     Net soil moisture tendency  
           ground_dWdt = (flxkup(i,j)-flxkdw(i,j)) / fieldCapac  
           
246  C-     Step forward soil moisture, level k :  C-     Step forward soil moisture, level k :
247            groundWnp1 = land_groundW(i,j,k,bi,bj)            groundWnp1 = land_groundW(i,j,k,bi,bj)
248       &               + land_deltaT*ground_dWdt       &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac
249            land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1)            land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1)
250    
251  C-     Run off: fraction 1-fractRunOff enters level below  C-     Run off: fraction 1-fractRunOff enters level below
252            grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )            grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )
253       &                *fieldCapac/land_deltaT       &                *fieldCapac/land_deltaT
           flxkdw(i,j) = flxkdw(i,j)  
      &                + (1. _d 0-fractRunOff)*grdWexcess  
254            land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)            land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
255       &                           + fractRunOff*grdWexcess       &                           + fractRunOff*grdWexcess
256    
257              IF ( land_calc_snow ) THEN
258    C-     account for water fluxes in energy budget:
259               enthalpGrdW = land_enthalp(i,j,k,bi,bj)
260         &                 - land_heatCs*land_groundT(i,j,k,bi,bj)
261               land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
262         &                            + fractRunOff*grdWexcess*enthalpGrdW
263               land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
264         &          + ( flxEng(i,j) - (flxkdw(i,j)+grdWexcess)*enthalpGrdW
265         &            )*land_deltaT/land_dzF(k)
266              ELSE
267               enthalpGrdW = 0. _d 0
268              ENDIF
269    C-     prepare fluxes for next level:
270              flxkdw(i,j) = flxkdw(i,j)
271         &                + (1. _d 0-fractRunOff)*grdWexcess
272              flxEng(i,j) = flxkdw(i,j)*enthalpGrdW
273    
274           ENDIF           ENDIF
275          ENDDO          ENDDO
276         ENDDO         ENDDO
# Line 171  C-     Run off: fraction 1-fractRunOff e Line 279  C-     Run off: fraction 1-fractRunOff e
279  C--   step forward ground Water: end  C--   step forward ground Water: end
280        ENDIF        ENDIF
281    
282    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
283    
284          IF (land_calc_grT ) THEN
285    C--   Compute ground temperature from enthalpy :
286    
287           DO k=1,land_nLev
288            DO j=1,sNy
289             DO i=1,sNx
290    C-     Ground Heat capacity, layer k:
291              mWater = land_rhoLiqW*land_waterCap
292         &            *land_groundW(i,j,k,bi,bj)
293              grd_HeatCp = land_heatCs + land_CpWater*mWater
294    C         temperature below freezing:
295              temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
296         &                                         / grd_HeatCp
297    C         temperature above freezing:
298              temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
299    #ifdef LAND_OLD_VERSION
300              land_enthalp(i,j,k,bi,bj) =
301         &          grd_HeatCp*land_groundT(i,j,k,bi,bj)
302    #else
303              land_groundT(i,j,k,bi,bj) =
304         &            MIN( temp_bf, MAX(temp_af, 0. _d 0) )
305    #endif
306             ENDDO
307            ENDDO
308           ENDDO
309    
310           IF ( land_impl_grT ) THEN
311            DO j=1,sNy
312             DO i=1,sNx
313              IF ( land_hSnow(i,j,bi,bj).GT.0. _d 0 ) THEN
314               land_skinT(i,j,bi,bj) = MIN(land_skinT(i,j,bi,bj), 0. _d 0)
315              ELSE
316               land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
317              ENDIF
318             ENDDO
319            ENDDO
320           ELSE
321            DO j=1,sNy
322             DO i=1,sNx
323               land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
324             ENDDO
325            ENDDO
326           ENDIF
327    
328    C--   Compute ground temperature: end
329          ENDIF
330    
331  #endif /* ALLOW_LAND */  #endif /* ALLOW_LAND */
332    
333        RETURN        RETURN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22