!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module operators 5,5
!BOP
! !MODULE: operators
!
! !DESCRIPTION:
! This module contains routines for various common
! mathematical operators, including gradient, divergence, curl,
! and a vertical velocity calculation. Note that these routines
! do {\em not} update ghost cells so results contain invalid
! entries in some ghost cells.
!
! !REVISION HISTORY:
! SVN:$Id: operators.F90 808 2006-04-28 17:06:38Z njn01 $
! !USES:
use kinds_mod
use blocks
use domain_size
use domain
use constants
use grid
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: div, &
grad, &
zcurl, &
wcalc
!EOP
!BOC
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: div
! !INTERFACE:
subroutine div(k,DIV_OUT,UX,UY,this_block) 2
! !DESCRIPTION:
! This routine returns the divergence at T points (times the cell
! area) of a vector field defined at U points. The divergence
! operator is defined as:
!
! \begin{equation}
! \nabla\cdot{\bf\rm u} = {1\over{\Delta_y}} \delta_x (\Delta_y u_x) +
! {1\over{\Delta_x}} \delta_y (\Delta_x u_y)
! \end{equation}
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS
integer (int_kind), intent(in) :: k ! vertical level
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UX,UY ! vector field defined at U-points
type (block), intent(in) :: &
this_block ! block information for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
DIV_OUT ! divergence at T-points times cell area
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j, &! dummy counters
bid ! local block id
!-----------------------------------------------------------------------
!
! compute divergence using a 4 point stencil
!
!-----------------------------------------------------------------------
bid = this_block%local_id
DIV_OUT = c0
do j=2,ny_block
do i=2,nx_block
if (k <= KMT(i,j,bid)) then
DIV_OUT(i,j) = p5*(UX(i ,j )*DYU(i ,j ,bid) + &
UX(i ,j-1)*DYU(i ,j-1,bid) - &
UX(i-1,j )*DYU(i-1,j ,bid) - &
UX(i-1,j-1)*DYU(i-1,j-1,bid) + &
UY(i ,j )*DXU(i ,j ,bid) + &
UY(i-1,j )*DXU(i-1,j ,bid) - &
UY(i ,j-1)*DXU(i ,j-1,bid) - &
UY(i-1,j-1)*DXU(i-1,j-1,bid))
endif
end do
end do
!-----------------------------------------------------------------------
!EOC
end subroutine div
!***********************************************************************
!BOP
! !IROUTINE: grad
! !INTERFACE:
subroutine grad(k, GRADX, GRADY, F, this_block) 3
! !DESCRIPTION:
! This routine computes the gradient in i,j directions at U points
! based on field defined at T-points.
!
! \begin{eqnarray}
! \nabla_x(F) &=& \delta_x F \\
! \nabla_y(F) &=& \delta_y F
! \end{eqnarray}
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! vertical level
real (r8), dimension(nx_block,ny_block), intent(in) :: &
F ! field defined at T points
type (block), intent(in) :: &
this_block ! block information for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
GRADX,GRADY ! gradient in (i,j) direction at U points
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j, &! dummy counters
bid ! local block id
!-----------------------------------------------------------------------
!
! compute gradient with a 4 point stencil
!
!-----------------------------------------------------------------------
bid = this_block%local_id
GRADX = c0
GRADY = c0
do j=1,ny_block-1
do i=1,nx_block-1
if (k <= KMU(i,j,bid)) then
GRADX(i,j) = DXUR(i,j,bid)*p5*(F(i+1,j+1) - F(i ,j) - &
F(i ,j+1) + F(i+1,j))
GRADY(i,j) = DYUR(i,j,bid)*p5*(F(i+1,j+1) - F(i ,j) + &
F(i ,j+1) - F(i+1,j))
endif
end do
end do
!-----------------------------------------------------------------------
!EOC
end subroutine grad
!***********************************************************************
!BOP
! !IROUTINE: zcurl
! !INTERFACE:
subroutine zcurl(k,CURL,UX,UY,this_block) 2
! !DESCRIPTION:
! This function returns the z-component of the curl of a vector
! field defined at U points.
!
! \begin{equation}
! \hat{\bf\rm z}\cdot\nabla\times{\bf\rm u} =
! {1\over{\Delta_y}} \delta_x(\Delta_y u_y) -
! {1\over{\Delta_x}} \delta_y(\Delta_x u_x)
! \end{equation}
!
! The result is actually multiplied by cell area and returned at
! T points.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! vertical level
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UX,UY ! vector field defined at U-points
type (block), intent(in) :: &
this_block ! block information for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
CURL ! z.curl(Ux,Uy) at T-points times cell area
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j, &! dummy counters
bid ! local block index
!-----------------------------------------------------------------------
!
! compute curl using stencil similar to divergence 4 point stencil
!
!-----------------------------------------------------------------------
bid = this_block%local_id
CURL = c0
do j=2,ny_block
do i=2,nx_block
if (k <= KMT(i,j,bid)) then
CURL(i,j) = p5*(UY(i ,j )*DYU(i ,j ,bid) + &
UY(i ,j-1)*DYU(i ,j-1,bid) - &
UY(i-1,j )*DYU(i-1,j ,bid) - &
UY(i-1,j-1)*DYU(i-1,j-1,bid) - &
UX(i ,j )*DXU(i ,j ,bid) - &
UX(i-1,j )*DXU(i-1,j ,bid) + &
UX(i ,j-1)*DXU(i ,j-1,bid) + &
UX(i-1,j-1)*DXU(i-1,j-1,bid))
endif
end do
end do
!-----------------------------------------------------------------------
!EOC
end subroutine zcurl
!***********************************************************************
!BOP
! !IROUTINE: wcalc
! !INTERFACE:
subroutine wcalc(WWW,UUU,VVV,this_block) 3
! !DESCRIPTION:
! This function calculates vertical velocity $w$ at tops of T-cells
! from horizontal velocity field (at U points) by integrating the
! continuity equation.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UUU, VVV ! horizontal velocity at all vertical levels
type (block), intent(in) :: &
this_block ! block information for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km), intent(out) :: &
WWW ! vertical velocity at T points for all levels
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j, &! dummy counters
k, &! vert level index
bid ! local block id
real (r8) :: &
fvn,fvs,fue,fuw, &! advective fluxes across sides of box
wtkb ! vertical velocity at bottom of box
!-----------------------------------------------------------------------
!
! calculate velocity by integrating continuity equation at
! each horizontal grid point
!
!-----------------------------------------------------------------------
bid = this_block%local_id
WWW = c0
!-----------------------------------------------------------------------
!
! partial bottom cell case
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
!-----------------------------------------------------------------------
!
! integrate from bottom up for every horizontal grid point
!
!-----------------------------------------------------------------------
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
wtkb = c0 ! vertical velocity zero at bottom of bottom cell
do k=KMT(i,j,bid),1,-1
!***
!*** advection fluxes
!***
fue = p5*(UUU(i ,j ,k)*DYU(i ,j ,bid)* &
DZU(i ,j ,k,bid) + &
UUU(i ,j-1,k)*DYU(i ,j-1 ,bid)* &
DZU(i ,j-1,k,bid))
fuw = p5*(UUU(i-1,j ,k)*DYU(i-1,j ,bid)* &
DZU(i-1,j ,k,bid) + &
UUU(i-1,j-1,k)*DYU(i-1,j-1 ,bid)* &
DZU(i-1,j-1,k,bid))
fvn = p5*(VVV(i ,j ,k)*DXU(i ,j ,bid)* &
DZU(i ,j ,k,bid) + &
VVV(i-1,j ,k)*DXU(i-1,j ,bid)* &
DZU(i-1,j ,k,bid))
fvs = p5*(VVV(i ,j-1,k)*DXU(i ,j-1 ,bid)* &
DZU(i ,j-1,k,bid) + &
VVV(i-1,j-1,k)*DXU(i-1,j-1 ,bid)* &
DZU(i-1,j-1,k,bid))
!***
!*** vertical velocity at top of box from continuity eq.
!***
WWW(i,j,k) = wtkb - &
(fvn - fvs + fue - fuw)*TAREA_R(i,j,bid)
wtkb = WWW(i,j,k) ! top value becomes bottom for next pass
enddo
end do ! horizontal loops
end do
!-----------------------------------------------------------------------
!
! no partial bottom cells
!
!-----------------------------------------------------------------------
else ! no partial bottom cells
!-----------------------------------------------------------------------
!
! integrate from bottom up for every horizontal grid point
!
!-----------------------------------------------------------------------
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
wtkb = c0 ! vertical velocity zero at bottom of bottom cell
do k=KMT(i,j,bid),1,-1
!***
!*** advection fluxes
!***
fue = p5*(UUU(i ,j ,k)*DYU(i ,j ,bid) + &
UUU(i ,j-1,k)*DYU(i ,j-1,bid))
fuw = p5*(UUU(i-1,j ,k)*DYU(i-1,j ,bid) + &
UUU(i-1,j-1,k)*DYU(i-1,j-1,bid))
fvn = p5*(VVV(i ,j ,k)*DXU(i ,j ,bid) + &
VVV(i-1,j ,k)*DXU(i-1,j ,bid))
fvs = p5*(VVV(i ,j-1,k)*DXU(i ,j-1,bid) + &
VVV(i-1,j-1,k)*DXU(i-1,j-1,bid))
!***
!*** vertical velocity at top of box from continuity eq.
!***
WWW(i,j,k) = wtkb - dz(k)* &
(fvn - fvs + fue - fuw)*TAREA_R(i,j,bid)
wtkb = WWW(i,j,k) ! top value becomes bottom for next pass
end do ! vertical loop
end do ! horizontal loops
end do
endif ! partial bottom cells
!-----------------------------------------------------------------------
!EOC
end subroutine wcalc
!***********************************************************************
end module operators
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||