oneleaf.f90

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

         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 
!-----------------------------------------------------------------------
        if(lai .gt. 0.001) then
           rbone = rb / lai
           eah = qaf * psrf / ( 0.622 + 0.378 * qaf )    ! pa

           CALL stomata (   vmax25 ,effcon ,slti  ,hlti   ,shti ,&
               hhti  ,trda ,trdm   ,trop   ,gradm ,binter ,thm  ,&
               psrf  ,po2m ,pco2m  ,pco2a  ,eah   ,ei     ,tl   ,par  ,&
               rbone ,raw  ,rstfac ,cint   ,assim ,respc  ,rs    )

        else
           rs = 1.e10; assim = 0.; respc = 0.
        endif

! above stomatal resistances are for the canopy, the stomatal rsistances 
! and the "rb" in the following calculations are the average for single leaf. thus,
        rs = rs * lai

!-----------------------------------------------------------------------
! dimensional and non-dimensional sensible and latent heat conductances
! for canopy and soil flux calculations.
!-----------------------------------------------------------------------
        delta = 0.0
        if( (fwet .lt. .99) .AND. (sigf*etrc .gt. 1.e-12) )then
           if(qsatl-qaf .gt. 0.) delta = 1.0
        endif
 
        cah = sigf / rah
        cgh = sigf / rd
        cfh = sigf * (lai + sai) / rb

        caw = sigf / raw
        cgw = sigf / rd
        cfw = sigf * ( fwet*(lai+sai)/rb + (1.-fwet)*delta*lai/(rb+rs) )

        wtshi = 1. / ( cah + cgh + cfh )
        wtsqi = 1. / ( caw + cgw + cfw )

        wta0 = cah * wtshi
        wtg0 = cgh * wtshi
        wtl0 = cfh * wtshi

        wtaq0 = caw * wtsqi
        wtgq0 = cgw * wtsqi
        wtlq0 = cfw * 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)

! longwave absorption and their derivatives 
        irab = (frl - 2. * stefnc * tl**4 + emg*stefnc*tg**4 ) * fac 
        dirab_dtl = - 8.* stefnc * tl**3                       * fac

! sensible heat fluxes and their derivatives
        fsenl = rhoair * cp * cfh * ( (wta0 + wtg0)*tl - wta0*thm - wtg0*tg )
        fsenl_dtl = rhoair * cp * cfh * (wta0 + wtg0)

! latent heat fluxes and their derivatives
!--->   fevpl = rhoair * cfw * ( (wtaq0 + wtgq0)*qsatl - wtaq0*qm - wtgq0*qg )
!--->   fevpl_dtl = rhoair * cfw * (wtaq0 + wtgq0) * qsatlDT 

        etr = sigf * rhoair * (1.-fwet) * delta * lai / (rb + rs) * ( (wtaq0 + wtgq0)*qsatl - wtaq0*qm - wtgq0*qg )
        etr_dtl =  sigf * rhoair * (1.-fwet) * delta * lai / (rb + rs) * (wtaq0 + wtgq0)*qsatlDT 
     
        if(etr.ge.sigf*etrc)then
           etr = sigf*etrc
           etr_dtl = 0.
        endif

        fevplwet = sigf * rhoair * fwet * (lai+sai) / rb * ( (wtaq0 + wtgq0)*qsatl - wtaq0*qm - wtgq0*qg )
        fevplwet_dtl = sigf * rhoair * fwet * (lai+sai) / rb * (wtaq0 + wtgq0)*qsatlDT 
        if(fevplwet.ge.ldew/dtime)then
           fevplwet = ldew/dtime
           fevplwet_dtl = 0.
        endif
 
        fevpl = etr + fevplwet
        fevpl_dtl = etr_dtl + fevplwet_dtl

!-----------------------------------------------------------------------
! difference of temperatures by quasi-newton-raphson method for the non-linear system equations
!-----------------------------------------------------------------------
        dtl(it) = (sabv + irab - fsenl - dlv*fevpl) &
            / ((lai+sai)*clai/dtime - dirab_dtl + fsenl_dtl + dlv*fevpl_dtl)

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

        if(abs(dtl(it)).gt.delmax) dtl(it) = delmax*dtl(it)/abs(dtl(it))
        if((it.ge.2) .and. (dtl(it-1)*dtl(it).lt.0.)) dtl(it) = 0.5*(dtl(it-1) + dtl(it))

        tl = tlbef + dtl(it)

!-----------------------------------------------------------------------
! square roots differences of temperatures and fluxes for use as the condition of convergences
!-----------------------------------------------------------------------
        del  = sqrt( dtl(it)*dtl(it) )
        dele = dtl(it) * dtl(it) * ( dirab_dtl**2 + fsenl_dtl**2 + dlv*fevpl_dtl**2 ) 
        dele = sqrt(dele)

