thermal.f90

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

  SUBROUTINE thermal ( ivt       ,lb        ,iub       ,dtime    ,&
             csol     ,porsl     ,phi0      ,bsw       ,dkmg     ,&
             dkdry    ,dksatu    ,lai       ,sai       ,z0m      ,&
             displa   ,sqrtdi    ,rootfr    ,effcon    ,vmax25   ,&
             slti     ,hlti      ,shti      ,hhti      ,trda     ,&
             trdm     ,trop      ,gradm     ,binter    ,extkn    ,&
             hu       ,ht        ,hq        ,us        ,vs       ,&
             tm       ,qm        ,rhoair    ,psrf      ,pco2m    ,&
             po2m     ,cosz      ,parsun    ,parsha    ,sabvsun  ,&
             sabvsha  ,sabg      ,frl       ,extkb     ,extkd    ,&
             thermk   ,fsno      ,sigf      ,dz        ,z        ,&
             zi       ,tlsun     ,tlsha     ,tss       ,wice     ,&
             wliq     ,ldew      ,scv       ,snowdp    ,imelt    ,&
             taux     ,tauy      ,fsena     ,fevpa     ,lfevpa   ,&
             fsenl    ,fevpl     ,etr       ,fseng     ,fevpg    ,&
             olrg     ,fgrnd     ,etrc      ,qseva     ,qsdew    ,&
             qsubl    ,qfros     ,sm        ,tref      ,trad     ,&
             rst      ,assim     ,respc     ,errore    ,solisb    ,solisd)

!=======================================================================
!      Source file: thermal.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! this is the main subroutine to execute the calculation of thermal processes
! and surface fluxes
! 
! (1) Leaf temperature
!     Foliage energy conservation is given by foliage energy budget equation
!                      Rnet - Hf - LEf = 0
!     The equation is solved by Newton-Raphson iteration, in which this iteration
!     includes the calculation of the photosynthesis and stomatal resistance, and the
!     integration of turbulent flux profiles. The sensible and latent heat
!     transfer between foliage and atmosphere and ground is linked by the equations:
!                      Ha = Hf + Hg and Ea = Ef + Eg
! 
! (2) Snow and soil temperatures
!     o The volumetric heat capacity is calculated as a linear combination
!       in terms of the volumetric fraction of the constituent phases.
!     o The thermal conductivity of soil is computed from
!       the algorithm of Johansen (as reported by Farouki 1981), and of snow is from
!       the formulation used in SNTHERM (Jordan 1991).
!     o Boundary conditions:
!       F = Rnet - Hg - LEg (top),  F= 0 (base of the soil column).
!     o Soil / snow temperature is predicted from heat conduction
!       in 10 soil layers and up to 5 snow layers.
!       The thermal conductivities at the interfaces between two neighbor layers
!       (j, j+1) are derived from an assumption that the flux across the interface
!       is equal to that from the node j to the interface and the flux from the
!       interface to the node j+1. The equation is solved using the Crank-Nicholson
!       method and resulted in a tridiagonal system equation.
! 
! (3) Phase change (see clm_meltfreeze.f90)
! 
! REVISION HISTORY:
!  the original code for soil temperatures calculation was provided by
!  Robert E. Dickinson and Z.-L. Yang, Nov. 1996
!  15 September 1999: Yongjiu Dai; Initial code
! 
! FLOW DIAGRAM FOR clm_thermal.f90
! 
! thermal ===> qsadv
!              obuini
!              obult
!              oneleaf    |===>  
!              twoleaf    |===>  fractdew
!              thermalk          qsadv
!              tridia            obuini
!              meltfreeze        obult
!                                stomata
! 
!=======================================================================

  USE PHYCON_MODULE ! physical constant
  IMPLICIT NONE
