twoleaf.f90

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

        evplwetsun_dtlsha = sigf * rhoair * fwet * laisun / rb * ( - wtlshaq0*qsatlshaDT )

        if(evplwetsun .ge. ldew/dtime*laisun/lai)then
           evplwetsun = ldew/dtime*laisun/lai
           evplwetsun_dtlsun = 0.
           evplwetsun_dtlsha = 0.
        endif

        evplsun = etrsun + evplwetsun
        evplsun_dtlsun = etrsun_dtlsun + evplwetsun_dtlsun
        evplsun_dtlsha = etrsun_dtlsha + evplwetsun_dtlsha

! ---------------
!_      evplsha = rhoair * cfshaw * ( (wtaq0 + wtgq0 + wtlsunq0)*qsatlsha - wtaq0*qm - wtgq0*qg - wtlsunq0*qsatlsun )
!_      evplsha_dtlsun = rhoair * cfshaw * ( - wtlsunq0 ) * qsatlsunDT 
!_      evplsha_dtlsha = rhoair * cfshaw * ( wtaq0 + wtgq0 + wtlsunq0 ) * qsatlshaDT 

        etrsha = sigf * rhoair * (1.-fwet) * delta2 * laisha / (rb + rssha) * ( (wtaq0 + wtgq0 + wtlsunq0)*qsatlsha &
               - wtaq0*qm - wtgq0*qg - wtlsunq0*qsatlsun )
        etrsha_dtlsun = sigf * rhoair * (1.-fwet) * delta2 * laisha / (rb + rssha) * ( - wtlsunq0*qsatlsunDT )
        etrsha_dtlsha = sigf * rhoair * (1.-fwet) * delta2 * laisha / (rb + rssha) * ( (wtaq0 + wtgq0 + wtlsunq0)*qsatlshaDT )

        if(etrsha .ge. sigf*etrc*laisha/lai)then
           etrsha = sigf*etrc*laisha/lai
           etrsha_dtlsun = 0.
           etrsha_dtlsha = 0.
        endif

        evplwetsha = sigf * rhoair * fwet * (laisha+sai) / (rb) * ( (wtaq0 + wtgq0 + wtlsunq0)*qsatlsha &
                    - wtaq0*qm - wtgq0*qg - wtlsunq0*qsatlsun )
        evplwetsha_dtlsun = sigf * rhoair * fwet * (laisha+sai) / rb * ( - wtlsunq0*qsatlsunDT )
        evplwetsha_dtlsha = sigf * rhoair * fwet * (laisha+sai) / rb * ( (wtaq0 + wtgq0 + wtlsunq0)*qsatlshaDT )

        if(evplwetsha .ge. ldew/dtime*(laisha+sai)/(lai+sai))then
           evplwetsha = ldew/dtime*(laisha+sai)/(lai+sai) 
           evplwetsha_dtlsun = 0.
           evplwetsha_dtlsha = 0
        endif

        evplsha = etrsha + evplwetsha

        evplsha_dtlsun = etrsha_dtlsun + evplwetsha_dtlsun
        evplsha_dtlsha = etrsha_dtlsha + evplwetsha_dtlsha

! functions and their derivatives with respect to temperatures
        ftsun = sabvsun + irabsun - senlsun - dlv*evplsun 
        ftsha = sabvsha + irabsha - senlsha - dlv*evplsha

        dftsunDTlsun = dirabsun_dtlsun - senlsun_dtlsun - dlv*evplsun_dtlsun
        dftsunDTlsha = dirabsun_dtlsha - senlsun_dtlsha - dlv*evplsun_dtlsha

        dftshaDTlsun = dirabsha_dtlsun - senlsha_dtlsun - dlv*evplsha_dtlsun
        dftshaDTlsha = dirabsha_dtlsha - senlsha_dtlsha - dlv*evplsha_dtlsha

!-----------------------------------------------------------------------
! difference of temperatures by quasi-newton-raphson method for the non-linear system equations
!-----------------------------------------------------------------------
        dtlsun(it) = - (ftsun * (dftshaDTlsha-(laisha+sai)*clai/dtime) - ftsha * dftsunDTlsha) &
               / ((dftsunDTlsun-laisun*clai/dtime) * (dftshaDTlsha-(laisha+sai)*clai/dtime) - dftsunDTlsha * dftshaDTlsun)

        dtlsha(it) = - (ftsun * dftshaDTlsun - ftsha * (dftsunDTlsun-laisun*clai/dtime)) &
               / (dftsunDTlsha * dftshaDTlsun - (dftsunDTlsun-laisun*clai/dtime) * (dftshaDTlsha-(laisha+sai)*clai/dtime))

        if(abs(dtlsun(it)).gt.delmax) dtlsun(it) = delmax*dtlsun(it)/abs(dtlsun(it))
        if(abs(dtlsha(it)).gt.delmax) dtlsha(it) = delmax*dtlsha(it)/abs(dtlsha(it))

