twoleaf.f90

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

        assimsha,   &! shaded leaf assimilation rate [umol co2 /m**2/ s] [+]
        respcsun,   &!
        respcsha,   &!
        rsoil,      &!
        gah2o,      &!
        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
        ndtlsgn1,   &! number of times dtlsun changes sign
        ndtlsgn2     ! number of times dtlsun changes sign
 
   real etrsun, etrsha, evplwetsun, evplwetsha
   real etrsun_dtlsun, etrsha_dtlsun, evplwetsun_dtlsun, evplwetsha_dtlsun
   real etrsun_dtlsha, etrsha_dtlsha, evplwetsun_dtlsha, evplwetsha_dtlsha

  real, dimension(:), allocatable :: &
        dtlsun,     &! difference of tlsun between two iterative step
        dtlsha       ! difference of tlsha 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 = 1.0  ! maximum change in  leaf temperature
       dtmin  = 0.01 ! max limit for temperature convergence
       dlemin = 0.1  ! max limit for energy flux convergence

       allocate (dtlsun(0:itmax+1))
       allocate (dtlsha(0:itmax+1))

       dtlsun(0) = 0.
       dtlsha(0) = 0.

       ndtlsgn1 = 0
       ndtlsgn2 = 0
!-----------------------------------------------------------------------
! leaf area index
!-----------------------------------------------------------------------
! 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

!-----------------------------------------------------------------------
! 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(tlsun,psrf,eisun,deisunDT,qsatlsun,qsatlsunDT)
       call qsadv(tlsha,psrf,eisha,deishaDT,qsatlsha,qsatlshaDT)

!-----------------------------------------------------------------------
! 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. !soil 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) 

         tlsunbef = tlsun
         tlshabef = tlsha

         del2 = del
         dele2 = dele
 
!-----------------------------------------------------------------------
! Aerodynamical resistances
!-----------------------------------------------------------------------
! 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./(ustar*ustar/um) 
        rah = 1./(temp1*ustar)
        raw = 1./(temp2*ustar)
 
! Bulk boundary layer resistance of leaves
        uaf = um*sqrt( 1./(ram*um) )
        cf = 0.01*sqrtdi/sqrt(uaf)
        rb = 1./(cf*uaf) 
        rd = 1./(csoilc*uaf) 
!-----------------------------------------------------------------------
! stomatal resistances for sunlit and shaded fractions of canopy.
! should do each iteration to account for differences in eah, leaf temperatures.
!-----------------------------------------------------------------------
        if(lai .gt. 0.001) then
           rbsun = rb / laisun
           rbsha = rb / laisha

           eah = qaf * psrf / ( 0.622 + 0.378 * qaf )    ! pa

! Sunlit leaves
           CALL stomata (   vmax25 ,effcon  ,slti     ,hlti     ,shti   ,&
               hhti  ,trda ,trdm   ,trop    ,gradm    ,binter   ,thm    ,&
               psrf  ,po2m ,pco2m  ,pco2a   ,eah      ,eisun    ,tlsun  ,parsun ,&
               rbsun ,raw  ,rstfac ,cintsun ,assimsun ,respcsun ,rssun   )

! Shaded leaves
           CALL stomata (   vmax25 ,effcon  ,slti     ,hlti     ,shti   ,&
               hhti  ,trda ,trdm   ,trop    ,gradm    ,binter   ,thm    ,&
               psrf  ,po2m ,pco2m  ,pco2a   ,eah      ,eisha    ,tlsha  ,parsha ,&
               rbsha ,raw  ,rstfac ,cintsha ,assimsha ,respcsha ,rssha   )

        else
           rssun = 1.e10 ; rssha = 1.e10
           assimsun = 0. ; assimsha = 0.
           respcsun = 0. ; respcsha = 0.
        endif

! above stomatal resistances are for the canopy, the stomatal resistance 
! the "rb" in the following calculations are the averages for single leaf. thus,
        rssun = rssun * laisun
        rssha = rssha * laisha