! 
!----------------------- Dummy argument --------------------------------
!
  integer, INTENT(in) :: &
        ivt,         &! land cover type
        lb,          &! lower bound of array 
        iub           ! upper bound of array
 
  real, INTENT(in) :: &
        dtime,       &! model time step [second]

        ! soil physical parameters
        csol(1:iub), &! heat capacity of soil solids [J/(m3 K)]
        porsl(1:iub),&! soil porosity [-]
        phi0(1:iub), &! soil water suction, negative potential [m]
        bsw(1:iub),  &! clapp and hornbereger "b" parameter [-]
        dkmg(1:iub), &! thermal conductivity of soil minerals [W/m-K]
        dkdry(1:iub),&! thermal conductivity of dry soil [W/m-K]
        dksatu(1:iub), &! thermal conductivity of saturated soil [W/m-K]

        ! vegetation parameters
        lai,         &! adjusted leaf area index for seasonal variation [-]
        sai,         &! stem area index  [-]
        z0m,         &! roughness length, momentum [m]
        displa,      &! displacement height [m]
        sqrtdi,      &! inverse sqrt of leaf dimension [m**-0.5]
        rootfr(1:iub),&!
       
        effcon,      &! quantum efficiency of RuBP regeneration (mol CO2 / mol quanta)
        vmax25,      &! maximum carboxylation rate at 25 C at canopy top
                      ! the range : 30.e-6 <-> 100.e-6 (mol co2 m-2 s-1)
        slti,        &! slope of low temperature inhibition function      (0.2)
        hlti,        &! 1/2 point of low temperature inhibition function  (288.16)
        shti,        &! slope of high temperature inhibition function     (0.3)
        hhti,        &! 1/2 point of high temperature inhibition function (313.16)
        trda,        &! temperature coefficient in gs-a model             (1.3)
        trdm,        &! temperature coefficient in gs-a model             (328.16)
        trop,        &! temperature coefficient in gs-a model             (298.16)
        gradm,       &! conductance-photosynthesis slope parameter
        binter,      &! conductance-photosynthesis intercept
        extkn,       &!

        ! atmospherical variables and observational height
        hu,          &! observational height of wind [m]
        ht,          &! observational height of temperature [m]
        hq,          &! observational height of humidity [m]
        us,          &! wind component in eastward direction [m/s]
        vs,          &! wind component in northward direction [m/s]
        tm,          &! temperature at agcm reference height [kelvin]
        qm,          &! specific humidity at agcm reference height [kg/kg]
        rhoair,      &! density air [kg/m3]
        psrf,        &! atmosphere pressure at the surface [pa]
        pco2m,       &!
        po2m,        &!

        ! radiative fluxes
        cosz,        &! cosine of the solar zenith angle
        solisb,      &!
        solisd,      &!
        parsun,      &!
        parsha,      &!
        sabvsun,     &! solar radiation absorbed by vegetation [W/m2]
        sabvsha,     &! solar radiation absorbed by vegetation [W/m2]
        sabg,        &! solar radiation absorbed by ground [W/m2]
        frl,         &! atmospheric infrared (longwave) radiation [W/m2]
        extkb,       &!
        extkd,       &!
        thermk,      &!

        ! state variable (1)
        etrc,        &! maximum possible transpiration rate [mm/s]
        fsno,        &! fraction of ground covered by snow
        sigf,        &! fraction of veg cover, excluding snow-covered veg [-]
        dz(lb : iub),&! layer thickiness [m]
        z (lb : iub),&! node depth [m]
        zi(lb - 1 : iub)     ! interface depth [m]

        ! state variables (2)
  real, INTENT(inout) :: &
        tlsun,        &! sunlit leaf temperature [K]
        tlsha,        &! shaded leaf temperature [K]
        tss (lb : iub),   &! soil temperature [K]
        wice(lb : iub),   &! ice lens [kg/m2]
        wliq(lb : iub),   &! liqui water [kg/m2]
        ldew,        &! depth of water on foliage [kg/(m2 s)] 
        scv,         &! snow cover, water equivalent [mm, kg/m2]
        snowdp        ! snow depth [m]

  integer, INTENT(out) :: & 
       imelt(lb : iub)  ! flag for melting or freezing [-]
 
        ! Output fluxes
  real, INTENT(out) :: &
        taux,        &! wind stress: E-W [kg/m/s**2]
        tauy,        &! wind stress: N-S [kg/m/s**2]
        fsena,       &! sensible heat from canopy height to atmosphere [W/m2]
        fevpa,       &! evapotranspiration from canopy height to atmosphere [mm/s]
        lfevpa,      &! latent heat flux from canopy height to atmosphere [W/m2]
        fsenl,       &! ensible heat from leaves [W/m2]
        fevpl,       &! evaporation+transpiration from leaves [mm/s]
        etr,         &! transpiration rate [mm/s]
        fseng,       &! sensible heat flux from ground [W/m2]
        fevpg,       &! evaporation heat flux from ground [mm/s]
        olrg,        &! outgoing long-wave radiation from ground+canopy
        fgrnd,       &! ground heat flux [W/m2]

        qseva,       &!ground surface evaporation rate (mm h2o/s)
        qsdew,       &!ground surface dew formation (mm h2o /s) [+]
        qsubl,       &!sublimation rate from snow pack (mm h2o /s) [+]
        qfros,       &!surface dew added to snow pack (mm h2o /s) [+]

        sm,          &! rate of snowmelt [kg/(m2 s)]
        tref,        &! 2 m height air temperature [kelvin]
        trad,        &! radiative temperature [K]

        rst,         &! stomatal resistance (s m-1)
        assim,       &! assimilation
        respc         ! respiration

