oneleaf.f90

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

  SUBROUTINE oneleaf ( dtime    ,htvp     ,lai      ,sai      ,displa   , & 
             sqrtdi   ,z0mv     ,z0hv     ,z0qv     ,effcon   ,vmax25   , &
             slti     ,hlti     ,shti     ,hhti     ,trda     ,trdm     , &
             trop     ,gradm    ,binter   , &

             hu       ,ht       ,hq       ,us       ,vs       ,thm      , &
             th       ,thv      ,qm       ,psrf     ,rhoair   ,cosz     , &
             par      ,sabv     ,frl      ,cint     ,thermk   ,rstfac   , &
             po2m     ,pco2m    ,sigf     ,etrc     ,tg       ,qg       , &
             dqgdT    ,emg      , &

             tl       ,ldew     ,taux     ,tauy     ,fseng    ,fevpg    , &
             cgrnd    ,cgrndl   ,cgrnds   ,tref     ,rst      ,assim    ,respc    , &
             fsenl    ,fevpl    ,etr      ,dlrad    ,ulrad    )
 
!=======================================================================
!      Source file: twoleaf.f90
! Original version: Yongjiu Dai, August 15, 2001
!
!=======================================================================

  USE PHYCON_MODULE ! physical constant
  IMPLICIT NONE
 
  real, INTENT(in) :: &
        dtime,      &! time step [second]
        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]

        effcon,     &! quantum efficiency of RuBP regeneration (mol CO2 / mol quanta)
        vmax25,     &! maximum carboxylation rate at 25 C at canopy top
                     ! the range : 30.e-6 <-> 100.e-6 (mol co2 m-2 s-1)
        slti,       &! slope of low temperature inhibition function      (0.2)
        hlti,       &! 1/2 point of low temperature inhibition function  (288.16)
        shti,       &! slope of high temperature inhibition function     (0.3)
        hhti,       &! 1/2 point of high temperature inhibition function (313.16)
        trda,       &! temperature coefficient in gs-a model             (1.3)
        trdm,       &! temperature coefficient in gs-a model             (328.16)
        trop,       &! temperature coefficient in gs-a model             (298.16)
        gradm,      &! conductance-photosynthesis slope parameter
        binter       ! conductance-photosynthesis intercept

! 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
        par,        &! par absorbed per unit lai [w/m**2]
        sabv,       &! solar radiation absorbed by vegetation [W/m2]
        frl,        &! atmospheric infrared (longwave) radiation [W/m2]

        cint(3),    &!
        thermk,     &!
        rstfac,     &!

        po2m,       &! atmospheric partial pressure  o2 (pa)
        pco2m,      &! atmospheric partial pressure co2 (pa)

        sigf,       &! fraction of veg cover, excluding snow-covered veg [-]
        etrc,       &! maximum possible transpiration rate (mm/s)
        tg,         &! ground surface temperature [K]
        qg,         &! specific humidity at ground surface [kg/kg]
        dqgdT,      &! temperature derivative of "qg"
        emg          ! vegetation emissivity

  real, INTENT(inout) :: &
        tl,         &! leaf temperature [K]
        ldew,       &! depth of water on foliage [mm]
        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) :: &
        rst,        &!
        assim,      &!
        respc,      &!
        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]
