baregroundfluxes.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine BareGroundFluxes (tg,     thm,   qg,    thv,   z0mg,   &
                            z0hg,   z0qg,  dqgdT, htvp,  beta,   &
                            zii,    ur,    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:
! Compute sensible and latent fluxes and their derivatives with respect
! to ground temperature using ground temperatures from previous time step.
!
! Method:
!
! 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: BareGroundFluxes.F90,v 1.6.6.3 2002/04/27 15:38:36 erik Exp $
!-----------------------------------------------------------------------

  use precision
  use clmtype
  use clm_varcon, only : cpair, vkc, grav
  use shr_const_mod, only : SHR_CONST_RGAS
  implicit none

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

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

  real(r8), intent(in) :: tg    ! ground surface temperature [K]
  real(r8), intent(in) :: thm   ! intermediate variable (forc_t+0.0098*forc_hgt_t) [K]
  real(r8), intent(in) :: qg    ! specific humidity at ground surface [kg/kg]
  real(r8), intent(in) :: thv   ! virtual potential temperature [K]
  real(r8), intent(in) :: z0mg  ! roughness length, momentum [m]
  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) :: beta  ! coefficient of convective velocity [-]
  real(r8), intent(in) :: zii   ! convective boundary height [m]
  real(r8), intent(in) :: ur    ! wind speed at reference height [m/s]

  real(r8), intent(inout) :: z0hg   ! roughness length, sensible heat [m]
  real(r8), intent(inout) :: z0qg   ! roughness length, latent heat [m]
  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/m**2/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----------------------------------------------------

  integer nmozsgn  ! number of times moz changes sign [-]
  integer niters   ! maximum number of iterations for surface temperature [-]
  integer iter     ! iteration index [-]
  real(r8) zldis   ! reference height "minus" zero displacement height [m]
  real(r8) displa  ! displacement height [m]
  real(r8) zeta    ! dimensionless height used in Monin-Obukhov theory [-]
  real(r8) wc      ! convective velocity [m/s]
  real(r8) dth     ! diff of virtual temp. between ref. height and surface [K]
  real(r8) dthv    ! diff of vir. poten. temp. between ref. height and surface [K]
  real(r8) dqh     ! diff of humidity between ref. height and surface [kg/kg]
  real(r8) obu     ! Monin-Obukhov length [m]
  real(r8) um      ! wind speed including the stability effect [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 [K]
  real(r8) qstar   ! moisture scaling parameter [kg/kg]
  real(r8) thvstar ! virtual potential temperature scaling parameter [K]
  real(r8) cf      ! heat transfer coefficient from leaves [-]
  real(r8) ram     ! aerodynamic resistance [s/m]
  real(r8) rah     ! thermal resistance [s/m]
  real(r8) raw     ! moisture resistance [s/m]
  real(r8) raih    ! temporary variable [kg/m2/s]
  real(r8) raiw    ! temporary variable [kg/m2/s]
  real(r8) obuold  ! monin-obukhov length from previous iteration [m]

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

!
! Compute sensible and latent fluxes and their derivatives with respect
! to ground temperature using ground temperatures from previous time step.
!

!
! Initialization variables
!

     dlrad  = 0.
     ulrad  = 0.

     nmozsgn = 0
     obuold = 0.
     dth   = thm-tg
     dqh   = clm%forc_q-qg
     dthv  = dth*(1.+0.61*clm%forc_q)+0.61*clm%forc_th*dqh
     zldis = clm%forc_hgt_u-0.

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

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

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

     niters=3
     do iter = 1, niters

        displa = 0.0_r8
        call FrictionVelocity(displa,z0mg,z0hg,z0qg,obu, &
                              iter,ur,um,ustar,temp1,temp2,clm)

        tstar = temp1*dth
        qstar = temp2*dqh
        z0hg = z0mg/exp(0.13 * (ustar*z0mg/1.5e-5)**0.45)
        z0qg = z0hg

        thvstar=tstar*(1.+0.61*clm%forc_q) + 0.61*clm%forc_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) EXIT

        obuold = obu

     enddo                       ! end stability iteration

!
! Determine aerodynamic resistances
!

     ram    = 1./(ustar*ustar/um)
     rah    = 1./(temp1*ustar)
     raw    = 1./(temp2*ustar)
     raih   = clm%forc_rho*cpair/rah
     raiw   = clm%forc_rho/raw

!
! Get derivative of fluxes with respect to ground temperature
!

     cgrnds = raih
     cgrndl = raiw*dqgdT
     cgrnd  = cgrnds + htvp*cgrndl

!
! Surface fluxes of momentum, sensible and latent heat
! using ground temperatures from previous time step
!

     clm%taux   = -clm%forc_rho*clm%forc_u/ram
     clm%tauy   = -clm%forc_rho*clm%forc_v/ram
     clm%eflx_sh_grnd  = -raih*dth
     clm%qflx_evap_soi  = -raiw*dqh
     clm%eflx_sh_tot  = clm%eflx_sh_grnd
     clm%qflx_evap_tot  = clm%qflx_evap_soi

!
! 2 m height air temperature
!

     clm%t_ref2m = (tg+temp1*dth * 1./vkc *log((2.+z0hg)/z0hg))

!
! Variables needed by history tape
!

     clm%t_veg = clm%forc_t
     clm%btran = 0.
     clm%rootr(:) = 0.
     cf = clm%forc_pbot/(SHR_CONST_RGAS*0.001*thm)*1.e06
     clm%rssun = 1./1.e15 * cf
     clm%rssha = 1./1.e15 * cf

end subroutine BareGroundFluxes