!------------------------ LOCAL VARIABLES ------------------------------
!
  integer i,j

  real at(lb : iub),     &!"a" vector for tridiagonal matrix
       bt(lb : iub),     &!"b" vector for tridiagonal matrix
       ct(lb : iub),     &!"c" vector for tridiagonal matrix
       rt(lb : iub),     &!"r" vector for tridiagonal solution
       tg,          &! ground surface temperature [K]

       cv(lb : iub),     &! heat capacity [J/(m2 K)]
       tk(lb : iub),     &! thermal conductivity [W/(m K)]
       tssbef(lb : iub), &! soil/snow temperature before update
       qred,        &! soil surface relative humidity

       z0mg,        &! roughness length over ground, momentum [m]
       z0hg,        &! roughness length over ground, sensible heat [m]
       z0qg,        &! roughness length over ground, latent heat [m]
       z0mv,        &! roughness length over vegetation, momentum [m]
       z0hv,        &! roughness length over vegetation, sensible heat [m]
       z0qv          ! roughness length over vegetation, latent heat [m]

  real htvp,        &! latent heat of vapor of water (or sublimation) [j/kg]
       fact(lb : iub),  &! used in computing tridiagonal matrix
       fn  (lb : iub),  &! heat diffusion through the layer interface [W/m2]
       fn1 (lb : iub),  &! heat diffusion through the layer interface [W/m2]
       dzm,         &! used in computing tridiagonal matrix
       dzp           ! used in computing tridiagonal matrix

  integer &
       niters,      &! maximum number of iterations for surface temperature
       iter,        &! iteration index
       nmozsgn       ! number of times moz changes sign

  real beta,        &! coefficient of conective velocity [-]
       zii,         &! convective boundary height [m]
       zldis,       &! reference height "minus" zero displacement heght [m]
       ur,          &! wind speed at reference height [m/s]
       thm,         &! intermediate variable (tm+0.0098*ht)
       th,          &! potential temperature (kelvin)
       thv,         &! virtual potential temperature (kelvin)
       dth,         &! diff of virtual temp. between ref. height and surface
       dqh,         &! diff of humidity between ref. height and surface
       dthv,        &! diff of vir. poten. temp. between ref. height and surface
       thvstar,     &! virtual potential temperature scaling parameter
       obu,         &! monin-obukhov length (m)
       zeta,        &! dimensionless height used in Monin-Obukhov theory
       wc,          &! convective velocity [m/s]
       um,          &! wind speed including the stablity effect [m/s]
       temp1,       &! relation for potential temperature profile
       temp2,       &! relation for specific humidity profile
       temp12m,     &! 
       ustar,       &! friction velocity [m/s]
       tstar,       &! temperature scaling parameter
       qstar,       &! moisture scaling parameter
       ram,         &! aerodynamical resistance [s/m]
       rah,         &! thermal resistance [s/m]
       raw,         &! moisture resistance [s/m]
       raih,        &! temporary variable [kg/m2/s]
       raiw,        &! temporary variable [kg/m2/s]
       emg           ! ground emissivity (0.97 for snow, 
                     ! glaciers and water surface; 0.96 for soil and wetland)

  real errore,      &! energy balnce error [w/m2]
       cgrnd,       &! deriv. of soil energy flux wrt to soil temp [w/m2/k]
       cgrndl,      &! deriv, of soil sensible heat flux wrt soil temp [w/m2/k]
       cgrnds,      &! deriv of soil latent heat flux wrt soil temp [w/m**2/k]
       hs,          &! net energy flux into the surface (w/m2)
       dhsdt,       &! d(hs)/dT

       eg,          &! water vapor pressure at temperature T [pa]
       qsatg,       &! saturated humidity [kg/kg]
       degdT,       &! d(eg)/dT
       qsatgdT,     &! d(qsatg)/dT

       fac,         &! soil wetness of surface layer
       psit,        &! negative potential of soil
       hr,          &! relative humidity
       phimax,      &!

       qg,          &! ground specific humidity [kg/kg]
       dqgdT,       &! d(qg)/dT
       wice0(lb : iub),  &! ice mass from previous time-step
       wliq0(lb : iub),  &! liquid mass from previous time-step
       wx,          &! patitial volume of ice and water of surface layer

       egsmax,      &! max. evaporation which soil can provide at one time step
       egidif,      &! the excess of evaporation over "egsmax"
       brr(lb : iub),    &! temporay set 
       xmf,         &! total latent heat of phase change of ground water
       dlrad,       &! downward longwave radiation blow the canopy [W/m2]
       ulrad,       &! upward longwave radiation above the canopy [W/m2]
       tinc,        &! temperature difference of two time step
       obuold        ! monin-obukhov length from previous iteration

  real fsun,        &!
       cint(3),     &!
       cintsun(3),  &!
       cintsha(3),  &!
       par,         &!
       sabv,        &!
       rstfac,      &!
       tl            !
 integer oneleafx
 real psnsun, psnsha
