thermal.f90

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

                                        ! tm*(pgcm/psrf)**(rgas/cp)
      th = tm*(100000./psrf)**(rgas/cp) ! potential T
      thv = th*(1.+0.61*qm)             ! virtual potential T
      ur = max(0.1,sqrt(us*us+vs*vs))   ! limit set to 0.1
      
!-----------------------------------------------------------------------
! 3.2 BARE PART
!     Compute sensible and latent fluxes and their derivatives with respect 
!     to ground temperature using ground temperatures from previous time step.
!-----------------------------------------------------------------------

      if(sigf <= 0.999) then

! Initialization variables
         nmozsgn = 0
         obuold = 0.

         dth   = thm-tg
         dqh   = qm-qg
         dthv  = dth*(1.+0.61*qm)+0.61*th*dqh
         zldis = hu-0.

         call obuini(ur,th,thm,thv,dth,dqh,dthv,zldis,z0mg,um,obu)
 
! Evaluated stability-dependent variables using moz from prior iteration
        niters=3
        do iter = 1, niters         ! begin stability iteration
           call obult(hu,ht,hq,0.,z0mg,z0hg,z0qg,obu,um,ustar,temp1,temp2,temp12m)

           tstar = temp1*dth
           qstar = temp2*dqh

           z0hg = z0mg/exp(0.13 * (ustar*z0mg/1.5e-5)**0.45)
           z0qg = z0hg

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

           obuold = obu

        enddo                       ! end stability iteration
!
! Get derivative of fluxes with repect to ground temperature
        ram    = 1./(ustar*ustar/um)
        rah    = 1./(temp1*ustar) 
        raw    = 1./(temp2*ustar) 

        raih   = (1.-sigf)*rhoair*cp/rah
        raiw   = (1.-sigf)*rhoair/raw          
        cgrnds = raih
        cgrndl = raiw*dqgdT
        cgrnd  = cgrnds + htvp*cgrndl

! surface fluxes of momentum, sensible and latent 
! using ground temperatures from previous time step
        taux   = -(1.-sigf)*rhoair*us/ram        
        tauy   = -(1.-sigf)*rhoair*vs/ram
        fseng  = -raih*dth
        fevpg  = -raiw*dqh 

        fsena  = fseng
        fevpa  = fevpg

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

! Equate canopy temperature to air over bareland.
!     Needed as sigf=0 carried over to next time step
        if(sigf < 0.001)then
           tlsun = tm; tlsha = tm; ldew = 0.
        endif
      end if
      print*,'sigf',sigf  
!-----------------------------------------------------------------------
! 3.3 VEGETATED PART
!     Calculate canopy temperature, latent and sensible fluxes from the canopy,
!     and leaf water change by evapotranspiration
!-----------------------------------------------------------------------
      print*,'vegetated part' 
      par = parsun + parsha; sabv = sabvsun + sabvsha
      print*,'sigf',sigf
      if(sigf > 0.001) then
         print*,'enter leaf'
       ! fraction of sunlit and shaded leaves of canopy
         fsun = ( 1. - exp(-extkb*lai) ) / max( extkb*lai, 1.e-6 )
         if(cosz<=0.0 .OR. sabv<1.) fsun = 0.
         print*,'fsun',fsun
       ! soil water strees factor on stomatal resistance
         phimax = -1.5e5
         rstfac = 0.
         do i = 1, iub
            wx = wliq(i)/(rhowat*dz(i))
            fac = min(1.0,wx/porsl(i))
            fac = max( fac, 0.001 )
            psit = -phi0(i) * fac ** (- bsw(i) )
            psit = max(phimax, psit)
            rstfac = rstfac + rootfr(i)*min(1.,(phimax-psit)/(phimax+phi0(i))) 
         enddo

oneleafx = 2
if(oneleafx == 0) then

            cint(1) = (1.-exp(-extkn*lai))/extkn
            cint(2) = (1.-exp(-extkd*lai))/extkd
            cint(3) = lai
            tl = tlsha

            call oneleaf ( dtime    ,htvp     ,lai      ,sai      ,displa   ,&
                 sqrtdi   ,z0mv     ,z0hv     ,z0qv     ,effcon   ,vmax25   ,&
                 slti     ,hlti     ,shti     ,hhti     ,trda     ,trdm     ,&
                 trop     ,gradm    ,binter   ,&

                 hu       ,ht       ,hq       ,us       ,vs       ,thm      ,&
                 th       ,thv      ,qm       ,psrf     ,rhoair   ,cosz     ,&
                 par      ,sabv     ,frl      ,cint     ,thermk   ,rstfac   ,&
                 po2m     ,pco2m    ,sigf     ,etrc     ,tg       ,qg       ,&
                 dqgdT    ,emg      ,&

                 tl       ,ldew     ,taux     ,tauy     ,fseng    ,&
                 fevpg    ,cgrnd    ,cgrndl   ,cgrnds   ,tref     ,rst      ,assim    ,&
                 respc    ,fsenl    ,fevpl    ,etr      ,dlrad    ,ulrad    )

            tlsun = tl; tlsha = tl

