canopyfluxes.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

     tlbef = clm%t_veg
     del2 = del

!
! Determine friction velocity, and potential temperature and humidity
! profiles of the surface boundary layer
!

     call FrictionVelocity (clm%displa, z0mv,  z0hv,  z0qv,  obu, &
                            itlef+1, ur, um, ustar, temp1, temp2, clm)

!
! Determine aerodynamic resistances
!

     ram(1)=1./(ustar*ustar/um)
     rah(1)=1./(temp1*ustar) 
     raw(1)=1./(temp2*ustar) 

!
! Bulk boundary layer resistance of leaves
!

     uaf = um*sqrt( 1./(ram(1)*um) )
     cf = 0.01/(sqrt(uaf)*sqrt(clm%dleaf))
     rb = 1./(cf*uaf)

!
! Aerodynamic resistances raw and rah between heights zpd+z0h and z0hg.
! If no vegetation, rah(2)=0 because zpd+z0h = z0hg.
! (Dickinson et al., 1993, pp.54)
!

     ram(2) = 0.               ! not used
     rah(2) = 1./(clm%csoilc*uaf)
     raw(2) = rah(2) 

!
! Stomatal resistances for sunlit and shaded fractions of canopy.
! Done each iteration to account for differences in eah, tv.
!

     svpts = el                        ! pa
     eah = clm%forc_pbot * qaf / 0.622 ! pa

     call Stomata(mpe,        clm%parsun, svpts,     eah,      thm,        &
                  o2,         co2,        clm%btran, rb,       clm%rssun,  &
                  clm%psnsun, clm%qe25,   clm%vcmx25,clm%mp,   clm%c3psn,  &
                  clm       )
                                                      
     call Stomata(mpe,        clm%parsha, svpts,     eah,      thm,        &
                  o2,         co2,        clm%btran, rb,       clm%rssha,  &
                  clm%psnsha, clm%qe25,   clm%vcmx25,clm%mp,   clm%c3psn,  &
                  clm       )

!
! Heat conductance for air, leaf and ground  
!

     call SensibleHCond(rah(1), rb,   rah(2), wta,  wtl,  &
                        wtg,    wta0, wtl0,   wtg0, wtal, &
                        wtga,   wtgl, clm     )

!
! Fraction of potential evaporation from leaf
!

     if (clm%fdry .gt. 0.0) then
        rppdry  = clm%fdry*rb*(clm%laisun/(rb+clm%rssun) + clm%laisha/(rb+clm%rssha))/ &
                  clm%elai
     else
        rppdry = 0.0
     endif
     efpot = clm%forc_rho*wtl*(qsatl-qaf)

     if (efpot > 0.) then
       if (clm%btran > 1.e-10) then

        clm%qflx_tran_veg = efpot*rppdry
        rpp = rppdry + clm%fwet

!
! No transpiration if btran below 1.e-10
!

       else
        rpp = clm%fwet
        clm%qflx_tran_veg = 0.
       endif

!
! Check total evapotranspiration from leaves
!

       rpp = min(rpp, (clm%qflx_tran_veg+clm%h2ocan/clm%dtime)/efpot)

     else

!
! No transpiration if potential evaporation less than zero
!

       rpp = 1.
       clm%qflx_tran_veg = 0.

     endif

!
! Update conductances for changes in rpp 
! Latent heat conductances for ground and leaf.
! Air has same conductance for both sensible and latent heat.
!

     call LatentHCond(raw(1), rb,    raw(2), rpp,   wtaq,  &
                      wtlq,   wtgq,  wtaq0,  wtlq0, wtgq0, &
                      wtalq,  wtgaq, wtglq,  clm    ) 

     dc1 = clm%forc_rho*cpair*wtl
     dc2 = hvap*clm%forc_rho*wtlq

     efsh = dc1*(wtga*clm%t_veg-wtg0*tg-wta0*thm)
     efe = dc2*(wtgaq*qsatl-wtgq0*qg-wtaq0*clm%forc_q)

!
! Evaporation flux from foliage
!

     erre = 0.
     if (efe*efeb < 0.0) then
        efeold = efe
        efe  = 0.1*efeold
        erre = efe - efeold
     endif
     clm%dt_veg = (clm%sabv + air + bir*clm%t_veg**4 + cir*tg**4 - efsh - efe) &
          / (- 4.*bir*clm%t_veg**3 +dc1*wtga +dc2*wtgaq*qsatldT)
     clm%t_veg = tlbef + clm%dt_veg
     dels = clm%t_veg-tlbef
     del  = abs(dels)
     err = 0.
     if (del > delmax) then
        clm%dt_veg = delmax*dels/del
        clm%t_veg = tlbef + clm%dt_veg
        err = clm%sabv + air + bir*tlbef**3*(tlbef + 4.*clm%dt_veg) &
             + cir*tg**4 - (efsh + dc1*wtga*clm%dt_veg)          &
             - (efe + dc2*wtgaq*qsatldT*clm%dt_veg)
     endif

!
! Fluxes from leaves to canopy space
!

     efpot = clm%forc_rho*wtl*(wtgaq*(qsatl+qsatldT*clm%dt_veg) &
          -wtgq0*qg-wtaq0*clm%forc_q)
     clm%qflx_evap_veg = rpp*efpot

