vertical_diffusion.f90

来自「CCSM Research Tools: Community Atmospher」· F90 代码 · 共 783 行 · 第 1/3 页

#include <misc.h>
#include <params.h>

module vertical_diffusion

!---------------------------------------------------------------------------------
! Module to compute vertical diffusion of momentum,  moisture, trace constituents
! and temperature. Uses turbulence module to get eddy diffusivities, including boundary
! layer diffusivities and nonlocal transport terms. Computes and applies molecular
! diffusion, if model top is high enough (above ~90 km).
!
! calling sequence:
!
!    vd_inti        initializes vertical diffustion constants
!      trbinti      initializes turblence/pbl constants
!     .
!     .
!    vd_intr            interface for vertical diffusion and pbl scheme
!      vdiff            performs vert diff and pbl
!        trbintr        compute turbulent diffusivities and boundary layer cg terms
!        vd_add_cg      add boudary layer cg term to profiles
!        vd_lu_decomp   perform lu decomposition of vertical diffusion matrices
!        vd_lu_qmolec   update decomposition for molecular diffusivity of a constituent
!        vd_lu_solve    solve the euler backward problem for vertical diffusion 
!
!---------------------------Code history--------------------------------
! Standardized:      J. Rosinski, June 1992
! Reviewed:          P. Rasch, B. Boville, August 1992
! Reviewed:          P. Rasch, April 1996
! Reviewed:          B. Boville, April 1996
! rewritten:         B. Boville, May 2000
! more changes       B. Boville, May-Aug 2001
!---------------------------------------------------------------------------------

  use precision,    only: r8
  use ppgrid,       only: pcols, pver, pverp
  use constituents, only: pcnst, pnats, cnst_name, cnst_add, qmincg
  use tracers,      only: ixcldw
  use pmgrid,       only: masterproc

  implicit none

  private          ! Make default type private to the module
  save
!
! Public interfaces
!
  public vd_inti   ! Initialization
  public vd_intr   ! Full routine
!
! Private data
!
  real(r8), private :: cpair      ! Specific heat of dry air
  real(r8), private :: gravit     ! Acceleration due to gravity
  real(r8), private :: rair       ! Gas const for dry air
  real(r8), private :: zvir       ! rh2o/rair - 1

  real(r8), parameter :: km_fac = 3.55E-7 ! molecular viscosity constant
  real(r8), parameter :: pr_num = 1.       ! Prandtl number
  real(r8), parameter :: d0     = 1.52E20  ! diffusion factor m-1 s-1 molec sqrt(kg/kmol/K)
                                           ! Aerononmy, Part B, Banks and Kockarts (1973), p39
                                           ! note text cites 1.52E18 cm-1 ...
  real(r8), parameter :: n_avog = 6.022E26 ! Avogadro's number (molec/kmol)

  real(r8), parameter :: kq_fac = 2.52E-13 ! molecular constituent diffusivity constant
  real(r8), parameter :: mw_dry = 29.      ! molecular weight of dry air ****REMOVE THIS***
  real(r8) :: mw_fac(pcnst+pnats)          ! sqrt(1/M_q + 1/M_d) in constituent diffusivity

  integer, private :: ntop        ! Top level to which vertical diffusion is applied (=1).
  integer, private :: nbot        ! Bottom level to which vertical diffusion is applied (=pver).
  integer, private :: ntop_eddy   ! Top level to which eddy vertical diffusion is applied.
  integer, private :: nbot_eddy   ! Bottom level to which eddy vertical diffusion is applied (=pver).
  integer, private :: ntop_molec  ! Top level to which molecular vertical diffusion is applied (=1).
  integer, private :: nbot_molec  ! Bottom level to which molecular vertical diffusion is applied.

  character(len=8), private :: vdiffnam(pcnst+pnats) ! names of v-diff tendencies

contains

!===============================================================================
  subroutine vd_inti(cpairx  ,cpwvx   ,gravx   ,rairx,  zvirx,    &
                     hypm    ,vkx)
!-----------------------------------------------------------------------
! Initialization of time independent fields for vertical diffusion.
! Calls initialization routine for boundary layer scheme.
!-----------------------------------------------------------------------
    use history,    only: addfld, add_default, phys_decomp
    use turbulence, only: trbinti
