vertical_diffusion.f90

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

    real(r8), intent(inout) :: q(pcols,pver,pcnst+pnats)  ! moisture and trace constituent input
    real(r8), intent(inout) :: dse(pcols,pver)     ! dry static energy

!---------------------------Local workspace-----------------------------

    real(r8) :: qtm(pcols,pver)                    ! temporary trace constituent input

    logical lqtst(pcols)                            ! adjust vertical profiles
    integer i                                      ! longitude index
    integer k                                      ! vertical index
    integer m                                      ! constituent index

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

! Add counter-gradient to input static energy profiles

    do k=ktopblmn-1,pver
       do i=1,ncol
          dse(i,k) = dse(i,k) + ztodt*rpdel(i,k)*gravit  &
               *(rhoi(i,k+1)*kvh(i,k+1)*cgh(i,k+1)       &
               - rhoi(i,k  )*kvh(i,k  )*cgh(i,k  ))
       end do
    end do

! Add counter-gradient to input mixing ratio profiles.
! Check for neg q's in each constituent and put the original vertical
! profile back if a neg value is found. A neg value implies that the
! quasi-equilibrium conditions assumed for the countergradient term are
! strongly violated.

    do m = 1, pcnst+pnats
       qtm(:ncol,ktopblmn-1:pver) = q(:ncol,ktopblmn-1:pver,m)
       do k = ktopblmn-1, pver
          do i = 1, ncol
             q(i,k,m) = q(i,k,m) + ztodt*rpdel(i,k)*gravit*kqfs(i,m)     &
                  *(rhoi(i,k+1)*kvh(i,k+1)*cgs(i,k+1)                    &
                  - rhoi(i,k  )*kvh(i,k  )*cgs(i,k  ))
          end do
       end do
       lqtst(:ncol) = all(q(:ncol,ktopblmn-1:pver,m) >= qmincg(m), 2)
       do k = ktopblmn-1, pver
          q(:ncol,k,m) = merge (q(:ncol,k,m), qtm(:ncol,k), lqtst(:ncol))
       end do
    end do

    return
  end subroutine vd_add_cg

!==============================================================================
  subroutine vd_lu_decomp(                                          &
       lchnk      ,ncol       ,                                     &
       kv         ,tmpi       ,rpdel      ,ztodt      ,ca         , &
       cc         ,dnom       ,ze         ,ntop       ,             &
       nbot       )
