/[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.3 by jmc, Fri May 14 16:14:48 2004 UTC revision 1.7 by jmc, Mon Oct 1 13:31:15 2007 UTC
# Line 21  C     !USES: Line 21  C     !USES:
21    
22  C     == Global variables ===  C     == Global variables ===
23  C-- size for MITgcm & Land package :  C-- size for MITgcm & Land package :
24  #include "LAND_SIZE.h"  #include "LAND_SIZE.h"
25    
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "LAND_PARAMS.h"  #include "LAND_PARAMS.h"
# Line 69  C     dhSnowMx     :: potential snow inc Line 69  C     dhSnowMx     :: potential snow inc
69  C     dhSnow       :: effective snow increase [m]  C     dhSnow       :: effective snow increase [m]
70  C     mIceDt       :: ground-ice growth rate (<- excess of snow) [kg/m2/s]  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, enthalpGrdW        _RL grd_HeatCp, enthalpGrdW
73        _RL fieldCapac, mWater        _RL fieldCapac, mWater
74        _RL groundWnp1, grdWexcess, fractRunOff        _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 flxEngU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL flxEngU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RL flxEngL, temp_af, temp_bf, mPmE, enWfx, enGr1        _RL flxEngL, temp_af, temp_bf, mPmE, enWfx, enGr1
79        _RL mSnow, dMsn, snowPrec          _RL mSnow, dMsn, snowPrec
80        _RL hNewSnow, dhSnowMx, dhSnow, mIceDt, ageFac        _RL hNewSnow, dhSnowMx, dhSnow, mIceDt, ageFac
81        INTEGER i,j,k,kp1        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        IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN
91  C--   Step forward ground temperature:  C--   Step forward ground temperature:
92    
# Line 107  C-     Thermal conductivity flux, lower Line 114  C-     Thermal conductivity flux, lower
114            flxkdw(i,j) = land_grdLambda*            flxkdw(i,j) = land_grdLambda*
115       &             ( land_groundT(i,j,k,bi,bj)       &             ( land_groundT(i,j,k,bi,bj)
116       &              -land_groundT(i,j,kp1,bi,bj) )       &              -land_groundT(i,j,kp1,bi,bj) )
117       &            *land_rec_dzC(kp1)       &            *land_rec_dzC(kp1)
118    
119  C-     Step forward ground enthalpy, level k :  C-     Step forward ground enthalpy, level k :
120            land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)            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)       &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j))/land_dzF(k)
# Line 123  C--   step forward ground temperature: e Line 130  C--   step forward ground temperature: e
130    
131  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
132    
133  #ifdef LAND_OLD_VERSION        IF ( land_calc_grW .OR. land_calc_snow ) THEN
       IF ( .TRUE. ) THEN  
 #else  
       IF ( land_calc_grW ) THEN  
 #endif  