!
!=======================================================================
! [1] Initial set 
!=======================================================================
! fluxes 
      taux   = 0.;  tauy   = 0.;  fsena  = 0.;  fevpa  = 0.; 
      lfevpa = 0.;  fsenl  = 0.;  fevpl  = 0.;  etr    = 0.; 
      fseng  = 0.;  fevpg  = 0.;  dlrad  = 0.;  ulrad  = 0.;
      cgrnds = 0.;  cgrndl = 0.;  cgrnd  = 0.;  tref   = 0.;

! temperature and water mass from previous time step
      tg = tss(lb)
      do i = lb, iub
        tssbef(i) = tss(i); wice0(i) = wice(i); wliq0(i) = wliq(i)
      enddo
   print*,'initial set'
!=======================================================================
! [2] Specific humidity and its derivative at ground surface
!=======================================================================

      qred = 1.
      call qsadv(tg,psrf,eg,degdT,qsatg,qsatgdT)
      print*,'qsadv'
      if(ivt/=11 .AND. ivt/=15)then            ! NOT wetland and ice land
        wx   = (wliq(1)/rhowat + wice(1)/dice)/dz(1)
        fac  = min(1.,wx/porsl(1))
        fac  = max( fac, 0.01 )

        psit = -phi0(1) * fac ** (- bsw(1) )   ! psit = max(smpmin, psit)
        hr   = exp(psit/roverg/tg)
        qred = (1.-fsno)*hr + fsno
      endif

      qg = qred*qsatg  
      dqgdT = qred*qsatgdT

      if(qsatg > qm .AND. qm > qred*qsatg)then
        qg = qm; dqgdT = 0.
      endif
      print*,'Specific humidity and its derivative at ground surface'  
!=======================================================================
! [3] Leaf and ground surface temperature and fluxes
!=======================================================================
! 3.1 Propositional variables

! emissivity
      if(scv>0. .OR. ivt==15)then
         emg = 0.97
      else
         emg = 0.96
      endif

! latent heat, we are arbitrarily assumed that the sublimation occured 
! only as wliq = 0
      htvp = dlv
      if(wliq(lb) <= 0. .AND. wice(lb) > 0.) htvp = dls

! roughness length
         z0mv = z0m;  z0hv = 0.1*z0mv; z0qv = 0.1*z0mv
      if(fsno > 0.)then
         z0mg = zsno; z0hg = z0mg; z0qg = z0mg            ! initial set
      else
         z0mg = zlnd; z0hg = z0mg; z0qg = z0mg
      endif

! potential temperatur at the reference height
      beta = 1.      ! -  (in computing W_*)
      zii = 1000.    ! m  (pbl height)

      thm = tm + 0.0098*ht              ! intermediate variable equivalent to