!
! Calculation of evaporative potentials (efpot) and
! interception losses; flux in kg m**-2 s-1.  ecidif 
! holds the excess energy if all intercepted water is evaporated
! during the timestep.  This energy is later added to the
! sensible heat flux.
!

     ecidif = 0.
     if (efpot > 0. .AND. clm%btran > 1.e-10) then
        clm%qflx_tran_veg = efpot*rppdry
     else
        clm%qflx_tran_veg = 0.
     endif
     ecidif = max(0._r8, clm%qflx_evap_veg-clm%qflx_tran_veg-clm%h2ocan/clm%dtime)
     clm%qflx_evap_veg = min(clm%qflx_evap_veg,clm%qflx_tran_veg+clm%h2ocan/clm%dtime)

!
! The energy loss due to above limits is added to 
! the sensible heat flux.
!

     clm%eflx_sh_veg = efsh + dc1*wtga*clm%dt_veg + err + erre +hvap*ecidif

!
! Re-calculate saturated vapor pressure, specific humidity, and their
! derivatives at the leaf surface
!

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

!
! Update vegetation/ground surface temperature, canopy air temperature, 
! canopy vapor pressure, aerodynamic temperature, and
! Monin-Obukhov stability parameter for next iteration. 
!

     taf = wtg0*tg + wta0*thm + wtl0*clm%t_veg
     qaf = wtlq0*qsatl+wtgq0*qg+clm%forc_q*wtaq0

!
! Update Monin-Obukhov length and wind speed including the stability effect
!

     dth = thm-taf       
     dqh = clm%forc_q-qaf

     tstar=temp1*dth
     qstar=temp2*dqh

     dthv=dth*(1.+0.61*clm%forc_q)+0.61*th*dqh

     thvstar=tstar*(1.+0.61*clm%forc_q) + 0.61*th*qstar
     zeta=zldis*vkc*grav*thvstar/(ustar**2*thv)

     if (zeta >= 0.) then     !stable
        zeta = min(2._r8,max(zeta,0.01_r8))
        um = max(ur,0.1_r8)
     else                     !unstable
        zeta = max(-100._r8,min(zeta,-0.01_r8))
        wc = beta*(-grav*ustar*thvstar*zii/thv)**0.333
        um = sqrt(ur*ur+wc*wc)
     endif
     obu = zldis/zeta

     if (obuold*obu < 0.) nmozsgn = nmozsgn+1
     if (nmozsgn >= 4) then 
        obu = zldis/(-0.01)
     endif

     obuold = obu

!
! Test for convergence
!

     itlef = itlef+1
     if (itlef > itmin) then
        dele = abs(efe-efeb)
        efeb = efe
        det  = max(del,del2)
        if (det < dtmin .AND. dele < dlemin) exit 
     endif

  enddo ITERATION     ! End stability iteration

!
! Energy balance check in canopy
!

  err = clm%sabv + air + bir*tlbef**3*(tlbef + 4.*clm%dt_veg) &
       + cir*tg**4 - clm%eflx_sh_veg - hvap*clm%qflx_evap_veg
  if (abs(err) > 0.1) then
     write(6,*) 'energy balance in canopy X',err
  endif

!
! Fluxes from ground to canopy space 
!

  delt  = wtal*tg-wtl0*clm%t_veg-wta0*thm
  delq  = wtalq*qg-wtlq0*qsatl-wtaq0*clm%forc_q
  clm%taux  = -clm%forc_rho*clm%forc_u/ram(1)
  clm%tauy  = -clm%forc_rho*clm%forc_v/ram(1)
  clm%eflx_sh_grnd = cpair*clm%forc_rho*wtg*delt
  clm%qflx_evap_soi = clm%forc_rho*wtgq*delq

!
! 2 m height air temperature
!

  clm%t_ref2m   = clm%t_ref2m + (taf + temp1*dth * &
       1./vkc *log((2.+z0hv)/z0hv))

!
! Downward longwave radiation below the canopy    
!

  dlrad = (1.-emv)*emg*clm%forc_lwrad &
       + emv*emg * sb * &
       tlbef**3*(tlbef + 4.*clm%dt_veg)

!
! Upward longwave radiation above the canopy    
!

  ulrad = ( (1.-emg)*(1.-emv)*(1.-emv)*clm%forc_lwrad &
       + emv*(1.+(1.-emg)*(1.-emv))*sb * tlbef**3 &
       *(tlbef + 4.*clm%dt_veg) + emg *(1.-emv) *sb * tg**4)

!
! Derivative of soil energy fluxes with respect to soil temperature
!

  cgrnds = cgrnds + cpair*clm%forc_rho*wtg*wtal
  cgrndl = cgrndl + clm%forc_rho*wtgq*wtalq*dqgdT
  cgrnd  = cgrnds  + cgrndl*htvp

!
! Update canopy water storage (kg/m2) 
!

  clm%h2ocan = max(0._r8,clm%h2ocan + (clm%qflx_tran_veg-clm%qflx_evap_veg)*clm%dtime)

end subroutine CanopyFluxes