leaftem.f90

CLM集合卡曼滤波数据同化算法

  SUBROUTINE leaftem ( dtime  ,sigf  ,z0mv   ,z0hv   ,z0qv   ,&
                       lai    ,sai   ,displa ,sqrtdi ,&
                       hu     ,ht    ,hq     ,us     ,vs     ,&
                       qm     ,thm   ,th     ,thv    ,psrf   ,&
                       rhoair ,cosz  ,sabv   ,solisb ,solisd ,&
                       frl    ,tg    ,tlef   ,ldew   ,qg     ,&
                       dqgdT  ,htvp  ,&
                       emg    ,etrc  ,fsenl  ,fevpl  ,etr    ,&
                       fseng  ,fevpg ,taux   ,tauy   ,tref   ,&
                       dlrad  ,ulrad ,cgrnds ,cgrndl ,cgrnd  ,rst, psnsun, psnsha )

!=======================================================================
!      Source file: leaf_tem.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! 1, calculation of leaf temperature:
!     Use newton-raphson iteration to solve for foliage temperature that balance
!     the surface energy budgets:
! 
!     f(tlef) = Net radiation - Sensible - Latent = 0
!     f(tlef) + d(f)/d(tlef) * dtlef = 0     (*)
! 
! notes:
!     (1) in solving for tlef, tg is given from the previous time-step
!     (2) the partial derivatives of aerodynamical resistances, which can not 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) two energy loss due to: numerical truncation of energy budget equation (*);
!          and "ecidif" (see the code) are drop into sensible heat
!     (6) the convergence criteria: the difference, del = tlef(n+1)-tlef(n) and
!         del2 = tlef(n)-tlef(n-1) less than 0.01 K, and the difference of water flux
!         from leaf between iteration step (n+1) and (n) less than 0.1 W/m2 ; or the
!         iterative steps over 40.
! 
! 2, calculation of leaf fluxes, transpiration, photosynthesis and
!    update the dew accumulation due to evaporation.
!
!=======================================================================

  USE PHYCON_MODULE ! physical constant
  IMPLICIT NONE
!
!----------------------- input/output variables ------------------------

  real, INTENT(in) :: &
        dtime        ! time step [second]

  real, INTENT(in) :: &
        htvp         ! latent heat of evaporation (/sublimation) [J/kg]

! vegetation parameters
  real, INTENT(in) :: &
        lai,        &! adjusted leaf area index for seasonal variation [-]
        sai,        &! stem area index  [-]
        displa,     &! displacement height [m]
        sqrtdi,     &! inverse sqrt of leaf dimension [m**-0.5]
        z0mv,       &! roughness length, momentum [m]
        z0hv,       &! roughness length, sensible heat [m]
        z0qv         ! roughness length, latent heat [m]

! input variables
  real, INTENT(in) :: &
        hu,         &! observational height of wind [m]
        ht,         &! observational height of temperature [m]
        hq,         &! observational height of humidity [m]
        us,         &! wind component in eastward direction [m/s]
        vs,         &! wind component in northward direction [m/s]
   
        thm,        &! intermediate variable (tm+0.0098*ht)
        th,         &! potential temperature (kelvin)
        thv,        &! virtual potential temperature (kelvin)
        qm,         &! specific humidity at reference height [kg/kg]
        psrf,       &! pressure at reference height [pa]
        rhoair,     &! density air [kg/m**3]

        cosz,       &! cosine of the solar zenith angle
        sabv,       &! solar radiation absorbed by vegetation [W/m2]
        solisb,     &!
        solisd,     &!
        frl          ! atmospheric infrared (longwave) radiation [W/m2]

  real, INTENT(in) :: &
        sigf,       &! fraction of veg cover, excluding snow-covered veg [-]
        tg,         &! ground surface temperature [K]
        qg,         &! specific humidity at ground surface [kg/kg]
        dqgdT,      &! temperature derivative of "qg"
        emg          ! vegetation emissivity

  real, INTENT(inout) :: &
        tlef,       &! leaf temperature [K]
        ldew,       &! depth of water on foliage [mm]
        etrc,       &! maximum possible transpiration rate (mm/s)
        taux,       &! wind stress: E-W [kg/m/s**2]
        tauy,       &! wind stress: N-S [kg/m/s**2]
        fseng,      &! sensible heat flux from ground [W/m2]
        fevpg,      &! evaporation heat flux from ground [mm/s]
        cgrnd,      &! deriv. of soil energy flux wrt to soil temp [w/m2/k]
        cgrndl,     &! deriv, of soil sensible heat flux wrt soil temp [w/m2/k]
        cgrnds,     &! deriv of soil latent heat flux wrt soil temp [w/m**2/k]
        tref         ! 2 m height air temperature (kelvin)

  real, INTENT(out) :: &
        fsenl,      &! sensible heat from leaves [W/m2]
        fevpl,      &! evaporation+transpiration from leaves [mm/s]
        etr,        &! transpiration rate [mm/s]
        dlrad,      &! downward longwave radiation blow the canopy [W/m2]
        ulrad,      &! upward longwave radiation above the canopy [W/m2]

        rst,        &!
        psnsun,     &!sunlit leaf photosynthesis [umol co2 /m**2/ s] [+]
        psnsha       !shaded leaf photosynthesis [umol co2 /m**2/ s] [+]