!-----------------------------------------------------------------------
!  saturated vapor pressures and canopy air temperature, canopy air humidity
!-----------------------------------------------------------------------
! Recalculate leaf saturated vapor pressure (ei_)for updated leaf temperature
! and adjust specific humidity (qsatl_) proportionately
        call qsadv(tl,psrf,ei,deiDT,qsatl,qsatlDT)

! update vegetation/ground surface temperature, canopy air temperature, 
! canopy air humidity
        taf = wta0*thm + wtg0*tg + wtl0*tl 

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

! update co2 partial pressure within canopy air
        gah2o  = 1.0/raw * tprcor/thm                     ! mol m-2 s-1
        gdh2o  = 1.0/rd  * tprcor/thm                     ! mol m-2 s-1
        pco2a = pco2m - 1.37*psrf/max(0.446,gah2o) * (assim - respc - rsoil)

!     if( pco2a < 30.) print*, pco2a, pco2m, assim-respc-rsoil, assim-respc, gah2o
!       
!       pco2g =  pco2m
!       if(sabv > 0.) then
!          pco2g = 2.* pco2m
!          gdh2o = gah2o
!       endif
!       pco2a = ( pco2m*gah2o + pco2g*gdh2o - 1.37*psrf*(assim-respc) ) / (gah2o + gdh2o)

!-----------------------------------------------------------------------
! 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 .ge. 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 .ge. 0.)then
          um = max(ur,.1)
        else
          wc = beta*(-grv*ustar*thvstar*zii/thv)**0.333
          um = sqrt(ur*ur+wc*wc)
        endif

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

        obuold = obu
 
!-----------------------------------------------------------------------
! Test for convergence
!-----------------------------------------------------------------------
      it = it+1

      if(it .gt. itmin) then
         det  = max(del,del2)
         del = max(dele,dele2)
         if(det .lt. dtmin .AND. del .lt. dlemin) exit 
      endif
 
      end do ITERATION

! ======================================================================
!     END stability iteration 
! ======================================================================

! canopy fluxes and total assimilation amd respiration

      rst = rs / lai
      respc = respc + rsoil

! canopy fluxes and total assimilation amd respiration
      fsenl = fsenl + fsenl_dtl*dtl(it-1)

!     if(fsenl .lt. -100) then
!        print*, fsenl, fsenl_dtl*dtl(it-1), dtl(it-1), tl, taf, thm, tg
!        print*,wta0, wtg0, wtl0, zeta 
!        stop
!     endif


      etr = etr + etr_dtl*dtl(it-1)
      fevpl = fevpl + fevpl_dtl*dtl(it-1)

!     write(6, 100)  it-1, tl, fsenl, dlv*etr, dlv*fevpl, rs/lai, sabv, irab, zeta, rb, thm
100   format(i5,f9.1,10f9.1)
      ! wind stresses 
      taux  = taux - sigf*rhoair*us/ram
      tauy  = tauy - sigf*rhoair*vs/ram

!-----------------------------------------------------------------------
! fluxes from ground to canopy space
!-----------------------------------------------------------------------

      fseng = fseng + cp*rhoair*cgh*(tg-taf)
      fevpg = fevpg +   rhoair*cgw*(qg-qaf)

!-----------------------------------------------------------------------
! downward (upward) longwave radiation below (above) the canopy
!-----------------------------------------------------------------------
      dlrad = sigf * thermk * frl + stefnc * fac*tl**4 

      ulrad = stefnc * ( fac*tl**4 + sigf*thermk*emg*tg**4 )  

!-----------------------------------------------------------------------
! Derivative of soil energy flux with respect to soil temperature (cgrnd)
!-----------------------------------------------------------------------
      cgrnds = cgrnds + cp*rhoair*cgh*(1.-wtg0)
      cgrndl = cgrndl + rhoair*cgw*(1.-wtgq0)*dqgdT
      cgrnd  = cgrnds + cgrndl*htvp

!-----------------------------------------------------------------------
! balance check
!-----------------------------------------------------------------------
!     err = sabv + irab - fsenl - dlv*fevpl
!     if(abs(err) .gt. .2) write(6,*) 'energy balance in oneleaf.f90',err

!     err = sabv + irab + dirab_dtl*dtl(it-1) - fsenl - dlv*fevpl
!     if(abs(err) .ne. 0.) fsenl = fsenl + err

!-----------------------------------------------------------------------
! Update dew accumulation (kg/m2)
!-----------------------------------------------------------------------
      ldew = max(0.,ldew + (etr-fevpl)*dtime)

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


  END SUBROUTINE oneleaf