!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module advection 2,26
!BOP
! !MODULE: advection
!
! !DESCRIPTION:
! This module contains arrays and variables necessary for performing
! advection of momentum and tracer quantities. Currently, the
! module supports leapfrog centered advection of momentum and
! both leapfrog centered advection and third-order upwinding of
! tracers.
!
! !REVISION HISTORY:
! SVN:$Id: advection.F90 22881 2010-05-11 04:23:39Z njn01 $
! !USES:
use POP_KindsMod
use POP_ErrorMod
use POP_CommMod
use POP_FieldMod
use POP_GridHorzMod
use POP_HaloMod
use kinds_mod
, only: r4, r8, int_kind, char_len, log_kind, rtavg
use constants
, only: c0, c1, p5, p125, p25, blank_fmt, delim_fmt, ndelim_fmt, c2, &
field_loc_center, field_type_scalar, &
field_loc_Eface, field_type_vector
use blocks
, only: nx_block, ny_block, block, get_block
use domain_size
use communicate
, only: my_task, master_task
use distribution
, only:
use grid
, only: dz, DXT, DYT, HUW, HUS, c2dz, KMT, HTE, UAREA_R, DZT, &
partial_bottom_cells, DYU, DZU, DXU, DZR, DZ2R, KMU, TAREA_R, HTN, &
sfc_layer_type, sfc_layer_varthick, FCORT, KMTE, &
KMTW, KMTEE, KMTN, KMTS, KMTNN, ugrid_to_tgrid
use domain
, only: nblocks_clinic, blocks_clinic, distrb_clinic, &
POP_haloClinic
use broadcast
, only: broadcast_scalar, broadcast_array
use diagnostics
, only: cfl_advect
use state_mod
, only: state
use operators
, only: zcurl
use tavg
, only: tavg_requested, ltavg_on, define_tavg_field, &
accumulate_tavg_field, tavg_in_which_stream
use io_types
, only: nml_in, nml_filename, stdout
use time_management
, only: max_blocks_clinic, km, nt, mix_pass, c2dtt
use timers
, only: timer_start, timer_stop, get_timer
use exit_mod
, only: sigAbort, exit_pop, flushm
use prognostic
, only: UVEL, VVEL, curtime, tracer_d
use passive_tracers
, only: tadvect_ctype_passive_tracers
use registry
use overflows
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_advection, &
comp_flux_vel_ghost, &
advt, &
advu
!EOP
!BOC
!-----------------------------------------------------------------------
!
! choices for tracer advection
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
tadvect_centered = 1, &! centered leap-frog tracer advection
tadvect_upwind3 = 2, &! 3rd-order upwind tracer advection
tadvect_lw_lim = 3 ! 1d Lax-Wendroff with 1d flux limiters
integer (int_kind), dimension(nt) :: &
tadvect_itype ! users tracer advection choice
!-----------------------------------------------------------------------
!
! coefficients for metric advection terms (KXU,KYU)
!
!-----------------------------------------------------------------------
real (r8), dimension (nx_block,ny_block,max_blocks_clinic) :: &
KXU,KYU
!-----------------------------------------------------------------------
!
! geometric arrays necessary for non-centered advection
!
!-----------------------------------------------------------------------
real (r8), dimension(:,:,:,:), allocatable :: &
AUX
!-----------------------------------------------------------------------
!
! geometric arrays necessary for upwind advection
!
!-----------------------------------------------------------------------
real (r8), dimension(:), allocatable :: &
talfzp,tbetzp,tgamzp, &
talfzm,tbetzm,tdelzm
real (r8), dimension(:,:,:), allocatable :: &
TALFXP,TBETXP,TGAMXP,TALFYP,TBETYP,TGAMYP, &
TALFXM,TBETXM,TDELXM,TALFYM,TBETYM,TDELYM
!-----------------------------------------------------------------------
!
! geometric arrays necessary for lw_lim advection
!
!-----------------------------------------------------------------------
real (r8), dimension(:), allocatable :: &
p5_dz_ph_r
real (r8), dimension(:,:,:), allocatable :: &
p5_DXT_ph_R, &! 1/(DXT(i,j)+DXT(i+1,j))
p5_DYT_ph_R, &! 1/(DYT(i,j)+DYT(i,j+1))
UTE_jbm2, &! UTE for j==jb-2
WTKB_jbm2, &! WTKB for j==jb-2
WTKB_jep2, &! WTKB for j==je+2
WTKB_ibm2, &! WTKB for i==ib-2
WTKB_iep2 ! WTKB for i==ie+2
real (r8), dimension(:,:,:,:), allocatable :: &
FLUX_VEL_prev, &! flux velocities from prev k iteration
UTE_to_UVEL_E, &! converts UTE to UVEL_E
VTN_to_VVEL_N ! converts VTN to VVEL_E
!-----------------------------------------------------------------------
!
! tavg ids for tavg diagnostics related to advection
! north, east, zonal, merid refer here to logical space only
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_WVEL, &! Vertical Velocity
tavg_WVEL2, &! Vertical Velocity Squared
tavg_UEU, &! flux of zonal momentum across east face
tavg_VNU, &! flux of zonal momentum across north face
tavg_WTU, &! flux of zonal momentum across top face
tavg_UEV, &! flux of merid momentum across east face
tavg_VNV, &! flux of merid momentum across north face
tavg_WTV, &! flux of merid momentum across top face
tavg_PV, &! potential vorticity
tavg_Q, &! z-derivative of pot density
tavg_PD, &! potential density
tavg_RHOU, &! pot density times U velocity
tavg_RHOV, &! pot density times V velocity
tavg_PVWM, &! pot vorticity flux through bottom
tavg_PVWP, &! pot vorticity flux through top
tavg_UPV, &! pot vorticity flux through east face
tavg_VPV, &! pot vorticity flux through north face
tavg_URHO, &! pot density flux through east face
tavg_VRHO, &! pot density flux through north face
tavg_WRHO, &! pot density flux through top face
tavg_UQ, &! advection of Q across east face
tavg_VQ ! advection of Q across north face
integer (int_kind), dimension(nt) :: &
tavg_ADV_TRACER, &! vertical average of tracer advective tendency
tavg_UE_TRACER, &! flux of tracer across east face
tavg_VN_TRACER, &! flux of tracer across north face
tavg_WT_TRACER ! flux of tracer across top face
!-----------------------------------------------------------------------
!
! advection timers
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
timer_advt, &! timer for tracer advection
timer_advu ! timer for momentum advection
!-----------------------------------------------------------------------
!
! advection flags
!
!-----------------------------------------------------------------------
logical (log_kind) :: &
luse_centered, &! is centered used by any tracer
luse_upwind3, &! is upwind3 used by any tracer
luse_lw_lim ! is lw_lim used by any tracer
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_advection
! !INTERFACE:
subroutine init_advection 1,51
! !DESCRIPTION:
! This subroutine initializes variables associated with advection
! schemes.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
!
! namelist input
!
!-----------------------------------------------------------------------
character (*), parameter :: &
docfmt = "( 3x, a16, 3x, a16 )"
character (char_len) :: &
tadvect_ctype ! character string for tracer advect choice
namelist /advect_nml/ tadvect_ctype
!-----------------------------------------------------------------------
!
! local variables for setting up upwind coefficients
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,k,n, &! dummy loop indices
iblock, &! local block number
nattempts, &! num of attempts to read namelist input
nml_error ! error flag for namelist read
real (r8) :: &
dxc,dxcw,dxce,dxce2, &
dyc,dycs,dycn,dycn2
real (r8), dimension(:), allocatable :: &
dzc
real (r8), dimension(nx_block,ny_block) :: &
WORK1 ! local temp space
type (block) :: &
this_block ! block information for current block
!-----------------------------------------------------------------------
!
! check that init_passive_tracers has been called
!
!-----------------------------------------------------------------------
if (.not. registry_match('init_passive_tracers')) then
call exit_POP
(sigAbort,'ERROR: init_passive_tracers must be ' /&
&/ 'called before init_advection')
end if
!-----------------------------------------------------------------------
!
! read input namelist for choice of tracer advection method
!
!-----------------------------------------------------------------------
tadvect_ctype = 'centered'
if (my_task == master_task) then
open (nml_in, file=nml_filename, status='old', iostat=nml_error)
if (nml_error /= 0) then
nml_error = -1
else
nml_error = 1
endif
!*** keep reading until find right namelist
do while (nml_error > 0)
read(nml_in, nml=advect_nml,iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
end if
call broadcast_scalar
(nml_error, master_task)
if (nml_error /= 0) then
call exit_POP
(sigAbort,'ERROR reading advection namelist')
endif
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,ndelim_fmt)
write(stdout,blank_fmt)
write(stdout,*) ' Advection:'
write(stdout,blank_fmt)
write(stdout,*) ' advect_nml namelist settings:'
write(stdout,blank_fmt)
write(stdout,advect_nml)
write(stdout,blank_fmt)
write(stdout,delim_fmt)
write(stdout,blank_fmt)
write(stdout,'(a17)') 'Advection options'
write(stdout,blank_fmt)
write(stdout,delim_fmt)
tadvect_itype = tadvect_ctype_to_tadvect_itype
(tadvect_ctype)
do n=3,nt
if (tadvect_ctype_passive_tracers(n) == 'base_model') then
tadvect_ctype_passive_tracers(n) = tadvect_ctype
else
tadvect_itype(n) = tadvect_ctype_to_tadvect_itype
( &
tadvect_ctype_passive_tracers(n))
endif
end do
if (all(tadvect_itype == tadvect_itype(1))) then
write(stdout,*) 'Using ' /&
&/ trim(tadvect_ctype) /&
&/ ' for all tracers.'
else
write(stdout,docfmt) 'Tracer Name', 'Advection Option'
write(stdout,docfmt) 'TEMP', trim(tadvect_ctype)
write(stdout,docfmt) 'SALT', trim(tadvect_ctype)
do n=3,nt
write(stdout,docfmt) trim(tracer_d(n)%short_name), &
trim(tadvect_ctype_passive_tracers(n))
end do
end if
endif
call broadcast_array
(tadvect_itype, master_task)
if (minval(tadvect_itype) < 0) then
call exit_POP
(sigAbort,'ERROR: unknown tracer advection method')
endif
!-----------------------------------------------------------------------
!
! advection type must be implemented by overflows if on and interactive
!
!-----------------------------------------------------------------------
if ( overflows_on .and. overflows_interactive ) then
do n=1,nt
if( tadvect_itype(n) == tadvect_centered ) then
write(stdout,*) &
'advection: ERROR: advection type not implemented ', &
'for interactive overflows'
write(stdout,*) ' advection type = ',tadvect_itype(n), &
' for tracer = ',n
write(stdout,*) ' only upwind3 and lw_lim implemented'
call exit_POP
(sigAbort,'ERROR non-implemented advection with overflows')
endif
enddo
endif
!-----------------------------------------------------------------------
!
! set use flags
!
!-----------------------------------------------------------------------
luse_centered = any(tadvect_itype == tadvect_centered)
luse_upwind3 = any(tadvect_itype == tadvect_upwind3)
luse_lw_lim = any(tadvect_itype == tadvect_lw_lim)
!-----------------------------------------------------------------------
!
! initialize metric advection coefficients
!
!-----------------------------------------------------------------------
do iblock = 1,nblocks_clinic
this_block = get_block
(blocks_clinic(iblock),iblock)
KXU(:,:,iblock) = (eoshift(HUW(:,:,iblock),dim=1,shift=1) - &
HUW(:,:,iblock))*UAREA_R(:,:,iblock)
KYU(:,:,iblock) = (eoshift(HUS(:,:,iblock),dim=2,shift=1) - &
HUS(:,:,iblock))*UAREA_R(:,:,iblock)
end do
!*** KXU,KYU only needed in physical domain
!*** no ghost cell update required assuming HUS,HUW,UAREA_R
!*** were defined correctly in ghost cells
!-----------------------------------------------------------------------
!
! allocate and initialize non-centered arrays if necessary
!
!-----------------------------------------------------------------------
if (any(tadvect_itype /= tadvect_centered)) then
allocate (AUX (nx_block,ny_block,nt,nblocks_clinic))
endif ! non-centered setup
!-----------------------------------------------------------------------
!
! allocate and initialize upwinding grid arrays if necessary
!
!-----------------------------------------------------------------------
if (luse_upwind3) then
allocate (talfzp(km), &
tbetzp(km), &
tgamzp(km), &
talfzm(km), &
tbetzm(km), &
tdelzm(km))
allocate (TALFXP(nx_block,ny_block,nblocks_clinic), &
TBETXP(nx_block,ny_block,nblocks_clinic), &
TGAMXP(nx_block,ny_block,nblocks_clinic), &
TALFYP(nx_block,ny_block,nblocks_clinic), &
TBETYP(nx_block,ny_block,nblocks_clinic), &
TGAMYP(nx_block,ny_block,nblocks_clinic), &
TALFXM(nx_block,ny_block,nblocks_clinic), &
TBETXM(nx_block,ny_block,nblocks_clinic), &
TDELXM(nx_block,ny_block,nblocks_clinic), &
TALFYM(nx_block,ny_block,nblocks_clinic), &
TBETYM(nx_block,ny_block,nblocks_clinic), &
TDELYM(nx_block,ny_block,nblocks_clinic))
allocate (dzc(0:km+1))
!***
!*** vertical grid coefficients
!***
dzc(0) = dz(1)
do k=1,km
dzc(k) = dz(k)
enddo
dzc(km+1) = dzc(km)
do k=1,km-1
talfzp(k) = dz(k)*(c2*dz(k)+dzc(k-1))/ &
((dz(k)+dz(k+1))* &
(dzc(k-1)+c2*dz(k)+dz(k+1)))
tbetzp(k) = dz(k+1)*(c2*dz(k)+dzc(k-1))/ &
((dz(k)+dz(k+1))* &
(dz(k)+dzc(k-1) ))
tgamzp(k) = -(dz(k)*dz(k+1))/ &
((dz(k)+dzc(k-1))* &
(dz(k+1)+dzc(k-1)+c2*dz(k)))
enddo
tbetzp(1) = tbetzp(1) + tgamzp(1)
tgamzp(1) = c0
talfzp(km) = c0
tbetzp(km) = c0
tgamzp(km) = c0
do k=1,km-1
talfzm(k) = dz(k)*(c2*dz(k+1)+dzc(k+2))/ &
((dz(k)+dz(k+1))* &
(dz(k+1)+dzc(k+2) ))
tbetzm(k) = dz(k+1)*(c2*dz(k+1)+dzc(k+2))/ &
((dz(k)+dz(k+1))* &
(dz(k)+dzc(k+2)+c2*dz(k+1)))
tdelzm(k) = -(dz(k)*dz(k+1))/ &
((dz(k+1)+dzc(k+2))* &
(dz(k)+dzc(k+2)+c2*dz(k+1)))
enddo
talfzm(km-1) = talfzm(km-1) + tdelzm(km-1)
tdelzm(km-1) = c0
talfzm(km) = c0
tbetzm(km) = c0
tdelzm(km) = c0
deallocate (dzc)
!***
!*** horizontal grid coeffs
!***
do iblock = 1,nblocks_clinic
this_block = get_block
(blocks_clinic(iblock),iblock)
!***
!*** zonal grid coefficients
!***
do j=this_block%jb,this_block%je
do i=this_block%ib-1,this_block%ie
dxc = DXT(i ,j,iblock)
dxcw = DXT(i-1,j,iblock)
dxce = DXT(i+1,j,iblock)
dxce2 = DXT(i+2,j,iblock)
TALFXP(i,j,iblock) = dxc*(c2*dxc+dxcw)/ &
((dxc+dxce)*(dxcw+c2*dxc+dxce))
TBETXP(i,j,iblock) = dxce*(c2*dxc+dxcw)/ &
((dxc+dxcw)*( dxc+dxce))
TGAMXP(i,j,iblock) = -( dxc*dxce)/ &
((dxc+dxcw)*(dxcw+c2*dxc+dxce))
TALFXM(i,j,iblock) = dxc *(c2*dxce+dxce2)/ &
((dxc +dxce)*( dxce+dxce2))
TBETXM(i,j,iblock) = dxce*(c2*dxce+dxce2)/ &
((dxc +dxce)*(dxc+c2*dxce+dxce2))
TDELXM(i,j,iblock) = -( dxc *dxce )/ &
((dxce2+dxce)*(dxc+c2*dxce+dxce2))
end do
end do
!***
!*** poloidal grid coefficients
!***
do j=this_block%jb-1,this_block%je
do i=this_block%ib,this_block%ie
dyc = DYT(i,j ,iblock)
dycs = DYT(i,j-1,iblock)
dycn = DYT(i,j+1,iblock)
dycn2 = DYT(i,j+2,iblock)
TALFYP(i,j,iblock) = dyc *(c2*dyc+dycs)/ &
((dyc+dycn)*(dycs+c2*dyc+dycn))
TBETYP(i,j,iblock) = dycn*(c2*dyc+dycs)/ &
((dyc+dycn)*(dycs+ dyc ))
TGAMYP(i,j,iblock) = -( dyc*dycn)/ &
((dyc+dycs)*(dycs+c2*dyc+dycn))
TALFYM(i,j,iblock) = dyc *(c2*dycn+dycn2)/ &
((dyc+dycn)*( dycn+dycn2))
TBETYM(i,j,iblock) = dycn*(c2*dycn+dycn2)/ &
((dyc+dycn)*(dyc+c2*dycn+dycn2))
TDELYM(i,j,iblock) = -( dyc *dycn )/ &
((dycn2+dycn)*(dyc+c2*dycn+dycn2))
end do
end do
end do
!*** assuming DXT,DYT were defined correctly in ghost cells
!*** these are valid from (ib-1:ie,jb-1:je) and are only
!*** accessed in that range so no halo update necessary
endif ! 3rd order upwind setup
!-----------------------------------------------------------------------
!
! allocate and initialize lw_lim grid arrays if necessary
!
!-----------------------------------------------------------------------
if (luse_lw_lim) then
allocate (p5_dz_ph_r(km), &
p5_DXT_ph_R(nx_block,ny_block,nblocks_clinic), &
p5_DYT_ph_R(nx_block,ny_block,nblocks_clinic), &
UTE_jbm2(nx_block,km,nblocks_clinic), &
WTKB_jbm2(nx_block,km,nblocks_clinic), &
WTKB_jep2(nx_block,km,nblocks_clinic), &
WTKB_ibm2(ny_block,km,nblocks_clinic), &
WTKB_iep2(ny_block,km,nblocks_clinic), &
FLUX_VEL_prev(nx_block,ny_block,5,nblocks_clinic) )
if (partial_bottom_cells) then
allocate (UTE_to_UVEL_E(nx_block,ny_block,km,nblocks_clinic), &
VTN_to_VVEL_N(nx_block,ny_block,km,nblocks_clinic))
else
allocate (UTE_to_UVEL_E(nx_block,ny_block,1,nblocks_clinic), &
VTN_to_VVEL_N(nx_block,ny_block,1,nblocks_clinic))
endif
p5_DXT_ph_R = c0
p5_DYT_ph_R = c0
UTE_jbm2 = c0
WTKB_jbm2 = c0
WTKB_jep2 = c0
WTKB_ibm2 = c0
WTKB_iep2 = c0
UTE_to_UVEL_E = c0
VTN_to_VVEL_N = c0
p5_dz_ph_r(1:km-1) = c1 / (dz(1:km-1) + dz(2:km))
p5_dz_ph_r(km) = p5 / dz(km)
!***
!*** horizontal grid coeffs
!***
do iblock = 1,nblocks_clinic
this_block = get_block
(blocks_clinic(iblock),iblock)
!***
!*** zonal grid coefficients
!*** compute in j halo because they are used
!***
do j=this_block%jb-2,this_block%je+2
do i=this_block%ib-2,this_block%ie+1
p5_DXT_ph_R(i,j,iblock) = c1 / (DXT(i,j,iblock) + &
DXT(i+1,j,iblock))
end do
end do
if (partial_bottom_cells) then
do k=1,km
do j=this_block%jb-2,this_block%je+2
do i=this_block%ib-2,this_block%ie+1
UTE_to_UVEL_E(i,j,k,iblock) = c1 / HTE(i,j,iblock) / &
min(DZT(i,j,k,iblock), DZT(i+1,j,k,iblock))
end do
end do
end do
else
do j=this_block%jb-2,this_block%je+2
do i=this_block%ib-2,this_block%ie+1
UTE_to_UVEL_E(i,j,1,iblock) = c1 / HTE(i,j,iblock)
end do
end do
endif
!***
!*** poloidal grid coefficients
!***
do j=this_block%jb-2,this_block%je+1
do i=this_block%ib,this_block%ie
p5_DYT_ph_R(i,j,iblock) = c1 / (DYT(i,j,iblock) + &
DYT(i,j+1,iblock))
end do
end do
if (partial_bottom_cells) then
do k=1,km
do j=this_block%jb-2,this_block%je+1
do i=this_block%ib,this_block%ie
VTN_to_VVEL_N(i,j,k,iblock) = c1 / HTN(i,j,iblock) / &
min(DZT(i,j,k,iblock), DZT(i,j+1,k,iblock))
end do
end do
end do
else
do j=this_block%jb-2,this_block%je+1
do i=this_block%ib,this_block%ie
VTN_to_VVEL_N(i,j,1,iblock) = c1 / HTN(i,j,iblock)
end do
end do
endif
end do
endif ! lw_lim setup
!-----------------------------------------------------------------------
!
! initialize timers
!
!-----------------------------------------------------------------------
call get_timer
(timer_advu,'ADVECTION_MOMENTUM', nblocks_clinic, &
distrb_clinic%nprocs)
if (all(tadvect_itype == tadvect_centered)) then
call get_timer
(timer_advt,'ADVECTION_TRACER_CENTERED', &
nblocks_clinic, distrb_clinic%nprocs)
else if (all(tadvect_itype == tadvect_upwind3)) then
call get_timer
(timer_advt,'ADVECTION_TRACER_UPWIND3', &
nblocks_clinic, distrb_clinic%nprocs)
else if (all(tadvect_itype == tadvect_lw_lim)) then
call get_timer
(timer_advt,'ADVECTION_TRACER_LW_LIM', &
nblocks_clinic, distrb_clinic%nprocs)
else
call get_timer
(timer_advt,'ADVECTION_TRACER', &
nblocks_clinic, distrb_clinic%nprocs)
endif
!-----------------------------------------------------------------------
!
! define tavg fields related to advection
!
!-----------------------------------------------------------------------
call define_tavg_field
(tavg_WVEL,'WVEL',3, &
long_name='Vertical Velocity', &
units='centimeter/s', grid_loc='3112', &
coordinates='TLONG TLAT z_w time')
call define_tavg_field
(tavg_WVEL2,'WVEL2',3, &
long_name='Vertical Velocity**2', &
units='centimeter^2/s^2', grid_loc='3112', &
coordinates='TLONG TLAT z_w time')
call define_tavg_field
(tavg_UEU,'UEU',3, &
long_name='East Flux of Zonal Momentum', &
units='cm/s^2', grid_loc='3321')
call define_tavg_field
(tavg_VNU,'VNU',3, &
long_name='North Flux of Zonal Momentum', &
units='cm/s^2', grid_loc='3231')
call define_tavg_field
(tavg_UEV,'UEV',3, &
long_name='East Flux of Meridional Momentum',&
units='cm/s^2', grid_loc='3321')
call define_tavg_field
(tavg_VNV,'VNV',3, &
long_name='North Flux of Meridional Momentum', &
units='cm/s^2', grid_loc='3231')
call define_tavg_field
(tavg_WTU,'WTU',3, &
long_name='Top flux of Zonal Momentum', &
units='cm/s^2', grid_loc='3222')
call define_tavg_field
(tavg_WTV,'WTV',3, &
long_name='Top flux of Meridional Momentum', &
units='cm/s^2', grid_loc='3222')
call define_tavg_field
(tavg_UE_TRACER(1),'UET',3, &
long_name='Flux of Heat in grid-x direction',&
units='degC/s', grid_loc='3211', &
coordinates='ULONG TLAT z_t time' )
call define_tavg_field
(tavg_VN_TRACER(1),'VNT',3, &
long_name='Flux of Heat in grid-y direction',&
units='degC/s', grid_loc='3121', &
coordinates='TLONG ULAT z_t time')
call define_tavg_field
(tavg_WT_TRACER(1),'WTT',3, &
long_name='Heat Flux Across Top Face', &
units='degC/s', grid_loc='3112', &
coordinates='TLONG TLAT z_w time' )
call define_tavg_field
(tavg_UE_TRACER(2),'UES',3, &
long_name='Salt Flux in grid-x direction', &
scale_factor=1000.0_rtavg, &
units='gram/kilogram/s', grid_loc='3211', &
coordinates='ULONG TLAT z_t time' )
call define_tavg_field
(tavg_VN_TRACER(2),'VNS',3, &
long_name='Salt Flux in grid-y direction', &
scale_factor=1000.0_rtavg, &
units='gram/kilogram/s', grid_loc='3121', &
coordinates='TLONG ULAT z_t time')
call define_tavg_field
(tavg_WT_TRACER(2),'WTS',3, &
long_name='Salt Flux Across Top Face', &
scale_factor=1000.0_rtavg, &
units='gram/kilogram/s', grid_loc='3112', &
coordinates='TLONG TLAT z_w time' )
call define_tavg_field
(tavg_ADV_TRACER(1),'ADVT',2, &
long_name='Vertically-Integrated T Advection Tendency',&
units='centimeter degC/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field
(tavg_ADV_TRACER(2),'ADVS',2, &
long_name='Vertically-Integrated S Advection Tendency', &
scale_factor=1000.0_rtavg, &
units='centimeter gram/kilogram/s', grid_loc='2110',&
coordinates='TLONG TLAT time')
do n=3,nt
call define_tavg_field
(tavg_UE_TRACER(n), &
'UE_' /&
&/ trim(tracer_d(n)%short_name),3, &
long_name=trim(tracer_d(n)%short_name) /&
&/ ' Flux in grid-x direction', &
units=trim(tracer_d(n)%tend_units), &
scale_factor=tracer_d(n)%scale_factor, &
grid_loc='3211', &
coordinates='ULONG TLAT z_t time' )
call define_tavg_field
(tavg_VN_TRACER(n), &
'VN_' /&
&/ trim(tracer_d(n)%short_name),3, &
long_name=trim(tracer_d(n)%short_name) /&
&/ ' Flux in grid-y direction', &
units=trim(tracer_d(n)%tend_units), &
scale_factor=tracer_d(n)%scale_factor, &
grid_loc='3121', &
coordinates='TLONG ULAT z_t time')
call define_tavg_field
(tavg_WT_TRACER(n), &
'WT_' /&
&/ trim(tracer_d(n)%short_name),3, &
long_name=trim(tracer_d(n)%short_name) /&
&/ ' Flux Across Top Face', &
units=trim(tracer_d(n)%tend_units), &
scale_factor=tracer_d(n)%scale_factor, &
grid_loc='3112', &
coordinates='TLONG TLAT z_w time' )
call define_tavg_field
(tavg_ADV_TRACER(n), &
'ADV_' /&
&/ trim(tracer_d(n)%short_name),2,&
long_name='Vertically-Integrated ' /&
&/trim(tracer_d(n)%short_name) /&
&/' Advection Tendency', &
units=trim(tracer_d(n)%flux_units), &
scale_factor=tracer_d(n)%scale_factor, &
grid_loc='2110', &
coordinates='TLONG TLAT time' )
end do
call define_tavg_field
(tavg_PV,'PV',3, &
long_name='Potential Vorticity', &
units='1/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field
(tavg_Q,'Q',3, &
long_name='Static Stability (d(rho(p_r))/dz)', &
units='gram/centimeter^4', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field
(tavg_PD,'PD',3, &
long_name='Potential Density Ref to Surface',&
units='gram/centimeter^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field
(tavg_RHOU,'RHOU',3, &
long_name='U times potential density', &
units='g/cm^2/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field
(tavg_RHOV,'RHOV',3, &
long_name='V times potential density', &
units='g/cm^2/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field
(tavg_URHO,'URHO',3, &
long_name='flux of pot density across east face', &
units='g/cm^2/s', grid_loc='3321')
call define_tavg_field
(tavg_VRHO,'VRHO',3, &
long_name='flux of pot density across north face', &
units='g/cm^2/s', grid_loc='3231')
call define_tavg_field
(tavg_WRHO,'WRHO',3, &
long_name='flux of pot density across top face', &
units='g/cm^2/s', grid_loc='3112', &
coordinates='TLONG TLAT z_w time' )
call define_tavg_field
(tavg_PVWM,'PVWM',3, &
long_name='flux of pot vorticity through top face', &
units='cm/s^2', grid_loc='3112', &
coordinates='TLONG TLAT z_w time' )
call define_tavg_field
(tavg_PVWP,'PVWP',3, &
long_name='flux of pot vorticity through bottom face', &
units='cm/s^2', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field
(tavg_UPV,'UPV',3, &
long_name='flux of pot vorticity through east face', &
units='cm/s^2', grid_loc='3321')
call define_tavg_field
(tavg_VPV,'VPV',3, &
long_name='flux of pot vorticity through north face', &
units='cm/s^2', grid_loc='3231')
call define_tavg_field
(tavg_UQ,'UQ',3, &
long_name='flux of Q through east face', &
units='g/cm^3/s', grid_loc='3321')
call define_tavg_field
(tavg_VQ,'VQ',3, &
long_name='flux of Q through north face', &
units='g/cm^3/s', grid_loc='3231')
call flushm
(stdout)
!-----------------------------------------------------------------------
!EOC
end subroutine init_advection
!***********************************************************************
!BOP
! !IROUTINE: tadvect_ctype_to_tadvect_itype
! !INTERFACE:
function tadvect_ctype_to_tadvect_itype(tadvect_ctype) 2
! !DESCRIPTION:
! map an advection ctype to an itype
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character(len=char_len), intent(in) :: &
tadvect_ctype ! character string for tracer advect choice
! !OUTPUT PARAMETERS:
integer(kind=int_kind) :: tadvect_ctype_to_tadvect_itype
!-----------------------------------------------------------------------
select case (tadvect_ctype(1:6))
case ('center')
tadvect_ctype_to_tadvect_itype = tadvect_centered
case ('upwind')
tadvect_ctype_to_tadvect_itype = tadvect_upwind3
case ('lw_lim')
tadvect_ctype_to_tadvect_itype = tadvect_lw_lim
case default
tadvect_ctype_to_tadvect_itype = -1000
end select
!-----------------------------------------------------------------------
!EOC
end function tadvect_ctype_to_tadvect_itype
!***********************************************************************
!BOP
! !IROUTINE: comp_flux_vel_ghost
! !INTERFACE:
subroutine comp_flux_vel_ghost(DH, errorCode) 1,7
! !DESCRIPTION:
! Compute and store tracer flux velocities in ghost cells.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
intent(in) :: &
DH ! change in surface height at T points
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: errorCode
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
type (block) :: &
this_block ! block information for current block
integer (int_kind) :: &
iblock, &! local block number
i,j, &! dummy indices for horizontal directions
k,n ! dummy indices for vertical level, tracer
real (r8), dimension(nx_block,ny_block,nblocks_clinic) :: &
UTE,UTW,VTN,VTS, &! tracer flux velocities across E,W,N,S faces
WTK, &! vertical velocity at bottom of T box
WTKB ! vertical velocity at bottom of T box
!-----------------------------------------------------------------------
errorCode = POP_Success
if (.not. luse_lw_lim) return
do k = 1,km
!$OMP PARALLEL DO PRIVATE(iblock,this_block)
do iblock = 1,nblocks_clinic
this_block = get_block
(blocks_clinic(iblock),iblock)
if (k == 1) then
WTK(:,:,iblock) = DH(:,:,iblock)
else
WTK(:,:,iblock) = WTKB(:,:,iblock)
endif
call comp_flux_vel
(k,UVEL(:,:,:,curtime,iblock),&
VVEL(:,:,:,curtime,iblock),&
WTK(:,:,iblock),&
UTE(:,:,iblock),UTW(:,:,iblock),&
VTN(:,:,iblock),VTS(:,:,iblock),&
WTKB(:,:,iblock),this_block)
end do
!$OMP END PARALLEL DO
call POP_HaloUpdate
(UTE, POP_haloClinic, &
POP_gridHorzLocEFace, &
POP_fieldKindVector, errorCode, &
fillValue = 0.0_POP_r8)
if (errorCode /= POP_Success) then
call POP_ErrorSet
(errorCode, &
'comp_flux_vel_ghost: error updating halo for UTE')
return
endif
call POP_HaloUpdate
(WTKB, POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode, &
fillValue = 0.0_POP_r8)
if (errorCode /= POP_Success) then
call POP_ErrorSet
(errorCode, &
'comp_flux_vel_ghost: error updating halo for WTKB')
return
endif
!$OMP PARALLEL DO PRIVATE(iblock,this_block)
do iblock = 1,nblocks_clinic
this_block = get_block
(blocks_clinic(iblock),iblock)
UTE_jbm2(:,k,iblock) = UTE(:,this_block%jb-2,iblock)
WTKB_jbm2(:,k,iblock) = WTKB(:,this_block%jb-2,iblock)
WTKB_jep2(:,k,iblock) = WTKB(:,this_block%je+2,iblock)
WTKB_ibm2(:,k,iblock) = WTKB(this_block%ib-2,:,iblock)
WTKB_iep2(:,k,iblock) = WTKB(this_block%ie+2,:,iblock)
end do
!$OMP END PARALLEL DO
end do
!-----------------------------------------------------------------------
!EOC
end subroutine comp_flux_vel_ghost
!***********************************************************************
!BOP
! !IROUTINE: advu
! !INTERFACE:
subroutine advu(k,LUK,LVK,WUK,UUU,VVV,this_block) 1,20
! !DESCRIPTION:
! This routine computes the $x,y$ components of the advection of
! momentum and metric terms given by:
! \begin{equation}
! L_U(u_x) + u_x u_y k_y - u_y^2 k_x
! \end{equation}
! \begin{equation}
! L_U(u_y) + u_x u_y k_x - u_x^2 k_y
! \end{equation}
! where
! \begin{equation}
! L_U(\alpha) = {1\over{\Delta_y}}
! \delta_x\left[
! \overline{\left(\overline{\Delta_y u_x}^y\right)}^{xy}
! \overline{\alpha}^x\right] +
! {1\over{\Delta_x}}
! \delta_y\left[
! \overline{\left(\overline{\Delta_x u_y}^x\right)}^{xy}
! \overline{\alpha}^y\right]
! + \delta_z(w^U\overline{\alpha}^z),
! \end{equation}
! \begin{equation}
! k_x = {1\over{\Delta_y}}\delta_x \Delta_y
! \end{equation}
! and
! \begin{equation}
! k_y = {1\over{\Delta_x}}\delta_y \Delta_x
! \end{equation}
!
! Comments:
! \begin{itemize}
! \item this routine must be called successively with k = 1,2,3,...
!
! \item the vertical velocity $w^U$ in U columns is determined by
! integrating the continuity equation $L(1) = 0$ from the surface
! down to level k. In the rigid-lid formulation, the integration
! starts with $w^U = 0$ at the surface. In the free-surface
! formulation, the integration starts with $w^U =$ the area-weighted
! average of $\partial\eta/\partial t$ at surrounding T points
! ($\partial\eta/\partial t$ is the time change of the surface
! height, and satisfies the barotropic continuity equation
! $\partial\eta/\partial t + \nabla\cdot H {\rm\bf U}=q_w$
! where ${\rm\bf U}$ is the barotropic velocity and $q_w$
! is the surface fresh water flux.)
! \end{itemize}
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! depth level index
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UUU,VVV ! U,V velocity for this block
! at the current time
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(inout) :: &
WUK ! on input: flux velocity at top of U box
! on output: flux velocity at bottom of U box
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
LUK, &! advection of U-momentum
LVK ! advection of V-momentum
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,n, &! loop indices
ibeg,iend, &! beginning and ending index of
jbeg,jend, &! physical domain
bid, &! local block index
ns1, &! dummy stream id number
ns2 ! dummy stream id number
integer (int_kind) :: &
tavg_UEU_stream, &! stream in which UEU is defined
tavg_UEV_stream, &! stream in which UEV is defined
tavg_VNU_stream, &! stream in which VNU is defined
tavg_VNV_stream, &! stream in which VNV is defined
tavg_WTU_stream, &! stream in which WTU is defined
tavg_WTV_stream ! stream in which WTV is defined
real (r8) :: &
cc ! scalar central weight
real (r8), dimension(nx_block,ny_block) :: &
UUE,UUW,VUN,VUS, &! momen flux velocities across E,W,N,S faces
FVN,FUE, &! stencil coeffs used by tavg
WORK, &! local temp space
WUKB ! vert velocity at bottom of level k U box
logical (log_kind) :: &
ltavg_requested_local1, &! temporary logical
ltavg_requested_local2, &! temporary logical
ltavg_on_local1, &! temporary logical
ltavg_on_local2 ! temporary logical
!-----------------------------------------------------------------------
!
! advection fluxes for U box (4-point averages of adv fluxes
! across T boxes).
!
!-----------------------------------------------------------------------
bid = this_block%local_id
call timer_start
(timer_advu, block_id=bid)
ibeg = this_block%ib
iend = this_block%ie
jbeg = this_block%jb
jend = this_block%je
UUW = c0
UUE = c0
VUN = c0
VUS = c0
if (partial_bottom_cells) then
do j=jbeg-1,jend+1
do i=ibeg-1,iend+1
UUW(i,j) = p25 *(UUU(i ,j ,k)*DYU(i ,j ,bid)* &
DZU(i ,j ,k,bid) + &
UUU(i-1,j ,k)*DYU(i-1,j ,bid)* &
DZU(i-1,j ,k,bid))+ &
p125*(UUU(i ,j-1,k)*DYU(i ,j-1 ,bid)* &
DZU(i ,j-1,k,bid) + &
UUU(i-1,j-1,k)*DYU(i-1,j-1 ,bid)* &
DZU(i-1,j-1,k,bid) + &
UUU(i ,j+1,k)*DYU(i ,j+1 ,bid)* &
DZU(i ,j+1,k,bid) + &
UUU(i-1,j+1,k)*DYU(i-1,j+1 ,bid)* &
DZU(i-1,j+1,k,bid))
UUE(i,j) = p25 *(UUU(i+1,j ,k)*DYU(i+1,j ,bid)* &
DZU(i+1,j ,k,bid) + &
UUU(i ,j ,k)*DYU(i ,j ,bid)* &
DZU(i ,j ,k,bid))+ &
p125*(UUU(i+1,j-1,k)*DYU(i+1,j-1 ,bid)* &
DZU(i+1,j-1,k,bid) + &
UUU(i ,j-1,k)*DYU(i ,j-1 ,bid)* &
DZU(i ,j-1,k,bid) + &
UUU(i+1,j+1,k)*DYU(i+1,j+1 ,bid)* &
DZU(i+1,j+1,k,bid) + &
UUU(i ,j+1,k)*DYU(i ,j+1 ,bid)* &
DZU(i ,j+1,k,bid))
VUS(i,j) = p25* (VVV(i ,j ,k)*DXU(i ,j ,bid)* &
DZU(i ,j ,k,bid) + &
VVV(i ,j-1,k)*DXU(i ,j-1 ,bid)* &
DZU(i ,j-1,k,bid))+ &
p125*(VVV(i-1,j ,k)*DXU(i-1,j ,bid)* &
DZU(i-1,j ,k,bid) + &
VVV(i-1,j-1,k)*DXU(i-1,j-1 ,bid)* &
DZU(i-1,j-1,k,bid) + &
VVV(i+1,j ,k)*DXU(i+1,j ,bid)* &
DZU(i+1,j ,k,bid) + &
VVV(i+1,j-1,k)*DXU(i+1,j-1 ,bid)* &
DZU(i+1,j-1,k,bid))
VUN(i,j) = p25* (VVV(i ,j+1,k)*DXU(i ,j+1 ,bid)* &
DZU(i ,j+1,k,bid) + &
VVV(i ,j ,k)*DXU(i ,j ,bid)* &
DZU(i ,j ,k,bid))+ &
p125*(VVV(i-1,j+1,k)*DXU(i-1,j+1 ,bid)* &
DZU(i-1,j+1,k,bid) + &
VVV(i-1,j ,k)*DXU(i-1,j ,bid)* &
DZU(i-1,j ,k,bid) + &
VVV(i+1,j+1,k)*DXU(i+1,j+1 ,bid)* &
DZU(i+1,j+1,k,bid) + &
VVV(i+1,j ,k)*DXU(i+1,j ,bid)* &
DZU(i+1,j ,k,bid))
end do
end do
else
do j=jbeg-1,jend+1
do i=ibeg-1,iend+1
UUW(i,j) = p25 *(UUU(i ,j ,k)*DYU(i ,j ,bid) + &
UUU(i-1,j ,k)*DYU(i-1,j ,bid))+ &
p125*(UUU(i ,j-1,k)*DYU(i ,j-1,bid) + &
UUU(i-1,j-1,k)*DYU(i-1,j-1,bid) + &
UUU(i ,j+1,k)*DYU(i ,j+1,bid) + &
UUU(i-1,j+1,k)*DYU(i-1,j+1,bid))
UUE(i,j) = p25 *(UUU(i+1,j ,k)*DYU(i+1,j ,bid) + &
UUU(i ,j ,k)*DYU(i ,j ,bid))+ &
p125*(UUU(i+1,j-1,k)*DYU(i+1,j-1,bid) + &
UUU(i ,j-1,k)*DYU(i ,j-1,bid) + &
UUU(i+1,j+1,k)*DYU(i+1,j+1,bid) + &
UUU(i ,j+1,k)*DYU(i ,j+1,bid))
VUS(i,j) = p25* (VVV(i ,j ,k)*DXU(i ,j ,bid) + &
VVV(i ,j-1,k)*DXU(i ,j-1,bid))+ &
p125*(VVV(i-1,j ,k)*DXU(i-1,j ,bid) + &
VVV(i-1,j-1,k)*DXU(i-1,j-1,bid) + &
VVV(i+1,j ,k)*DXU(i+1,j ,bid) + &
VVV(i+1,j-1,k)*DXU(i+1,j-1,bid))
VUN(i,j) = p25* (VVV(i ,j+1,k)*DXU(i ,j+1,bid) + &
VVV(i ,j ,k)*DXU(i ,j ,bid))+ &
p125*(VVV(i-1,j+1,k)*DXU(i-1,j+1,bid) + &
VVV(i-1,j ,k)*DXU(i-1,j ,bid) + &
VVV(i+1,j+1,k)*DXU(i+1,j+1,bid) + &
VVV(i+1,j ,k)*DXU(i+1,j ,bid))
end do
end do
endif ! partial bottom cells
!-----------------------------------------------------------------------
!
! calculate vertical velocity at bottom of kth level
! (vertical velocity is nonzero at bottom of U columns
! if topography is varying)
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
WUKB = WUK + (VUN - VUS + UUE - UUW)*UAREA_R(:,:,bid)
else
WUKB = WUK + c2dz(k)*p5*(VUN - VUS + UUE - UUW)* &
UAREA_R(:,:,bid)
endif
!#if drifter_particles
!#if drifter_particles != 2
!-----------------------------------------------------------------------
!
! save the vertical velocity for advecting drifters
!
!-----------------------------------------------------------------------
! WVEL(:,:,k) = WUKB
!#endif
!#endif
!-----------------------------------------------------------------------
!
! advect momentum
!
! horizontal advection
!
!-----------------------------------------------------------------------
LUK = c0
LVK = c0
if (partial_bottom_cells) then
do j=jbeg,jend
do i=ibeg,iend
cc = VUS(i,j+1) - VUS(i,j) + UUW(i+1,j) - UUW(i,j)
LUK(i,j) = p5*( cc*UUU(i ,j ,k) + &
VUS(i,j+1)*UUU(i ,j+1,k) - &
VUS(i,j )*UUU(i ,j-1,k) + &
UUW(i+1,j)*UUU(i+1,j ,k) - &
UUW(i ,j)*UUU(i-1,j ,k))* &
UAREA_R(i,j,bid)/ &
DZU(i,j,k,bid)
LVK(i,j) = p5*( cc*VVV(i ,j ,k) + &
VUS(i,j+1)*VVV(i ,j+1,k) - &
VUS(i,j )*VVV(i ,j-1,k) + &
UUW(i+1,j)*VVV(i+1,j ,k) - &
UUW(i ,j)*VVV(i-1,j ,k))* &
UAREA_R(i,j,bid)/ &
DZU(i,j,k,bid)
end do
end do
else ! no partial bottom cells
do j=jbeg,jend
do i=ibeg,iend
cc = VUS(i,j+1) - VUS(i,j) + UUW(i+1,j) - UUW(i,j)
LUK(i,j) = p5*( cc*UUU(i ,j ,k) + &
VUS(i,j+1)*UUU(i ,j+1,k) - &
VUS(i,j )*UUU(i ,j-1,k) + &
UUW(i+1,j)*UUU(i+1,j ,k) - &
UUW(i ,j)*UUU(i-1,j ,k))* &
UAREA_R(i,j,bid)
LVK(i,j) = p5*( cc*VVV(i ,j ,k) + &
VUS(i,j+1)*VVV(i ,j+1,k) - &
VUS(i,j )*VVV(i ,j-1,k) + &
UUW(i+1,j)*VVV(i+1,j ,k) - &
UUW(i ,j)*VVV(i-1,j ,k))* &
UAREA_R(i,j,bid)
end do
end do
endif ! partial bottom cells
!-----------------------------------------------------------------------
!
! vertical advection through top of U box
!
!-----------------------------------------------------------------------
if (k == 1) then
LUK = LUK + dzr(k)*WUK*UUU(:,:,k)
LVK = LVK + dzr(k)*WUK*VVV(:,:,k)
else
if (partial_bottom_cells) then
LUK = LUK + p5/DZU(:,:,k,bid)*WUK*(UUU(:,:,k-1) + &
UUU(:,:,k))
LVK = LVK + p5/DZU(:,:,k,bid)*WUK*(VVV(:,:,k-1) + &
VVV(:,:,k))
else
LUK = LUK + dz2r(k)*WUK*(UUU(:,:,k-1) + UUU(:,:,k))
LVK = LVK + dz2r(k)*WUK*(VVV(:,:,k-1) + VVV(:,:,k))
endif
endif
!-----------------------------------------------------------------------
!
! vertical advection through bottom of U box
! for k=km, UUU(k+1) is not defined, but WUKB=0
!
!-----------------------------------------------------------------------
if (k < km) then
if (partial_bottom_cells) then
LUK = LUK - p5/DZU(:,:,k,bid)*WUKB*(UUU(:,:,k) + &
UUU(:,:,k+1))
LVK = LVK - p5/DZU(:,:,k,bid)*WUKB*(VVV(:,:,k) + &
VVV(:,:,k+1))
else
LUK = LUK - dz2r(k)*WUKB*(UUU(:,:,k) + UUU(:,:,k+1))
LVK = LVK - dz2r(k)*WUKB*(VVV(:,:,k) + VVV(:,:,k+1))
endif
endif
!-----------------------------------------------------------------------
!
! add metric terms, and zero fields at land points
!
!-----------------------------------------------------------------------
do j=jbeg,jend
do i=ibeg,iend
if (k <= KMU(i,j,bid)) then
LUK(i,j) = LUK(i,j) + UUU(i,j,k)*VVV(i,j,k)*KYU(i,j,bid) - &
VVV(i,j,k)**2*KXU(i,j,bid)
LVK(i,j) = LVK(i,j) + UUU(i,j,k)*VVV(i,j,k)*KXU(i,j,bid) - &
UUU(i,j,k)**2*KYU(i,j,bid)
else
LUK(i,j) = c0
LVK(i,j) = c0
endif
end do
end do
!-----------------------------------------------------------------------
!
! accumulate time average if necessary
!
!-----------------------------------------------------------------------
if (mix_pass /= 1) then
if (tavg_requested
(tavg_UEU)) then
tavg_UEU_stream = tavg_in_which_stream
(tavg_UEU)
ltavg_on_local1 = ltavg_on(tavg_UEU_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if (tavg_requested
(tavg_UEV)) then
tavg_UEV_stream = tavg_in_which_stream
(tavg_UEV)
ltavg_on_local2 = ltavg_on(tavg_UEV_stream)
ltavg_requested_local2 = .true.
else
ltavg_on_local2 = .false.
ltavg_requested_local2 = .false.
endif
if ( (ltavg_on_local1 .and. ltavg_requested_local1 ) .or. &
(ltavg_on_local2 .and. ltavg_requested_local2 ) ) then
if (partial_bottom_cells) then
FUE = UUE*p5*UAREA_R(:,:,bid)/DZU(:,:,k,bid)
else
FUE = UUE*p5*UAREA_R(:,:,bid)
endif
endif ! ltavg_on
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
do j=jbeg,jend
do i=ibeg,iend
WORK(i,j) = FUE(i,j)*(UUU(i,j,k) + UUU(i+1,j,k))
end do
end do
call accumulate_tavg_field
(WORK,tavg_UEU,bid,k)
endif ! ltavg_on
if ( ltavg_on_local2 .and. ltavg_requested_local2 ) then
do j=jbeg,jend
do i=ibeg,iend
WORK(i,j) = FUE(i,j)*(VVV(i ,j,k) + VVV(i+1,j,k))
end do
end do
call accumulate_tavg_field
(WORK,tavg_UEV,bid,k)
endif ! ltavg_on
if (tavg_requested
(tavg_VNU)) then
tavg_VNU_stream = tavg_in_which_stream
(tavg_VNU)
ltavg_on_local1 = ltavg_on(tavg_VNU_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if (tavg_requested
(tavg_VNV)) then
tavg_VNV_stream = tavg_in_which_stream
(tavg_VNV)
ltavg_on_local2 = ltavg_on(tavg_VNV_stream)
ltavg_requested_local2 = .true.
else
ltavg_on_local2 = .false.
ltavg_requested_local2 = .false.
endif
if ( (ltavg_on_local1 .and. ltavg_requested_local1 ) .or. &
(ltavg_on_local2 .and. ltavg_requested_local2 ) ) then
if (partial_bottom_cells) then
FVN = VUN*p5*UAREA_R(:,:,bid)/DZU(:,:,k,bid)
else
FVN = VUN*p5*UAREA_R(:,:,bid)
endif
endif ! ltavg_on
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
do j=jbeg,jend
do i=ibeg,iend
WORK(i,j) = FVN(i,j)*(UUU(i,j ,k) + UUU(i,j+1,k))
end do
end do
call accumulate_tavg_field
(WORK,tavg_VNU,bid,k)
endif ! ltavg_on
if ( ltavg_on_local2 .and. ltavg_requested_local2 ) then
do j=jbeg,jend
do i=ibeg,iend
WORK(i,j) = FVN(i,j)*(VVV(i,j ,k) + VVV(i,j+1,k))
end do
end do
call accumulate_tavg_field
(WORK,tavg_VNV,bid,k)
endif ! ltavg_on
if (tavg_requested
(tavg_WTU)) then
tavg_WTU_stream = tavg_in_which_stream
(tavg_WTU)
ltavg_on_local1 = ltavg_on(tavg_WTU_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
if (k == 1) then
WORK = dzr(k)*WUK*UUU(:,:,k)
else
WORK = dz2r(k)*WUK*(UUU(:,:,k) + UUU(:,:,k-1))
endif
call accumulate_tavg_field
(WORK,tavg_WTU,bid,k)
endif ! ltavg_on
if (tavg_requested
(tavg_WTV)) then
tavg_WTV_stream = tavg_in_which_stream
(tavg_WTV)
ltavg_on_local1 = ltavg_on(tavg_WTV_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
if (k == 1) then
WORK = dzr(k)*WUK*VVV(:,:,k)
else
WORK = dz2r(k)*WUK*(VVV(:,:,k) + VVV(:,:,k-1))
endif
call accumulate_tavg_field
(WORK,tavg_WTV,bid,k)
endif ! ltavg_on
endif ! mix_pass
!-----------------------------------------------------------------------
!
! reset W for next k level
!
!-----------------------------------------------------------------------
WUK = WUKB ! bottom value becomes top value for next pass
call timer_stop
(timer_advu, block_id=bid)
!-----------------------------------------------------------------------
!EOC
end subroutine advu
!***********************************************************************
!BOP
! !IROUTINE: advt
! !INTERFACE:
subroutine advt(k,LTK,WTK,TMIX,TRCR,UUU,VVV,this_block) 1,100
! !DESCRIPTION:
! Advection of tracers - this routine actually computes only
! vertical velocities and then calls individual routines based
! on the type of advection chosen
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k ! depth level index
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TMIX, &! tracers for this block at mix time
TRCR ! tracers for this block at current time
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UUU, VVV ! U,V for this block at current time
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(inout) :: &
WTK ! on input flux velocity at top of T box
! on output flux velocity at bottom of T box
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(out) :: &
LTK ! returned as L(T) for nth tracer at level k
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local or common variables:
!
!-----------------------------------------------------------------------
logical (log_kind), dimension(nt) :: &
tr_mask ! which tracers are using a particular advection scheme
logical (log_kind) :: &
ltavg_requested_local1, &! temporary logical
ltavg_requested_local2, &! temporary logical
ltavg_on_local1, &! temporary logical
ltavg_on_local2 ! temporary logical
integer (int_kind) :: &
i,j,n, &! dummy loop indices
ib, ie, jb, je, &! domain limits
bid, &! local block address
tavg_stream ! dummy stream id
real (r8), dimension(nx_block,ny_block) :: &
UTE,UTW,VTN,VTS, &! tracer flux velocities across E,W,N,S faces
FVN,FUE, &! north and east advective stencil weights
WTKB, &! vertical velocity at bottom of T box
RHOK1, &! pot density at k relative to k=1
RHOK1M, &! pot density at k-1
WORK, &! local temp space
WORK1,WORK2, &! local temp space
WORK3,WORK4
real (r8), dimension(nx_block,ny_block,nt) :: &
TRACER_E, &! tracer value on east face
TRACER_N, &! tracer value on north face
FLUX_T ! tracer tendency due to flux through top face, volume normalized
!-----------------------------------------------------------------------
!
! advection fluxes for T box.
!
!-----------------------------------------------------------------------
bid = this_block%local_id
jb = this_block%jb
je = this_block%je
ib = this_block%ib
ie = this_block%ie
call timer_start
(timer_advt, block_id=bid)
!-----------------------------------------------------------------------
!
! calculate vertical velocity at bottom of kth level
! (vertical velocity is zero at bottom of T columns)
!
! use flux velocities from k+1 of previous iteration if available
!
!-----------------------------------------------------------------------
if (k == 1 .or. .not. luse_lw_lim) then
call comp_flux_vel
(k,UUU,VVV,WTK, &
UTE,UTW,VTN,VTS,WTKB,this_block)
else
UTE = FLUX_VEL_prev(:,:,1,bid)
UTW = FLUX_VEL_prev(:,:,2,bid)
VTN = FLUX_VEL_prev(:,:,3,bid)
VTS = FLUX_VEL_prev(:,:,4,bid)
WTKB = FLUX_VEL_prev(:,:,5,bid)
end if
!-----------------------------------------------------------------------
!
! advect tracers
!
!-----------------------------------------------------------------------
LTK = c0
!*** lw_lim advection
if (luse_lw_lim) then
tr_mask = tadvect_itype == tadvect_lw_lim
!-----------------------------------------------------------------------
! compute WTKB at bottom of level k+1 T box and store flux velocities
! WTKBp1 is stored in FLUX_VEL_prev(:,:,5,bid)
!-----------------------------------------------------------------------
call comp_flux_vel
(k+1,UUU,VVV,WTKB, &
FLUX_VEL_prev(:,:,1,bid),&
FLUX_VEL_prev(:,:,2,bid),&
FLUX_VEL_prev(:,:,3,bid),&
FLUX_VEL_prev(:,:,4,bid),&
FLUX_VEL_prev(:,:,5,bid),this_block)
call advt_lw_lim
(k,LTK,TMIX,WTK,WTKB,FLUX_VEL_prev(:,:,5,bid),&
UTE,UTW,VTN,VTS,TRACER_E,TRACER_N,FLUX_T,&
tr_mask,this_block)
end if !luse_lw_lim
!*** upwind3 advection
if (luse_upwind3) then
tr_mask = tadvect_itype == tadvect_upwind3
call advt_upwind3
(k,LTK,TRCR,WTK,WTKB,UTE,UTW,VTN,VTS, &
TRACER_E,TRACER_N,FLUX_T,tr_mask,this_block)
end if
!*** centered advection
if (luse_centered) then
tr_mask = tadvect_itype == tadvect_centered
call advt_centered
(k,LTK,TRCR,WTK,WTKB,UTE,VTN,tr_mask,this_block)
end if
call timer_stop
(timer_advt, block_id=bid)
!-----------------------------------------------------------------------
!
! compute diagnostics if necessary
!
!-----------------------------------------------------------------------
if (mix_pass /= 1 ) then
if (partial_bottom_cells) then
FVN = p5*VTN*TAREA_R(:,:,bid)/DZT(:,:,k,bid)
FUE = p5*UTE*TAREA_R(:,:,bid)/DZT(:,:,k,bid)
else
FVN = p5*VTN*TAREA_R(:,:,bid)
FUE = p5*UTE*TAREA_R(:,:,bid)
endif
if (tavg_requested
(tavg_WVEL)) then
tavg_stream = tavg_in_which_stream
(tavg_WVEL)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(WTK,tavg_WVEL,bid,k)
endif
if (tavg_requested
(tavg_WVEL2)) then
tavg_stream = tavg_in_which_stream
(tavg_WVEL2)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(WTK**2,tavg_WVEL2,bid,k)
endif
do n=1,nt
if (tadvect_itype(n) == tadvect_centered) then
if (tavg_requested
(tavg_UE_TRACER(n))) then
tavg_stream = tavg_in_which_stream
(tavg_UE_TRACER(n))
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FUE*( TRCR(:,:,k,n) + &
eoshift(TRCR(:,:,k,n),dim=1,shift=1))
call accumulate_tavg_field
(WORK,tavg_UE_TRACER(n),bid,k)
endif
if (tavg_requested
(tavg_VN_TRACER(n))) then
tavg_stream = tavg_in_which_stream
(tavg_VN_TRACER(n))
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FVN*( TRCR(:,:,k,n) + &
eoshift(TRCR(:,:,k,n),dim=2,shift=1))
call accumulate_tavg_field
(WORK,tavg_VN_TRACER(n),bid,k)
endif
if (tavg_requested
(tavg_WT_TRACER(n))) then
tavg_stream = tavg_in_which_stream
(tavg_WT_TRACER(n))
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
if (k == 1) then
if (sfc_layer_type /= sfc_layer_varthick) then
WORK = dzr(k)*WTK*TRCR(:,:,k,n)
else
WORK = c0
endif
else
WORK = dz2r(k)*WTK*(TRCR(:,:,k ,n) &
+ TRCR(:,:,k-1,n))
endif
call accumulate_tavg_field
(WORK,tavg_WT_TRACER(n),bid,k)
endif
else
if (tavg_requested
(tavg_UE_TRACER(n))) then
tavg_stream = tavg_in_which_stream
(tavg_UE_TRACER(n))
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = c2*FUE*TRACER_E(:,:,n)
call accumulate_tavg_field
(WORK,tavg_UE_TRACER(n),bid,k)
endif
if (tavg_requested
(tavg_VN_TRACER(n))) then
tavg_stream = tavg_in_which_stream
(tavg_VN_TRACER(n))
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = c2*FVN*TRACER_N(:,:,n)
call accumulate_tavg_field
(WORK,tavg_VN_TRACER(n),bid,k)
endif
if (tavg_requested
(tavg_WT_TRACER(n))) then
tavg_stream = tavg_in_which_stream
(tavg_WT_TRACER(n))
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(FLUX_T(:,:,n),tavg_WT_TRACER(n),bid,k)
endif
endif
if (tavg_requested
(tavg_ADV_TRACER(n))) then
tavg_stream = tavg_in_which_stream
(tavg_ADV_TRACER(n))
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = c0
do j=jb,je
do i=ib,ie
if (k <= KMT(i,j,bid)) then
WORK(i,j) = -dz(k)*LTK(i,j,n)
endif
end do
end do
call accumulate_tavg_field
(WORK,tavg_ADV_TRACER(n),bid,k)
endif
enddo
!***
!*** potential density referenced to k=1 needed for a variety of
!*** tavg fields
!***
if (tavg_requested
(tavg_PD) .or. &
tavg_requested
(tavg_RHOU) .or. &
tavg_requested
(tavg_RHOV) .or. &
tavg_requested
(tavg_URHO) .or. &
tavg_requested
(tavg_VRHO) .or. &
tavg_requested
(tavg_WRHO)) then
call state
(k,1,TRCR(:,:,k,1), &
TRCR(:,:,k,2), this_block, &
RHOFULL=RHOK1)
if (k == 1) then
RHOK1M = RHOK1
else
call state
(k-1,1,TRCR(:,:,k,1), &
TRCR(:,:,k,2), this_block, &
RHOFULL=RHOK1M)
endif
endif
if (tavg_requested
(tavg_PD)) then
tavg_stream = tavg_in_which_stream
(tavg_PD)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(RHOK1,tavg_PD,bid,k)
endif
if (tavg_requested
(tavg_URHO)) then
tavg_stream = tavg_in_which_stream
(tavg_URHO)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FUE*(RHOK1 + eoshift(RHOK1,dim=1,shift=1))
call accumulate_tavg_field
(WORK,tavg_URHO,bid,k)
endif
if (tavg_requested
(tavg_VRHO)) then
tavg_stream = tavg_in_which_stream
(tavg_VRHO)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FVN*(RHOK1 + eoshift(RHOK1,dim=2,shift=1))
call accumulate_tavg_field
(WORK,tavg_VRHO,bid,k)
endif
if (tavg_requested
(tavg_WRHO)) then
tavg_stream = tavg_in_which_stream
(tavg_WRHO)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = dz2r(k)*WTK*(RHOK1 + RHOK1M)
call accumulate_tavg_field
(WORK,tavg_WRHO,bid,k)
endif
if (tavg_requested
(tavg_RHOU)) then
tavg_stream = tavg_in_which_stream
(tavg_RHOU)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call ugrid_to_tgrid
(WORK1,UUU(:,:,k),bid)
WORK = RHOK1*WORK1
call accumulate_tavg_field
(WORK,tavg_RHOU,bid,k)
endif
if (tavg_requested
(tavg_RHOV)) then
tavg_stream = tavg_in_which_stream
(tavg_RHOV)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call ugrid_to_tgrid
(WORK1,VVV(:,:,k),bid)
WORK = RHOK1*WORK1
call accumulate_tavg_field
(WORK,tavg_RHOV,bid,k)
endif
!***
!*** vertical density gradient and potential vorticity
!***
if (tavg_requested
(tavg_Q) .or. &
tavg_requested
(tavg_PV) .or. &
tavg_requested
(tavg_PVWM) .or. &
tavg_requested
(tavg_PVWP) .or. &
tavg_requested
(tavg_UPV) .or. &
tavg_requested
(tavg_VPV) .or. &
tavg_requested
(tavg_UQ) .or. &
tavg_requested
(tavg_VQ) ) then
call state
(k,k,TRCR(:,:,k,1) &
,TRCR(:,:,k,2), this_block, &
RHOOUT=WORK)
if (k == 1 ) then
WORK3 = WORK
else
call state
(k-1,k,TRCR(:,:,k-1,1) &
,TRCR(:,:,k-1,2), this_block, &
RHOOUT=WORK3)
WORK3 = p5*(WORK3 + WORK)
endif
if (k == km) then
WORK4 = WORK
else
call state
(k+1,k,TRCR(:,:,k+1,1) &
,TRCR(:,:,k+1,2), this_block, &
RHOOUT=WORK4)
do j=jb,je
do i=ib,ie
if (k /= KMT(i,j,bid)) then
WORK4(i,j) = p5*(WORK4(i,j) + WORK(i,j))
else
WORK4(i,j) = WORK(i,j)
endif
end do
end do
endif
call zero_ghost_cells
(this_block,WORK1)
do j=jb,je
do i=ib,ie
if (k <= KMT(i,j,bid)) then
WORK1(i,j) = (WORK3(i,j) - WORK4(i,j))*dzr(k) ! drho/dz
else
WORK1(i,j) = c0
endif
end do
end do
call zcurl
(k,WORK3,UUU(:,:,k),VVV(:,:,k),this_block)
WORK2 = WORK1*(WORK3*TAREA_R(:,:,bid) + &
FCORT(:,:,bid)) ! PV = pot vorticity
endif
if (tavg_requested
(tavg_Q)) then
tavg_stream = tavg_in_which_stream
(tavg_Q)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(WORK1,tavg_Q,bid,k)
endif
if (tavg_requested
(tavg_PV)) then
tavg_stream = tavg_in_which_stream
(tavg_PV)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(WORK2,tavg_PV,bid,k)
endif
!***
!*** advection of potential vorticity, Q
!***
if (tavg_requested
(tavg_PVWM)) then
tavg_stream = tavg_in_which_stream
(tavg_PVWM)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(WORK2*WTK,tavg_PVWM,bid,k)
endif
if (tavg_requested
(tavg_PVWP)) then
tavg_stream = tavg_in_which_stream
(tavg_PVWP)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
call accumulate_tavg_field
(WORK2*WTKB,tavg_PVWP,bid,k)
endif
if (tavg_requested
(tavg_UPV)) then
tavg_stream = tavg_in_which_stream
(tavg_UPV)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FUE*(WORK2 + eoshift(WORK2,dim=1,shift=1))
call accumulate_tavg_field
(WORK,tavg_UPV,bid,k)
endif
if (tavg_requested
(tavg_VPV)) then
tavg_stream = tavg_in_which_stream
(tavg_VPV)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FVN*(WORK2 + eoshift(WORK2,dim=2,shift=1))
call accumulate_tavg_field
(WORK,tavg_VPV,bid,k)
endif
if (tavg_requested
(tavg_UQ)) then
tavg_stream = tavg_in_which_stream
(tavg_UQ)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FUE*(WORK1 + eoshift(WORK1,dim=1,shift=1))
call accumulate_tavg_field
(WORK,tavg_UQ,bid,k)
endif
if (tavg_requested
(tavg_VQ)) then
tavg_stream = tavg_in_which_stream
(tavg_VQ)
ltavg_on_local1 = ltavg_on(tavg_stream)
ltavg_requested_local1 = .true.
else
ltavg_on_local1 = .false.
ltavg_requested_local1 = .false.
endif
if ( ltavg_on_local1 .and. ltavg_requested_local1 ) then
WORK = FVN*(WORK1 + eoshift(WORK1,dim=2,shift=1))
call accumulate_tavg_field
(WORK,tavg_VQ,bid,k)
endif
endif ! mix_pass
call cfl_advect
(k,bid,UTE,UTW,VTN,VTS,WTK,WTKB,this_block)
!-----------------------------------------------------------------------
!
! set bottom vertical velocity to top for next k-level
!
!-----------------------------------------------------------------------
WTK = WTKB
!-----------------------------------------------------------------------
!EOC
end subroutine advt
!***********************************************************************
!BOP
! !IROUTINE: comp_flux_vel
! !INTERFACE:
subroutine comp_flux_vel(k,UUU,VVV,WTK, & 3,5
UTE,UTW,VTN,VTS,WTKB,this_block)
! !DESCRIPTION:
! This routine computes the vertical velocity at bottom of a layer.
! It also returns the tracer flux velocities across the lateral faces.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k ! depth level index
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UUU, VVV ! U,V for this block at current time
real (r8), dimension(nx_block,ny_block), intent(in) :: &
WTK ! vertical velocity at top of T box
type (block), intent(in) :: &
this_block ! block information for this block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
UTE,UTW,VTN,VTS, &! tracer flux velocities across E,W,N,S faces
WTKB ! vertical velocity at bottom of T box
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j, &! tracer loop index
ibeg,iend, &! loop limits
jbeg,jend, &! loop limits
bid ! local block index
real (r8), dimension(nx_block,ny_block) :: &
FC ! local temp space
!-----------------------------------------------------------------------
call zero_ghost_cells
(this_block,UTE)
call zero_ghost_cells
(this_block,UTW)
call zero_ghost_cells
(this_block,VTN)
call zero_ghost_cells
(this_block,VTS)
if (k > km) then
WTKB = c0
return
endif
jbeg = this_block%jb-1
jend = this_block%je+1
ibeg = this_block%ib-1
iend = this_block%ie+1
if (luse_lw_lim) jend = this_block%je+2
bid = this_block%local_id
if (partial_bottom_cells) then
do j=jbeg,jend
do i=ibeg,iend
UTE(i,j) = 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))
UTW(i,j) = 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))
VTN(i,j) = 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))
VTS(i,j) = 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))
end do
end do
else
do j=jbeg,jend
do i=ibeg,iend
UTE(i,j) = p5*(UUU(i ,j ,k)*DYU(i ,j ,bid) + &
UUU(i ,j-1,k)*DYU(i ,j-1,bid))
UTW(i,j) = p5*(UUU(i-1,j ,k)*DYU(i-1,j ,bid) + &
UUU(i-1,j-1,k)*DYU(i-1,j-1,bid))
VTN(i,j) = p5*(VVV(i ,j ,k)*DXU(i ,j ,bid) + &
VVV(i-1,j ,k)*DXU(i-1,j ,bid))
VTS(i,j) = p5*(VVV(i ,j-1,k)*DXU(i ,j-1,bid) + &
VVV(i-1,j-1,k)*DXU(i-1,j-1,bid))
end do
end do
endif ! partial bottom cells
if (luse_lw_lim) then
UTE(:,this_block%jb-2) = UTE_jbm2(:,k,bid)
UTW(2:nx_block,this_block%jb-2) = &
UTE(1:nx_block-1,this_block%jb-2)
UTE(this_block%ib-2,:) = UTW(this_block%ib-1,:)
VTN(:,this_block%jb-2) = VTS(:,this_block%jb-1)
endif
!-----------------------------------------------------------------------
!
! calculate vertical velocity at bottom of kth level
! (vertical velocity is zero at bottom of T columns)
!
!-----------------------------------------------------------------------
if (k < km) then
FC = (VTN - VTS + UTE - UTW)*TAREA_R(:,:,bid)
if ( overflows_on ) then
WTKB = WTK+dz(k)*FC
call ovf_wtkb_check
(k,WTKB,this_block)
endif
if (partial_bottom_cells) then
WTKB = merge(WTK+FC, c0, k < KMT(:,:,bid))
else
WTKB = merge(WTK+dz(k)*FC, c0, k < KMT(:,:,bid))
endif
if (luse_lw_lim) then
WTKB(:,this_block%jb-2) = WTKB_jbm2(:,k,bid)
WTKB(:,this_block%je+2) = WTKB_jep2(:,k,bid)
WTKB(this_block%ib-2,:) = WTKB_ibm2(:,k,bid)
WTKB(this_block%ie+2,:) = WTKB_iep2(:,k,bid)
endif
else
WTKB = c0
endif
!-----------------------------------------------------------------------
!EOC
end subroutine comp_flux_vel
!***********************************************************************
!BOP
! !IROUTINE: advt_centered
! !INTERFACE:
subroutine advt_centered(k,LTK,TRCR,WTK,WTKB,UTE,VTN,tr_mask,this_block) 1
! !DESCRIPTION:
! This routine computes the tracer advection tendencies using
! a centered differencing for leapfrog advection given by:
! \begin{equation}
! L_T(\phi) = {1\over{\Delta_y}}
! \delta_x\left[
! \overline{\left(\overline{\Delta_y u_x}^y\right)}^{xy}
! \overline{\phi}^x\right] +
! {1\over{\Delta_x}}
! \delta_y\left[
! \overline{\left(\overline{\Delta_x u_y}^x\right)}^{xy}
! \overline{\phi}^y\right]
! + \delta_z(w\overline{\phi}^z),
! \end{equation}
! where $\phi$ is the tracer concentration.
!
! Comments:
! \begin{itemize}
! \item this routine must be called successively with k = 1,2,3,...
!
! \item the vertical velocity $w$ in T columns is determined by
! integrating the continuity equation $L(1) = 0$ from the surface
! down to level k. In the rigid-lid formulation, the integration
! starts with $w = 0$ at the surface. In the free-surface
! formulation, the integration starts with
! $w = \partial\eta/\partial t$
! ($\partial\eta/\partial t$ is the time change of the surface
! height, and satisfies the barotropic continuity equation
! $\partial\eta/\partial t + \nabla\cdot H {\rm\bf U}=q_w$
! where ${\rm\bf U}$ is the barotropic velocity and $q_w$
! is the surface fresh water flux.)
! \end{itemize}
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! depth level index
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TRCR ! tracers for this block at current time
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UTE,VTN, &! tracer flux velocities across E,N faces
WTK, &! vert velocity at top of level k T box
WTKB ! vert velocity at bottom of level k T box
logical (log_kind), dimension(nt), intent(in) :: &
tr_mask ! true for tracers using centered
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(inout) :: &
LTK ! returned as L(T) for nth tracer at level k
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,n, &! tracer loop index
bid ! local block index
!-----------------------------------------------------------------------
!
! advect tracers
! horizontal advection
!
!-----------------------------------------------------------------------
bid = this_block%local_id
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
do n=1,nt
if (.not. tr_mask(n)) cycle
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
LTK(i,j,n) = p5*((VTN(i,j)-VTN(i,j-1)+UTE(i,j)-UTE(i-1,j)) &
*TRCR(i ,j ,k,n) + &
VTN(i ,j )*TRCR(i ,j+1,k,n) - &
VTN(i ,j-1)*TRCR(i ,j-1,k,n) + &
UTE(i ,j )*TRCR(i+1,j ,k,n) - &
UTE(i-1,j )*TRCR(i-1,j ,k,n))* &
TAREA_R(i,j,bid)/ &
DZT(i,j,k,bid)
end do
end do
end do
else ! no partial bottom cells
do n=1,nt
if (.not. tr_mask(n)) cycle
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
LTK(i,j,n) = p5*((VTN(i,j)-VTN(i,j-1)+UTE(i,j)-UTE(i-1,j)) &
*TRCR(i ,j ,k,n) + &
VTN(i ,j )*TRCR(i ,j+1,k,n) - &
VTN(i ,j-1)*TRCR(i ,j-1,k,n) + &
UTE(i ,j )*TRCR(i+1,j ,k,n) - &
UTE(i-1,j )*TRCR(i-1,j ,k,n))* &
TAREA_R(i,j,bid)
end do
end do
end do
endif ! partial bottom cells
!-----------------------------------------------------------------------
!
! vertical advection through top and bottom of T box
!
!-----------------------------------------------------------------------
do n = 1,nt
if (.not. tr_mask(n)) cycle
!*** no advection thru surface (k=1) in vs model
if (k == 1) then
if (sfc_layer_type /= sfc_layer_varthick) then
LTK(:,:,n) = LTK(:,:,n) + &
dzr(k)*WTK*TRCR(:,:,k,n)
endif
else
if (partial_bottom_cells) then
LTK(:,:,n) = LTK(:,:,n) + p5/DZT(:,:,k,bid)*WTK* &
(TRCR(:,:,k-1,n) + TRCR(:,:,k ,n))
else
LTK(:,:,n) = LTK(:,:,n) + dz2r(k)*WTK* &
(TRCR(:,:,k-1,n) + TRCR(:,:,k ,n))
endif
endif
!***
!*** for k=km, TRACER(k+1) is not defined, but WTKB=0
!***
if (k < km) then
if (partial_bottom_cells) then
LTK(:,:,n) = LTK(:,:,n) - p5/DZT(:,:,k,bid)*WTKB* &
(TRCR(:,:,k,n) + TRCR(:,:,k+1,n))
else
LTK(:,:,n) = LTK(:,:,n) - dz2r(k)*WTKB* &
(TRCR(:,:,k,n) + TRCR(:,:,k+1,n))
endif
endif
enddo
!-----------------------------------------------------------------------
!EOC
end subroutine advt_centered
!***********************************************************************
!BOP
! !IROUTINE: advt_upwind3
! !INTERFACE:
subroutine advt_upwind3(k,LTK,TRCR,WTK,WTKB,UTE,UTW,VTN,VTS, & 1,1
TRACER_E,TRACER_N,FLUX_T,tr_mask,this_block)
! !DESCRIPTION:
! This routine computes the advection of tracers using a 3rd-order
! upwinding defined by the QUICKEST scheme in
! Leonard, B.P. 1979, {\em Comp. Meth. Applied Math. and Engineering},
! {\bf 19}, 59-98.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! depth level index
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TRCR ! tracers for this block at current time
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UTE,UTW,VTN,VTS, &! tracer flux velocities across E,W,N,S faces
WTK, &! vert velocity at top of level k T box
WTKB ! vert velocity at bottom of level k T box
logical (log_kind), dimension(nt), intent(in) :: &
tr_mask ! true for tracers using upwind3
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(inout) :: &
LTK, &! returned as L(T) for nth tracer at level k
TRACER_E, &! tracer value on east face
TRACER_N, &! tracer value on north face
FLUX_T ! tracer tendency due to flux through top face, volume normalized
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,n, &! loop index
bid ! local block address for this block
real (r8), dimension(nx_block,ny_block) :: &
WORK1, &! local temp space
FVN,FVS,FUE,FUW, &! flux velocities
TPLUS,TMINUS,AZMINUS,DZMINUS
real (r8), dimension(nx_block,ny_block,nt) :: &
AUXB
!-----------------------------------------------------------------------
!
! advect tracers
! horizontal advection
!
!-----------------------------------------------------------------------
bid = this_block%local_id
if (partial_bottom_cells) then
WORK1 = TAREA_R(:,:,bid) / DZT(:,:,k,bid)
FVN = VTN*WORK1
FVS = -VTS*WORK1 ! note sign
FUE = UTE*WORK1
FUW = -UTW*WORK1 ! note sign
else
FVN = VTN*TAREA_R(:,:,bid)
FVS = -VTS*TAREA_R(:,:,bid) ! note sign change
FUE = UTE*TAREA_R(:,:,bid)
FUW = -UTW*TAREA_R(:,:,bid) ! note sign change
endif
call hupw3
(k,LTK,TRCR,FVN,FVS,FUE,FUW,TRACER_E,TRACER_N, &
bid, tr_mask, this_block)
!-----------------------------------------------------------------------
!
! vertical advection through top and bottom of T box
!
!-----------------------------------------------------------------------
do j=1,ny_block
do i=1,nx_block
if ( k < KMT(i,j,bid)-1 ) then
AZMINUS(i,j) = talfzm(k)
DZMINUS(i,j) = tdelzm(k)
else
AZMINUS(i,j) = talfzm(k) + tdelzm(k)
DZMINUS(i,j) = c0
endif
end do
end do
do n = 1,nt
if (.not. tr_mask(n)) cycle
if ( k < km-1 .and. k > 1) then
TPLUS = talfzp(k)*TRCR(:,:,k+1,n) &
+ tbetzp(k)*TRCR(:,:,k ,n) &
+ tgamzp(k)*TRCR(:,:,k-1,n)
TMINUS = AZMINUS *TRCR(:,:,k+1,n) &
+ tbetzm(k)*TRCR(:,:,k ,n) &
+ DZMINUS *TRCR(:,:,k+2,n)
AUXB(:,:,n) = (WTKB-abs(WTKB))*TPLUS + &
(WTKB+abs(WTKB))*TMINUS
else if ( k == 1) then
TPLUS = talfzp(k)*TRCR(:,:,k+1,n) &
+ tbetzp(k)*TRCR(:,:,k ,n)
TMINUS = AZMINUS *TRCR(:,:,k+1,n) &
+ tbetzm(k)*TRCR(:,:,k ,n) &
+ DZMINUS *TRCR(:,:,k+2,n)
AUXB(:,:,n) = (WTKB-abs(WTKB))*TPLUS + &
(WTKB+abs(WTKB))*TMINUS
else if ( k == km-1) then
TPLUS = talfzp(k)*TRCR(:,:,k+1,n) &
+ tbetzp(k)*TRCR(:,:,k ,n) &
+ tgamzp(k)*TRCR(:,:,k-1,n)
TMINUS = AZMINUS *TRCR(:,:,k+1,n) &
+ tbetzm(k)*TRCR(:,:,k ,n)
AUXB(:,:,n) = (WTKB-abs(WTKB))*TPLUS + &
(WTKB+abs(WTKB))*TMINUS
else
AUXB(:,:,n) = c0
endif
!*** no advection thru surface (k=1) in vs model
if (k == 1) then
if (sfc_layer_type /= sfc_layer_varthick) then
FLUX_T(:,:,n) = dzr(k)*WTK*TRCR(:,:,k,n)
LTK(:,:,n) = LTK(:,:,n) &
+ FLUX_T(:,:,n) &
- dz2r(k)*AUXB(:,:,n)
else
FLUX_T(:,:,n) = c0
LTK(:,:,n) = LTK(:,:,n) &
- dz2r(k)*AUXB(:,:,n)
endif
else
if (partial_bottom_cells) then
WORK1 = p5 / DZT(:,:,k,bid)
FLUX_T(:,:,n) = WORK1*AUX(:,:,n,bid)
LTK(:,:,n) = LTK(:,:,n) + WORK1* &
(AUX(:,:,n,bid) - AUXB(:,:,n))
else
FLUX_T(:,:,n) = dz2r(k)*AUX(:,:,n,bid)
LTK(:,:,n) = LTK(:,:,n) + dz2r(k)* &
(AUX(:,:,n,bid) - AUXB(:,:,n))
endif
endif
AUX(:,:,n,bid) = AUXB(:,:,n) ! set bottom to top for next k level
enddo
!-----------------------------------------------------------------------
!EOC
end subroutine advt_upwind3
!***********************************************************************
!BOP
! !IROUTINE: hupw3
! !INTERFACE:
subroutine hupw3(k,XOUT,X,CN,CS,CE,CW,TRACER_E,TRACER_N, & 1,4
bid,tr_mask,this_block)
! !DESCRIPTION:
! This routine is the horizontal stencil operator for 3rd-order upwind
! advection.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k, &! vertical level index
bid ! local block address for this block
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
X ! input tracer array
real (r8), dimension(nx_block,ny_block), intent(in) :: &
CN,CS,CE,CW ! stencil weights based on flux velocities
logical (log_kind), dimension(nt), intent(in) :: &
tr_mask ! true for tracers using upwind3
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(inout) :: &
XOUT, &! output tracer advection
TRACER_E, &! tracer value on east face
TRACER_N ! tracer value on north face
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,n ! dummy loop indices
real (r8) :: &
am, bm, dm, ap, bp, gp, work
!-----------------------------------------------------------------------
!
! 3rd order upwinding stencil
!
!-----------------------------------------------------------------------
do n=1,nt
if (.not. tr_mask(n)) cycle
call zero_ghost_cells
(this_block,TRACER_E(:,:,n))
call zero_ghost_cells
(this_block,TRACER_N(:,:,n))
!-----------------------------------------------------------------------
! Compute grid-x direction contribution
!-----------------------------------------------------------------------
do j=this_block%jb,this_block%je
do i=this_block%ib-1,this_block%ie
if (k <= KMTE(i,j,bid)) then
work = TBETXP(i,j,bid)
ap = TALFXP(i,j,bid)
else
work = TBETXP(i,j,bid) + TALFXP(i,j,bid)
ap = c0
endif
if (k <= KMTW(i,j,bid)) then
bp = work
gp = TGAMXP(i,j,bid)
else
bp = work + TGAMXP(i,j,bid)
gp = c0
endif
if (k <= KMTEE(i,j,bid)) then
am = TALFXM(i,j,bid)
dm = TDELXM(i,j,bid)
else
am = TALFXM(i,j,bid) + TDELXM(i,j,bid)
dm = c0
endif
bm = TBETXM(i,j,bid)
if (CE(i,j) > c0) then
TRACER_E(i,j,n) = ap*X(i+1,j,k,n) + &
bp*X(i,j,k,n) + gp*X(i-1,j,k,n)
else
TRACER_E(i,j,n) = am*X(i+1,j,k,n) + &
bm*X(i,j,k,n) + dm*X(i+2,j,k,n)
endif
end do ! i loop for grid-x contribution
!*** i loop broken up to avoid dependency on previous i iteration
!*** and to allow different lower loop bounds
end do ! j loop for grid-x TRACER_E evaluation
!*** j loop broken up for overflow TRACER_E modification
if ( overflows_on .and. overflows_interactive ) then
call ovf_advt
(k,TRACER_E(:,:,n),TRACER_N(:,:,n),n,this_block, &
CE,CW,CN,CS)
endif
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
XOUT(i,j,n) = CE(i,j)*TRACER_E(i,j,n) + CW(i,j)*TRACER_E(i-1,j,n)
end do ! i loop for grid-x contribution
end do ! j loop for grid-x contribution
!-----------------------------------------------------------------------
! Compute grid-y direction contribution
!-----------------------------------------------------------------------
do j=this_block%jb-1,this_block%je
do i=this_block%ib,this_block%ie
if (k <= KMTN(i,j,bid)) then
work = TBETYP(i,j,bid)
ap = TALFYP(i,j,bid)
else
work = TBETYP(i,j,bid) + TALFYP(i,j,bid)
ap = c0
endif
if (k <= KMTS(i,j,bid)) then
bp = work
gp = TGAMYP(i,j,bid)
else
bp = work + TGAMYP(i,j,bid)
gp = c0
endif
if (k <= KMTNN(i,j,bid)) then
am = TALFYM(i,j,bid)
dm = TDELYM(i,j,bid)
else
am = TALFYM(i,j,bid) + TDELYM(i,j,bid)
dm = c0
endif
bm = TBETYM(i,j,bid)
if (CN(i,j) > c0) then
TRACER_N(i,j,n) = ap*X(i,j+1,k,n) + &
bp*X(i,j,k,n) + gp*X(i,j-1,k,n)
else
TRACER_N(i,j,n) = am*X(i,j+1,k,n) + &
bm*X(i,j,k,n) + dm*X(i,j+2,k,n)
endif
end do ! i loop for grid-y TRACER_N evaluation
end do ! j loop for grid-y TRACER_N evaluation
!*** i,j loops broken for overflow TRACER_N modification
if ( overflows_on .and. overflows_interactive ) then
call ovf_advt
(k,TRACER_E(:,:,n),TRACER_N(:,:,n),n,this_block, &
CE,CW,CN,CS)
endif
do j=this_block%jb-1,this_block%je
do i=this_block%ib,this_block%ie
!*** The formula below is not correct for j==this_block%jb-1
!*** because TRACER_N(:,this_block%jb-2,n) is not computed.
!*** This is OK because XOUT(:,this_block%jb-1,n) is not used.
XOUT(i,j,n) = XOUT(i,j,n) + &
CN(i,j)*TRACER_N(i,j,n) + CS(i,j)*TRACER_N(i,j-1,n)
end do ! i loop for grid-y contribution
end do ! j loop for grid-y contribution
end do ! tracer loop
!-----------------------------------------------------------------------
!EOC
end subroutine hupw3
!***********************************************************************
!BOP
! !IROUTINE: advt_lw_lim
! !INTERFACE:
subroutine advt_lw_lim(k,LTK,TMIX,WTK,WTKB,WTKBp1,UTE,UTW,VTN,VTS, & 1,1
TRACER_E,TRACER_N,FLUX_T,tr_mask,this_block)
! !DESCRIPTION:
! This routine computes the advection of tracers using a 2nd-order
! forward in time scheme with one-dimensional flux limiters.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! depth level index
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TMIX ! tracers for this block at mix time
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UTE,UTW,VTN,VTS, &! tracer flux velocities across E,W,N,S faces
WTK, &! vert velocity at top of level k T box
WTKB, &! vert velocity at bottom of level k T box
WTKBp1 ! vert velocity at bottom of level k+1 T box
logical (log_kind), dimension(nt), intent(in) :: &
tr_mask ! true for tracers using lw_lim
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(inout) :: &
LTK, &! returned as L(T) for nth tracer at level k
TRACER_E, &! tracer value on east face
TRACER_N, &! tracer value on north face
FLUX_T ! tracer tendency due to flux through top face, volume normalized
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,n, &! loop index
bid ! local block address for this block
real (r8) :: &
adv_dt ! timestep for advection
real (r8), dimension(nx_block,ny_block) :: &
WORK1, &! local temp space
WTK_EFF, &! effective velocity at top of cell
FVN,FVS,FUE,FUW, &! flux velocities
UVEL_E_dt, &! dt * velocity on east face
VVEL_N_dt, &! dt * velocity on north face
KMASKE,KMASKN ! land mask on east and north faces
real (r8), dimension(nx_block,ny_block,nt) :: &
AUXB
!-----------------------------------------------------------------------
!
! set up scalars & arrays for call to lw_lim
!
!-----------------------------------------------------------------------
bid = this_block%local_id
adv_dt = c2dtt(k)
if (partial_bottom_cells) then
UVEL_E_dt = adv_dt * UTE * UTE_to_UVEL_E(:,:,k,bid)
VVEL_N_dt = adv_dt * VTN * VTN_to_VVEL_N(:,:,k,bid)
WORK1 = TAREA_R(:,:,bid) / DZT(:,:,k,bid)
FVN = VTN*WORK1
FVS = -VTS*WORK1 ! note sign
FUE = UTE*WORK1
FUW = -UTW*WORK1 ! note sign
else
UVEL_E_dt = adv_dt * UTE * UTE_to_UVEL_E(:,:,1,bid)
VVEL_N_dt = adv_dt * VTN * VTN_to_VVEL_N(:,:,1,bid)
FVN = VTN*TAREA_R(:,:,bid)
FVS = -VTS*TAREA_R(:,:,bid) ! note sign change
FUE = UTE*TAREA_R(:,:,bid)
FUW = -UTW*TAREA_R(:,:,bid) ! note sign change
endif
KMASKE = merge(c1, c0, k <= KMT(:,:,bid) .and. k <= KMTE(:,:,bid))
KMASKN = merge(c1, c0, k <= KMT(:,:,bid) .and. k <= KMTN(:,:,bid))
!-----------------------------------------------------------------------
!
! set flux at surface (k=1)
! store top of cell volume averaged flux
!
!-----------------------------------------------------------------------
if (k == 1 .and. sfc_layer_type == sfc_layer_varthick) then
WTK_EFF = c0
else
WTK_EFF = WTK
endif
if (partial_bottom_cells) then
do n = 1,nt
if (.not. tr_mask(n)) cycle
if (k == 1) AUX(:,:,n,bid) = WTK_EFF * TMIX(:,:,k,n)
FLUX_T(:,:,n) = AUX(:,:,n,bid) / DZT(:,:,k,bid)
enddo
else
do n = 1,nt
if (.not. tr_mask(n)) cycle
if (k == 1) AUX(:,:,n,bid) = WTK_EFF * TMIX(:,:,k,n)
FLUX_T(:,:,n) = AUX(:,:,n,bid)*dzr(k)
enddo
endif
call lw_lim
(k,adv_dt,UVEL_E_dt,VVEL_N_dt,LTK,TMIX, &
WTK_EFF,WTKB,WTKBp1,AUXB, &
FVN,FVS,FUE,FUW,KMASKE,KMASKN, &
TRACER_E,TRACER_N, &
bid, tr_mask, this_block)
!-----------------------------------------------------------------------
!
! copy bottom flux to top for next k level
!
!-----------------------------------------------------------------------
do n = 1,nt
if (tr_mask(n)) AUX(:,:,n,bid) = AUXB(:,:,n)
enddo
!-----------------------------------------------------------------------
!EOC
end subroutine advt_lw_lim
!***********************************************************************
!BOP
! !IROUTINE: lw_lim
! !INTERFACE:
subroutine lw_lim(k,adv_dt,UVEL_E_dt,VVEL_N_dt,XOUT,X, & 1,5
WTK,WTKB,WTKBp1,AUXB, &
CN,CS,CE,CW,KMASKE,KMASKN, &
TRACER_E,TRACER_N,bid,tr_mask,this_block)
! !DESCRIPTION:
! This routine computes the advection of tracers using a 2nd-order
! forward in time scheme with one-dimensional flux limiters.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k, &! vertical level index
bid ! local block address for this block
real (r8), intent(in) :: &
adv_dt ! advection timestep
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
X ! input tracer array
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UVEL_E_dt, &! dt * velocity on east face
VVEL_N_dt, &! dt * velocity on north face
WTK,WTKB,WTKBp1, &! vert velocity at top of T box for levels k,k+1,k+2
CN,CS,CE,CW, &! stencil weights based on flux velocities
KMASKE,KMASKN ! land mask on east and north faces
logical (log_kind), dimension(nt), intent(in) :: &
tr_mask ! true for tracers using lw_lim
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(inout) :: &
AUXB, &! tracer flux through bottom of cell
XOUT, &! output tracer advection
TRACER_E, &! tracer value on east face
TRACER_N ! tracer value on north face
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
nu ,&! i/o unit number
reclength ,&! record length
i,j,n ! dummy loop indices
real (r8) :: &
adv_dt_r, &
dzk_div_dt, &
dzkp1_div_dt, &
dTRm1, dTR, dTRp1, &! tracer differences at cell faces
psi, &! fraction of tracer difference appearing in TRACER_E
psi_dTR, &! psi * dTR
work1, work2, work3 ! temporary work space
real (r8), dimension(nx_block,ny_block) :: &
DIV, &! divergence of velocity field
XSTAR, &! tracer w/ advective tendency terms added
MU_z,LW_z, &! coefficients for z contribution
MU_x,LW_x, &! coefficients for x contribution
MU_y,LW_y ! coefficients for y contribution
!-----------------------------------------------------------------------
! compute tracer independent coefficients
!-----------------------------------------------------------------------
adv_dt_r = c1 / adv_dt
if (partial_bottom_cells) then
do j=this_block%jb-2,this_block%je+2
do i=this_block%ib-2,this_block%ie+2
if (WTKB(i,j) > c0) then
LW_z(i,j) = (DZT(i,j,k+1,bid) - adv_dt * WTKB(i,j)) / &
(DZT(i,j,k,bid) + DZT(i,j,k+1,bid))
if (WTKBp1(i,j) > c0) then
MU_z(i,j) = (DZT(i,j,k+1,bid) * adv_dt_r - WTKBp1(i,j)) / WTKB(i,j)
else if (WTKBp1(i,j) < c0) then
MU_z(i,j) = -WTKBp1(i,j) / WTKB(i,j) * &
(DZT(i,j,k+1,bid) + adv_dt * WTKBp1(i,j)) / &
(DZT(i,j,k+1,bid) + DZT(i,j,k+2,bid))
else
MU_z(i,j) = c0
end if
else if (WTKB(i,j) < c0) then
LW_z(i,j) = (DZT(i,j,k,bid) + adv_dt * WTKB(i,j)) / &
(DZT(i,j,k,bid) + DZT(i,j,k+1,bid))
if (WTK(i,j) < c0) then
MU_z(i,j) = -(DZT(i,j,k,bid) * adv_dt_r + WTK(i,j)) / WTKB(i,j)
else if (WTK(i,j) > c0) then
MU_z(i,j) = -WTK(i,j) / WTKB(i,j) * &
(DZT(i,j,k,bid) - adv_dt * WTK(i,j)) / &
(DZT(i,j,k-1,bid) + DZT(i,j,k,bid))
else
MU_z(i,j) = c0
end if
if (WTK(i,j) > c0) then
end if
end if
end do
end do
else
dzk_div_dt = dz(k) * adv_dt_r
if (k < km) dzkp1_div_dt = dz(k+1) * adv_dt_r
work1 = dz(k) * p5_dz_ph_r(k)
if (k < km) work2 = dz(k+1) * p5_dz_ph_r(k)
work3 = adv_dt * p5_dz_ph_r(k)
do j=this_block%jb-2,this_block%je+2
do i=this_block%ib-2,this_block%ie+2
if (WTKB(i,j) > c0) then
LW_z(i,j) = work2 - work3 * WTKB(i,j)
if (WTKBp1(i,j) > c0) then
MU_z(i,j) = (dzkp1_div_dt - WTKBp1(i,j)) / WTKB(i,j)
else if (WTKBp1(i,j) < c0) then
MU_z(i,j) = -WTKBp1(i,j) / WTKB(i,j) * &
(dz(k+1) + adv_dt * WTKBp1(i,j)) * p5_dz_ph_r(k+1)
else
MU_z(i,j) = c0
end if
else if (WTKB(i,j) < c0) then
LW_z(i,j) = work1 + work3 * WTKB(i,j)
if (WTK(i,j) < c0) then
MU_z(i,j) = -(dzk_div_dt + WTK(i,j)) / WTKB(i,j)
else if (WTK(i,j) > c0) then
MU_z(i,j) = -WTK(i,j) / WTKB(i,j) * &
(dz(k) - adv_dt * WTK(i,j)) * p5_dz_ph_r(k-1)
else
MU_z(i,j) = c0
end if
end if
end do
end do
endif
do j=this_block%jb-2,this_block%je+2
i=this_block%ib-2
if (UVEL_E_dt(i,j) < c0) then
LW_x(i,j) = (DXT(i+1,j,bid) + UVEL_E_dt(i,j)) * p5_DXT_ph_R(i,j,bid)
end if
do i=this_block%ib-1,this_block%ie
if (UVEL_E_dt(i,j) > c0) then
LW_x(i,j) = (DXT(i,j,bid) - UVEL_E_dt(i,j)) * p5_DXT_ph_R(i,j,bid)
if (UVEL_E_dt(i-1,j) > c0) then
MU_x(i,j) = (DXT(i,j,bid) - UVEL_E_dt(i-1,j)) / UVEL_E_dt(i,j)
else
MU_x(i,j) = c0
end if
else if (UVEL_E_dt(i,j) < c0) then
LW_x(i,j) = (DXT(i+1,j,bid) + UVEL_E_dt(i,j)) * p5_DXT_ph_R(i,j,bid)
if (UVEL_E_dt(i+1,j) < c0) then
MU_x(i,j) = -(DXT(i+1,j,bid) + UVEL_E_dt(i+1,j)) / UVEL_E_dt(i,j)
else
MU_x(i,j) = c0
end if
else
LW_x(i,j) = DXT(i,j,bid) * p5_DXT_ph_R(i,j,bid)
end if
end do
i=this_block%ie+1
if (UVEL_E_dt(i,j) > c0) then
LW_x(i,j) = (DXT(i,j,bid) - UVEL_E_dt(i,j)) * p5_DXT_ph_R(i,j,bid)
end if
do i=this_block%ib-1,this_block%ie
if ((UVEL_E_dt(i,j) > c0) .and. (UVEL_E_dt(i-1,j) < c0)) then
MU_x(i,j) = -UVEL_E_dt(i-1,j)/UVEL_E_dt(i,j)*LW_x(i-1,j)
end if
if ((UVEL_E_dt(i,j) < c0) .and. (UVEL_E_dt(i+1,j) > c0)) then
MU_x(i,j) = -UVEL_E_dt(i+1,j)/UVEL_E_dt(i,j)*LW_x(i+1,j)
end if
end do
end do
j=this_block%jb-2
do i=this_block%ib,this_block%ie
if (VVEL_N_dt(i,j) < c0) then
LW_y(i,j) = (DYT(i,j+1,bid) + VVEL_N_dt(i,j)) * p5_DYT_ph_R(i,j,bid)
end if
end do
do j=this_block%jb-1,this_block%je
do i=this_block%ib,this_block%ie
if (VVEL_N_dt(i,j) > c0) then
LW_y(i,j) = (DYT(i,j,bid) - VVEL_N_dt(i,j)) * p5_DYT_ph_R(i,j,bid)
if (VVEL_N_dt(i,j-1) > c0) then
MU_y(i,j) = (DYT(i,j,bid) - VVEL_N_dt(i,j-1)) / VVEL_N_dt(i,j)
else
MU_y(i,j) = c0
end if
else if (VVEL_N_dt(i,j) < c0) then
LW_y(i,j) = (DYT(i,j+1,bid) + VVEL_N_dt(i,j)) * p5_DYT_ph_R(i,j,bid)
if (VVEL_N_dt(i,j+1) < c0) then
MU_y(i,j) = -(DYT(i,j+1,bid) + VVEL_N_dt(i,j+1)) / VVEL_N_dt(i,j)
else
MU_y(i,j) = c0
end if
else
LW_y(i,j) = DYT(i,j,bid) * p5_DYT_ph_R(i,j,bid)
end if
end do
end do
j=this_block%je+1
do i=this_block%ib,this_block%ie
if (VVEL_N_dt(i,j) > c0) then
LW_y(i,j) = (DYT(i,j,bid) - VVEL_N_dt(i,j)) * p5_DYT_ph_R(i,j,bid)
end if
end do
do j=this_block%jb-1,this_block%je
do i=this_block%ib,this_block%ie
if ((VVEL_N_dt(i,j) > c0) .and. (VVEL_N_dt(i,j-1) < c0)) then
MU_y(i,j) = -VVEL_N_dt(i,j-1)/VVEL_N_dt(i,j)*LW_y(i,j-1)
end if
if ((VVEL_N_dt(i,j) < c0) .and. (VVEL_N_dt(i,j+1) > c0)) then
MU_y(i,j) = -VVEL_N_dt(i,j+1)/VVEL_N_dt(i,j)*LW_y(i,j+1)
end if
end do
end do
if (partial_bottom_cells) then
DIV = (WTK - WTKB) / DZT(:,:,k,bid) + CE + CW + CN + CS
else
DIV = (WTK - WTKB) * dzr(k) + CE + CW + CN + CS
end if
!-----------------------------------------------------------------------
! loop over tracers
!-----------------------------------------------------------------------
do n=1,nt
if (.not. tr_mask(n)) cycle
call zero_ghost_cells
(this_block,TRACER_E(:,:,n))
call zero_ghost_cells
(this_block,TRACER_N(:,:,n))
call zero_ghost_cells
(this_block,AUXB(:,:,n))
!-----------------------------------------------------------------------
! Compute vertical contribution
! Computations are done in the grix-x & grid-y halo regions
!
! Compute grid-x direction contribution
! Computations are done in the grid-y halo regions because the
! grid-y computations, which use the grid-y halo region,
! rely on the grid-x results.
!
! j loops for z and x contributions fused to improve performance
!
! j loop for y contribution cannot easily be fused
! because of data dependencies
!-----------------------------------------------------------------------
do j=this_block%jb-2,this_block%je+2
do i=this_block%ib-2,this_block%ie+2
if (k+1 <= KMT(i,j,bid)) then
dTR = (X(i,j,k+1,n) - X(i,j,k,n))
if (WTKB(i,j) > c0) then
AUXB(i,j,n) = WTKB(i,j) * X(i,j,k+1,n)
if (k+2 <= KMT(i,j,bid)) then
dTRp1 = (X(i,j,k+2,n) - X(i,j,k+1,n))
!!! if (dTR * dTRp1 > c0) then
!!! psi = min(LW_z(i,j), MU_z(i,j) * dTRp1 / dTR)
!!! AUXB(i,j,n) = WTKB(i,j) * (X(i,j,k+1,n) - psi * dTR)
!!! end if
if (dTR > c0 .and. dTRp1 > c0) then
psi_dTR = min(LW_z(i,j) * dTR, MU_z(i,j) * dTRp1)
AUXB(i,j,n) = WTKB(i,j) * (X(i,j,k+1,n) - psi_dTR)
else if (dTR < c0 .and. dTRp1 < c0) then
psi_dTR = max(LW_z(i,j) * dTR, MU_z(i,j) * dTRp1)
AUXB(i,j,n) = WTKB(i,j) * (X(i,j,k+1,n) - psi_dTR)
end if
end if
else if (WTKB(i,j) < c0) then
AUXB(i,j,n) = WTKB(i,j) * X(i,j,k,n)
if (k > 1) then
dTRm1 = (X(i,j,k,n) - X(i,j,k-1,n))
!!! if (dTR * dTRm1 > c0) then
!!! psi = min(LW_z(i,j), MU_z(i,j) * dTRm1 / dTR)
!!! AUXB(i,j,n) = WTKB(i,j) * (X(i,j,k,n) + psi * dTR)
!!! end if
if (dTR > c0 .and. dTRm1 > c0) then
psi_dTR = min(LW_z(i,j) * dTR, MU_z(i,j) * dTRm1)
AUXB(i,j,n) = WTKB(i,j) * (X(i,j,k,n) + psi_dTR)
else if (dTR < c0 .and. dTRm1 < c0) then
psi_dTR = max(LW_z(i,j) * dTR, MU_z(i,j) * dTRm1)
AUXB(i,j,n) = WTKB(i,j) * (X(i,j,k,n) + psi_dTR)
end if
end if
else
AUXB(i,j,n) = c0
end if
else ! (k+1>KMT)
AUXB(i,j,n) = c0
end if
if (partial_bottom_cells) then
XOUT(i,j,n) = ( AUX(i,j,n,bid) - AUXB(i,j,n) - &
(WTK(i,j) - WTKB(i,j)) * X(i,j,k,n) ) / DZT(i,j,k,bid)
else
XOUT(i,j,n) = ( AUX(i,j,n,bid) - AUXB(i,j,n) - &
(WTK(i,j) - WTKB(i,j)) * X(i,j,k,n) ) * dzr(k)
end if
XSTAR(i,j) = X(i,j,k,n) - adv_dt * XOUT(i,j,n)
end do ! i loop for z contribution
do i=this_block%ib-1,this_block%ie
dTR = (XSTAR(i+1,j) - XSTAR(i,j)) * KMASKE(i,j)
if (CE(i,j) > c0) then
dTRm1 = (XSTAR(i,j) - XSTAR(i-1,j)) * KMASKE(i-1,j)
!!! if (dTR * dTRm1 > c0) then
!!! psi = min(LW_x(i,j), MU_x(i,j) * dTRm1 / dTR)
!!! TRACER_E(i,j,n) = XSTAR(i,j) + psi * dTR
if (dTR > c0 .and. dTRm1 > c0) then
psi_dTR = min(LW_x(i,j) * dTR, MU_x(i,j) * dTRm1)
TRACER_E(i,j,n) = XSTAR(i,j) + psi_dTR
else if (dTR < c0 .and. dTRm1 < c0) then
psi_dTR = max(LW_x(i,j) * dTR, MU_x(i,j) * dTRm1)
TRACER_E(i,j,n) = XSTAR(i,j) + psi_dTR
else
TRACER_E(i,j,n) = XSTAR(i,j)
end if
else if (CE(i,j) < c0) then
dTRp1 = (XSTAR(i+2,j) - XSTAR(i+1,j)) * KMASKE(i+1,j)
!!! if (dTR * dTRp1 > c0) then
!!! psi = min(LW_x(i,j), MU_x(i,j) * dTRp1 / dTR)
!!! TRACER_E(i,j,n) = XSTAR(i+1,j) - psi * dTR
if (dTR > c0 .and. dTRp1 > c0) then
psi_dTR = min(LW_x(i,j) * dTR, MU_x(i,j) * dTRp1)
TRACER_E(i,j,n) = XSTAR(i+1,j) - psi_dTR
else if (dTR < c0 .and. dTRp1 < c0) then
psi_dTR = max(LW_x(i,j) * dTR, MU_x(i,j) * dTRp1)
TRACER_E(i,j,n) = XSTAR(i+1,j) - psi_dTR
else
TRACER_E(i,j,n) = XSTAR(i+1,j)
end if
else
TRACER_E(i,j,n) = XSTAR(i,j) + LW_x(i,j) * dTR
end if
end do ! i loop for grid-x contribution
!*** i loop broken up to avoid dependency on previous i
!*** iteration and to avoid using an XSTAR value that
!*** contains grid-x direction contribution
end do ! j loop for z and grid-x contributions
!*** j loop broken up for overflow TRACER_E modification
if ( overflows_on .and. overflows_interactive ) then
call ovf_advt
(k,TRACER_E(:,:,n),TRACER_N(:,:,n),n,this_block, &
CE,CW,CN,CS)
endif
do j=this_block%jb-2,this_block%je+2
do i=this_block%ib,this_block%ie
work1 = &
CE(i,j) * TRACER_E(i,j,n) + CW(i,j) * TRACER_E(i-1,j,n) - &
(CE(i,j) + CW(i,j)) * X(i,j,k,n)
XOUT(i,j,n) = XOUT(i,j,n) + work1
XSTAR(i,j) = XSTAR(i,j) - adv_dt * work1
end do ! i loop for grid-x contribution
end do ! j loop for z and grid-x contributions
!-----------------------------------------------------------------------
! Compute grid-y direction contribution
! Add divergence term
!-----------------------------------------------------------------------
do j=this_block%jb-1,this_block%je
do i=this_block%ib,this_block%ie
dTR = (XSTAR(i,j+1) - XSTAR(i,j)) * KMASKN(i,j)
if (CN(i,j) > c0) then
dTRm1 = (XSTAR(i,j) - XSTAR(i,j-1)) * KMASKN(i,j-1)
!!! if (dTR * dTRm1 > c0) then
!!! psi = min(LW_y(i,j), MU_y(i,j) * dTRm1 / dTR)
!!! TRACER_N(i,j,n) = XSTAR(i,j) + psi * dTR
if (dTR > c0 .and. dTRm1 > c0) then
psi_dTR = min(LW_y(i,j) * dTR, MU_y(i,j) * dTRm1)
TRACER_N(i,j,n) = XSTAR(i,j) + psi_dTR
else if (dTR < c0 .and. dTRm1 < c0) then
psi_dTR = max(LW_y(i,j) * dTR, MU_y(i,j) * dTRm1)
TRACER_N(i,j,n) = XSTAR(i,j) + psi_dTR
else
TRACER_N(i,j,n) = XSTAR(i,j)
end if
else if (CN(i,j) < c0) then
dTRp1 = (XSTAR(i,j+2) - XSTAR(i,j+1)) * KMASKN(i,j+1)
!!! if (dTR * dTRp1 > c0) then
!!! psi = min(LW_y(i,j), MU_y(i,j) * dTRp1 / dTR)
!!! TRACER_N(i,j,n) = XSTAR(i,j+1) - psi * dTR
if (dTR > c0 .and. dTRp1 > c0) then
psi_dTR = min(LW_y(i,j) * dTR, MU_y(i,j) * dTRp1)
TRACER_N(i,j,n) = XSTAR(i,j+1) - psi_dTR
else if (dTR < c0 .and. dTRp1 < c0) then
psi_dTR = max(LW_y(i,j) * dTR, MU_y(i,j) * dTRp1)
TRACER_N(i,j,n) = XSTAR(i,j+1) - psi_dTR
else
TRACER_N(i,j,n) = XSTAR(i,j+1)
end if
else
TRACER_N(i,j,n) = XSTAR(i,j) + LW_y(i,j) * dTR
end if
end do ! i loop for grid-y contribution
end do ! j loop for grid-y contribution
!*** i,j loops broken for overflow TRACER_N modification
if ( overflows_on .and. overflows_interactive ) then
call ovf_advt
(k,TRACER_E(:,:,n),TRACER_N(:,:,n),n,this_block, &
CE,CW,CN,CS)
endif
do j=this_block%jb-1,this_block%je
do i=this_block%ib,this_block%ie
!*** The formula below is not correct for j==this_block%jb-1
!*** because TRACER_N(:,this_block%jb-2,n) is not computed.
!*** This is OK because XOUT(:,this_block%jb-1,n) is not used.
XOUT(i,j,n) = XOUT(i,j,n) + &
CN(i,j) * TRACER_N(i,j,n) + CS(i,j) * TRACER_N(i,j-1,n) - &
(CN(i,j) + CS(i,j) - DIV(i,j)) * X(i,j,k,n)
!*** do not update XSTAR in order to avoid using an XSTAR
!*** value that contains grid-y direction contribution
end do ! i loop for grid-y contribution
end do ! j loop for grid-y contribution
end do ! tracer loop
!-----------------------------------------------------------------------
end subroutine lw_lim
!-----------------------------------------------------------------------
!EOC
!***********************************************************************
!BOP
! !IROUTINE: zero_ghost_cells
! !INTERFACE:
subroutine zero_ghost_cells(this_block,FIELD) 10
! !DESCRIPTION:
! This routine sets ghost cell values to zero.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block information for this block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(inout) :: &
FIELD ! FIELD to be masked
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j ! dummy loop indices
!-----------------------------------------------------------------------
do j=1,this_block%jb-1
do i=1,nx_block
FIELD(i,j) = c0
end do
end do
do j=this_block%jb,this_block%je
do i=1,this_block%ib-1
FIELD(i,j) = c0
end do
do i=this_block%ie+1,nx_block
FIELD(i,j) = c0
end do
end do
do j=this_block%je+1,ny_block
do i=1,nx_block
FIELD(i,j) = c0
end do
end do
end subroutine zero_ghost_cells
!-----------------------------------------------------------------------
!EOC
!***********************************************************************
end module advection
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||