--- MITgcm/pkg/land/land_stepfwd.F 2004/03/11 14:42:00 1.2 +++ MITgcm/pkg/land/land_stepfwd.F 2004/06/03 16:43:14 1.6 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/land/land_stepfwd.F,v 1.2 2004/03/11 14:42:00 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/land/land_stepfwd.F,v 1.6 2004/06/03 16:43:14 jmc Exp $ C $Name: $ #include "LAND_OPTIONS.h" @@ -44,15 +44,18 @@ C i,j,k :: loop counters C kp1 :: k+1 C grd_HeatCp :: Heat capacity of the ground [J/m3/K] +C enthalpGrdW :: enthalpy of ground water [J/m3] C fieldCapac :: field capacity (of water) [m] C mWater :: water content of the ground [kg/m3] -C fractRunOff :: fraction of water in excess which leaves as runoff +C groundWnp1 :: hold temporary future soil moisture [] C grdWexcess :: ground water in excess [m/s] -C groundWnp1 :: hold temporary future soil moisture -C enthalpGrdW :: enthalpy of ground water [J/m3] +C fractRunOff :: fraction of water in excess which leaves as runoff C flxkup :: downward flux of water, upper interface (k-1,k) C flxdwn :: downward flux of water, lower interface (k,k+1) -C flxEng :: downward energy flux associated with water flux +C flxEngU :: downward energy flux associated with water flux (W/m2) +C upper interface (k-1,k) +C flxEngL :: downward energy flux associated with water flux (W/m2) +C lower interface (k,k+1) C temp_af :: ground temperature if above freezing C temp_bf :: ground temperature if below freezing C mPmE :: hold temporary (liquid) Precip minus Evap [kg/m2/s] @@ -62,16 +65,28 @@ C dMsn :: mass of melting snow [kg/m2] C snowPrec :: snow precipitation [kg/m2/s] C hNewSnow :: fresh snow accumulation [m] +C dhSnowMx :: potential snow increase [m] +C dhSnow :: effective snow increase [m] +C mIceDt :: ground-ice growth rate (<- excess of snow) [kg/m2/s] C ageFac :: snow aging factor [1] - _RL grd_HeatCp, fieldCapac, mWater - _RL fractRunOff, grdWexcess, groundWnp1, enthalpGrdW + _RL grd_HeatCp, enthalpGrdW + _RL fieldCapac, mWater + _RL groundWnp1, grdWexcess, fractRunOff _RL flxkup(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL flxkdw(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL flxEng(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL temp_af, temp_bf, mPmE, enWfx, enGr1 - _RL mSnow, dMsn, snowPrec, hNewSnow, ageFac + _RL flxEngU(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL flxEngL, temp_af, temp_bf, mPmE, enWfx, enGr1 + _RL mSnow, dMsn, snowPrec + _RL hNewSnow, dhSnowMx, dhSnow, mIceDt, ageFac INTEGER i,j,k,kp1 +#ifdef LAND_DEBUG + LOGICAL dBug + INTEGER iprt,jprt,lprt + DATA iprt, jprt , lprt / 19 , 20 , 6 / + 1010 FORMAT(A,I3,1P4E11.3) +#endif + IF (land_calc_grT .AND. .NOT.land_impl_grT ) THEN C-- Step forward ground temperature: @@ -115,10 +130,20 @@ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + IF ( land_calc_grW .OR. land_calc_snow ) THEN +C-- Initialize run-off arrays. + 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 + ENDIF + #ifdef LAND_OLD_VERSION IF ( .TRUE. ) THEN #else - IF ( land_calc_snow ) THEN + IF ( land_calc_grW ) THEN #endif C-- need (later on) ground temp. to be consistent with updated enthalpy: DO k=1,land_nLev @@ -127,12 +152,19 @@ IF ( land_frc(i,j,bi,bj).GT.0. ) THEN mWater = land_rhoLiqW*land_waterCap & *land_groundW(i,j,k,bi,bj) + mWater = MAX( mWater, 0. _d 0 ) grd_HeatCp = land_heatCs + land_CpWater*mWater temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater) & / grd_HeatCp temp_af = land_enthalp(i,j,k,bi,bj) / grd_HeatCp land_groundT(i,j,k,bi,bj) = & MIN( temp_bf, MAX(temp_af, 0. _d 0) ) +#ifdef LAND_DEBUG + dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt + IF (dBug) write(6,1010) + & 'LAND_STEPFWD: k,temp,af,bf=', + & k,land_groundT(i,j,k,bi,bj),temp_af,temp_bf +#endif ENDIF ENDDO ENDDO @@ -148,23 +180,45 @@ mPmE = land_Pr_m_Ev(i,j,bi,bj) enWfx = land_EnWFlux(i,j,bi,bj) enGr1 = land_enthalp(i,j,1,bi,bj)*land_dzF(1) +#ifdef LAND_DEBUG + dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt + IF (dBug) write(6,1010) + & 'LAND_STEPFWD:mPmE,enWfx,enGr1/dt,hSnow=',0, + & mPmE,enWfx,enGr1/land_deltaT,land_hSnow(i,j,bi,bj) +#endif C- snow aging: land_snowAge(i,j,bi,bj) = & ( land_deltaT + land_snowAge(i,j,bi,bj)*ageFac ) IF ( enWfx.LT.0. ) THEN -C- snow precip in excess (Snow > Evap) : +C- snow precip in excess ( > Evap of snow) or snow prec & Evap of Liq.Water: C => start to melt (until ground at freezing point) and then accumulate snowPrec = -enWfx -MAX( enGr1/land_deltaT, 0. _d 0 ) - snowPrec = MAX( snowPrec*recip_Lfreez , 0. _d 0 ) +C- snow accumulation cannot be larger that net precip + snowPrec = MAX( 0. _d 0 , + & MIN( snowPrec*recip_Lfreez, mPmE ) ) mPmE = mPmE - snowPrec - flxEng(i,j) = enWfx + land_Lfreez*snowPrec + flxEngU(i,j) = enWfx + land_Lfreez*snowPrec hNewSnow = land_deltaT * snowPrec / land_rhoSnow - land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow C- refresh snow age: land_snowAge(i,j,bi,bj) = land_snowAge(i,j,bi,bj) & *EXP( -hNewSnow/hNewSnowAge ) +C- update snow thickness: +c land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + hNewSnow +C glacier & ice-sheet missing: excess of snow put directly into run-off + dhSnowMx = MAX( 0. _d 0, + & land_hMaxSnow - land_hSnow(i,j,bi,bj) ) + dhSnow = MIN( hNewSnow, dhSnowMx ) + land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) + dhSnow + mIceDt = land_rhoSnow * (hNewSnow-dhSnow) / land_deltaT + land_runOff(i,j,bi,bj) = mIceDt/land_rhoLiqW + land_enRnOf(i,j,bi,bj) = -mIceDt*land_Lfreez +#ifdef LAND_DEBUG + IF (dBug) write(6,1010) + & 'LAND_STEPFWD: 3,snP,mPmE,hNsnw,hSnw=', + & 3,snowPrec,mPmE,hNewSnow,land_hSnow(i,j,bi,bj) +#endif ELSE -C- rain precip (whatever Evap is) or Evap exceeds snow precip : +C- rain precip (whatever Evap is) or Evap of snow exceeds snow precip: C => snow melts or sublimates c snowMelt = MIN( enWfx*recip_Lfreez , c & land_hSnow(i,j,bi,bj)*land_rhoSnow/land_deltaT ) @@ -173,14 +227,19 @@ IF ( dMsn .GE. mSnow ) THEN dMsn = mSnow land_hSnow(i,j,bi,bj) = 0. _d 0 - flxEng(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT + flxEngU(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT ELSE - flxEng(i,j) = 0. _d 0 + flxEngU(i,j) = 0. _d 0 land_hSnow(i,j,bi,bj) = land_hSnow(i,j,bi,bj) & - dMsn / land_rhoSnow ENDIF c IF (mPmE.GT.0.) land_snowAge(i,j,bi,bj) = timeSnowAge mPmE = mPmE + dMsn/land_deltaT +#ifdef LAND_DEBUG + IF (dBug) write(6,1010) + & 'LAND_STEPFWD: 4,dMsn,mPmE,hSnw,enWfx=', + & 4,dMsn,mPmE,land_hSnow(i,j,bi,bj),flxEngU(i,j) +#endif ENDIF flxkup(i,j) = mPmE/land_rhoLiqW c land_Pr_m_Ev(i,j,bi,bj) = mPmE @@ -199,7 +258,7 @@ DO j=1,sNy DO i=1,sNx flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW - flxEng(i,j) = 0. _d 0 + flxEngU(i,j) = 0. _d 0 ENDDO ENDDO ENDIF @@ -219,57 +278,108 @@ ENDIF fieldCapac = land_waterCap*land_dzF(k) - 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 - DO j=1,sNy DO i=1,sNx IF ( land_frc(i,j,bi,bj).GT.0. ) THEN +#ifdef LAND_DEBUG + dBug = bi.eq.lprt .AND. i.EQ.iprt .AND. j.EQ.jprt +#endif + +#ifdef LAND_OLD_VERSION + IF ( .TRUE. ) THEN + IF ( k.EQ.land_nLev ) THEN +#else + IF ( land_groundT(i,j,k,bi,bj).LT.0. _d 0 ) THEN +C- Frozen level: only account for upper level fluxes + IF ( flxkup(i,j) .LT. 0. _d 0 ) THEN +C- Step forward soil moisture (& enthapy), level k : + land_groundW(i,j,k,bi,bj) = land_groundW(i,j,k,bi,bj) + & + land_deltaT * flxkup(i,j) / fieldCapac + IF ( land_calc_snow ) + & land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj) + & + land_deltaT * flxEngU(i,j) / land_dzF(k) + ELSE +C- Frozen level: incoming water flux goes directly into run-off + land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj) + & + flxkup(i,j) + land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj) + & + flxEngU(i,j) + ENDIF +C- prepare fluxes for next level: + flxkup(i,j) = 0. _d 0 + flxEngU(i,j) = 0. _d 0 + + ELSE + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C- Diffusion flux of water, lower interface (k,k+1): - flxkdw(i,j) = fieldCapac* - & ( land_groundW(i,j,k,bi,bj) - & -land_groundW(i,j,kp1,bi,bj) ) - & / land_wTauDiff + IF ( k.EQ.land_nLev .OR. + & land_groundT(i,j,kp1,bi,bj).LT.0. _d 0 ) THEN +#endif /* LAND_OLD_VERSION */ +C- no Diffusion of water if one level is frozen : + flxkdw(i,j) = 0. _d 0 + flxEngL = 0. _d 0 + ELSE + flxkdw(i,j) = fieldCapac* + & ( land_groundW(i,j,k,bi,bj) + & -land_groundW(i,j,kp1,bi,bj) ) + & / land_wTauDiff +C- energy flux associated with water flux: take upwind Temp + IF ( flxkdw(i,j).GE.0. ) THEN + flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater + & *land_groundT(i,j,k,bi,bj) + ELSE + flxEngL = flxkdw(i,j)*land_rhoLiqW*land_CpWater + & *land_groundT(i,j,kp1,bi,bj) + ENDIF + ENDIF C- Step forward soil moisture, level k : - groundWnp1 = land_groundW(i,j,k,bi,bj) + groundWnp1 = land_groundW(i,j,k,bi,bj) & + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac - land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1) + +#ifdef LAND_DEBUG + IF(dBug)write(6,1010)'LAND_STEPFWD: grdW-1,fx_ku,kd,grdW-1=' + & ,5,land_groundW(i,j,k,bi,bj)-1., + & flxkup(i,j),flxkdw(i,j),groundWnp1-1. +#endif + +C- Water in excess will leave as run-off or go to level below + land_groundW(i,j,k,bi,bj) = MIN(1. _d 0, groundWnp1) + grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) ) + & *fieldCapac/land_deltaT C- Run off: fraction 1-fractRunOff enters level below - grdWexcess = ( groundWnp1 - MIN(1. _d 0, groundWnp1) ) - & *fieldCapac/land_deltaT - land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj) - & + fractRunOff*grdWexcess - - IF ( land_calc_snow ) THEN -C- account for water fluxes in energy budget: - enthalpGrdW = land_enthalp(i,j,k,bi,bj) - & - land_heatCs*land_groundT(i,j,k,bi,bj) - land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj) - & + fractRunOff*grdWexcess*enthalpGrdW - land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj) - & + ( flxEng(i,j) - (flxkdw(i,j)+grdWexcess)*enthalpGrdW - & )*land_deltaT/land_dzF(k) - ELSE - enthalpGrdW = 0. _d 0 - ENDIF + land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj) + & + fractRunOff*grdWexcess C- prepare fluxes for next level: - flxkdw(i,j) = flxkdw(i,j) - & + (1. _d 0-fractRunOff)*grdWexcess - flxEng(i,j) = flxkdw(i,j)*enthalpGrdW + flxkup(i,j) = flxkdw(i,j) + & + (1. _d 0-fractRunOff)*grdWexcess + + IF ( land_calc_snow ) THEN + enthalpGrdW = land_rhoLiqW*land_CpWater + & *land_groundT(i,j,k,bi,bj) +C-- Account for water fluxes in energy budget: update ground Enthalpy + land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj) + & + ( flxEngU(i,j) - flxEngL - grdWexcess*enthalpGrdW + & )*land_deltaT/land_dzF(k) + + land_enRnOf(i,j,bi,bj) = land_enRnOf(i,j,bi,bj) + & + fractRunOff*grdWexcess*enthalpGrdW +C- prepare fluxes for next level: + flxEngU(i,j) = flxEngL + & + (1. _d 0-fractRunOff)*grdWexcess*enthalpGrdW + ENDIF +#ifdef LAND_DEBUG + IF (dBug) write(6,1010) 'LAND_STEPFWD: Temp,FlxE,FlxW=', + & 7, land_groundT(i,j,k,bi,bj), flxEngU(i,j), flxkup(i,j) +#endif + ENDIF +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +#ifdef LAND_DEBUG + IF (dBug) write(6,1010) 'LAND_STEPFWD: RO,enRO=', + & 8, land_runOff(i,j,bi,bj),land_enRnOf(i,j,bi,bj) +#endif ENDIF ENDDO @@ -281,8 +391,8 @@ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - IF (land_calc_grT ) THEN -C-- Compute ground temperature from enthalpy : + IF ( land_calc_grT ) THEN +C-- Compute ground temperature from enthalpy (if not already done): DO k=1,land_nLev DO j=1,sNy @@ -290,6 +400,7 @@ C- Ground Heat capacity, layer k: mWater = land_rhoLiqW*land_waterCap & *land_groundW(i,j,k,bi,bj) + mWater = MAX( mWater, 0. _d 0 ) grd_HeatCp = land_heatCs + land_CpWater*mWater C temperature below freezing: temp_bf = (land_enthalp(i,j,k,bi,bj)+land_Lfreez*mWater)