canopyfluxes.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <preproc.h>

subroutine CanopyFluxes (z0mv,   z0hv,  z0qv,  thm,   th,     &
                         thv,    tg,    qg,    dqgdT, htvp,   &
                         emv,    emg,   dlrad, ulrad, cgrnds, &
                         cgrndl, cgrnd, clm    )

!-----------------------------------------------------------------------
!
!  CLMCLMCLMCLMCLMCLMCLMCLMCLMCL  A community developed and sponsored, freely
!  L                           M  available land surface process model.
!  M --COMMUNITY LAND MODEL--  C
!  C                           L
!  LMCLMCLMCLMCLMCLMCLMCLMCLMCLM
!
!-----------------------------------------------------------------------
! Purpose:
! This subroutine:
! 1. Calculates the leaf temperature: 
! 2. Calculates the leaf fluxes, transpiration, photosynthesis and 
!    updates the canopy water storage.
!
! Method:
! Use Newton-Raphson iteration to solve for the foliage 
! temperature that balances the surface energy budget:
!
! f(t_veg) = Net radiation - Sensible - Latent = 0
! f(t_veg) + d(f)/d(t_veg) * dt_veg = 0     (*)
!
! Note:
! (1) In solving for t_veg, t_grnd is given from the previous timestep.
! (2) The partial derivatives of aerodynamical resistances, which cannot 
!     be determined analytically, are ignored for d(H)/dT and d(LE)/dT
! (3) The weighted stomatal resistance of sunlit and shaded foliage is used 
! (4) Canopy air temperature and humidity are derived from => Hc + Hg = Ha
!                                                          => Ec + Eg = Ea
! (5) Energy loss is due to: numerical truncation of energy budget equation
!     (*); and "ecidif" (see the code) which is put into the sensible heat 
! (6) The convergence criteria: the difference, del = t_veg(n+1)-t_veg(n) and 
!     del2 = t_veg(n)-t_veg(n-1) less than 0.01 K, and the difference of 
!     water flux from the leaf between the iteration step (n+1) and (n) 
!     less than 0.1 W/m2; or max iterations of 40.
!
! Author:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision 
! April 2002: Vertenstein/Oleson/Levis; Final form
!
!-----------------------------------------------------------------------
! $Id: CanopyFluxes.F90,v 1.6.6.2 2002/04/27 15:38:37 erik Exp $
!-----------------------------------------------------------------------

  use precision
  use clmtype
  use clm_varcon, only : sb, cpair, hvap, vkc, grav, denice, denh2o, tfrz
  implicit none

!----Arguments----------------------------------------------------------

  type (clm1d), intent(inout) :: clm	 !CLM 1-D Module

  real(r8), intent(in) :: z0mv  ! roughness length, momentum [m]
  real(r8), intent(in) :: z0hv  ! roughness length, sensible heat [m]
  real(r8), intent(in) :: z0qv  ! roughness length, latent heat [m]
  real(r8), intent(in) :: thm   ! intermediate variable (forc_t+0.0098*forc_hgt_t) [K]
  real(r8), intent(in) :: th    ! potential temperature [K]
  real(r8), intent(in) :: thv   ! virtual potential temperature [K]
  real(r8), intent(in) :: tg    ! ground surface temperature [K]
  real(r8), intent(in) :: qg    ! specific humidity at ground surface [kg/kg]
  real(r8), intent(in) :: dqgdT ! temperature derivative of "qg"
  real(r8), intent(in) :: htvp  ! latent heat of evaporation (/sublimation) [J/kg]
  real(r8), intent(in) :: emv   ! ground emissivity [-]
  real(r8), intent(in) :: emg   ! vegetation emissivity [-]

  real(r8), intent(inout) :: cgrnd  ! deriv. of soil energy flux wrt to soil temp [W/m2/K]
  real(r8), intent(inout) :: cgrndl ! deriv. of soil sensible heat flux wrt soil temp [W/m2/K]
  real(r8), intent(inout) :: cgrnds ! deriv. of soil latent heat flux wrt soil temp [W/m2/K]

  real(r8), intent(out) :: dlrad    ! downward longwave radiation below the canopy [W/m2]
  real(r8), intent(out) :: ulrad    ! upward longwave radiation above the canopy [W/m2]