!--->   if(dtlsun(it-1)*dtlsun(it) .lt. 0.) ndtlsgn1 = ndtlsgn1+1
!--->   if(dtlsha(it-1)*dtlsha(it) .lt. 0.) ndtlsgn2 = ndtlsgn2+1

        if(it.ge.2)then
        if(dtlsun(it-1)*dtlsun(it) .lt. 0.) dtlsun(it) = 0.5*(dtlsun(it-1) + dtlsun(it))
        if(dtlsha(it-1)*dtlsha(it) .lt. 0.) dtlsha(it) = 0.5*(dtlsha(it-1) + dtlsha(it))
        endif

        tlsun = tlsunbef + dtlsun(it)
        tlsha = tlshabef + dtlsha(it)

!-----------------------------------------------------------------------
! square roots differences of temperatures and fluxes for use as the condition of convergences
!-----------------------------------------------------------------------
        del  = sqrt( dtlsun(it)*dtlsun(it) + dtlsha(it)*dtlsha(it) )
        dele = (    dirabsun_dtlsun * dtlsun(it) +     dirabsun_dtlsha * dtlsha(it))**2 &
             + (    dirabsha_dtlsun * dtlsun(it) +     dirabsha_dtlsha * dtlsha(it))**2 &
             + (     senlsun_dtlsun * dtlsun(it) +      senlsun_dtlsha * dtlsha(it))**2 &
             + (     senlsha_dtlsun * dtlsun(it) +      senlsha_dtlsha * dtlsha(it))**2 &
             + ( dlv*evplsun_dtlsun * dtlsun(it) +  dlv*evplsun_dtlsha * dtlsha(it))**2 &
             + ( dlv*evplsha_dtlsun * dtlsun(it) +  dlv*evplsha_dtlsha * dtlsha(it))**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(tlsun,psrf,eisun,deisunDT,qsatlsun,qsatlsunDT)

        call qsadv(tlsha,psrf,eisha,deishaDT,qsatlsha,qsatlshaDT)
 
! update vegetation/ground surface temperature, canopy air temperature, 
! canopy air humidity
        taf = wta0*thm + wtg0*tg + wtlsun0*tlsun +  wtlsha0*tlsha

        qaf = wtaq0*qm + wtgq0*qg + wtlsunq0*qsatlsun +  wtlshaq0*qsatlsha

! update co2 partial pressure within canopy air
        gah2o  = 1.0/raw * tprcor/thm                     ! mol m-2 s-1
        pco2a = pco2m - 1.37*psrf/max(0.446,gah2o) &
              * (assimsun + assimsha - respcsun - respcsha - rsoil)

!-----------------------------------------------------------------------
! 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,0.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) obu = zldis/(-0.01)

        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

      assim = assimsun + assimsha
      respc = respcsun + respcsha + rsoil
      rst = 1./ (laisun/rssun + laisha/rssha)

      fsenl = senlsun + senlsha + senlsun_dtlsun*dtlsun(it-1) + senlsun_dtlsha*dtlsha(it-1) &
                                + senlsha_dtlsun*dtlsun(it-1) + senlsha_dtlsha*dtlsha(it-1)
      etr = etrsun + etrsha
      fevpl = evplsun + evplsha + evplsun_dtlsun*dtlsun(it-1) + evplsun_dtlsha*dtlsha(it-1) &
                                  + evplsha_dtlsun*dtlsun(it-1) + evplsha_dtlsha*dtlsha(it-1)

!     write(6,100) it-1, tlsun,tlsha,dlv*etr,rssun/laisun,rssha/laisha,rstfac,thm,fwet,laisun,laisha,&
!     sabvsun, sabvsha
100   format(i5, 2f9.1,3f9.1,7f9.2)

!     print 200, wta0*100. , wtg0*100. , wtlsun0*100. ,  wtlsha0*100., raw, zeta, taf-thm, fsenl

200   format(8f10.2)
      ! 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 * ( fac1*tlsun**4 + fac2*tlsha**4 )

      ulrad = stefnc * ( fac1*tlsun**4 + fac2*tlsha**4 + sigf*emg*thermk*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 = sabvsun+sabvsha &
!         + irabsun+dirabsun_dtlsun*dtlsun(it-1)+dirabsun_dtlsha*dtlsha(it-1) &
!         + irabsha+dirabsha_dtlsun*dtlsun(it-1)+dirabsha_dtlsha*dtlsha(it-1) &
!         - fsenl-dlv*fevpl
!     if(abs(err) .ne. 0.) fsenl = fsenl + err
!
!     print*, 'sabvsun+sabvsha, irabsun+irabsha, fsenl, dlv*fevpl'
!     print*, sabvsun+sabvsha + irabsun+irabsha - fsenl - dlv*fevpl, &
!             sabvsun+sabvsha, irabsun+irabsha, fsenl, dlv*fevpl
!
!     if(abs(err) .gt. .2) write(6,*) 'energy balance in canopy X',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))

      if(dth<0..and. tref<thm) write(6,130)  'theta_star, taf, t2m, thm', temp1*dth, taf, tref, thm
 
130   format(1x,4f12.3)
  END SUBROUTINE twoleaf