!------------------------------Arguments--------------------------------
    real(r8), intent(in) :: cpairx                 ! specific heat of dry air
    real(r8), intent(in) :: cpwvx                  ! spec. heat of water vapor at const. pressure
    real(r8), intent(in) :: gravx                  ! acceleration due to gravity
    real(r8), intent(in) :: rairx                  ! gas constant for dry air
    real(r8), intent(in) :: zvirx                  ! rh2o/rair - 1
    real(r8), intent(in) :: hypm(pver)             ! reference pressures at midpoints
    real(r8), intent(in) :: vkx                    ! Von Karman's constant

!---------------------------Local workspace-----------------------------
    integer :: k                                   ! vertical loop index

!-----------------------------------------------------------------------
! Set physical constants for vertical diffusion and pbl:
    cpair  = cpairx
    gravit = gravx
    rair   = rairx
    zvir   = zvirx

! Range of levels to which v-diff is applied.
    ntop_eddy  = 1       ! no reason not to make this 1, if >1, must be <= nbot_molec
    nbot_eddy  = pver    ! should always be pver
    ntop_molec = 1       ! should always be 1
    nbot_molec = 0       ! should be set below about 70 km

! Molecular diffusion turned on above ~60 km (50 Pa) if model top is above ~90 km (.1 Pa).
! Note that computing molecular diffusivities is a trivial expense, but constituent
! diffusivities depend on their molecular weights. Decomposing the diffusion matric
! for each constituent is a needless expense unless the diffusivity is significant.
    if (hypm(1) .lt. 0.1) then
       do k = 1, pver
          if (hypm(k) .lt. 50.) nbot_molec  = k
       end do
       if (masterproc) then
          write (6,fmt='(a,i3,5x,a,i3)') 'NTOP_MOLEC =',ntop_molec, 'NBOT_MOLEC =',nbot_molec
          write (6,fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY  =',ntop_eddy,  'NBOT_EDDY  =',nbot_eddy
       endif
    end if

! The vertical diffusion solver must operate over the full range of molecular and eddy diffusion
    ntop = 1
    nbot = pver

! Molecular weight factor in constitutent diffusivity
! ***** FAKE THIS FOR NOW USING MOLECULAR WEIGHT OF DRY AIR FOR ALL TRACERS ****
    do k = 1, pcnst+pnats
       mw_fac(k) = d0 * mw_dry * sqrt(1./mw_dry + 1./mw_dry) / n_avog
    end do

! Initialize turbulence variables
    call trbinti(gravit, cpair, rair, zvir, ntop_eddy, nbot_eddy,   hypm, vkx)

! Set names of diffused variable tendencies and declare them as history variables
    do k = 1, pcnst+pnats
       vdiffnam(k) = 'VD'//cnst_name(k)
       if(k==1)vdiffnam(k) = 'VD01'    !**** compatibility with old code ****
       call addfld (vdiffnam(k),'kg/kg/s ',pver, 'A','Vertical diffusion of '//cnst_name(k),phys_decomp)
    end do

! Only tracer 1 (Q) is output by default
    call add_default (vdiffnam(1), 1, ' ')

    return
  end subroutine vd_inti

!===============================================================================
  subroutine vd_intr(                                     &
       ztodt    ,state    ,                               &
       taux     ,tauy     ,shflx    ,cflx     ,pblh     , &
       tpert    ,qpert    ,ustar    ,obklen   ,ptend    ) 
!-----------------------------------------------------------------------
! interface routine for vertical diffusion and pbl scheme
!-----------------------------------------------------------------------
    use physics_types,  only: physics_state, physics_ptend
    use history,        only: outfld
!!$    use geopotential, only: geopotential_dse
!------------------------------Arguments--------------------------------
    real(r8), intent(in) :: taux(pcols)            ! x surface stress (n)
    real(r8), intent(in) :: tauy(pcols)            ! y surface stress (n)
    real(r8), intent(in) :: shflx(pcols)           ! surface sensible heat flux (w/m2)
    real(r8), intent(in) :: cflx(pcols,pcnst+pnats)! surface constituent flux (kg/m2/s)
    real(r8), intent(in) :: ztodt                  ! 2 delta-t

    type(physics_state), intent(in)  :: state      ! Physics state variables
    type(physics_ptend), intent(inout)  :: ptend   ! indivdual parameterization tendencies

    real(r8), intent(out) :: pblh(pcols)           ! planetary boundary layer height
    real(r8), intent(out) :: tpert(pcols)          ! convective temperature excess
    real(r8), intent(out) :: qpert(pcols,pcnst+pnats) ! convective humidity and constituent excess
    real(r8), intent(out) :: ustar(pcols)          ! surface friction velocity
    real(r8), intent(out) :: obklen(pcols)         ! Obukhov length
!
!---------------------------Local storage-------------------------------
    integer :: lchnk                               ! chunk identifier
    integer :: ncol                                ! number of atmospheric columns
    integer :: i,k,m                               ! longitude,level,constituent indices

    real(r8) :: dtk(pcols,pver)                    ! T tendency from KE dissipation
    real(r8) :: kvh(pcols,pverp)                   ! diffusion coefficient for heat
    real(r8) :: kvm(pcols,pverp)                   ! diffusion coefficient for momentum
    real(r8) :: cgs(pcols,pverp)                   ! counter-gradient star (cg/flux)
    real(r8) :: rztodt                             ! 1./ztodt

!-----------------------------------------------------------------------
! local constants
    rztodt = 1./ztodt

    lchnk = state%lchnk
    ncol  = state%ncol

! Operate on copies of the input states, convert to tendencies at end.
! The input values of u,v are required for computing the kinetic energy dissipation.

    ptend%q(:ncol,:,:) = state%q(:ncol,:,:)
    ptend%s(:ncol,:) = state%s(:ncol,:)
    ptend%u(:ncol,:) = state%u(:ncol,:)
    ptend%v(:ncol,:) = state%v(:ncol,:)

! All variables are modified by vertical diffusion

    ptend%name  = "vertical diffusion"
    ptend%lq(:) = .TRUE.
    ptend%ls    = .TRUE.
    ptend%lu    = .TRUE.
    ptend%lv    = .TRUE.

! Call the real vertical diffusion code.
    call vdiff(                                                            &
         lchnk       ,ncol        ,                                        &
         ptend%u     ,ptend%v     ,state%t     ,ptend%q     ,ptend%s     , &
         state%pmid  ,state%pint  ,state%lnpint,state%lnpmid,state%pdel  , &
         state%rpdel ,state%zm    ,state%zi    ,ztodt       ,taux        , &
         tauy        ,shflx       ,cflx        ,pblh        ,ustar       , &
         kvh         ,kvm         ,tpert       ,qpert(1,1)  ,cgs         , &
         obklen      ,state%exner ,dtk         )

! Convert the new profiles into vertical diffusion tendencies.
! Convert KE dissipative heat change into "temperature" tendency.

    do k = 1, pver
       do i = 1, ncol
          ptend%s(i,k) = (ptend%s(i,k) - state%s(i,k)) * rztodt
          ptend%u(i,k) = (ptend%u(i,k) - state%u(i,k)) * rztodt
          ptend%v(i,k) = (ptend%v(i,k) - state%v(i,k)) * rztodt
       end do
       do m = 1, pcnst+pnats
          do i = 1, ncol
             ptend%q(i,k,m) = (ptend%q(i,k,m) - state%q(i,k,m))*rztodt
          end do
       end do
    end do
#ifdef PERGRO
! For pergro case, do not diffuse cloud water: replace with input values
    ptend%lq(ixcldw) = .FALSE.
    ptend%q(:ncol,:,ixcldw) = 0.
#endif

! Save the vertical diffusion variables on the history file
    dtk(:ncol,:) = dtk(:ncol,:)/cpair                ! normalize heating for history
    call outfld ('DTVKE   ',dtk,pcols,lchnk)
    dtk(:ncol,:) = ptend%s(:ncol,:)/cpair            ! normalize heating for history using dtk
    call outfld ('DTV     ',dtk    ,pcols,lchnk)
    call outfld ('DUV     ',ptend%u,pcols,lchnk)
    call outfld ('DVV     ',ptend%v,pcols,lchnk)
    do m = 1, pcnst+pnats
       call outfld(vdiffnam(m),ptend%q(1,1,m),pcols,lchnk)
    end do

    return
  end subroutine vd_intr

!===============================================================================
  subroutine vdiff (lchnk      ,ncol       ,                                     &