lake.f90

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

      SUBROUTINE lake (msl     ,istep  ,idlak  ,dlat   ,dtime   ,&
                       zlak    ,dzlak  ,zilak  ,hu     ,ht      ,&
                       hq      ,us     ,vs     ,tm     ,qm      ,&
                       snowrate,rhoair ,psrf   ,sabg   ,frl     ,&
                       tg      ,tlak   ,wliq   ,wice   ,scv     ,&
                       snowdp  ,trad   ,tref   ,taux   ,tauy    ,&
                       fsena   ,fevpa  ,lfevpa ,fseng  ,fevpg   ,&
                       olrg    ,fgrnd  )

! ------------------------ code history ---------------------------
! source file:       lake.f90
! purpose:           lake temperature and snow on frozen lake
! author:            Gordon Bonan
! re-packed:         Yongjiu Dai, September 15, 1999
!
! ------------------------ notes ----------------------------------
! calculate lake temperatures from one-dimensional thermal
! stratification model based on eddy diffusion concepts to 
! represent vertical mixing of heat
!
! d ts    d            d ts     1 ds
! ---- = -- [(km + ke) ----] + -- --
!  dt    dz             dz     cw dz   
!
! where: ts = temperature (kelvin)
!         t = time (s)
!         z = depth (m)
!        km = molecular diffusion coefficient (m**2/s)
!        ke = eddy diffusion coefficient (m**2/s)
!        cw = heat capacity (j/m**3/kelvin)
!         s = heat source term (w/m**2)

! there are two types of lakes: 
!    deep lakes are 50 m. shallow lakes are 10 m deep.
!    for unfrozen deep lakes:    ke > 0 and    convective mixing
!    for unfrozen shallow lakes: ke = 0 and no convective mixing

! use crank-nicholson method to set up tridiagonal system of equations to
! solve for ts at time n+1, where the temperature equation for layer i is
! r_i = a_i [ts_i-1] n+1 + b_i [ts_i] n+1 + c_i [ts_i+1] n+1

! the solution conserves energy as

! cw*([ts(  1)] n+1 - [ts(  1)] n)*dz(  1)/dt + ... +
! cw*([ts(msl)] n+1 - [ts(msl)] n)*dz(msl)/dt = fin

! where 
! [ts] n   = old temperature (kelvin)
! [ts] n+1 = new temperature (kelvin)
! fin      = heat flux into lake (w/m**2)
!          = beta*sabg+frl-olrg-fsena-lfevpa-hm + phi(1) + ... + phi(msl) 
!
! -----------------------------------------------------------------
      USE phycon_module

      IMPLICIT NONE
      
