C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/land/land_stepfwd.F,v 1.7 2007/10/01 13:31:15 jmc Exp $ C $Name: checkpoint63a $ #include "LAND_OPTIONS.h" CBOP C !ROUTINE: LAND_STEPFWD C !INTERFACE: SUBROUTINE LAND_STEPFWD( I land_frc, bi, bj, myTime, myIter, myThid) C !DESCRIPTION: \bv C *==========================================================* C | S/R LAND_STEPFWD C | o Land model main S/R: step forward land variables C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables === C-- size for MITgcm & Land package : #include "LAND_SIZE.h" #include "EEPARAMS.h" #include "LAND_PARAMS.h" #include "LAND_VARS.h" C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C land_frc :: land fraction [0-1] C bi,bj :: Tile index C myTime :: Current time of simulation ( s ) C myIter :: Current iteration number in simulation C myThid :: Number of this instance of the routine _RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) INTEGER bi, bj, myIter, myThid _RL myTime CEOP #ifdef ALLOW_LAND C == Local variables == 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 groundWnp1 :: hold temporary future soil moisture [] C grdWexcess :: ground water in excess [m/s] 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 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] C enWfx :: hold temporary energy flux of Precip [W/m2] C enGr1 :: ground enthalpy of level 1 [J/m2] C mSnow :: mass of snow [kg/m2] 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, 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 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: DO k=1,land_nLev kp1 = MIN(k+1,land_nLev) IF (k.EQ.1) THEN DO j=1,sNy DO i=1,sNx flxkup(i,j) = land_HeatFlx(i,j,bi,bj) 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 C- Thermal conductivity flux, lower interface (k,k+1): flxkdw(i,j) = land_grdLambda* & ( land_groundT(i,j,k,bi,bj) & -land_groundT(i,j,kp1,bi,bj) ) & *land_rec_dzC(kp1) C- Step forward ground enthalpy, level k : land_enthalp(i,j,k,bi,bj) = land_enthalp(i,j,k,bi,bj) & + land_deltaT * (flxkup(i,j)-flxkdw(i,j))/land_dzF(k) ENDIF ENDDO ENDDO ENDDO C-- step forward ground temperature: end ENDIF 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_grW ) THEN #endif C-- need (later on) ground temp. to be consistent with updated enthalpy: DO k=1,land_nLev DO j=1,sNy DO i=1,sNx 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 ENDDO ENDIF IF ( land_calc_snow ) THEN C-- Step forward Snow thickness (also account for rain temperature) ageFac = 1. _d 0 - land_deltaT/timeSnowAge DO j=1,sNy DO i=1,sNx IF ( land_frc(i,j,bi,bj).GT.0. ) THEN 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 ( > 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 ) C- snow accumulation cannot be larger that net precip snowPrec = MAX( 0. _d 0 , & MIN( snowPrec*recip_Lfreez, mPmE ) ) mPmE = mPmE - snowPrec flxEngU(i,j) = enWfx + land_Lfreez*snowPrec hNewSnow = land_deltaT * snowPrec / land_rhoSnow 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_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 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 ) mSnow = land_hSnow(i,j,bi,bj)*land_rhoSnow dMsn = enWfx*recip_Lfreez*land_deltaT IF ( dMsn .GE. mSnow ) THEN dMsn = mSnow land_hSnow(i,j,bi,bj) = 0. _d 0 flxEngU(i,j) = enWfx - land_Lfreez*mSnow/land_deltaT ELSE 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 IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) & land_snowAge(i,j,bi,bj) = 0. _d 0 C- avoid negative (but very small, < 1.e-34) hSnow that occurs because C of truncation error. Might need to rewrite this part. c IF ( land_hSnow(i,j,bi,bj).LE. 0. _d 0 ) THEN c land_hSnow(i,j,bi,bj) = 0. _d 0 c land_snowAge(i,j,bi,bj) = 0. _d 0 c ENDIF ENDIF ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx flxkup(i,j) = land_Pr_m_Ev(i,j,bi,bj)/land_rhoLiqW flxEngU(i,j) = 0. _d 0 ENDDO ENDDO ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF (land_calc_grW) THEN C-- Step forward ground Water: DO k=1,land_nLev IF (k.EQ.land_nLev) THEN kp1 = k fractRunOff = 1. _d 0 ELSE kp1 = k+1 fractRunOff = land_fractRunOff ENDIF fieldCapac = land_waterCap*land_dzF(k) 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_rhoLiqW 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): 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) & + land_deltaT * (flxkup(i,j)-flxkdw(i,j)) / fieldCapac #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 land_runOff(i,j,bi,bj) = land_runOff(i,j,bi,bj) & + fractRunOff*grdWexcess*land_rhoLiqW C- prepare fluxes for next level: 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 ENDDO ENDDO C-- step forward ground Water: end ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF ( land_calc_grT ) THEN C-- Compute ground temperature from enthalpy (if not already done): DO k=1,land_nLev DO j=1,sNy DO i=1,sNx 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) & / grd_HeatCp C temperature above freezing: temp_af = land_enthalp(i,j,k,bi,bj) / grd_HeatCp #ifdef LAND_OLD_VERSION land_enthalp(i,j,k,bi,bj) = & grd_HeatCp*land_groundT(i,j,k,bi,bj) #else land_groundT(i,j,k,bi,bj) = & MIN( temp_bf, MAX(temp_af, 0. _d 0) ) #endif ENDDO ENDDO ENDDO IF ( land_impl_grT ) THEN DO j=1,sNy DO i=1,sNx IF ( land_hSnow(i,j,bi,bj).GT.0. _d 0 ) THEN land_skinT(i,j,bi,bj) = MIN(land_skinT(i,j,bi,bj), 0. _d 0) ELSE land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj) ENDIF ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx land_skinT(i,j,bi,bj) = land_groundT(i,j,1,bi,bj) ENDDO ENDDO ENDIF C-- Compute ground temperature: end ENDIF #endif /* ALLOW_LAND */ RETURN END