vertical_diffusion.f90

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

                    u          ,v          ,t          ,q          ,dse        , &
                    pmid       ,pint       ,piln       ,pmln       ,pdel       , &
                    rpdel      ,zm         ,zi         ,ztodt      ,taux       , &
                    tauy       ,shflx      ,cflx       ,pblh       ,ustar      , &
                    kvh        ,kvm        ,tpert      ,qpert      ,cgs        , &
                    obklen     ,exner      ,dtk        )
!-----------------------------------------------------------------------
! Driver routine to compute vertical diffusion of momentum, moisture, trace 
! constituents and dry static energy. The new temperature is computed from
! the diffused dry static energy.

! Turbulent diffusivities and boundary layer nonlocal transport terms are 
! obtained from the turbulence module.
!-----------------------------------------------------------------------
    use turbulence,   only: trbintr
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                   ! chunk identifier
    integer, intent(in) :: ncol                    ! number of atmospheric columns

    real(r8), intent(in) :: exner(pcols,pver)      ! (ps/pmid)**cappa
    real(r8), intent(in) :: pmid(pcols,pver)       ! midpoint pressures
    real(r8), intent(in) :: pint(pcols,pverp)      ! interface pressures
    real(r8), intent(in) :: piln (pcols,pverp)     ! Log interface pressures
    real(r8), intent(in) :: pmln (pcols,pver)      ! Log midpoint pressures
    real(r8), intent(in) :: pdel(pcols,pver)       ! thickness between interfaces
    real(r8), intent(in) :: rpdel(pcols,pver)      ! 1./pdel
    real(r8), intent(in) :: t(pcols,pver)          ! temperature input
    real(r8), intent(in) :: ztodt                  ! 2 delta-t
    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(inout) :: u(pcols,pver)       ! u wind
    real(r8), intent(inout) :: v(pcols,pver)       ! v wind
    real(r8), intent(inout) :: q(pcols,pver,pcnst+pnats)! moisture and trace constituents
    real(r8), intent(inout) :: dse(pcols,pver)     ! dry static energy [J/kg]
    real(r8), intent(in) :: zm(pcols,pver)         ! midpoint geoptl height above sfc
    real(r8), intent(in) :: zi(pcols,pverp)        ! interface geoptl height above sfc

    real(r8), intent(out) :: pblh(pcols)           ! planetary boundary layer height
    real(r8), intent(out) :: ustar(pcols)          ! surface friction velocity
    real(r8), intent(out) :: kvh(pcols,pverp)      ! diffusivity for heat
    real(r8), intent(out) :: kvm(pcols,pverp)      ! viscosity (diffusivity for momentum)
    real(r8), intent(out) :: tpert(pcols)          ! convective temperature excess
    real(r8), intent(out) :: qpert(pcols)          ! convective humidity excess
    real(r8), intent(out) :: cgs(pcols,pverp)      ! counter-grad star (cg/flux)
    real(r8), intent(out) :: obklen(pcols)         ! Obukhov length
    real(r8), intent(out) :: dtk(pcols,pver)       ! T tendency from KE dissipation
!
!---------------------------Local workspace-----------------------------
    real(r8) :: tmpm(pcols,pver)                   ! potential temperature, ze term in tri-diag sol'n

    real(r8) :: ca(pcols,pver)                     ! -upper diag of tri-diag matrix
    real(r8) :: cc(pcols,pver)                     ! -lower diag of tri-diag matrix
    real(r8) :: dnom(pcols,pver)                   ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))

    real(r8) :: kqfs(pcols,pcnst+pnats)            ! kinematic surf constituent flux (kg/m2/s)
    real(r8) :: kvq(pcols,pverp)                   ! diffusivity for constituents
    real(r8) :: kq_scal(pcols,pverp)               ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
    real(r8) :: mkvisc                             ! molecular kinematic viscosity c*tint**(2/3)/rho
    real(r8) :: tint                               ! interface temperature
    real(r8) :: tmp1(pcols)                        ! temporary storage
    real(r8) :: tmpi1(pcols,pverp)                 ! heat cg term [J/kg/m], or interface KE dissipation
    real(r8) :: tmpi2(pcols,pverp)                 ! rho or dt*(g*rho)**2/dp at interfaces

    real(r8) :: dinp_u(pcols,pverp)                ! vertical difference at interfaces, input u
    real(r8) :: dinp_v(pcols,pverp)                ! vertical difference at interfaces, input v
    real(r8) :: dout_u                             ! vertical difference at interfaces, output u
    real(r8) :: dout_v                             ! vertical difference at interfaces, output v

    integer :: i,k,m                               ! longitude,level,constituent indices
    integer :: ktopbl(pcols)                       ! index of first midpoint inside pbl
    integer :: ktopblmn                            ! min value of ktopbl

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