else if (oneleafx == 1) then
         if(fsun.le.0.01)then
            fsun = 0.
            cint(1) = (1.-exp(-extkn*lai))/extkn 
            cint(2) = (1.-exp(-extkd*lai))/extkd 
            cint(3) = lai
            tl = tlsha

            call oneleaf ( dtime    ,htvp     ,lai      ,sai      ,displa   ,&
                 sqrtdi   ,z0mv     ,z0hv     ,z0qv     ,effcon   ,vmax25   ,&
                 slti     ,hlti     ,shti     ,hhti     ,trda     ,trdm     ,&
                 trop     ,gradm    ,binter   ,&

                 hu       ,ht       ,hq       ,us       ,vs       ,thm      ,&
                 th       ,thv      ,qm       ,psrf     ,rhoair   ,cosz     ,&
                 par      ,sabv     ,frl      ,cint     ,thermk   ,rstfac   ,&
                 po2m     ,pco2m    ,sigf     ,etrc     ,tg       ,qg       ,&
                 dqgdT    ,emg      ,&

                 tl       ,ldew     ,taux     ,tauy     ,fseng    ,&
                 fevpg    ,cgrnd    ,cgrndl   ,cgrnds   ,tref     ,rst      ,assim    ,&
                 respc    ,fsenl    ,fevpl    ,etr      ,dlrad    ,ulrad    )

            tlsun = tl; tlsha = tl

         else
 
! scaling-up coefficients from leaf to canopy
            cintsun(1) = (1.-exp(-(extkn+extkb)*lai))/(extkn+extkb)
            cintsun(2) = (1.-exp(-(extkb+extkd)*lai))/(extkb+extkd)
            cintsun(3) = (1.-exp(-extkb*lai))/extkb
  
            cintsha(1) = (1.-exp(-extkn*lai))/extkn - cintsun(1)
            cintsha(2) = (1.-exp(-extkd*lai))/extkd - cintsun(2)
            cintsha(3) = lai - cintsun(3)

            call twoleaf ( dtime    ,htvp     ,lai      ,sai      ,displa   ,&
                 sqrtdi   ,z0mv     ,z0hv     ,z0qv     ,effcon   ,vmax25   ,&
                 slti     ,hlti     ,shti     ,hhti     ,trda     ,trdm     ,&
                 trop     ,gradm    ,binter   ,&

                 hu       ,ht       ,hq       ,us       ,vs       ,thm      ,&
                 th       ,thv      ,qm       ,psrf     ,rhoair   ,cosz     ,&
                 parsun   ,parsha   ,sabvsun  ,sabvsha  ,frl      ,fsun     ,&
                 cintsun  ,cintsha  ,thermk   ,rstfac   ,po2m     ,pco2m    ,&
                 sigf     ,etrc     ,tg       ,qg       ,dqgdT    ,emg      ,&

                 tlsun    ,tlsha    ,ldew     ,taux     ,tauy     ,fseng    ,&
                 fevpg    ,cgrnd    ,cgrndl   ,cgrnds   ,tref     ,rst      ,assim    ,&
                 respc    ,fsenl    ,fevpl    ,etr      ,dlrad    ,ulrad    )
  
         endif

else if (oneleafx == 2) then
            tl = tlsha
            call leaftem ( dtime  ,sigf  ,z0mv   ,z0hv   ,z0qv   ,&
                           lai    ,sai   ,displa ,sqrtdi ,&
                           hu     ,ht    ,hq     ,us     ,vs     ,&
                           qm     ,thm   ,th     ,thv    ,psrf   ,&
                           rhoair ,cosz  ,sabv   ,solisb ,solisd ,&
                           frl    ,tg    ,tl     ,ldew   ,qg     ,&
                           dqgdT  ,htvp  ,&
                           emg    ,etrc  ,fsenl  ,fevpl  ,etr    ,&
                           fseng  ,fevpg ,taux   ,tauy   ,tref   ,&
                           dlrad  ,ulrad ,cgrnds ,cgrndl ,cgrnd  ,rst, psnsun ,psnsha ) 
            assim = (psnsun + psnsha)*1.e-6 
            respc = 0.
            tlsun = tl; tlsha = tl
endif

      endif

!=======================================================================
! [4] Gound temperature
!=======================================================================
print*,'ground temperature'
! 4.1 Thermal conductivity and Heat capacity
      call thermalk ( ivt  ,lb   ,iub   ,&
                      tss  ,wice ,wliq  ,scv    ,&
                      csol ,dkmg ,dkdry ,dksatu ,porsl ,&
                      dz   ,z    ,zi    ,tk     ,cv )

! 4.2 net ground heat flux into the surface and its temperature derivative
      hs    = sabg + dlrad &
            + (1.-sigf)*emg*frl - emg*stefnc*tg**4 &
            - (fseng+fevpg*htvp) 

      dhsdT = - cgrnd - 4.*emg * stefnc * tg**3

      j       = lb
      fact(j) = dtime / cv(j) &
              * dz(j) / (0.5*(z(j)-zi(j-1)+capr*(z(j+1)-zi(j-1))))

      do j = lb + 1, iub
        fact(j) = dtime/cv(j)
      enddo

      do j = lb, iub - 1
        fn(j) = tk(j)*(tss(j+1)-tss(j))/(z(j+1)-z(j))
      enddo
      fn(iub) = 0.