!
!------------------------ local variables ------------------------------
!
  integer           &!
       iter          ! iteration index

  real zldis,       &! reference height "minus" zero displacement heght [m]
       zii,         &! convective boundary layer height [m]
       lsai,        &! leaf + stem area index  [-]
       emv,         &! ground emissivity
       zeta,        &! dimensionless height used in Monin-Obukhov theory
       beta,        &! coefficient of conective velocity [-]
       wc,          &! convective velocity [m/s]
       dth,         &! diff of virtual temp. between ref. height and surface 
       dthv,        &! diff of vir. poten. temp. between ref. height and surface
       dqh,         &! diff of humidity between ref. height and surface
       rib,         &! bulk Richardson number
       obu,         &! monin-obukhov length (m)
       um,          &! wind speed including the stablity effect [m/s]
       ur,          &! wind speed at reference height [m/s]
       uaf,         &! velocity of air within foliage [m/s]
       temp1,       &! relation for potential temperature profile
       temp2,       &! relation for specific humidity profile
       temp12m,      &! relation for specific humidity profile
       ustar,       &! friction velocity [m/s]
       tstar,       &! temperature scaling parameter
       qstar,       &! moisture scaling parameter
       thvstar,     &! virtual potential temperature scaling parameter
       taf,         &! air temperature within canopy space [K]
       qaf           ! humidity of canopy air [kg/kg]

  real fdry,        &! fraction of foliage that is green and dry [-]
       fwet,        &! fraction of foliage covered by water [-]
       rpp,         &! fraction of potential evaporation from leaf [-]
       rppq,        &! fraction of potential evaporation from leaf [-]
       rppdry,      &! fraction of potential evaporation through transp [-]

       resp,        &! leaf respiration from stomat.F
       cf,          &! heat transfer coefficient from leaves [-]

       rb,          &! leaf boundary layer resistance [s/m]
       ram(2),      &! aerodynamical resistance [s/m]
       rah(2),      &! thermal resistance [s/m]
       raw(2),      &! moisture resistance [s/m]

       wta,         &! heat conduactance for air [m/s]
       wtg,         &! heat conduactance for ground [m/s]
       wtl,         &! heat conduactance for leaf [m/s]
       wta0,        &! normalized heat conduactance for air [-]
       wtl0,        &! normalized heat conduactance for leaf [-]
       wtg0,        &! normalized heat conduactance for ground [-]
       wtal,        &! normalized heat conductance for air and leaf [-]
       wtgl,        &! normalized heat conductance for leaf and ground [-]
       wtga,        &! normalized heat cond. for air and ground  [-]

       wtaq,        &! latent heat conduactance for air [m/s]
       wtlq,        &! latent heat conduactance for leaf [m/s]
       wtgq,        &! latent heat conduactance for ground [m/s]
       wtaq0,       &! normalized latent heat conduactance for air [-]
       wtlq0,       &! normalized latent heat conduactance for leaf [-]
       wtgq0,       &! normalized heat conduactance for ground [-]
       wtalq,       &! normalized latent heat cond. for air and leaf [-]
       wtglq,       &! normalized latent heat cond. for leaf and ground [-]
       wtgaq,       &! normalized latent heat cond. for air and ground [-]

       el,          &! vapor pressure on leaf surface [pa]
       deldT,       &! derivative of "el" on "tlef" [pa/K]
       qsatl,       &! leaf specific humidity [kg/kg]
       qsatldT,     &! derivative of "qsatl" on "tlef"

       dtlef,       &! difference of tlef between two iterative step
       air,bir,cir, &! atmos. radiation temporay set
       dc1,dc2       ! derivative of energy flux [W/m2/K]

  real delt,        &! temporary
       delq          ! temporary

  integer           &!
       itlef,       &! counter for leaf temperature iteration [-]
       itmax,       &! maximum number of iteration [-]
       itmin         ! minimum number of iteration [-]

  real del,         &! absolute change in leaf temp in current iteration [K]
       del2,        &! change in leaf temperature in previous iteration [K]
       dele,        &! change in latent heat flux from leaf [K]
       delmax,      &! maximum change in  leaf temperature [K]
       dels,        &! change in leaf temperature in current iteration [K]
       det,         &! maximum leaf temp. change in two consecutive iter [K]
       dlemin,      &! max limit for energy flux convergence [w/m2]
       dtmin,       &! max limit for temperature convergence [K]
       efeb          ! latent heat flux from leaf (previous iter) [mm/s]

  real efpot,       &! potential latent energy flux [kg/m2/s]
       efe,         &! water flux from leaf [mm/s]
       efsh,        &! sensible heat from leaf [mm/s]
       epss,        &! minimum canopy weyness factor [-]
       ratd          ! ratio of transpiration to potential evaporation [-]

  integer nmozsgn    ! number of times moz changes sign

  real obuold,      &! monin-obukhov length from previous iteration
       tlbef,       &! leaf temperature from previous iteration [K]
       ecidif,      &! excess energies [W/m2]
       erre,        &! balance error
       err           ! balance error

