leaftem.f90

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

         else
            parsun = 0.
            parsha = (solisb+solisd)         *laifra/max(laisha,mpe)
         endif

! define growing season

         if(tlef > tfrz)then
            igs = 1.
         else
            igs = 0.
         end if

! soil water transpiration factor, it is from the BATS formulation
         btran = etrc / 2.e-4

! maximum possible transpiration rate
         etrc   = sigf * etrc

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

         taf = (tg + thm)/2.
         qaf = (qm+qg)/2.

         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 (itlef <= itmax) 

         tlbef = tlef
         del2 = del
 
! Evaluate stability-dependent variables using moz from prior iteration
 
        call obult ( hu,ht,hq,                &
                    displa,z0mv,z0hv,z0qv,&
                    obu,um,ustar,temp1,temp2, temp12m )
 
! Aerodynamic resistance
 
        ram(1)=1./(ustar*ustar/um)
        rah(1)=1./(temp1*ustar) 
        raw(1)=1./(temp2*ustar) 
 
! Bulk boundary layer resistance of leaves
 
        uaf = um*sqrt( 1./(ram(1)*um) )
        cf = 0.01*sqrtdi/sqrt(uaf)
        rb = 1./(cf*uaf)

! Aerodynamic resistances raw and rah between heights zpd+z0h and z0hg.
! if no vegetation, rah(2)=0 because zpd+z0h = z0hg.
! (Dickinson et al., 1993, pp.54)

        ram(2) = 0.                  ! not used
        rah(2) = 1./(csoilc*uaf)
        raw(2) = rah(2) 
  
! Heat conductance for air, leaf and ground
  
        call condch(sigf,lsai,rah(1),rb,rah(2), &
                    wta,wtl,wtg,wta0,wtl0,wtg0,wtal,wtga,wtgl)

! -----------------------------------------------------------------
! stomatal resistances for sunlit and shaded fractions of canopy.
! should do each iteration to account for differences in eah, tv.
! -----------------------------------------------------------------

        if(lai > 0.001) then
           svpts = el               ! pa
           eah = psrf * qaf / 0.622 ! = q*p/(0.622+0.378q) (pa)
           foln = folnvt

           call stomataB (tfrz ,mpe    ,parsun, &
                          tlef ,svpts  ,eah  ,thm    ,psrf  , &
                          o2   ,co2    ,igs  ,btran  ,rb    , &
                          rssun,psnsun ,qe25 ,aqe    ,kc25  , &
                          ko25 ,vcmx25 ,akc  ,ako    ,avcmx , &
                          bp   ,mp     ,foln ,folnmx ,c3psn )

           call stomataB (tfrz ,mpe    ,parsha, &
                          tlef ,svpts  ,eah  ,thm    ,psrf  , &
                          o2   ,co2    ,igs  ,btran  ,rb    , &
                          rssha,psnsha ,qe25 ,aqe    ,kc25  , &
                          ko25 ,vcmx25 ,akc  ,ako    ,avcmx , &
                          bp   ,mp     ,foln ,folnmx ,c3psn )
 
! Fraction of potential evaporation from leaf
 
             rppdry = fdry*rb*(laisun/(rb+rssun)+laisha/(rb+rssha))/lai

        else
           rssun = 2.e4; rssha=2.e4; psnsun = 0.; psnsha = 0.; rppdry = 0.
        endif
! -----------------------------------------------------------------

        efpot = rhoair*wtl*(qsatl-qaf)

        if(efpot > 0. .AND. etrc > 1.e-12) then
           etr = efpot*rppdry
           rpp = rppdry + fwet
 
! ratd is ratio of transpiration (from dry leaves) to potential
! evaporation for given soil conditions. If transpiration rate exceeds the
! soil potential evaporation (ratd >1) , adjust fraction of dry canopy
! evaporation and stomatal resistance. (Note. Evaporation from wet leaves
! remains unaffected by these changes)
 
           ratd = etr/etrc   ! max transp for the given soil moisture [mm/s]
           if (ratd >= 1) then
              rppdry = rppdry/ratd
              etr = etr/ratd
              rpp = rppdry+fwet
           end if
           epss = 1.e-10
 
! Check total evapotranspiration from leaves
 
           rpp = min(rpp,(etr+ldew/dtime)/efpot - epss)
 
! No transpiration, if potential evaporation from foliage is zero
 
        else
           rpp = 1.
           etr = 0.
        end if
 
! Update conductances for changes in rpp 
! Latent heat conductances for ground and leaf
! Air has same conductance for both sensible and latent heat
 
        call condcq(sigf,lsai,raw(1),rb,raw(2),rpp, &
                    wtaq,wtlq,wtgq,wtaq0,wtlq0,wtgq0,wtalq,wtgaq,wtglq)
 
! the partial derivatives of areodynamical resistance are ignored 
! which cannot be determined analtically
 
        dc1 = rhoair*cp*wtl
        dc2 = dlv*rhoair*wtlq

        efsh = dc1*(wtga*tlef-wtg0*tg-wta0*thm)
        efe = dc2*(wtgaq*qsatl-wtgq0*qg-wtaq0*qm)
 
! "efe" was limited as its sign changes frequently, this limit may
! result in imbalance in "dlv*fevpl" and "efe + dc2*wtgaq*qsatldT*dtlef" 
        if (efe*efeb < 0.0)then
!***        erre = 0.9*efe 
            efe = 0.1*efe
        endif

        dtlef = (sabv + air + bir*tlef**4 + cir*tg**4 - efsh - efe) &
              / (- 4.*bir*tlef**3 +dc1*wtga +dc2*wtgaq*qsatldT)

        tlef = tlbef + dtlef

        dels = tlef-tlbef
        del  = abs(dels)
        err = 0.