!----Local Variables----------------------------------------------------

  real(r8) zldis       ! reference height "minus" zero displacement height [m]
  real(r8) zii         ! convective boundary layer height [m]
  real(r8) zeta        ! dimensionless height used in Monin-Obukhov theory
  real(r8) beta        ! coefficient of convective velocity [-]
  real(r8) wc          ! convective velocity [m/s]
  real(r8) dth         ! diff of virtual temp. between ref. height and surface 
  real(r8) dthv        ! diff of vir. poten. temp. between ref. height and sfc
  real(r8) dqh         ! diff of humidity between ref. height and surface
  real(r8) obu         ! Monin-Obukhov length [m]
  real(r8) um          ! wind speed including the stability effect [m/s]
  real(r8) ur          ! wind speed at reference height [m/s]
  real(r8) uaf         ! velocity of air within foliage [m/s]
  real(r8) temp1       ! relation for potential temperature profile
  real(r8) temp2       ! relation for specific humidity profile
  real(r8) ustar       ! friction velocity [m/s]
  real(r8) tstar       ! temperature scaling parameter
  real(r8) qstar       ! moisture scaling parameter
  real(r8) thvstar     ! virtual potential temperature scaling parameter
  real(r8) taf         ! air temperature within canopy space [K]
  real(r8) qaf         ! humidity of canopy air [kg/kg]
  real(r8) rpp         ! fraction of potential evaporation from leaf [-]
  real(r8) rppdry      ! fraction of potential evaporation through transp [-]
  real(r8) cf          ! heat transfer coefficient from leaves [-]
  real(r8) rb          ! leaf boundary layer resistance [s/m]
  real(r8) ram(2)      ! aerodynamic resistance [s/m]
  real(r8) rah(2)      ! thermal resistance [s/m]
  real(r8) raw(2)      ! moisture resistance [s/m]
  real(r8) wta         ! heat conductance for air [m/s]
  real(r8) wtg         ! heat conductance for ground [m/s]
  real(r8) wtl         ! heat conductance for leaf [m/s]
  real(r8) wta0        ! normalized heat conductance for air [-]
  real(r8) wtl0        ! normalized heat conductance for leaf [-]
  real(r8) wtg0        ! normalized heat conductance for ground [-]
  real(r8) wtal        ! normalized heat conductance for air and leaf [-]
  real(r8) wtgl        ! normalized heat conductance for leaf and ground [-]
  real(r8) wtga        ! normalized heat conductance for air and ground  [-]
  real(r8) wtaq        ! latent heat conductance for air [m/s]
  real(r8) wtlq        ! latent heat conductance for leaf [m/s]
  real(r8) wtgq        ! latent heat conductance for ground [m/s]
  real(r8) wtaq0       ! normalized latent heat conductance for air [-]
  real(r8) wtlq0       ! normalized latent heat conductance for leaf [-]
  real(r8) wtgq0       ! normalized heat conductance for ground [-]
  real(r8) wtalq       ! normalized latent heat cond. for air and leaf [-]
  real(r8) wtglq       ! normalized latent heat cond. for leaf and ground [-]
  real(r8) wtgaq       ! normalized latent heat cond. for air and ground [-]
  real(r8) el          ! vapor pressure on leaf surface [pa]
  real(r8) deldT       ! derivative of "el" wrt. "t_veg" [pa/K]
  real(r8) qsatl       ! leaf specific humidity [kg/kg]
  real(r8) qsatldT     ! derivative of "qsatl" wrt "t_veg"
  real(r8) air,bir,cir ! temporary variables for longwave radiation
  real(r8) dc1,dc2     ! derivative of energy flux [W/m2/K]
  real(r8) delt        ! temporary variable
  real(r8) delq        ! temporary variable
  real(r8) del         ! absolute change in leaf temp in current iteration [K]
  real(r8) del2        ! change in leaf temperature in previous iteration [K]
  real(r8) dele        ! change in latent heat flux from leaf [K]
  real(r8) delmax      ! maximum change in leaf temperature [K]
  real(r8) dels        ! change in leaf temperature in current iteration [K]
  real(r8) det         ! maximum leaf temp. change in two consecutive iter [K]
  real(r8) dlemin      ! max limit for energy flux convergence [w/m2]
  real(r8) dtmin       ! max limit for temperature convergence [K]
  real(r8) efeb        ! latent heat flux from leaf (previous iter) [mm/s]
  real(r8) efeold      ! latent heat flux from leaf (previous iter) [mm/s]
  real(r8) efpot       ! potential latent energy flux [kg/m2/s]
  real(r8) efe         ! water flux from leaf [mm/s]
  real(r8) efsh        ! sensible heat from leaf [mm/s]
  real(r8) obuold      ! monin-obukhov length from previous iteration [m]
  real(r8) tlbef       ! leaf temperature from previous iteration [K]
  real(r8) ecidif      ! excess energy [W/m2]
  real(r8) err         ! balance error
  real(r8) erre        ! balance error
  real(r8) co2         ! atmospheric co2 concentration (pa)
  real(r8) o2          ! atmospheric o2 concentration (pa)
  real(r8) svpts       ! saturation vapor pressure at t_veg (pa)
  real(r8) eah         ! canopy air vapor pressure (pa)
  real(r8) s_node      ! vol_liq/eff_porosity
  real(r8) smp_node    ! matrix potential
  real(r8) vol_ice(1:nlevsoi)      ! partial volume of ice lens in layer
  real(r8) eff_porosity(1:nlevsoi) ! effective porosity in layer
  real(r8) vol_liq(1:nlevsoi)      ! partial volume of liquid water in layer
  real(r8) rresis(1:nlevsoi)       ! soil water contribution to root resistance

