/[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.6 by jmc, Thu Jun 3 16:43:14 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    
# Line 114  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 157  C--   need (later on) ground temp. to be Line 157  C--   need (later on) ground temp. to be
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  #ifdef LAND_DEBUG
163              dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt              dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt
# Line 187  C--   Step forward Snow thickness (also Line 187  C--   Step forward Snow thickness (also
187       &       mPmE,enWfx,enGr1/land_deltaT,land_hSnow(i,j,bi,bj)       &       mPmE,enWfx,enGr1/land_deltaT,land_hSnow(i,j,bi,bj)
188  #endif  #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 ( > Evap of snow) or snow prec & Evap of Liq.Water:  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  C-    snow accumulation cannot be larger that net precip  C-    snow accumulation cannot be larger that net precip
197              snowPrec = MAX( 0. _d 0 ,              snowPrec = MAX( 0. _d 0 ,
# Line 205  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  #ifdef LAND_DEBUG
216              IF (dBug) write(6,1010)              IF (dBug) write(6,1010)
# Line 231  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
# Line 245  c           IF (mPmE.GT.0.) land_snowAge Line 245  c           IF (mPmE.GT.0.) land_snowAge
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 301  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 333  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
# Line 351  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 361  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 408  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 428  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.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22