! 4.3 Set up vector r and vectors a, b, c that define tridiagonal matrix
      j     = lb
      dzp   = z(j+1)-z(j)
      at(j) = 0.
      bt(j) = 1+(1.-cn)*fact(j)*tk(j)/dzp-fact(j)*dhsdT
      ct(j) =  -(1.-cn)*fact(j)*tk(j)/dzp
      rt(j) = tss(j) +  fact(j)*( hs - dhsdT*tss(j) + cn*fn(j) )


      do j = lb + 1, iub - 1
         dzm   = (z(j)-z(j-1))
         dzp   = (z(j+1)-z(j))
         at(j) =   - (1.-cn)*fact(j)* tk(j-1)/dzm
         bt(j) = 1.+ (1.-cn)*fact(j)*(tk(j)/dzp + tk(j-1)/dzm)
         ct(j) =   - (1.-cn)*fact(j)* tk(j)/dzp
         rt(j) = tss(j) + cn*fact(j)*( fn(j) - fn(j-1) )
      end do

      j     =  iub
      dzm   = (z(j)-z(j-1))
      at(j) =   - (1.-cn)*fact(j)*tk(j-1)/dzm
      bt(j) = 1.+ (1.-cn)*fact(j)*tk(j-1)/dzm
      ct(j) = 0.
      rt(j) = tss(j) - cn*fact(j)*fn(j-1)


! 4.4 Solve for tss
      i = size(at)
      call tridia (i ,at ,bt ,ct ,rt ,tss) 

!=======================================================================
! [5] Melting or Freezing 
!=======================================================================

      do j = lb, iub - 1
         fn1(j) = tk(j)*(tss(j+1)-tss(j))/(z(j+1)-z(j))
      enddo
      fn1(iub) = 0.

      j = lb
      brr(j) = cn*fn(j) + (1.-cn)*fn1(j)

      do j = lb + 1, iub
         brr(j) = cn*(fn(j)-fn(j-1)) + (1.-cn)*(fn1(j)-fn1(j-1))
      enddo

      call meltf ( lb, iub, dtime, &
                   fact(lb), brr(lb), hs, dhsdT, &
                   tssbef(lb), tss(lb), wliq(lb), wice(lb), imelt(lb), &
                   scv, snowdp, sm, xmf )
      tg = tss(lb)

!=======================================================================
! [6] Correct fluxes to present soil temperature
!=======================================================================
 
      tinc = tss(lb) - tssbef(lb)
      fseng =  fseng + tinc*cgrnds 
      fevpg =  fevpg + tinc*cgrndl

!     calculation of evaporative potential; flux in kg m**-2 s-1.  
!     egidif holds the excess energy if all water is evaporated
!     during the timestep.  this energy is later added to the
!     sensible heat flux.

      egsmax = (wice(lb)+wliq(lb)) / dtime

      egidif = max( 0., fevpg - egsmax )
      fevpg = min ( fevpg, egsmax )
      fseng = fseng + htvp*egidif

! ground heat flux
      fgrnd = sabg + dlrad + (1.-sigf)*emg*frl &
            - emg*stefnc*tssbef(lb)**3*(tssbef(lb) + 4.*tinc) &
            - (fseng+fevpg*htvp)

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

      fsena = fsenl + fseng
      fevpa = fevpl + fevpg
      lfevpa= dlv*fevpl + htvp*fevpg   ! W/m2 (accouting for sublimation)
      
      qseva = 0.
      qsubl = 0.
      qfros = 0.
      qsdew = 0.

      if(fevpg >= 0.)then
! not allow for sublimation in melting (melting ==> evap. ==> sublimation)
         qseva = min(wliq(lb)/dtime, fevpg)
         qsubl = fevpg - qseva
      else
         if(tg < tfrz)then
            qfros = abs(fevpg)
         else
            qsdew = abs(fevpg)
         endif
      endif

! outgoing long-wave radiation from canopy + ground
      olrg = ulrad &
           + (1.-sigf)*(1.-emg)*frl &
           + (1.-sigf)*emg*stefnc * tssbef(lb)**4 &
! for conservation we put the increase of ground longwave to outgoing
           + 4.*emg*stefnc*tssbef(lb)**3*tinc

! radiative temperature
      trad = (olrg/stefnc)**0.25

!=======================================================================
! [7] energy balance error
!=======================================================================

      errore = sabv + sabg + frl - olrg &
          - fsena - dlv*fevpl - htvp*fevpg - xmf
      do j = lb, iub
         errore = errore - (tss(j)-tssbef(j))/fact(j)
      enddo

!     if(abs(errore)>.2) &
!     write(6,*) 'energy  balance violation in thermal.f90',errore,sabv,sabg,frl,olrg,fsena,dlv*fevpl,htvp*fevpg,xmf

!     if(tlsun > 350.) stop

  END SUBROUTINE thermal