134  C--   Initialize run-off arrays.  C--   Initialize run-off arrays.
135          DO j=1,sNy          DO j=1,sNy
136           DO i=1,sNx           DO i=1,sNx
# Line 135  C--   Initialize run-off arrays. Line 138  C--   Initialize run-off arrays.
138             land_enRnOf(i,j,bi,bj) = 0. _d 0             land_enRnOf(i,j,bi,bj) = 0. _d 0
139           ENDDO           ENDDO
140          ENDDO          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:  C--   need (later on) ground temp. to be consistent with updated enthalpy:
149          DO k=1,land_nLev          DO k=1,land_nLev
150           DO j=1,sNy           DO j=1,sNy
# Line 142  C--   need (later on) ground temp. to be Line 152  C--   need (later on) ground temp. to be
152             IF ( land_frc(i,j,bi,bj).GT.0. ) THEN             IF ( land_frc(i,j,bi,bj).GT.0. ) THEN
153              mWater = land_rhoLiqW*land_waterCap              mWater = land_rhoLiqW*land_waterCap
154       &              *land_groundW(i,j,k,bi,bj)       &              *land_groundW(i,j,k,bi,bj)
155                mWater = MAX( mWater, 0. _d 0 )
156              grd_HeatCp = land_heatCs + land_CpWater*mWater              grd_HeatCp = land_heatCs + land_CpWater*mWater
157              temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)              temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
158       &                                           / grd_HeatCp       &                                           / grd_HeatCp
159              temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp              temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
160              land_groundT(i,j,k,bi,bj) =              land_groundT(i,j,k,bi,bj) =
161       &              MIN( temp_bf, MAX(temp_af, 0. _d 0) )       &              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             ENDIF
169            ENDDO            ENDDO
170           ENDDO           ENDDO
# Line 163  C--   Step forward Snow thickness (also Line 180  C--   Step forward Snow thickness (also
180             mPmE  = land_Pr_m_Ev(i,j,bi,bj)             mPmE  = land_Pr_m_Ev(i,j,bi,bj)
181             enWfx = land_EnWFlux(i,j,bi,bj)             enWfx = land_EnWFlux(i,j,bi,bj)
182             enGr1 = land_enthalp(i,j,1,bi,bj)*land_dzF(1)             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:  C-    snow aging:
190             land_snowAge(i,j,bi,bj) =             land_snowAge(i,j,bi,bj) =
191       &         ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac )       &         ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac )
192             IF ( enWfx.LT.0. ) THEN             IF ( enWfx.LT.0. ) THEN
193  C-    snow precip in excess (Snow > Evap) :  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  C     => start to melt (until ground at freezing point) and then accumulate
195              snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )              snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 )
196              snowPrec = MAX( snowPrec*recip_Lfreez , 0. _d 0 )  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              mPmE = mPmE - snowPrec
200              flxEngU(i,j) = enWfx + land_Lfreez*snowPrec              flxEngU(i,j) = enWfx + land_Lfreez*snowPrec
201              hNewSnow = land_deltaT * snowPrec / land_rhoSnow              hNewSnow = land_deltaT * snowPrec / land_rhoSnow
# Line 180  C-    refresh snow age: Line 205  C-    refresh snow age:
205  C-    update snow thickness:  C-    update snow thickness:
206  c           land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow  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  C     glacier & ice-sheet missing: excess of snow put directly into run-off
208              dhSnowMx = MAX( 0. _d 0,              dhSnowMx = MAX( 0. _d 0,
209       &                      land_hMaxSnow - land_hSnow(i,j,bi,bj) )       &                      land_hMaxSnow - land_hSnow(i,j,bi,bj) )
210              dhSnow = MIN( hNewSnow, dhSnowMx )              dhSnow = MIN( hNewSnow, dhSnowMx )
211              land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow              land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow
212              mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT              mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT
213              land_runOff(i,j,bi,bj) = mIceDt/land_rhoLiqW              land_runOff(i,j,bi,bj) =  mIceDt
214              land_enRnOf(i,j,bi,bj) = -mIceDt*land_Lfreez              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             ELSE
221  C-    rain precip (whatever Evap is) or Evap exceeds snow precip :  C-    rain precip (whatever Evap is) or Evap of snow exceeds snow precip:
222  C     => snow melts or sublimates  C     => snow melts or sublimates
223  c           snowMelt = MIN( enWfx*recip_Lfreez ,  c           snowMelt = MIN( enWfx*recip_Lfreez ,
224  c    &                 land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT )  c    &                 land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT )
# Line 201  c    &                 land_hSnow(i,j,bi Line 231  c    &                 land_hSnow(i,j,bi
231              ELSE              ELSE
232                flxEngU(i,j) = 0. _d 0                flxEngU(i,j) = 0. _d 0
233                land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj)                land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj)
234       &                              - dMsn / land_rhoSnow       &                              - dMsn / land_rhoSnow
235              ENDIF              ENDIF
236  c           IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge  c           IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge
237              mPmE = mPmE + dMsn/land_deltaT              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             ENDIF
244             flxkup(i,j) = mPmE/land_rhoLiqW             flxkup(i,j) = mPmE/land_rhoLiqW
245  c          land_Pr_m_Ev(i,j,bi,bj) = mPmE  c          land_Pr_m_Ev(i,j,bi,bj) = mPmE
246             IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 )             IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 )
247       &          land_snowAge(i,j,bi,bj) = 0. _d 0       &          land_snowAge(i,j,bi,bj) = 0. _d 0
248  C-    avoid negative (but very small, < 1.e-34) hSnow that occurs because  C-    avoid negative (but very small, < 1.e-34) hSnow that occurs because
249  C      of truncation error. Might need to rewrite this part.  C      of truncation error. Might need to rewrite this part.
250  c          IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN  c          IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN
251  c             land_hSnow(i,j,bi,bj)   = 0. _d 0  c             land_hSnow(i,j,bi,bj)   = 0. _d 0
# Line 246  C--   Step forward ground Water: Line 281  C--   Step forward ground Water:
281         DO j=1,sNy         DO j=1,sNy
282          DO i=1,sNx          DO i=1,sNx
283           IF ( land_frc(i,j,bi,bj).GT.0. ) THEN           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  #ifdef LAND_OLD_VERSION
289            IF ( .TRUE. ) THEN            IF ( .TRUE. ) THEN
# Line 263  C-     Step forward soil moisture (& ent Line 301  C-     Step forward soil moisture (& ent
301             ELSE             ELSE
302  C-     Frozen level: incoming water flux goes directly into run-off  C-     Frozen level: incoming water flux goes directly into run-off
303              land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)              land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
304       &                             + flxkup(i,j)       &                             + flxkup(i,j)*land_rhoLiqW
305              land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)              land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
306       &                             + flxEngU(i,j)       &                             + flxEngU(i,j)
307             ENDIF             ENDIF
# Line 295  C-     energy flux associated with water Line 333  C-     energy flux associated with water
333       &                 *land_groundT(i,j,kp1,bi,bj)       &                 *land_groundT(i,j,kp1,bi,bj)
334               ENDIF               ENDIF
335             ENDIF             ENDIF
336    
337  C-     Step forward soil moisture, level k :  C-     Step forward soil moisture, level k :
338             groundWnp1 = land_groundW(i,j,k,bi,bj)             groundWnp1 = land_groundW(i,j,k,bi,bj)
339       &       + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac       &       + 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  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)             land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1)
349             grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )             grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) )
# Line 307  C-     Water in excess will leave as run Line 351  C-     Water in excess will leave as run
351    
352  C-     Run off: fraction 1-fractRunOff enters level below  C-     Run off: fraction 1-fractRunOff enters level below
353             land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)             land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj)
354       &                        + fractRunOff*grdWexcess       &                        + fractRunOff*grdWexcess*land_rhoLiqW
355  C-     prepare fluxes for next level:  C-     prepare fluxes for next level:
356             flxkup(i,j) = flxkdw(i,j)             flxkup(i,j) = flxkdw(i,j)
357       &              + (1. _d 0-fractRunOff)*grdWexcess       &              + (1. _d 0-fractRunOff)*grdWexcess
# Line 317  C-     prepare fluxes for next level: Line 361  C-     prepare fluxes for next level:
361       &                   *land_groundT(i,j,k,bi,bj)       &                   *land_groundT(i,j,k,bi,bj)
362  C--    Account for water fluxes in energy budget: update ground Enthalpy  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)              land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj)
364       &         + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW       &         + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW
365       &           )*land_deltaT/land_dzF(k)       &           )*land_deltaT/land_dzF(k)
366    
367              land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)              land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj)
# Line 326  C-     prepare fluxes for next level: Line 370  C-     prepare fluxes for next level:
370              flxEngU(i,j) = flxEngL              flxEngU(i,j) = flxEngL
371       &              + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW       &              + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW
372             ENDIF             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            ENDIF
378  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  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           ENDIF
385          ENDDO          ENDDO
# Line 348  C--   Compute ground temperature from en Line 400  C--   Compute ground temperature from en
400  C-     Ground Heat capacity, layer k:  C-     Ground Heat capacity, layer k:
401            mWater = land_rhoLiqW*land_waterCap            mWater = land_rhoLiqW*land_waterCap
402       &            *land_groundW(i,j,k,bi,bj)       &            *land_groundW(i,j,k,bi,bj)
403              mWater = MAX( mWater, 0. _d 0 )
404            grd_HeatCp = land_heatCs + land_CpWater*mWater            grd_HeatCp = land_heatCs + land_CpWater*mWater
405  C         temperature below freezing:  C         temperature below freezing:
406            temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)            temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)
# Line 355  C         temperature below freezing: Line 408  C         temperature below freezing:
408  C         temperature above freezing:  C         temperature above freezing:
409            temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp            temp_af =  land_enthalp(i,j,k,bi,bj) / grd_HeatCp
410  #ifdef LAND_OLD_VERSION  #ifdef LAND_OLD_VERSION
411            land_enthalp(i,j,k,bi,bj) =            land_enthalp(i,j,k,bi,bj) =
412       &          grd_HeatCp*land_groundT(i,j,k,bi,bj)       &          grd_HeatCp*land_groundT(i,j,k,bi,bj)
413  #else  #else
414            land_groundT(i,j,k,bi,bj) =            land_groundT(i,j,k,bi,bj) =
415       &            MIN( temp_bf, MAX(temp_af, 0. _d 0) )       &            MIN( temp_bf, MAX(temp_af, 0. _d 0) )
416  #endif  #endif
417           ENDDO           ENDDO
# Line 375  C         temperature above freezing: Line 428  C         temperature above freezing:
428            ENDIF            ENDIF
429           ENDDO           ENDDO
430          ENDDO          ENDDO
431         ELSE         ELSE
432          DO j=1,sNy          DO j=1,sNy
433           DO i=1,sNx           DO i=1,sNx
434             land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)             land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj)
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22