biogeophysics2.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine Biogeophysics2 (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 is the main subroutine to execute the calculation of soil/snow and
! ground temperatures and update surface fluxes based on the new ground
! temperature 
! 
! Method:
! Calling sequence is:
! Biogeophysics2:                    surface biogeophysics driver
!    -> SoilTemperature:             soil/snow and ground temperatures      
!          -> SoilTermProp           thermal conductivities and heat 
!                                     capacities        
!          -> Tridiagonal            tridiagonal matrix solution            
!          -> PhaseChange            phase change of liquid/ice contents        
!
! (1) Snow and soil temperatures
!     o The volumetric heat capacity is calculated as a linear combination 
!       in terms of the volumetric fraction of the constituent phases. 
!     o The thermal conductivity of soil is computed from 
!       the algorithm of Johansen (as reported by Farouki 1981), and the 
!       conductivity of snow is from the formulation used in
!       SNTHERM (Jordan 1991).
!     o Boundary conditions:  
!       F = Rnet - Hg - LEg (top),  F = 0 (base of the soil column).
!     o Soil / snow temperature is predicted from heat conduction 
!       in 10 soil layers and up to 5 snow layers. 
!       The thermal conductivities at the interfaces between two 
!       neighboring layers (j, j+1) are derived from an assumption that 
!       the flux across the interface is equal to that from the node j 
!       to the interface and the flux from the interface to the node j+1. 
!       The equation is solved using the Crank-Nicholson method and 
!       results in a tridiagonal system of equations.
!
! (2) Phase change (see PhaseChange.F90)
!
! 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: Biogeophysics2.F90,v 1.2.10.2 2002/04/27 15:38:36 erik Exp $
!-----------------------------------------------------------------------

  use precision
  use clmtype
  use clm_varcon, only : hvap, cpair, grav, vkc, tfrz, sb
  use clm_varpar, only : nlevsoi
  implicit none

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

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

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

  integer j                          ! do loop index
  real(r8) fact(clm%snl+1 : nlevsoi) ! used in computing tridiagonal matrix
  real(r8) egsmax ! max. evaporation which soil can provide at one time step [mm/s]
  real(r8) egidif ! the excess of evaporation over "egsmax" [mm/s]
  real(r8) xmf    ! total latent heat of phase change of ground water [W/m2]
  real(r8) tinc   ! temperature difference between two time steps [K]

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

!
! Determine soil temperatures including surface soil temperature
!

  call SoilTemperature(clm      , clm%tssbef, clm%htvp, clm%emg, clm%cgrnd, &
                       clm%dlrad, clm%tg    , xmf     , fact )

!
! Correct fluxes to present soil temperature
!

  tinc = clm%t_soisno(clm%snl+1) - clm%tssbef(clm%snl+1)
  clm%eflx_sh_grnd =  clm%eflx_sh_grnd + tinc*clm%cgrnds 
  clm%qflx_evap_soi =  clm%qflx_evap_soi + tinc*clm%cgrndl

!
! egidif holds the excess energy if all water is evaporated from
! the top soil layer during the timestep.  This energy is added to
! the sensible heat flux.
!

  egsmax = (clm%h2osoi_ice(clm%snl+1)+clm%h2osoi_liq(clm%snl+1)) / clm%dtime
  egidif = max( 0._r8, clm%qflx_evap_soi - egsmax )
  clm%qflx_evap_soi = min ( clm%qflx_evap_soi, egsmax )
  clm%eflx_sh_grnd = clm%eflx_sh_grnd + clm%htvp*egidif

!
! Ground heat flux
!

  clm%eflx_soil_grnd = clm%sabg + clm%dlrad + (1-clm%frac_veg_nosno)*clm%emg*clm%forc_lwrad &
       - clm%emg*sb*clm%tssbef(clm%snl+1)**3*(clm%tssbef(clm%snl+1) + 4.*tinc) &
       - (clm%eflx_sh_grnd+clm%qflx_evap_soi*clm%htvp)

!
! Total fluxes (vegetation + ground)
!

  clm%eflx_sh_tot = clm%eflx_sh_veg + clm%eflx_sh_grnd
  clm%qflx_evap_tot = clm%qflx_evap_veg + clm%qflx_evap_soi
  clm%eflx_lh_tot= hvap*clm%qflx_evap_veg + clm%htvp*clm%qflx_evap_soi   ! (account for sublimation)


!
! Assign ground evaporation to sublimation from soil ice or to dew
! on snow or ground 
!

  clm%qflx_evap_grnd = 0.
  clm%qflx_sub_snow = 0.
  clm%qflx_dew_snow = 0.
  clm%qflx_dew_grnd = 0.

  if (clm%qflx_evap_soi >= 0.) then
     ! Do not allow for sublimation in melting (melting ==> evap. ==> sublimation)
     clm%qflx_evap_grnd = min(clm%h2osoi_liq(clm%snl+1)/clm%dtime, clm%qflx_evap_soi)
     clm%qflx_sub_snow = clm%qflx_evap_soi - clm%qflx_evap_grnd
  else
     if (clm%tg < tfrz) then
        clm%qflx_dew_snow = abs(clm%qflx_evap_soi)
     else
        clm%qflx_dew_grnd = abs(clm%qflx_evap_soi)
     endif
  endif

!
! Outgoing long-wave radiation from vegetation + ground
!

  clm%eflx_lwrad_out = clm%ulrad &
       + (1-clm%frac_veg_nosno)*(1.-clm%emg)*clm%forc_lwrad &
       + (1-clm%frac_veg_nosno)*clm%emg*sb * clm%tssbef(clm%snl+1)**4 &
       ! For conservation we put the increase of ground longwave to outgoing
  + 4.*clm%emg*sb*clm%tssbef(clm%snl+1)**3*tinc

!
! Radiative temperature
!

  clm%t_rad = (clm%eflx_lwrad_out/sb)**0.25

!
! Soil energy balance check
!

  clm%errsoi = 0. 
  do j = clm%snl+1, nlevsoi
     clm%errsoi = clm%errsoi - (clm%t_soisno(j)-clm%tssbef(j))/fact(j) 
  enddo
  clm%errsoi = clm%errsoi + clm%eflx_soil_grnd - xmf

!
! Variables needed by history tape
!

 clm%dt_grnd        = tinc
 clm%eflx_lh_vege   = (clm%qflx_evap_veg - clm%qflx_tran_veg) * hvap
 clm%eflx_lh_vegt   = clm%qflx_tran_veg * hvap       
 clm%eflx_lh_grnd   = clm%qflx_evap_soi * clm%htvp
 clm%eflx_lwrad_net = clm%eflx_lwrad_out -  clm%forc_lwrad  

end subroutine Biogeophysics2