!-----------------------------------------------------------------------
! dimensional and non-dimensional sensible and latent heat conductances
! for canopy and soil flux calculations.
!-----------------------------------------------------------------------
        delta1 = 0.0
        delta2 = 0.0
        if( (fwet .lt. .99) .AND. (sigf*etrc .gt. 1.e-12) )then
           if(qsatlsun-qaf .gt. 0.) delta1 = 1.0
           if(qsatlsha-qaf .gt. 0.) delta2 = 1.0
        endif
 
        cah = sigf / rah
        cgh = sigf / rd
        cfsunh = sigf * laisun / rb
        cfshah = sigf * (laisha + sai) / rb

        caw = sigf / raw
        cgw = sigf / rd
        cfsunw = sigf * ( fwet * laisun / rb + (1. - fwet) * delta1 * laisun / (rb + rssun) )
        cfshaw = sigf * ( fwet * (laisha + sai) / rb + (1. - fwet) * delta2 * laisha / (rb + rssha) )

        wtshi = 1. / ( cah + cgh + cfsunh + cfshah )
        wtsqi = 1. / ( caw + cgw + cfsunw + cfshaw )

        wta0 = cah * wtshi
        wtg0 = cgh * wtshi
        wtlsun0 = cfsunh * wtshi
        wtlsha0 = cfshah * wtshi

        wtaq0 = caw * wtsqi
        wtgq0 = cgw * wtsqi
        wtlsunq0 = cfsunw * wtsqi
        wtlshaq0 = cfshaw * wtsqi
 
!-----------------------------------------------------------------------
! IR radiation, sensible and latent heat fluxes and their derivatives
!-----------------------------------------------------------------------
! the partial derivatives of areodynamical resistance are ignored 
! which cannot be determined analtically
        fac = sigf * (1. - thermk)
        fac1 = fac * fsun
        fac2 = fac * fsha

! longwave absorption and their derivatives 
        irabsun = (frl - 2. * stefnc * tlsun**4 + emg*stefnc*tg**4 ) * fac1 
        dirabsun_dtlsun = - 8.* stefnc * tlsun**3                    * fac1
        dirabsun_dtlsha = 0.

        irabsha = (frl - 2. * stefnc * tlsha**4 + emg*stefnc*tg**4 ) * fac2 
        dirabsha_dtlsha = - 8.* stefnc * tlsha**3                    * fac2 
        dirabsha_dtlsun = 0.

! sensible heat fluxes and their derivatives
        senlsun = rhoair * cp * cfsunh &
                * ( (wta0 + wtg0 + wtlsha0)*tlsun &
                   - wta0*thm - wtg0*tg - wtlsha0*tlsha )
        senlsun_dtlsun = rhoair * cp * cfsunh * (wta0 + wtg0 + wtlsha0)
        senlsun_dtlsha = rhoair * cp * cfsunh * ( - wtlsha0 )

        senlsha = rhoair * cp * cfshah &
                * ( (wta0 + wtg0 + wtlsun0)*tlsha &
                   - wta0*thm - wtg0*tg - wtlsun0*tlsun )
        senlsha_dtlsun = rhoair * cp * cfshah * ( - wtlsun0 ) 
        senlsha_dtlsha = rhoair * cp * cfshah * ( wta0 + wtg0 + wtlsun0 )

! latent heat fluxes and their derivatives
!-      evplsun = rhoair * cfsunw * ( (wtaq0 + wtgq0 + wtlshaq0)*qsatlsun - wtaq0*qm - wtgq0*qg - wtlshaq0*qsatlsha )
!-      evplsun_dtlsun = rhoair * cfsunw * ( wtaq0 + wtgq0 + wtlshaq0 ) * qsatlsunDT 
!-      evplsun_dtlsha = rhoair * cfsunw * ( - wtlshaq0 ) * qsatlshaDT 

        etrsun = sigf * rhoair * (1.-fwet) * delta1 * laisun / (rb + rssun) * ( (wtaq0 + wtgq0 + wtlshaq0)*qsatlsun &
               - wtaq0*qm - wtgq0*qg - wtlshaq0*qsatlsha )
        etrsun_dtlsun = sigf * rhoair * (1.-fwet) * delta1 * laisun / (rb + rssun) * (wtaq0 + wtgq0 + wtlshaq0)*qsatlsunDT 
        etrsun_dtlsha = sigf * rhoair * (1.-fwet) * delta1 * laisun / (rb + rssun) * ( - wtlshaq0*qsatlshaDT )

        if(etrsun .ge.  sigf*etrc*laisun/lai)then
           etrsun =  sigf*etrc*laisun/lai
           etrsun_dtlsun = 0.
           etrsun_dtlsha = 0.
        end if

        evplwetsun = sigf * rhoair * fwet * laisun / rb * ( (wtaq0 + wtgq0 + wtlshaq0)*qsatlsun &
                    - wtaq0*qm - wtgq0*qg - wtlshaq0*qsatlsha )
        evplwetsun_dtlsun = sigf * rhoair * fwet * laisun / rb * (wtaq0 + wtgq0 + wtlshaq0)*qsatlsunDT