! Compute the vertical differences of the input u,v for KE dissipation, interface k-
! Note, velocity=0 at surface, so then difference at the bottom interface is -u,v(pver)
    do i = 1, ncol
       dinp_u(i,1) = 0.
       dinp_v(i,1) = 0.
       dinp_u(i,pverp) = -u(i,pver)
       dinp_v(i,pverp) = -v(i,pver)
    end do
    do k = 2, pver
       do i = 1, ncol
          dinp_u(i,k) = u(i,k) - u(i,k-1)
          dinp_v(i,k) = v(i,k) - v(i,k-1)
       end do
    end do

! Compute potential temperature
    tmpm(:ncol,:pver) = t(:ncol,:pver) * exner(:ncol,:pver)

! Get diffusivities and c-g terms from turbulence module
    call trbintr(lchnk   ,ncol    ,                            &
                 tmpm    ,t       ,q        ,zm     ,zi      , &
                 pmid    ,u       ,v       ,taux    ,tauy    , &
                 shflx   ,cflx    ,obklen  ,ustar   ,pblh    , &
                 kvm     ,kvh     ,tmpi1   ,cgs     ,kqfs    , &
                 tpert   ,qpert   ,ktopbl  ,ktopblmn)

! Compute rho at interfaces p/RT,  Ti = (Tm_k + Tm_k-1)/2,  interface k-
    do k = 2, pver
       do i = 1, ncol
          tmpi2(i,k) = pint(i,k) * 2. / (rair*(t(i,k) + t(i,k-1)))
       end do
    end do
    do i = 1, ncol
       tmpi2(i,pverp) = pint(i,pverp) / (rair*t(i,pver))
    end do

! Copy turbulent diffusivity for constituents from heat diffusivity

    kvq(:ncol,:) = kvh(:ncol,:)

! Compute molecular kinematic viscosity, heat diffusivity and factor for constituent diffusivity
! For the moment, keep molecular diffusivities all the same 

    kq_scal(:ncol,ntop_molec) = 0.
    do k = ntop_molec+1, nbot_molec
       do i = 1, ncol
          tint     = 0.5*(t(i,k) + t(i,k-1))
          mkvisc   = km_fac * tint**(2./3.) / tmpi2(i,k)
          kvm(i,k) = kvm(i,k) + mkvisc
          kvh(i,k) = kvh(i,k) + mkvisc*pr_num
          kq_scal(i,k) = sqrt(tint) / tmpi2(i,k)
       end do
    end do

! Call gravity wave drag to get gw tendency and diffusivities
! ************ ADD GW CALL HERE *****************************

! Add gw momentum forcing and KE dissipation
! ************ ADD CODE HERE ********************************

! Add the nonlocal transport terms to dry static energy, specific humidity and 
! other constituents in the boundary layer.

    call vd_add_cg(                                                   &
         lchnk      ,ncol       ,                                     &
         ktopblmn   ,q          ,dse        ,tmpi2      ,kvh        , &
         tmpi1      ,cgs        ,kqfs       ,rpdel      ,ztodt      )

! Add the (explicit) surface fluxes to the lowest layer

    do i = 1, ncol
       tmp1(i)  = ztodt*gravit*rpdel(i,pver)
       u(i,pver) = u(i,pver) + tmp1(i) * taux(i)
       v(i,pver) = v(i,pver) + tmp1(i) * tauy(i)
       dse(i,pver) = dse(i,pver) + tmp1(i) * shflx(i)
    end do
    do m = 1, pcnst+pnats
       do i = 1, ncol
          q(i,pver,m) = q(i,pver,m) + tmp1(i) * cflx(i,m)
       end do
    end do

! Calculate dt * (g*rho)^2/dp at interior interfaces,  interface k-
    do k = 2, pver
       do i = 1, ncol
          tmpi2(i,k) = ztodt * (gravit*tmpi2(i,k))**2 / (pmid(i,k) - pmid(i,k-1))
       end do
    end do

! Diffuse momentum
    call vd_lu_decomp(                                        &
         lchnk    ,ncol     ,                                 &
         kvm      ,tmpi2    ,rpdel    ,ztodt    ,ca         , &
         cc       ,dnom     ,tmpm     ,ntop     ,nbot       )
    call vd_lu_solve(                                                &
         lchnk    ,ncol     ,                                        &
         u        ,ca       ,tmpm     ,dnom     ,ntop     ,nbot     )
    call vd_lu_solve(                                                &
         lchnk    ,ncol     ,                                        &
         v        ,ca       ,tmpm     ,dnom     ,ntop     ,nbot     )