! for use Bonan's stomatal resistance
!-----------------------------------------------------------------------
  real fsun          !sunlit fraction of canopy
  real fsha          !shaded fraction of canopy
  real laisun        !sunlit leaf area index, one-sided
  real laisha        !shaded leaf area index, one-sided
  real laifra        !leaf area fraction of canopy
  real btran         !soil water transpiration factor (0 - 1)
  real igs           !growing season index (0=off, 1=on)

! constant atmospheric co2 and o2
  real po2           !partial pressure  o2 (mol/mol)
  real pco2          !partial pressure co2 (mol/mol)

  data po2,pco2 /0.209,355.e-06/

  real co2           !atmospheric co2 concentration (pa)
  real o2            !atmospheric o2 concentration (pa)

  real svpts         !saturation vapor pressure at tlef (pa)
  real eah           !canopy air vapor pressure (pa)
  real foln          !foliage nitrogen (%)
  real rssun         !sunlit leaf stomatal resistance [s/m]
  real rssha         !shaded leaf stomatal resistance [s/m]

  real gdir          !relative projected leaf+stem area in solar direction
  real ext           !optical depth of direct beam per unit leaf area
  real mpe           !prevents overflow error if division by zero

  real parsun        ! visible direct incident flux
  real parsha        ! visible diffuse incident flux
  real solisba       !
  real solisda       !

!-----------------------------------------------------------------------
  real, parameter :: qe25 = 0.08   !quantum efficiency at 25c (umol co2 / umol photon)
  real, parameter :: ko25 = 30000. !o2 michaelis-menten constant at 25c (pa)
  real, parameter :: kc25 = 30.    !co2 michaelis-menten constant at 25c (pa)
  real, parameter :: vcmx25 = 29.3 !maximum rate of carboxylation at 25c (umol co2/m**2/s)
  real, parameter :: aqe = 1.      !q10 for qe25
  real, parameter :: ako = 1.2     !q10 for ko25
  real, parameter :: akc = 2.1     !q10 for kc25
  real, parameter :: avcmx = 2.4   !q10 for vcmx25
  real, parameter :: bp = 10000.   !minimum leaf conductance (umol/m**2/s)
  real, parameter :: mp = 9.       !slope for conductance-to-photosynthesis relationship
  real, parameter :: folnmx = 1.5  !foliage nitrogen concentration when f(n)=1 (%)
  real, parameter :: folnvt = 2.   !foliage nitrogen concentration (%)
  real, parameter :: c3psn = 1.    !photosynthetic pathway: 0. = c4, 1. = c3
!-----------------------------------------------------------------------
! initialization
         del    = 0.0  ! change in leaf temperature from previous iteration
         itlef  = 0    ! counter for leaf temperature iteration
         efeb   = 0.0  ! latent head flux from leaf for previous iteration
 
! assign iteration parameters
         delmax = 1.0  ! maximum change in  leaf temperature
         itmax  = 40   ! maximum number of iteration
         itmin  = 2    ! minimum number of iteration
         dtmin  = 0.01 ! max limit for temperature convergence
         dlemin = 0.1  ! max limit for energy flux convergence

! net absorbed longwave radiation by canopy and ground
! =air+bir*tlef**4+cir*tg**4
         emv = 1.
         air =   sigf * emv * (1.+(1.-emv)*(1.-emg)) * frl
         bir = - sigf * (2.-emv*(1.-emg)) * emv * stefnc
         cir =   sigf * emv*emg*stefnc

! get fraction of wet and dry canopy surface (fwet & fdry)
         lsai = lai + sai
         call fractdew (sigf,lai,sai,dewmx,ldew,fwet,fdry)

! saturated vapor pressure and humidity and their derivation
         call qsadv(tlef,psrf,el,deldT,qsatl,qsatldT)

! for use Bonan's stomatal resistance scheme
! atmospheric co2 and o2 are currently constants

         co2 = pco2*psrf
         o2  = po2*psrf

         mpe = 1.e-6
         gdir = 0.5
         ext = gdir/max(0.001,cosz)*sqrt(0.85)
         fsun = (1.-exp(-ext*lai))/max(ext*lai,mpe)
         if(fsun < 0.01) fsun = 0.  

! partion visible canopy absorption to sunlit and shaded fractions
! to get average absorbed par for sunlit and shaded leaves

         fsha = 1.-fsun
         laisun = lai*fsun
         laisha = lai*fsha
         laifra = 1.                             ! = lai / max(lsai,mpe)
         solisba = solisb*(1.-exp(-ext*lsai))    !define direct PAR absorbed
         solisda = solisd*(1.-exp(-sqrt(0.85)*lsai))   !diffuse PAR absorbed
         if(fsun > 0.)then                                 ! PAR flux on leaves
            parsun = (solisba + fsun*solisda)*laifra/max(laisun,mpe)
            parsha =            fsha*solisda *laifra/max(laisha,mpe)