!
!------------------------ local variables ------------------------------
!
   real zldis,      &! reference height "minus" zero displacement heght [m]
        zii,        &! convective boundary layer height [m]
        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,    &! 
        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]
        eah,        &! canopy air vapor pressure (pa)
        pco2g,      &! co2 pressure (pa) at ground surface (pa)
        pco2a        ! canopy air co2 pressure (pa)

   real fdry,       &! fraction of foliage that is green and dry [-]
        fwet,       &! fraction of foliage covered by water [-]
        cf,         &! heat transfer coefficient from leaves [-]
        rb,         &! leaf boundary layer resistance [s/m]
        rbone,      &! canopy bulk boundary layer resistance 
        rd,         &! aerodynamical resistance between ground and canopy air
        ram,        &! aerodynamical resistance [s/m]
        rah,        &! thermal resistance [s/m]
        raw,        &! moisture resistance [s/m]
        clai,       &! canopy heat capacity [Jm-2K-1]

        cah,        &! heat conduactance for air [m/s]
        cgh,        &! heat conduactance for ground [m/s]
        cfh,        &! heat conduactance for leaf [m/s]

        caw,        &! latent heat conduactance for air [m/s]
        cgw,        &! latent heat conduactance for ground [m/s]
        cfw,        &! latent heat conduactance for leaf [m/s]

        wtshi,      &! sensible heat resistance for air, grd and leaf [-]
        wtsqi,      &! latent heat resistance for air, grd and leaf [-]

        wta0,       &! normalized heat conduactance for air [-]
        wtg0,       &! normalized heat conduactance for ground [-]
        wtl0,       &! normalized heat conductance for air and leaf [-]

        wtaq0,      &! normalized latent heat conduactance for air [-]
        wtgq0,      &! normalized heat conduactance for ground [-]
        wtlq0,      &! normalized latent heat cond. for air and leaf [-]

        ei,         &! vapor pressure on leaf surface [pa]
        deidT,      &! derivative of "ei" on "tl" [pa/K]
        qsatl,      &! leaf specific humidity [kg/kg]
        qsatldT,    &! derivative of "qsatl" on "tlef"
 
        tl,         &! leaf temperature
 
        delta,      &!
        fac,        &!

        irab,       &!
        dirab_dtl,  &!
        fsenl_dtl,  &!
        fevpl_dtl    !
 
   real del,        &! absolute change in leaf temp in current iteration [K]
        del2,       &! change in leaf temperature in previous iteration [K]
        dele,       &! change in heat fluxes from leaf [K]
        dele2,      &! change in heat fluxes from leaf [K]
        det,        &! maximum leaf temp. change in two consecutive iter [K]
        del,        &! maximum heat fluxes change in two consecutive iter [K]
        delmax,     &! maximum change in  leaf temperature [K]
        dtmin,      &! max limit for temperature convergence [K]
        dlemin       ! max limit for energy flux convergence [w/m2]
 
   real obuold,     &! monin-obukhov length from previous iteration
        tlbef,      &! leaf temperature from previous iteration [K]
        ecidif,     &! excess energies [W/m2]
        err          ! balance error

   real rs,         &! leaf stomatal resistance [s/m]
        rsoil,      &!
        gah2o,      &!
        gdh2o,      &!
        tprcor       !

   integer          &!
        it,         &! counter for leaf temperature iteration [-]
        itmax,      &! maximum number of iteration [-]
        itmin,      &! minimum number of iteration [-]
        nmozsgn      ! number of times moz changes sign

   real fevplwet, fevplwet_dtl, etr_dtl

   real, dimension(:), allocatable :: &
        dtl          ! difference of tl between two iterative step
!-----------------------------------------------------------------------
! initialization of errors and  iteration parameters
!-----------------------------------------------------------------------
       it     = 1    ! counter for leaf temperature iteration
       del    = 0.0  ! change in leaf temperature from previous iteration
       dele   = 0.0  ! latent head flux from leaf for previous iteration
 
! assign iteration parameters
       itmax  = 40   ! maximum number of iteration
       itmin  = 3    ! minimum number of iteration
       delmax = 3.0  ! maximum change in  leaf temperature
       dtmin  = 0.01 ! max limit for temperature convergence
       dlemin = 0.1  ! max limit for energy flux convergence

       allocate (dtl(0:itmax+1))
       dtl(0) = 0.

!-----------------------------------------------------------------------
! get fraction of wet and dry canopy surface (fwet & fdry)
! initial saturated vapor pressure and humidity and their derivation
!-----------------------------------------------------------------------
       clai = 4.2 * 1000. * 0.2

       call fractdew (sigf,lai,sai,dewmx,ldew,fwet,fdry)

       call qsadv(tl,psrf,ei,deiDT,qsatl,qsatlDT)

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

       taf = 0.5 * (tg + thm)
       qaf = 0.5 * (qm + qg)

       pco2a = pco2m
       tprcor = 44.6*273.16*psrf/1.013e5
       rsoil = 0.  !respiration (mol m-2 s-1)
!      rsoil = 1.22e-6*exp(308.56*(1./56.02-1./(tg-227.13)))
!      rsoil = rstfac * 0.23 * 15. * 2.**((tg-273.16-10.)/10.) * 1.e-6
!      rsoil = 5.22 * 1.e-6
       rsoil = 0.22 * 1.e-6

       ur = max(0.1, sqrt(us*us+vs*vs))    ! limit set to 0.1
       dth = thm - taf
       dqh = qm - qaf
       dthv = dth*(1.+0.61*qm) + 0.61*th*dqh
       zldis = hu - displa

       call obuini(ur,th,thm,thv,dth,dqh,dthv,zldis,z0mv,um,obu)


! ======================================================================
!      BEGIN stability iteration 
! ======================================================================

      ITERATION : do while (it .le. itmax) 

         tlbef = tl