! Calculate kinetic energy dissipation
! 1. compute dissipation term at interfaces
    k = pverp
    do i = 1, ncol
       tmpi1(i,1) = 0.
       tmpi1(i,k) = 0.5 * ztodt * gravit * &
            ( (-u(i,k-1) + dinp_u(i,k))*taux(i) + (-v(i,k-1)+dinp_v(i,k))*tauy(i) )
    end do
    do k = 2, pver
       do i = 1, ncol
          dout_u = u(i,k) - u(i,k-1)
          dout_v = v(i,k) - v(i,k-1)
          tmpi1(i,k) = 0.25 * tmpi2(i,k) * kvm(i,k) *   &
               (dout_u**2 + dout_v**2 + dout_u*dinp_u(i,k) + dout_v*dinp_v(i,k))
       end do
    end do
! 2. Compute dissipation term at midpoints, add to dry static energy
    do k = 1, pver
       do i = 1, ncol
          dtk(i,k) = (tmpi1(i,k+1) + tmpi1(i,k)) * rpdel(i,k)
          dse(i,k) = dse(i,k) + dtk(i,k)
       end do
    end do

! Diffuse dry static energy
    call vd_lu_decomp(                                      &
         lchnk    ,ncol     ,                               &
         kvh      ,tmpi2    ,rpdel    ,ztodt    ,ca       , &
         cc       ,dnom     ,tmpm     ,ntop     ,nbot     )
    call vd_lu_solve(                                               &
         lchnk    ,ncol     ,                                       &
         dse      ,ca       ,tmpm     ,dnom     ,ntop     ,nbot    )

! Diffuse constituents.
! New decomposition required only if kvq != kvh (gw or molec diffusion)
    call vd_lu_decomp(                                                &
         lchnk      ,ncol       ,                                     &
         kvq        ,tmpi2      ,rpdel      ,ztodt      ,ca         , &
         cc         ,dnom       ,tmpm       ,ntop_eddy  ,nbot_eddy  )
    do m = 1, pcnst+pnats
! update decomposition in molecular diffusion range
       if (ntop_molec < nbot_molec) then
          call vd_lu_qmolec(                                                &
               ncol       ,                                                 &
               kvq        ,kq_scal    ,mw_fac(m)  ,tmpi2      ,rpdel      , &
               ztodt      ,ca         ,cc         ,dnom       ,tmpm       , &
               ntop_molec ,nbot_molec )
       end if
       call vd_lu_solve(                                              &
            lchnk   ,ncol     ,                                       &
            q(1,1,m),ca       ,tmpm     ,dnom     ,ntop    ,nbot     )
    end do

    return
  end subroutine vdiff

!==============================================================================
  subroutine vd_add_cg (                                &
       lchnk    ,ncol     ,                             &
       ktopblmn ,q        ,dse      ,rhoi    ,kvh     , &
       cgh      ,cgs      ,kqfs     ,rpdel   ,ztodt   )
!-----------------------------------------------------------------------
! Add the "counter-gradient" term to the dry static energy and tracer fields
! in the boundary layer.
! Note, ktopblmn gives the minimum vertical index of the first midpoint
! within the boundary layer.

! This is really a boundary layer routine, but, since the pbl also supplies the 
! diffusivities in the boundary layer, there is not a clear divsion between
! pbl and eddy vertical diffusion code.
!------------------------------Arguments--------------------------------

    integer, intent(in) :: lchnk                   ! chunk identifier
    integer, intent(in) :: ncol                    ! number of atmospheric columns
    integer, intent(in) :: ktopblmn                ! min value of ktopbl (highest bl top)

    real(r8), intent(in) :: rhoi(pcols,pverp)      ! density (p/RT) at interfaces
    real(r8), intent(in) :: kvh(pcols,pverp)       ! coefficient for heat and tracers
    real(r8), intent(in) :: cgh(pcols,pverp)       ! countergradient term for heat [J/kg/m]
    real(r8), intent(in) :: cgs(pcols,pverp)       ! unscaled cg term for constituents
    real(r8), intent(in) :: kqfs(pcols,pcnst+pnats)! kinematic surf constituent flux (kg/m2/s)
    real(r8), intent(in) :: rpdel(pcols,pver)      ! 1./pdel  (thickness bet interfaces)
    real(r8), intent(in) :: ztodt                  ! 2 delta-t