! Constant atmospheric co2 and o2
  real(r8) po2                   ! partial pressure  o2 (mol/mol)
  real(r8) pco2                  ! partial pressure co2 (mol/mol)
  data po2,pco2 /0.209,355.e-06/

  real(r8) :: mpe = 1.e-6        ! prevents overflow error if division by zero

  integer i       ! loop index
  integer itlef   ! counter for leaf temperature iteration [-]
  integer itmax   ! maximum number of iteration [-]
  integer itmin   ! minimum number of iteration [-]
  integer nmozsgn ! number of times stability changes sign

!----End Variable List--------------------------------------------------

!
! Initialization
!

  del   = 0.0  ! change in leaf temperature from previous iteration
  itlef = 0    ! counter for leaf temperature iteration
  efeb  = 0.0  ! latent heat flux from leaf for previous iteration

  wtlq = 0.0
  wtlq0 = 0.0
  wtgq0 = 0.0
  wtalq = 0.0
  wtgaq = 0.0
  wtglq = 0.0
  wtaq = 0.0
  wtgq = 0.0
  wtaq0 = 0.0
  wtlq0 = 0.0
  wtgq0 = 0.0
  wtalq = 0.0
  wtgaq = 0.0
  wtglq = 0.0

!
! Assign iteration parameters
!

  delmax = 1.0  ! maximum change in leaf temperature
  itmax  = 40   ! maximum number of iterations
  itmin  = 2    ! minimum number of iterations
  dtmin  = 0.01 ! max limit for temperature convergence
  dlemin = 0.1  ! max limit for energy flux convergence

!
! Effective porosity of soil, partial volume of ice and liquid (needed
! for btran)
!

  do i = 1,nlevsoi
     vol_ice(i) = min(clm%watsat(i), clm%h2osoi_ice(i)/(clm%dz(i)*denice))
     eff_porosity(i) = clm%watsat(i)-vol_ice(i)
     vol_liq(i) = min(eff_porosity(i), clm%h2osoi_liq(i)/(clm%dz(i)*denh2o))
  enddo

!
! Root resistance factors
!

  clm%btran = 1.e-10
  do i = 1,nlevsoi
     if (clm%t_soisno(i) > tfrz) then
        s_node = max(vol_liq(i)/eff_porosity(i),0.01_r8)
        smp_node = max(clm%smpmax, -clm%sucsat(i)*s_node**(-clm%bsw(i)))
        rresis(i) = (1.-smp_node/clm%smpmax)/(1.+clm%sucsat(i)/clm%smpmax)
        clm%rootr(i) = clm%rootfr(i)*rresis(i)
        clm%btran = clm%btran + clm%rootr(i)
     else
        clm%rootr(i) = 0.
     endif
  enddo

!
! Normalize root resistances to get layer contribution to ET
!

  do i = 1,nlevsoi
     clm%rootr(i)  = clm%rootr(i)/clm%btran
  enddo

!
! Net absorbed longwave radiation by canopy and ground
! =air+bir*t_veg**4+cir*t_grnd**4
!

  air =   emv * (1.+(1.-emv)*(1.-emg)) * clm%forc_lwrad
  bir = - (2.-emv*(1.-emg)) * emv * sb
  cir =   emv*emg*sb

!
! Saturated vapor pressure, specific humidity, and their derivatives
! at the leaf surface
!

  call QSat (clm%t_veg, clm%forc_pbot, el, deldT, qsatl, &
             qsatldT)

!
! Determine atmospheric co2 and o2
!

  co2 = pco2*clm%forc_pbot
  o2  = po2*clm%forc_pbot

!
! Initialize flux profile
!

  nmozsgn = 0
  obuold = 0.
  zii=1000.         ! m  (pbl height)
  beta=1.           ! -  (in computing W_*)

  taf = (tg + thm)/2.
  qaf = (clm%forc_q+qg)/2.

  ur = max(1.0_r8,sqrt(clm%forc_u*clm%forc_u+clm%forc_v*clm%forc_v))
  dth = thm-taf
  dqh = clm%forc_q-qaf
  dthv = dth*(1.+0.61*clm%forc_q)+0.61*th*dqh
  zldis = clm%forc_hgt_u - clm%displa

!
! Initialize Monin-Obukhov length and wind speed including stability effect
!

  call MoninObukIni(ur, thv, dthv, zldis, z0mv, &
                    um, obu  )

!
! Begin stability iteration
!

  ITERATION : do while (itlef <= itmax)