!-----------------------------------------------------------------------
! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the 
! tridiagonal diffusion matrix. 
! The diagonal elements (1+ca(k)+cc(k)) are not required by the solver.
! Also determine ze factor and denominator for ze and zf (see solver).
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                   ! chunk identifier
    integer, intent(in) :: ncol                    ! number of atmospheric columns
    integer, intent(in) :: ntop                    ! top level to operate on
    integer, intent(in) :: nbot                    ! bottom level to operate on

    real(r8), intent(in) :: kv(pcols,pverp)        ! vertical diffusion coefficients
    real(r8), intent(in) :: tmpi(pcols,pverp)      ! dt*(g/R)**2/dp*pi(k+1)/(.5*(tm(k+1)+tm(k))**2
    real(r8), intent(in) :: rpdel(pcols,pver)      ! 1./pdel  (thickness bet interfaces)
    real(r8), intent(in) :: ztodt                  ! 2 delta-t

    real(r8), intent(out) :: ca(pcols,pver)        ! -upper diagonal
    real(r8), intent(out) :: cc(pcols,pver)        ! -lower diagonal
    real(r8), intent(out) :: dnom(pcols,pver)      ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                                                   ! 1./(b(k) - c(k)*e(k-1))
    real(r8), intent(out) :: ze(pcols,pver)        ! term in tri-diag. matrix system

!---------------------------Local workspace-----------------------------
    integer :: i                                   ! longitude index
    integer :: k                                   ! vertical index

!-----------------------------------------------------------------------
! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the 
! tridiagonal diffusion matrix. The diagonal elements  (cb=1+ca+cc) are
! a combination of ca and cc; they are not required by the solver.

    do k = nbot-1, ntop, -1
       do i = 1, ncol
          ca(i,k  ) = kv(i,k+1)*tmpi(i,k+1)*rpdel(i,k  )
          cc(i,k+1) = kv(i,k+1)*tmpi(i,k+1)*rpdel(i,k+1)
       end do
    end do

! The bottom element of the upper diagonal (ca) is zero (element not used).
! The subdiagonal (cc) is not needed in the solver.

    do i=1,ncol
       ca(i,nbot) = 0.
    end do

! Calculate e(k).  This term is 
! required in solution of tridiagonal matrix defined by implicit diffusion eqn.

    do i = 1, ncol
       dnom(i,nbot) = 1./(1. + cc(i,nbot))
       ze(i,nbot)   = cc(i,nbot)*dnom(i,nbot)
    end do
    do k = nbot-1, ntop+1, -1
       do i=1,ncol
          dnom(i,k) = 1./ (1. + ca(i,k) + cc(i,k) - ca(i,k)*ze(i,k+1))
          ze(i,k) = cc(i,k)*dnom(i,k)
       end do
    end do
    do i=1,ncol
       dnom(i,ntop) = 1./ (1. + ca(i,ntop) - ca(i,ntop)*ze(i,ntop+1))
    end do

    return
  end subroutine vd_lu_decomp

!==============================================================================
  subroutine vd_lu_qmolec(                                          &
       ncol       ,                                                 &
       kv         ,kq_scal    ,mw_facm    ,tmpi       ,rpdel      , &
       ztodt      ,ca         ,cc         ,dnom       ,ze         , &
       ntop       ,nbot       )
!-----------------------------------------------------------------------
! Add the molecular diffusivity to the turbulent and gw diffusivity for a 
! consitutent.
! Update the superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the 
! tridiagonal diffusion matrix, also ze and denominator.
!------------------------------Arguments--------------------------------
    integer, intent(in)     :: ncol                 ! number of atmospheric columns
    integer, intent(in)     :: ntop                 ! top level to operate on
    integer, intent(in)     :: nbot                 ! bottom level to operate on

    real(r8), intent(in)    :: kv(pcols,pverp)      ! vertical diffusion coefficients
    real(r8), intent(in)    :: kq_scal(pcols,pverp) ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
    real(r8), intent(in)    :: mw_facm              ! sqrt(1/M_q + 1/M_d) for this constituent
    real(r8), intent(in)    :: tmpi(pcols,pverp)    ! dt*(g/R)**2/dp*pi(k+1)/(.5*(tm(k+1)+tm(k))**2
    real(r8), intent(in)    :: rpdel(pcols,pver)    ! 1./pdel  (thickness bet interfaces)
    real(r8), intent(in)    :: ztodt                ! 2 delta-t

    real(r8), intent(inout) :: ca(pcols,pver)       ! -upper diagonal
    real(r8), intent(inout) :: cc(pcols,pver)       ! -lower diagonal
    real(r8), intent(inout) :: dnom(pcols,pver)     ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                                                    ! 1./(b(k) - c(k)*e(k-1))
    real(r8), intent(inout) :: ze(pcols,pver)       ! term in tri-diag. matrix system

!---------------------------Local workspace-----------------------------
    integer :: i                                    ! longitude index
    integer :: k                                    ! vertical index

    real(r8) :: kmq                                 ! molecular diffusivity for constituent

!-----------------------------------------------------------------------
! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the 
! tridiagonal diffusion matrix. The diagonal elements  (cb=1+ca+cc) are
! a combination of ca and cc; they are not required by the solver.

    do k = nbot-1, ntop, -1
       do i = 1, ncol
          kmq = kq_scal(i,k+1) * mw_facm
          ca(i,k  ) = (kv(i,k+1)+kmq) * tmpi(i,k+1) * rpdel(i,k  )
          cc(i,k+1) = (kv(i,k+1)+kmq) * tmpi(i,k+1) * rpdel(i,k+1)
       end do
    end do

! Calculate e(k).  This term is 
! required in solution of tridiagonal matrix defined by implicit diffusion eqn.

    do k = nbot, ntop+1, -1
       do i=1,ncol
          dnom(i,k) = 1./ (1. + ca(i,k) + cc(i,k) - ca(i,k)*ze(i,k+1))
          ze(i,k) = cc(i,k)*dnom(i,k)
       end do
    end do
    do i=1,ncol
       dnom(i,ntop) = 1./ (1. + ca(i,ntop) - ca(i,ntop)*ze(i,ntop+1))
    end do

    return
  end subroutine vd_lu_qmolec

!==============================================================================
  subroutine vd_lu_solve(                                           &
       lchnk      ,ncol       ,                                     &
       q          ,ca         ,ze         ,dnom       ,             &
       ntop       ,nbot       )
!-----------------------------------------------------------------------
! Solve the implicit vertical diffusion equation with zero flux boundary conditions.
! Actual surface fluxes are explicit (rather than implicit) and are applied separately. 
! Procedure for solution of the implicit equation follows Richtmyer and 
! Morton (1967,pp 198-200).

! The equation solved is
!
! -ca(k)*q(k+1) + cb(k)*q(k) - cc(k)*q(k-1) = d(k), 
!
! where d(k) is the input profile and q(k) is the output profile
!
! The solution has the form
!
! q(k)  = ze(k)*q(k-1) + zf(k)
!
! ze(k) = cc(k) * dnom(k)
!
! zf(k) = [d(k) + ca(k)*zf(k+1)] * dnom(k)
!
! dnom(k) = 1/[cb(k) - ca(k)*ze(k+1)] =  1/[1 + ca(k) + cc(k) - ca(k)*ze(k+1)]

! Note that the same routine is used for temperature, momentum and tracers,
! and that input variables are replaced.
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                   ! chunk identifier
    integer, intent(in) :: ncol                    ! number of atmospheric columns
    integer, intent(in) :: ntop                    ! top level to operate on
    integer, intent(in) :: nbot                    ! bottom level to operate on

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

    real(r8), intent(inout) :: q(pcols,pver)       ! constituent field

!---------------------------Local workspace-----------------------------
    real(r8) :: zf(pcols,pver)                     ! term in tri-diag solution

    integer i,k                                    ! longitude, vertical indices

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

! Calculate zf(k).  Terms zf(k) and ze(k) are required in solution of 
! tridiagonal matrix defined by implicit diffusion eqn.
! Note that only levels ntop through nbot need be solved for.

    do i = 1, ncol
       zf(i,nbot) = q(i,nbot)*dnom(i,nbot)
    end do
    do k = nbot-1, ntop, -1
       do i=1,ncol
          zf(i,k) = (q(i,k) + ca(i,k)*zf(i,k+1))*dnom(i,k)
       end do
    end do

! Perform back substitution

    do i=1,ncol
       q(i,ntop) = zf(i,ntop)
    end do
    do k=ntop+1, nbot, +1
       do i=1,ncol
          q(i,k) = zf(i,k) + ze(i,k)*q(i,k-1)
       end do
    end do

    return
  end subroutine vd_lu_solve

end module vertical_diffusion