stomatab.f90

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

  subroutine stomataB (tfrz ,mpe    ,apar , &
                       tv   ,ei     ,ea   ,tgcm   ,pgcm  , &
                       o2   ,co2    ,igs  ,btran  ,rb    , &
                       rs   ,psn    ,qe25 ,aqe    ,kc25  , &
                       ko25 ,vcmx25 ,akc  ,ako    ,avcmx , &
                       bp   ,mp     ,foln ,folnmx ,c3psn )

   IMPLICIT NONE
! ------------------------ code history ---------------------------
! source file:       stomata.f90
! purpose:           leaf stomatal resistance and leaf photosynthesis
! date last revised: March 1996 - lsm version 1
! author:            Gordon Bonan

! ------------------------ input/output variables -----------------
  real, INTENT(in) :: &
        tfrz,       &!freezing point (kelvin)
        mpe,        &!prevents division by zero errors
        tv,         &!foliage temperature (kelvin)
        ei,         &!vapor pressure inside leaf (sat vapor press at tv) (pa)
        ea,         &!vapor pressure of canopy air (pa)
        apar,       &!par absorbed per unit lai (w/m**2)
        o2,         &!atmospheric o2 concentration (pa)
        co2,        &!atmospheric co2 concentration (pa)
        pgcm,       &!air pressure at agcm reference height (pa)
        tgcm,       &!air temperature at agcm reference height (kelvin)
        igs,        &!growing season index (0=off, 1=on)
        btran,      &!soil water transpiration factor (0 to 1)
        foln,       &!foliage nitrogen concentration (%)

        qe25,       &!quantum efficiency at 25c (umol co2 / umol photon)
        ko25,       &!o2 michaelis-menten constant at 25c (pa)
        kc25,       &!co2 michaelis-menten constant at 25c (pa)
        vcmx25,     &!maximum rate of carboxylation at 25c (umol co2/m**2/s)
        aqe,        &!q10 for qe25
        ako,        &!q10 for ko25
        akc,        &!q10 for kc25
        avcmx,      &!q10 for vcmx25
        bp,         &!minimum leaf conductance (umol/m**2/s)
        mp,         &!slope for conductance-to-photosynthesis relationship 
        folnmx,     &!foliage nitrogen concentration when f(n)=1 (%)
        c3psn        !photosynthetic pathway: 0. = c4, 1. = c3

  real, INTENT(inout) :: &
        rb           !boundary layer resistance (s/m)

! output
  real, INTENT(out) :: &
        rs,         &!leaf stomatal resistance (s/m)
        psn          !foliage photosynthesis (umol co2 /m**2/ s) [always +]

! ------------------------ local variables ------------------------

  integer, parameter :: niter = 3  ! number of iterations

  integer  iter      !iteration index

  real ab,          &!used in statement functions
       bc,          &!used in statement functions
       f1,          &!generic temperature response (statement function)
       f2,          &!generic temperature inhibition (statement function)
       tc,          &!foliage temperature (degree celsius)
       cs,          &!co2 concentration at leaf surface (pa)
       kc,          &!co2 michaelis-menten constant (pa)
       ko,          &!o2 michaelis-menten constant (pa)
       a,b,c,q,     &!intermediate calculations for rs
       r1,r2,       &!roots for rs
       fnf,         &!foliage nitrogen adjustment factor (0 to 1)
       ppf,         &!absorb photosynthetic photon flux (umol photons/m**2/s)
       wc,          &!rubisco limited photosynthesis (umol co2/m**2/s)
       wj,          &!light limited photosynthesis (umol co2/m**2/s)
       we,          &!export limited photosynthesis (umol co2/m**2/s)
       cp,          &!co2 compensation point (pa)
       ci,          &!internal co2 (pa)
       awc,         &!intermediate calcuation for wc
       vcmx,        &!maximum rate of carboxylation (umol co2/m**2/s)
       j,           &!electron transport (umol co2/m**2/s)
       cea,         &!constrain ea or else model blows up
       cf,          &!s m**2/umol -> s/m
       rsmax0        ! maximum stomatal resistance [s/m]

! -----------------------------------------------------------------

      f1(ab,bc) = ab**((bc-25.)/10.)
      f2(ab) = 1. + exp((-2.2e05+710.*(ab+273.16))/(8.314*(ab+273.16)))

! initialize rs=rsmax and psn=0 because will only do calculations
! for apar > 0, in which case rs <= rsmax and psn >= 0

! set constants
      rsmax0 = 2.e4
      cf = pgcm/(8.314*tgcm)*1.e06 

      if(apar <= 0.) then           ! night time
!        rs = min(rsmax0, 1./bp * cf)
         rs = 1./(bp*btran) * cf
         psn = 0.

      else                          ! day time
         fnf = min(foln/max(mpe,folnmx), 1.0)
         tc = tv-tfrz                            
         ppf = 4.6*apar                  
         j = ppf*qe25
         kc = kc25 * f1(akc,tc)       
         ko = ko25 * f1(ako,tc)
         awc = kc * (1.+o2/ko)
         cp = 0.5*kc/ko*o2*0.21
         vcmx = vcmx25 * f1(avcmx,tc) / f2(tc) * fnf * btran

! first guess ci

         ci = 0.7*co2*c3psn + 0.4*co2*(1.-c3psn)  

! rb: s/m -> s m**2 / umol

         rb = rb/cf 

! constrain ea

         cea = max(0.25*ei*c3psn+0.40*ei*(1.-c3psn), min(ea,ei) ) 

! ci iteration

         do iter = 1, niter
            wj = max(ci-cp,0.)*j/(ci+2.*cp)*c3psn + j*(1.-c3psn)
            wc = max(ci-cp,0.)*vcmx/(ci+awc)*c3psn + vcmx*(1.-c3psn)
            we = 0.5*vcmx*c3psn + 4000.*vcmx*ci/pgcm*(1.-c3psn) 
            psn = min(wj,wc,we) * igs

            cs = max( co2-1.37*rb*pgcm*psn, mpe )
            a = mp*psn*pgcm*cea / (cs*ei) + bp*btran
            b = ( mp*psn*pgcm/cs + bp*btran ) * rb - 1.
            c = -rb
            if (b >= 0.) then
               q = -0.5*( b + sqrt(b*b-4.*a*c) )
            else
               q = -0.5*( b - sqrt(b*b-4.*a*c) )
            end if
            r1 = q/a
            r2 = c/q
            rs = max(r1,r2)
            ci = max( cs-psn*pgcm*1.65*rs, 0. )
         end do

! rs, rb:  s m**2 / umol -> s/m 

!        rs = min(rsmax0, rs*cf)
         rs = rs*cf
         rb = rb*cf 

      endif

  END SUBROUTINE stomataB