! ------------------------ input/output variables -----------------

  integer, INTENT(in) :: &
        msl,       &!number of soil layers
        istep,     &!time step
        idlak       !index of lake, 1 = deep lake, 2 = shallow lake

  real, INTENT(in) :: &
        dlat,      &!latitude (radians)
        dtime,     &!time step (s)
        dzlak(msl),&!soil layer thickness (m)
        zlak(msl), &!depth (m)
        zilak(0:msl), &!
        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]
        snowrate,  &!rate of snowfall [mm/s]
        rhoair,    &!density air [kg/m3]
        psrf,      &!atmosphere pressure at the surface [pa]
        sabg,      &!solar radiation absorbed by ground [W/m2]
        frl         !atmospheric infrared (longwave) radiation [W/m2]

  real, INTENT(inout) :: &
        tg,        &!surface temperature (kelvin)
        tlak(msl), &!lake temperature (kelvin)
        wliq(msl), &!
        wice(msl), &!
        scv,       &!snow water equivalent [mm]
        snowdp      !snow depth [mm]

  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]
        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]

        tref,      &!2 m height air temperature [kelvin]
        trad        !radiative temperature [K]

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

  real  ax,        &!
        bx,        &!
        beta1,     &!coefficient of conective velocity [-]
        degdT,     &!d(eg)/dT
        dqh,       &!diff of humidity between ref. height and surface
        dth,       &!diff of virtual temp. between ref. height and surface
        dthv,      &!diff of vir. poten. temp. between ref. height and surface
        dzsur,     &!
        eg,        &!water vapor pressure at temperature T [pa]
        emg,       &!ground emissivity (0.97 for snow,
        errore,    &!lake temperature energy conservation error (w/m**2)
        hm,        &!energy residual [W/m2]
        htvp,      &!latent heat of vapor of water (or sublimation) [j/kg]
        obu,       &!monin-obukhov length (m)
        obuold,    &!monin-obukhov length of previous iteration
        qsatg,     &!saturated humidity [kg/kg]
        qsatgdT,   &!d(qsatg)/dT
        qstar,     &!moisture scaling parameter
        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) [+]
        qmelt,     &!snow melt [mm/s]
        ram,       &!aerodynamical resistance [s/m]
        rah,       &!thermal resistance [s/m]
        raw,       &!moisture resistance [s/m]
        stftg3,    &!
        temp1,     &!relation for potential temperature profile
        temp2,     &!relation for specific humidity profile
        temp12m,   &!
        tgbef,     &!
        thm,       &!intermediate variable (tm+0.0098*ht)
        th,        &!potential temperature (kelvin)
        thv,       &!virtual potential temperature (kelvin)
        thvstar,   &!virtual potential temperature scaling parameter
        tksur,     &!thermal conductivity of snow/soil (w/m/kelvin)
        tstar,     &!temperature scaling parameter
        um,        &!wind speed including the stablity effect [m/s]
        ur,        &!wind speed at reference height [m/s]
        ustar,     &!friction velocity [m/s]
        wc,        &!convective velocity [m/s]
        zeta,      &!dimensionless height used in Monin-Obukhov theory
        zii,       &!convective boundary height [m]
        zldis,     &!reference height "minus" zero displacement heght [m]
        z0mg,      &!roughness length over ground, momentum [m]
        z0hg,      &!roughness length over ground, sensible heat [m]
        z0qg        !roughness length over ground, latent heat [m]

  real  beta(2),   &!fraction solar rad absorbed at surface: depends on lake type
        za(2),     &!base of surface absorption layer (m): depends on lake type
        eta(2),    &!light extinction coefficient (/m): depends on lake type
        p0          !neutral value of turbulent prandtl number
     
  real  a(msl),    &!"a" vector for tridiagonal matrix
        b(msl),    &!"b" vector for tridiagonal matrix
        c(msl),    &!"c" vector for tridiagonal matrix
        r(msl),    &!"r" vector for tridiagonal solution
        rhow(msl), &!density of water (kg/m**3)
        phi(msl),  &!solar radiation absorbed by layer (w/m**2)
        kme(msl),  &!molecular + eddy diffusion coefficient (m**2/s)

        cwat,      &!specific heat capacity of water (j/m**3/kelvin)
        ws,        &!surface friction velocity (m/s)
        ks,        &!coefficient
        in,        &!relative flux of solar radiation into layer
        out,       &!relative flux of solar radiation out of layer
        ri,        &!richardson number
        fin,       &!heat flux into lake - flux out of lake (w/m**2)
        ocvts,     &!(cwat*(tlak[n  ])*dzlak
        ncvts,     &!(cwat*(tlak[n+1])*dzlak

        m1,        &!intermediate variable for calculating r, a, b, c
        m2,        &!intermediate variable for calculating r, a, b, c
        m3,        &!intermediate variable for calculating r, a, b, c
        ke,        &!eddy diffusion coefficient (m**2/s)
        km,        &!molecular diffusion coefficient (m**2/s)
        zin,       &!depth at top of layer (m)
        zout,      &!depth at bottom of layer (m)
        drhodz,    &!d [rhow] /dz (kg/m**4)
        n2,        &!brunt-vaisala frequency (/s**2)
        num,       &!used in calculating ri
        den,       &!used in calculating ri
        tav,       &!used in aver temp for convectively mixed layers
        nav,       &!used in aver temp for convectively mixed layers
        phidum,    &!temporary value of phi
        u2m         !2 m wind speed (m/s)

  integer i,j       !do loop or array index

! -----------------------------------------------------------------
!*[1] constants and model parameters
! -----------------------------------------------------------------
! constants for lake temperature model
      beta = (/0.4, 0.4/)                              ! (deep lake, shallow lake)
      za   = (/0.6, 0.5/)    
      eta  = (/0.1, 0.5/)  
      p0   = 1.  

! aerodynamical roughness 
      if(tg >= tfrz)then                               ! for unfrozen lake
         z0mg = 0.001
      else                                             ! for frozen lake
         z0mg = 0.04
      end if
      z0hg = z0mg
      z0qg = z0mg

! latent heat 
      if (tm > tfrz) then
         htvp = dlv
      else
         htvp = dlm
      end if

! surface emissivity
      emg = 0.97

! ----------------------------------------------------------------------
!*[2] SURFACE TEMPERATURE and FLUXES
! ----------------------------------------------------------------------

      dzsur = dzlak(1) + snowdp

      call qsadv(tg,psrf,eg,degdT,qsatg,qsatgdT)

! potential temperatur at the reference height

      beta1=1.       ! -  (in computing W_*)
      zii = 1000.    ! m  (pbl height)

      thm = tm + 0.0098*ht              ! intermediate variable equivalent to
                                        ! 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 1

! Initialization variables

      nmozsgn = 0
      obuold = 0.

      dth   = thm-tg
      dqh   = qm-qsatg
      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)

! ----------------------------------------------------------------------
      niters = 3

      do iter = 1, niters         ! begin stability iteration
         tgbef = tg
         if(tg > tfrz) then
            tksur = tkwat
         else
            tksur = tkice
         end if