! check magnitude of change in leaf temperature
! limit to maximum allowed value

!* 11/29/00
        if(del > delmax)then
           dtlef = delmax*dels/del
           tlef = tlbef + dtlef
           err = sabv + air + bir*tlbef**3*(tlbef + 4.*dtlef) &
               + cir*tg**4 - (efsh + dc1*wtga*dtlef)          &
               - (efe + dc2*wtgaq*qsatldT*dtlef)
        endif

        efpot = rhoair*wtl*(wtgaq*(qsatl+qsatldT*dtlef) &
                           -wtgq0*qg-wtaq0*qm)
        fevpl = rpp*efpot

! calculation of evaporative potentials (efpot) and interception losses; 
! flux in kg m**-2 s-1.  ecidif holds the excess energy if all intercepted water 
! is evaporated during the timestep.  this energies is later added to the
! sensible heat flux.

        ecidif = 0.
        if(efpot > 0. .AND. etrc > 1.e-12) then
           etr = efpot*rppdry
        else
           etr = 0.
        endif
        ecidif = max(0., fevpl-etr-ldew/dtime)
        fevpl = min(fevpl,etr+ldew/dtime)

! the energy loss due to above three limited are added to 
! the snesible heat flux

!***    fsenl = efsh + dc1*wtga*dtlef + err + erre + dlv*ecidif
        fsenl = efsh + dc1*wtga*dtlef + err + dlv*ecidif

! Recalculate leaf saturated vapor pressure (eg) for updated leaf temperature
! and adjust specific humidity (qsatl) proportionately
 
        call qsadv(tlef,psrf,el,deldT,qsatl,qsatldT)
 
! update vegetation/ground surface temperature, canopy air temperature, 
! canopy vapor pressure, aerodynamic temperature, and
! monin-obukhov stability parameter moz for next iteration
 
        taf = wtg0*tg + wta0*thm + wtl0*tlef

        qaf = wtlq0*qsatl+wtgq0*qg+qm*wtaq0

! Update monin-obukhov length and wind speed including the stability effect
        dth = thm-taf       
        dqh = qm-qaf

        tstar=temp1*dth
        qstar=temp2*dqh

        thvstar=tstar+0.61*th*qstar
        zeta=zldis*vonkar*grv*thvstar/(ustar**2*thv)
        if(zeta >= 0.)then      !stable
           zeta = min(2.,max(zeta,1.e-6))
        else                    !unstable
           zeta = max(-100.,min(zeta,-1.e-6))
        endif
        obu = zldis/zeta

        if (zeta >= 0.) then
          um=max(ur,0.1)
        else
          wc=beta*(-grv*ustar*thvstar*zii/thv)**0.333
          um=sqrt(ur*ur+wc*wc)
        endif

        if(obuold*obu < 0.) nmozsgn = nmozsgn+1
        if(nmozsgn >= 4)then 
           obu = zldis/(-0.01)
        end if

        obuold = obu
 
! Test for convergence
 
      itlef = itlef+1
      if(itlef > itmin) then
         dele = abs(efe-efeb)
         efeb = efe
         det  = max(del,del2)
         if(det < dtmin .AND. dele < dlemin) exit 
      endif
 
! Repeat iteration
 
      end do ITERATION
!                   +--------------------------+
!>------------------| END stability iteration: |-----------------------<
!                   +--------------------------+

      rst = 1./(laisun/rssun+laisha/rssha)
      psnsun = psnsun*laisun
      psnsha = psnsha*laisha

! balance check
      err = sabv + air + bir*tlbef**3*(tlbef + 4.*dtlef) &
          + cir*tg**4 - fsenl - dlv*fevpl
      if(abs(err) > 0.2) then
         write(6,*) 'energy balance in canopy X',err
      endif

! fluxes from ground to canopy space
 
      delt  = wtal*tg-wtl0*tlef-wta0*thm
      delq  = wtalq*qg-wtlq0*qsatl-wtaq0*qm

      taux  = taux - sigf*rhoair*us/ram(1)
      tauy  = tauy - sigf*rhoair*vs/ram(1)
      fseng = fseng + cp*rhoair*wtg*delt
      fevpg = fevpg +   rhoair*wtgq*delq

! 2 m height air temperature
!     tref   = tref + sigf*(taf + temp1*dth * 1./vonkar *alog((2.+z0hv)/z0hv))
      tref   = tref + sigf*(thm + temp1*dth * (1./temp12m - 1./temp1))

! downward longwave radiation below the canopy
    
      dlrad = sigf*(1.-emv)*emg*frl &
            + sigf * emv*emg * stefnc * tlbef**3*(tlbef + 4.*dtlef)

! upward longwave radiation above the canopy
    
      ulrad = sigf* (  &
                     (1.-emg)*(1.-emv)*(1.-emv)*frl &
           + emv*(1.+(1.-emg)*(1.-emv))*stefnc * tlbef**3 &
                                               *(tlbef + 4.*dtlef) &
                       + emg *(1.-emv) *stefnc * tg**4 )

! Derivative of soil energy flux with respect to soil temperature (cgrnd)
 
      cgrnds = cgrnds + cp*rhoair*wtg*wtal
      cgrndl = cgrndl + rhoair*wtgq*wtalq*dqgdT
      cgrnd  = cgrnds + cgrndl*htvp
 
! Update dew accumulation (kg/m2)
 
      ldew = max(0.,ldew + (etr-fevpl)*dtime)

  END SUBROUTINE leaftem