!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module vertical_mix 5,18
!BOP
! !MODULE: vertical_mix
! !DESCRIPTION:
! This module contains the routines for computing vertical mixing
! tendencies, implicit vertical mixing and convection. Routines for
! computing the diffusion coefficients themselves are in separate
! modules for the specific parameterizations.
! Currently, the options for vertical mixing include:
! \begin{itemize}
! \item constant: a constant (in space and time) vertical diffusion
! \item Rich no: vertical diffusion based on local Richardson number
! (e.g. Pacanowski and Philander)
! \item KPP: k-profile vertical mixing (from Large et al.)
! \end{itemize}
!
! !REVISION HISTORY:
! SVN:$Id: vertical_mix.F90 22881 2010-05-11 04:23:39Z njn01 $
! !USES:
use kinds_mod
use blocks
use distribution
use domain
use constants
use grid
use state_mod
use broadcast
use io
use timers
use tavg
use time_management
use diagnostics
use vmix_const
use vmix_rich
use vmix_kpp
use exit_mod
use prognostic
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_vertical_mix, &
vmix_coeffs, &
vdifft, vdiffu, &
impvmixt, impvmixt_correct, impvmixu, &
convad, impvmixt_tavg
! !PUBLIC DATA MEMBERS:
real (r8), dimension(:,:,:,:,:), allocatable, public, target :: &
VDC ! tracer diffusivity - public to allow
! possible modification by Gent-McWilliams
! horizontal mixing parameterization
real (r8), dimension(:,:,:,:), allocatable, public, target :: &
VDC_GM ! Gent-McWilliams contribution to VDC
integer (int_kind), parameter, public :: &
vmix_type_const = 1, & ! integer identifiers for desired
vmix_type_rich = 2, & ! mixing parameterization
vmix_type_kpp = 3
integer (int_kind), public :: &
vmix_itype ! users choice for vmix parameterization
logical (log_kind), public :: &
implicit_vertical_mix ! flag for computing vertical mixing
! implicitly in time
!EOP
!BOC
!-----------------------------------------------------------------------
!
! vertical mixing quantities
!
!-----------------------------------------------------------------------
real (r8), dimension(:,:,:,:), allocatable, target :: &
VVC ! momentum viscosity
!-----------------------------------------------------------------------
!
! convection variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
convection_itype, &! user choice of method to use for convection
nconvad ! number of times to convectively adjust
real (r8) :: &
convect_diff, &! tracer diffusion to use for convection
convect_visc ! viscosity to use for convection (momentum)
integer (int_kind), parameter :: &
convect_type_adjust = 1, &! ids for convection choice
convect_type_diff = 2
!-----------------------------------------------------------------------
!
! vectors needed for implicit vertical mixing
!
!-----------------------------------------------------------------------
real (r8) :: &
aidif ! time-centering parameter for implicit vmix
real (r8), dimension(:), allocatable :: &
afac_u, afac_t
!-----------------------------------------------------------------------
!
! these arrays store flux information to avoid re-computing fluxes
! for each k level
!
!-----------------------------------------------------------------------
real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic) :: &
VTF ! vertical tracer flux at top of T-box
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
VUF,VVF ! vertical momentum fluxes at top of box
!-----------------------------------------------------------------------
!
! bottom drag and bottom (geothermal) heat flux
!
!-----------------------------------------------------------------------
logical (log_kind) :: &
lbottom_heat_flx ! true if bottom heat flux /= 0
real (r8) :: &
bottom_drag, &! drag coefficient for bottom drag
bottom_heat_flx, &! bottom (geothermal) heat flux (W/m2)
bottom_heat_flx_depth ! depth below which bottom heat applied
!-----------------------------------------------------------------------
!
! timer ids and tavg ids for time-averaged vertical mix diagnostics
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
timer_vmix_coeffs, &! timer number for vertical mix coeffs
timer_vdiffu, &! timer number for explicit momentum vmix
timer_vdifft, &! timer number for explicit tracer vmix
timer_impvmixt, &! timer number for implicit tracer vmix
timer_impvmixu ! timer number for implicit momentum vmix
integer (int_kind) :: &
tavg_VDC_T, &! tavg id for total vertical TEMP diffusivity
tavg_VDC_S, &! tavg id for total vertical SALT diffusivity
tavg_VVC, &! tavg id for total vertical momentum viscosity
tavg_VUF, &! tavg id for vertical flux of U momentum
tavg_VVF, &! tavg id for vertical flux of V momentum
tavg_PEC, &! tavg id for pot energy release convection
tavg_NCNV ! tavg id for number of convective adjustments
integer (int_kind), dimension(nt) :: &
tavg_DIA_IMPVF_TRACER ! tavg id for diabatic implicit vertical flux of tracer
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_vertical_mix
! !INTERFACE:
subroutine init_vertical_mix 1,36
! !DESCRIPTION:
! Initializes various mixing quantities and calls initialization
! routines for specific parameterizations.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, &! vertical level index
n, &! dummy loop index
nu, &! i/o unit
nml_error ! namelist i/o error flag
character (char_len) :: &
vmix_choice ! input choice for desired parameterization
character (char_len) :: &
convection_type ! input choice for method for convection
namelist /vertical_mix_nml/ vmix_choice, aidif, bottom_drag, &
bottom_heat_flx, bottom_heat_flx_depth, &
implicit_vertical_mix, &
convection_type, nconvad, &
convect_diff, convect_visc
!-----------------------------------------------------------------------
!
! read input namelist and set mixing type
!
!-----------------------------------------------------------------------
aidif = c1
bottom_drag = 1.0e-3_r8
bottom_heat_flx = c0
bottom_heat_flx_depth = 1000.00e2_r8
implicit_vertical_mix = .true.
convection_type = 'diffusion'
nconvad = 2
convect_diff = 1000.0_r8
convect_visc = 1000.0_r8
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
do while (nml_error > 0)
read(nml_in, nml=vertical_mix_nml, iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
endif
call broadcast_scalar(nml_error, master_task)
if (nml_error /= 0) then
call exit_POP(sigAbort, &
'ERROR reading vertical_mix_nml')
endif
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,ndelim_fmt)
write(stdout,blank_fmt)
write(stdout,'(a27)') 'Vertical mixing options'
write(stdout,blank_fmt)
write(stdout,delim_fmt)
if (implicit_vertical_mix) then
write(stdout,'(a33)') ' Implicit vertical mixing enabled'
else
write(stdout,'(a34)') ' Implicit vertical mixing disabled'
endif
select case (vmix_choice)
case ('const')
vmix_itype = vmix_type_const
write(stdout,'(a32)') ' Using constant vertical mixing'
case ('rich')
vmix_itype = vmix_type_rich
write(stdout,'(a32)') ' Using Richardson number mixing'
case ('kpp')
vmix_itype = vmix_type_kpp
write(stdout,'(a32)') ' Using KPP vertical mixing '
case default
vmix_itype = -1000
end select
write(stdout,'(a17,2x,1pe12.5)') ' aidif =',aidif
write(stdout,'(a17,2x,1pe12.5)') ' bottom_drag =',bottom_drag
select case (convection_type)
case ('adjustment')
convection_itype = convect_type_adjust
write(stdout,'(a30)') ' Convective adjustment enabled'
write(stdout,'(a7,i4,a7)') ' Using ',nconvad,' passes'
case ('diffusion')
convection_itype = convect_type_diff
write(stdout,'(a41)') ' Convection simulated by strong diffusion'
write(stdout,'(a18,1pe10.3)') ' convect_diff = ',convect_diff
write(stdout,'(a18,1pe10.3)') ' convect_visc = ',convect_visc
case default
convection_itype = -1000
end select
if (bottom_heat_flx /= c0) then
write(stdout,'(a36)') ' Bottom geothermal heat flux enabled'
write(stdout,'(a17,1pe10.3)') ' Flux (W/m2) = ', &
bottom_heat_flx
write(stdout,'(a17,1pe10.3)') ' Depth (cm) = ', &
bottom_heat_flx_depth
else
write(stdout,'(a37)') ' Bottom geothermal heat flux disabled'
endif
endif
call broadcast_scalar(vmix_itype, master_task)
call broadcast_scalar(aidif, master_task)
call broadcast_scalar(bottom_drag, master_task)
call broadcast_scalar(bottom_heat_flx, master_task)
call broadcast_scalar(bottom_heat_flx_depth, master_task)
call broadcast_scalar(implicit_vertical_mix, master_task)
call broadcast_scalar(convection_itype, master_task)
call broadcast_scalar(nconvad, master_task)
call broadcast_scalar(convect_diff, master_task)
call broadcast_scalar(convect_visc, master_task)
bottom_heat_flx = bottom_heat_flx * hflux_factor
if (bottom_heat_flx /= c0) then
lbottom_heat_flx = .true.
else
lbottom_heat_flx = .false.
endif
!-----------------------------------------------------------------------
!
! do some error and consistency checks
!
!-----------------------------------------------------------------------
if (vmix_itype == -1000) then
call exit_POP(sigAbort,'Unknown vertical mixing choice')
endif
if (convection_itype == -1000) then
call exit_POP(sigAbort, &
'ERROR: Unknown option for convection type')
endif
if (implicit_vertical_mix .and. &
convection_itype == convect_type_adjust) then
call exit_POP(sigAbort, &
'Convective adjustment should not be used with implicit vmix')
endif
if (.not. implicit_vertical_mix .and. &
convection_itype == convect_type_diff) then
call exit_POP(sigAbort, &
'Implicit vertical mixing required for diffusive convection')
endif
if (.not. implicit_vertical_mix .and. &
vmix_itype == vmix_type_kpp) then
call exit_POP(sigAbort, &
'Implicit vertical mixing required for KPP')
endif
!-----------------------------------------------------------------------
!
! allocate VDC, VVC arrays and define options for chosen
! parameterization
!
!-----------------------------------------------------------------------
select case (vmix_itype)
case(vmix_type_const)
if (implicit_vertical_mix) then
allocate (VDC(nx_block,ny_block,km,1,nblocks_clinic), &
VVC(nx_block,ny_block,km, nblocks_clinic))
else
allocate (VDC(nx_block,ny_block,1,1,nblocks_clinic), &
VVC(nx_block,ny_block,1 ,nblocks_clinic))
endif
call init_vmix_const(VDC,VVC)
call get_timer(timer_vmix_coeffs,'VMIX_COEFFICIENTS_CONSTANT', &
nblocks_clinic, distrb_clinic%nprocs)
case(vmix_type_rich)
if (implicit_vertical_mix) then
allocate (VDC(nx_block,ny_block,km,1,nblocks_clinic), &
VVC(nx_block,ny_block,km ,nblocks_clinic))
else
!*** note that VVC must be 3-d because it is used
!*** in a different k-loop from where it is defined
allocate (VDC(nx_block,ny_block,1,1,nblocks_clinic), &
VVC(nx_block,ny_block,km ,nblocks_clinic))
endif
call init_vmix_rich(VDC,VVC)
call get_timer(timer_vmix_coeffs,'VMIX_COEFFICIENTS_RICH', &
nblocks_clinic, distrb_clinic%nprocs)
case(vmix_type_kpp)
allocate (VDC(nx_block,ny_block,0:km+1,2,nblocks_clinic), &
VVC(nx_block,ny_block,km, nblocks_clinic))
call init_vmix_kpp(VDC,VVC)
call get_timer(timer_vmix_coeffs,'VMIX_COEFFICIENTS_KPP', &
nblocks_clinic, distrb_clinic%nprocs)
end select
call get_timer(timer_vdifft,'VMIX_EXPLICIT_TRACER', &
nblocks_clinic, distrb_clinic%nprocs)
call get_timer(timer_vdiffu,'VMIX_EXPLICIT_MOMENTUM', &
nblocks_clinic, distrb_clinic%nprocs)
!-----------------------------------------------------------------------
!
! set up implicit coefficients if implicit vertical mixing chosen
!
!-----------------------------------------------------------------------
if (implicit_vertical_mix) then
call get_timer(timer_impvmixt,'VMIX_IMPLICIT_TRACER', &
nblocks_clinic, distrb_clinic%nprocs)
call get_timer(timer_impvmixu,'VMIX_IMPLICIT_MOMENTUM', &
nblocks_clinic, distrb_clinic%nprocs)
allocate(afac_u(km), afac_t(km))
do k=1,km
afac_u(k) = aidif*dzwr(k)
afac_t(k) = aidif*dzwr(k)
enddo
endif
!-----------------------------------------------------------------------
!
! define fields for accumulating tavg diagnostics
!
!-----------------------------------------------------------------------
if (vmix_itype == vmix_type_kpp) then
call define_tavg_field(tavg_VDC_T,'VDC_T',3, &
long_name='total diabatic vertical TEMP diffusivity', &
units='cm^2/s', grid_loc='3113', &
coordinates='TLONG TLAT z_w_bot time')
call define_tavg_field(tavg_VDC_S,'VDC_S',3, &
long_name='total diabatic vertical SALT diffusivity', &
units='cm^2/s', grid_loc='3113', &
coordinates='TLONG TLAT z_w_bot time')
call define_tavg_field(tavg_VVC,'VVC',3, &
long_name='total vertical momentum viscosity', &
units='cm^2/s', grid_loc='3113', &
coordinates='TLONG TLAT z_w_bot time')
endif
call define_tavg_field(tavg_VUF,'VUF',3, &
long_name='Zonal viscous stress', &
units=' ', grid_loc='3222')
call define_tavg_field(tavg_VVF,'VVF',3, &
long_name='Meridional viscous stress', &
units=' ', grid_loc='3222')
if (convection_itype == convect_type_adjust) then
call define_tavg_field(tavg_PEC,'PEC',3, &
long_name='Potential energy release convection', &
units='g/cm^3', grid_loc='3111')
call define_tavg_field(tavg_NCNV,'NCNV',3, &
long_name='Convective adjustments per second', &
units='adj/s', grid_loc='3111')
endif
if (implicit_vertical_mix) then
do n = 1,nt
call define_tavg_field(tavg_DIA_IMPVF_TRACER(n), 'DIA_IMPVF_' /&
&/ trim(tracer_d(n)%short_name),3, &
long_name=trim(tracer_d(n)%short_name) /&
&/ ' Flux Across Bottom Face from Diabatic Implicit Vertical Mixing', &
units=trim(tracer_d(n)%flux_units), grid_loc='3113', &
scale_factor=tracer_d(n)%scale_factor, &
coordinates='TLONG TLAT z_w_bot time')
enddo
endif
!-----------------------------------------------------------------------
!EOC
call flushm(stdout)
end subroutine init_vertical_mix
!***********************************************************************
!BOP
! !IROUTINE: vmix_coeffs
! !INTERFACE:
subroutine vmix_coeffs(k, TMIX, UMIX, VMIX, RHOMIX, & 2,16
STF, SHF_QSW, &
this_block, &
SMF, SMFT)
! !DESCRIPTION:
! This is a driver routine which calls the appropriate
! parameterization to compute the mixing coefficients VDC, VVC.
! This routine must be called successively for $k=1,2,3$...
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k ! vertical level index
real (r8), dimension(nx_block,ny_block,km,nt) :: TMIX
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UMIX, VMIX, &! U,V at mix time level
RHOMIX ! density at mix time level
real (r8), dimension(nx_block,ny_block,nt), intent(in) :: &
STF ! surface forcing for all tracers
real (r8), dimension(nx_block,ny_block), intent(in) :: &
SHF_QSW ! short-wave forcing
real (r8), dimension(nx_block,ny_block,2), intent(in), optional :: &
SMF, &! surface momentum forcing at U points
SMFT ! surface momentum forcing at T points
! *** must pass either one or the other
type (block), intent(in) :: &
this_block ! block info for current block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
bid, &! local block address
kk ! vertical level index
!-----------------------------------------------------------------------
!
! determine which block we are working on
!
!-----------------------------------------------------------------------
bid = this_block%local_id
call timer_start(timer_vmix_coeffs, block_id=bid)
!-----------------------------------------------------------------------
!
! call proper routine based on vmix choice
!
!-----------------------------------------------------------------------
select case(vmix_itype)
case(vmix_type_const)
if (convection_itype == convect_type_diff) then
call vmix_coeffs_const(k,VDC(:,:,:,:,bid), &
VVC(:,:,:, bid),TMIX, &
this_block, convect_diff, convect_visc)
else
call vmix_coeffs_const(k,VDC(:,:,:,:,bid), &
VVC(:,:,:, bid),TMIX, &
this_block)
endif
case(vmix_type_rich)
if (convection_itype == convect_type_diff) then
call vmix_coeffs_rich(k,VDC(:,:,:,:,bid), &
VVC(:,:,:, bid), &
TMIX,UMIX,VMIX,RHOMIX, &
this_block, convect_diff, convect_visc)
else
call vmix_coeffs_rich(k,VDC(:,:,:,:,bid), &
VVC(:,:,:, bid), &
TMIX,UMIX,VMIX,RHOMIX, &
this_block)
endif
case(vmix_type_kpp)
if (k == 1) then ! for KPP, compute coeffs for all levels
if (.not. present(SMFT) .and. .not. present(SMF)) &
call exit_POP(sigAbort, &
'vmix_coeffs: must supply either SMF,SMFT')
if (present(SMFT)) then
call vmix_coeffs_kpp(VDC(:,:,:,:,bid), &
VVC(:,:,:, bid), &
TMIX,UMIX,VMIX,RHOMIX, &
STF,SHF_QSW, &
this_block, &
convect_diff, convect_visc, &
SMFT=SMFT)
else
call vmix_coeffs_kpp(VDC(:,:,:,:,bid), &
VVC(:,:,:, bid), &
TMIX,UMIX,VMIX,RHOMIX, &
STF,SHF_QSW, &
this_block, &
convect_diff, convect_visc, &
SMF=SMF)
endif
if (tavg_requested(tavg_VDC_T) .and. mix_pass /= 1) then
do kk=1,km
! kk index is no longer shifted because z_w_bot is now an output coordinate
! VDC is at cell bottom
call accumulate_tavg_field(VDC(:,:,kk,1,bid),tavg_VDC_T,bid,kk)
end do
endif
if (tavg_requested(tavg_VDC_S) .and. mix_pass /= 1) then
do kk=1,km
! kk index is no longer shifted because z_w_bot is now an output coordinate
! VDC is at cell bottom
call accumulate_tavg_field(VDC(:,:,kk,2,bid),tavg_VDC_S,bid,kk)
end do
endif
if (tavg_requested(tavg_VVC) .and. mix_pass /= 1) then
do kk=1,km
! kk index is no longer shifted because z_w_bot is now an output coordinate
! VVC is at cell bottom
call accumulate_tavg_field(VVC(:,:,kk,bid),tavg_VVC,bid,kk)
end do
endif
endif
end select
call timer_stop(timer_vmix_coeffs, block_id=bid)
!-----------------------------------------------------------------------
!
! call cfl diagnostics if not implicit mixing
!
!-----------------------------------------------------------------------
if (.not. implicit_vertical_mix) then
call cfl_vdiff(k, VDC(:,:,:,:,bid), VVC(:,:,:,bid), this_block)
endif
!-----------------------------------------------------------------------
!EOC
end subroutine vmix_coeffs
!***********************************************************************
!BOP
! !IROUTINE: vdifft
! !INTERFACE:
subroutine vdifft(k, VDTK, TOLD, STF, this_block) 1,2
! !DESCRIPTION:
! Computes vertical diffusion of tracers
! \begin{equation}
! D_V(\varphi) = \delta_z(\kappa\delta_z\varphi)
! \end{equation}
! \begin{eqnarray}
! \kappa\delta_z\varphi = Q_\varphi &{\rm at }& z=0 \\
! \kappa\delta_z\varphi = 0 &{\rm at }& z=-H_T
! \end{eqnarray}
!
! This routine must be called successively with $k = 1,2,3,$...
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! vertical level
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TOLD ! tracers at old time level
real (r8), dimension(nx_block,ny_block,nt), intent(in) :: &
STF ! surface forcing for all tracers
type (block), intent(in) :: &
this_block ! block info for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(out) :: &
VDTK ! returns VDIFF(Tracer(:,:,k,n))
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
n, &! dummy loop counter for tracer number
bid, &! local block index
kp1, &! k+1
kvdc, &! k index into mixing coefficient array
mt2 ! index for when temperature VDC is different
real (r8), dimension(nx_block,ny_block) :: &
VTFB ! vertical tracer flux across bottom of T-box
!-----------------------------------------------------------------------
!
! set up vertical indices
!
!-----------------------------------------------------------------------
bid = this_block%local_id
call timer_start(timer_vdifft, block_id=bid)
if (k < km) then
kp1 = k + 1
else
kp1 = km
endif
kvdc = min(k,size(VDC,DIM=3)) !*** reduce to 1 if VDC 2-d array
!-----------------------------------------------------------------------
!
! start loop over tracers
! set mt2 to be either 1 (if VDC same for all tracers) or n
! if coefficients are different for tracers
!
!-----------------------------------------------------------------------
do n = 1,nt
mt2 = min(n,size(VDC,DIM=4))
!-----------------------------------------------------------------------
!
! surface tracer fluxes
!
!-----------------------------------------------------------------------
if (k == 1) then
VTF(:,:,n,bid) = merge(STF(:,:,n), c0, KMT(:,:,bid) >= 1)
endif
!-----------------------------------------------------------------------
!
! vertical tracer flux vdc*Dz(T) at bottom of T box
! calculate Dz(vdc*Dz(T))
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
!CDIR COLLAPSE
VTFB = merge(VDC(:,:,kvdc,mt2,bid)* &
(TOLD(:,:,k ,n) - TOLD(:,:,kp1,n))/ &
(p5*(DZT(:,:,k,bid) + DZT(:,:,kp1,bid))) &
,c0, KMT(:,:,bid) > k)
if (lbottom_heat_flx .and. n == 1) then
!CDIR COLLAPSE
VTFB = merge( -bottom_heat_flx, VTFB, &
k == KMT(:,:,bid) .and. &
(zt(k) + p5*DZT(:,:,k,bid)) >= &
bottom_heat_flx_depth)
endif
!CDIR COLLAPSE
VDTK(:,:,n) = merge((VTF(:,:,n,bid) - VTFB)/DZT(:,:,k,bid), &
c0, k <= KMT(:,:,bid))
else
!CDIR COLLAPSE
VTFB = merge(VDC(:,:,kvdc,mt2,bid)* &
(TOLD(:,:,k ,n) - TOLD(:,:,kp1,n))*dzwr(k) &
,c0, KMT(:,:,bid) > k)
if (lbottom_heat_flx .and. n == 1) then
!CDIR COLLAPSE
VTFB = merge( -bottom_heat_flx, VTFB, &
k == KMT(:,:,bid) .and. &
zw(k) >= bottom_heat_flx_depth)
endif
!CDIR COLLAPSE
VDTK(:,:,n) = merge((VTF(:,:,n,bid) - VTFB)*dzr(k), &
c0, k <= KMT(:,:,bid))
endif
!-----------------------------------------------------------------------
!
! set top value of VTF to bottom value for next pass at level k+1
!
!-----------------------------------------------------------------------
!CDIR COLLAPSE
VTF(:,:,n,bid) = VTFB
enddo
call timer_stop(timer_vdifft, block_id=bid)
!-----------------------------------------------------------------------
!EOC
end subroutine vdifft
!***********************************************************************
!BOP
! !IROUTINE: vdiffu
! !INTERFACE:
subroutine vdiffu(k, VDUK, VDVK, UOLD, VOLD, SMF, this_block) 1,6
! !DESCRIPTION:
! Computes vertical diffusion of momentum:
! \begin{equation}
! D_V(u,v) = \delta_z(\mu\delta_z(u,v))
! \end{equation}
! \begin{eqnarray}
! \mu\delta_z(u,v) = (\tau_x,\tau_y) &{\rm at }& z=0 \\
! \mu\delta_z(u,v) = c|{\bf\rm u}|(u,v) &{\rm at }& z=-H_U
! \end{eqnarray}
!
! This routine must be called successively with $k = 1,2,3,$...
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! vertical level index
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UOLD, &! U velocity at old time level
VOLD ! V velocity at old time level
real (r8), dimension(nx_block,ny_block,2), intent(in) :: &
SMF ! surface momentum fluxes
type (block), intent(in) :: &
this_block ! block info for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
VDUK, &! returns Vdiff(U)
VDVK ! returns Vdiff(V)
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
bid, &! local block index
i,j, &! loop indices
kp1, &! k+1
kvvc ! index into viscosity array
real (r8) :: &
vmag ! temp for bottom drag
real (r8), dimension(nx_block,ny_block) :: &
VUFB,VVFB, &! vertical momentum fluxes at bottom of box
WORK ! local work array
!-----------------------------------------------------------------------
!
! set vertical level indices
!
!-----------------------------------------------------------------------
bid = this_block%local_id
call timer_start(timer_vdiffu, block_id=bid)
if (k < km) then
kp1 = k + 1
else
kp1 = km
endif
kvvc = min(k,size(VVC,DIM=3)) !*** reduce to 1 if VVC 2-d array
!-----------------------------------------------------------------------
!
! surface momentum fluxes (wind stress)
!
!-----------------------------------------------------------------------
if (k == 1) then
VUF(:,:,bid) = merge(SMF(:,:,1), c0, KMU(:,:,bid) >= 1)
VVF(:,:,bid) = merge(SMF(:,:,2), c0, KMU(:,:,bid) >= 1)
endif
!-----------------------------------------------------------------------
!
! vertical momentum flux vvc*Dz{(Ub,Vb)} at bottom of U box
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
if (k < km) then
WORK = p5*(DZU(:,:,k,bid)+DZU(:,:,kp1,bid))
else
WORK = p5*DZU(:,:,kp1,bid)
endif
!CDIR COLLAPSE
VUFB = VVC(:,:,kvvc,bid)*(UOLD(:,:,k ) - &
UOLD(:,:,kp1))/WORK
!CDIR COLLAPSE
VVFB = VVC(:,:,kvvc,bid)*(VOLD(:,:,k ) - &
VOLD(:,:,kp1))/WORK
else
!CDIR COLLAPSE
VUFB = VVC(:,:,kvvc,bid)*(UOLD(:,:,k ) - &
UOLD(:,:,kp1))*dzwr(k)
!CDIR COLLAPSE
VVFB = VVC(:,:,kvvc,bid)*(VOLD(:,:,k ) - &
VOLD(:,:,kp1))*dzwr(k)
endif
!-----------------------------------------------------------------------
!
! bottom momentum fluxes from quadratic drag law
!
!-----------------------------------------------------------------------
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
if (k == KMU(i,j,bid)) then
vmag = bottom_drag*sqrt(UOLD(i,j,k)**2 + VOLD(i,j,k)**2)
VUFB(i,j) = vmag*UOLD(i,j,k)
VVFB(i,j) = vmag*VOLD(i,j,k)
endif
end do
end do
!-----------------------------------------------------------------------
!
! calculate Dz(vdc*Dz{(U,V)})
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
VDUK = merge((VUF(:,:,bid) - VUFB)/DZU(:,:,k,bid), &
c0, k <= KMU(:,:,bid))
VDVK = merge((VVF(:,:,bid) - VVFB)/DZU(:,:,k,bid), &
c0, k <= KMU(:,:,bid))
else
VDUK = merge((VUF(:,:,bid) - VUFB)*dzr(k), &
c0, k <= KMU(:,:,bid))
VDVK = merge((VVF(:,:,bid) - VVFB)*dzr(k), &
c0, k <= KMU(:,:,bid))
endif
!-----------------------------------------------------------------------
!
! set top value of (VUF,VVF) to bottom value for next pass (level k+1)
!
!-----------------------------------------------------------------------
VUF(:,:,bid) = VUFB
VVF(:,:,bid) = VVFB
call timer_stop(timer_vdiffu, block_id=bid)
!-----------------------------------------------------------------------
!
! accumulate time-average of vertical diffusion fluxes if reqd
!
!-----------------------------------------------------------------------
if (tavg_requested(tavg_VUF) .and. mix_pass /= 1) then
call accumulate_tavg_field(VUF(:,:,bid),tavg_VUF,bid,k)
endif
if (tavg_requested(tavg_VVF) .and. mix_pass /= 1) then
call accumulate_tavg_field(VVF(:,:,bid),tavg_VVF,bid,k)
endif
!-----------------------------------------------------------------------
!EOC
end subroutine vdiffu
!***********************************************************************
!BOP
! !IROUTINE: impvmixt
! !INTERFACE:
subroutine impvmixt(TNEW, TOLD, PSFC, nfirst, nlast, this_block) 4,2
! !DESCRIPTION:
! Computes the implicit vertical mixing of tracers
! using tridiagonal solver in each vertical column.
! Since the top boundary condition (given tracer flux) does not depend
! on the tracer value, its influence is accounted for in the explicit
! part (subroutine vdifft). The bottom b.c. is zero flux.
!
! !REVISION HISTORY:
! same as module
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(inout) :: &
TNEW ! on input, contains right hand side
! on output, contains updated tracers at new time
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TOLD ! old tracer to update with del(Tracer)
real (r8), dimension(nx_block,ny_block), intent(in) :: &
PSFC ! surface pressure for use in determining
! variable thickness surface layer
integer (int_kind), intent(in) :: &
nfirst, nlast ! only update tracers from nfirst to nlast
type (block), intent(in) :: &
this_block ! block info for current block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,k,n, &! dummy loop indices
ib,ie,jb,je, &! beg,end indices for physical domain
mt2, &! index for selecting tracer coefficient
bid ! local block address
real (r8), dimension(nx_block,ny_block) :: &
A,B,C,D ! various temporaries
real (r8), dimension(nx_block,ny_block,km) :: &
E,F ! various work arrays
real (r8), dimension(km) :: &
hfac_t
real (r8), dimension(nx_block,ny_block) :: &
H1 ! factor containing full thickness of sfc layer
!-----------------------------------------------------------------------
!
! check and initialize some necessary quantities
!
!-----------------------------------------------------------------------
if (nfirst > nlast .or. nfirst > nt) return
bid = this_block%local_id
ib = this_block%ib
ie = this_block%ie
jb = this_block%jb
je = this_block%je
call timer_start(timer_impvmixt,block_id=bid)
!*** note that this array is overwritten for partial bottom cells
do k=1,km
hfac_t(k) = dz(k)/c2dtt(k)
end do
!*** correct hfac for full surface layer thickness
if (sfc_layer_type == sfc_layer_varthick) then
H1 = hfac_t(1) + PSFC/(grav*c2dtt(1))
else
H1 = hfac_t(1)
endif
!-----------------------------------------------------------------------
!
! loop over tracers
!
!-----------------------------------------------------------------------
do n = nfirst,nlast
mt2 = min(n,size(VDC,DIM=4))
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
!*** perform tridiagonal solve in the vertical for every
!*** horizontal grid point
A(i,j) = afac_t(1)*VDC(i,j,1,mt2,bid)
D(i,j) = H1(i,j) + A(i,j)
E(i,j,1) = A(i,j)/D(i,j)
B(i,j) = H1(i,j)*E(i,j,1)
F(i,j,1) = hfac_t(1)*TNEW(i,j,1,n)/D(i,j)
end do
end do
do k=2,km
if (partial_bottom_cells) then
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
C(i,j) = A(i,j)
A(i,j) = aidif*VDC(i,j,k,mt2,bid)/ &
(p5*(DZT(i,j,k ,bid) + &
DZT(i,j,k+1,bid)))
hfac_t(k) = DZT(i,j,k,bid)/c2dtt(k)
end do
end do
else
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
C(i,j) = A(i,j)
A(i,j) = afac_t(k)*VDC(i,j,k,mt2,bid)
end do
end do
endif ! partial_bottom_cells
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
if (k > KMT(i,j,bid)) then
F(i,j,k) = c0
else
if (k == KMT(i,j,bid)) then
D(i,j) = hfac_t(k)+B(i,j)
else
D(i,j) = hfac_t(k)+A(i,j)+B(i,j)
endif
E(i,j,k) = A(i,j)/D(i,j)
B(i,j) = (hfac_t(k) + B(i,j))*E(i,j,k)
F(i,j,k) = (hfac_t(k)*TNEW(i,j,k,n) + &
C(i,j)*F(i,j,k-1))/D(i,j)
endif
end do ! i
end do ! j
end do ! k
do k=km-1,1,-1
do j=jb,je
do i=ib,ie
!*** back substitution
if (k < KMT(i,j,bid)) &
F(i,j,k) = F(i,j,k) + E(i,j,k)*F(i,j,k+1)
end do
end do
end do
!-----------------------------------------------------------------------
!
! The solution of the system (UZ) is DeltaTracer -- it has
! already been multiplied by dtt so we can advance the tracers
! directly
!
!-----------------------------------------------------------------------
!CDIR COLLAPSE
do k=1,km
do j=jb,je
do i=ib,ie
TNEW(i,j,k,n) = TOLD(i,j,k,n) + F(i,j,k)
end do
end do
end do
!-----------------------------------------------------------------------
!
! end tracer loop and stop timer
!
!-----------------------------------------------------------------------
end do
call timer_stop(timer_impvmixt,block_id=bid)
!-----------------------------------------------------------------------
!EOC
end subroutine impvmixt
!***********************************************************************
subroutine impvmixt_tavg(TNEW, bid) 1,2
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TNEW ! on input, contains tracer to update from
! on output, contains updated tracers at new time
integer (int_kind), intent(in) :: &
bid ! local block address
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
real (r8), dimension(nx_block,ny_block) :: &
WORK1, WORK2
integer (int_kind) :: &
k,n, &! dummy loop indices
mt2 ! index for selecting tracer coefficient
!-----------------------------------------------------------------------
do n = 1,nt
mt2 = min(n,size(VDC,DIM=4))
if (tavg_requested(tavg_DIA_IMPVF_TRACER(n))) then
do k=1,km-1
if (allocated(VDC_GM)) then
WORK1 = VDC(:,:,k,mt2,bid) - VDC_GM(:,:,k,bid)
else
WORK1 = VDC(:,:,k,mt2,bid)
endif
if (partial_bottom_cells) then
WORK2 = merge(WORK1*(TNEW(:,:,k,n) - TNEW(:,:,k+1,n))/ &
(p5*(DZT(:,:,k,bid) + DZT(:,:,k+1,bid))) &
,c0, k < KMT(:,:,bid))
if (lbottom_heat_flx .and. n == 1) then
WORK2 = merge( -bottom_heat_flx, WORK2, &
k == KMT(:,:,bid) .and. &
(zt(k) + p5*DZT(:,:,k,bid)) >= &
bottom_heat_flx_depth)
endif
else
WORK2 = merge(WORK1*(TNEW(:,:,k,n) - TNEW(:,:,k+1,n))*dzwr(k) &
,c0, k < KMT(:,:,bid))
if (lbottom_heat_flx .and. n == 1) then
WORK2 = merge( -bottom_heat_flx, WORK2, &
k == KMT(:,:,bid) .and. &
zw(k) >= bottom_heat_flx_depth)
endif
endif
call accumulate_tavg_field(WORK2,tavg_DIA_IMPVF_TRACER(n),bid,k)
end do
endif
end do
end subroutine impvmixt_tavg
!***********************************************************************
!BOP
! !IROUTINE: impvmixt_correct
! !INTERFACE:
subroutine impvmixt_correct(TNEW, PSFC, RHS, nfirst, nlast, this_block) 1,2
! !DESCRIPTION:
!
! Computes implicit vertical mixing of tracers for a corrector step
! using tridiagonal solver in each vertical column.
! Since the top boundary condition (given tracer flux) does not depend
! on the tracer value, its influence is accounted for in the explicit
! part (subroutine vdifft). The bottom b.c. is zero flux.
!
! !REVISION HISTORY:
! same as module
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(inout) :: &
TNEW ! on input, contains tracer to update from
! on output, contains updated tracers at new time
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(in) :: &
PSFC ! surface pressure for use in determining
! variable thickness surface layer
real (r8), dimension(nx_block,ny_block,nt), intent(in) :: &
RHS ! right hand side of linear system
integer (int_kind), intent(in) :: &
nfirst, nlast ! only update tracers from nfirst to nlast
type (block), intent(in) :: &
this_block ! block info for current block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,k,n, &! dummy loop indices
ib,ie,jb,je, &! beg,end indices for physical domain
mt2, &! index for selecting tracer coefficient
bid ! local block address
real (r8), dimension(nx_block,ny_block) :: &
A,B,C,D ! various temporaries
real (r8), dimension(nx_block,ny_block,km) :: &
E,F ! various work arrays
real (r8), dimension(km) :: &
hfac_t
real (r8), dimension(nx_block,ny_block) :: &
H1 ! factor containing full thickness of sfc layer
!-----------------------------------------------------------------------
!
! check and initialize some necessary quantities
!
!-----------------------------------------------------------------------
if (nfirst > nlast .or. nfirst > nt) return
bid = this_block%local_id
ib = this_block%ib
ie = this_block%ie
jb = this_block%jb
je = this_block%je
call timer_start(timer_impvmixt,block_id=bid)
!*** note that this array is overwritten for partial bottom cells
do k=1,km
hfac_t(k) = dz(k)/c2dtt(k)
end do
!*** correct hfac for full surface layer thickness
if (sfc_layer_type == sfc_layer_varthick) then
H1 = hfac_t(1) + PSFC/(grav*c2dtt(1))
else
H1 = hfac_t(1)
endif
!-----------------------------------------------------------------------
!
! loop over tracers
!
!-----------------------------------------------------------------------
do n = nfirst,nlast
mt2 = min(n,size(VDC,DIM=4))
!*** perform tridiagonal solve in the vertical for every
!*** horizontal grid point
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
A(i,j) = afac_t(1)*VDC(i,j,1,mt2,bid)
D(i,j) = H1(i,j) + A(i,j)
E(i,j,1) = A(i,j)/D(i,j)
B(i,j) = H1(i,j)*E(i,j,1)
F(i,j,1) = hfac_t(1)*RHS(i,j,n)/D(i,j)
end do
end do
do k=2,km
if (partial_bottom_cells) then
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
C(i,j) = A(i,j)
A(i,j) = aidif*VDC(i,j,k,mt2,bid)/ &
(p5*(DZT(i,j,k ,bid) + DZT(i,j,k+1,bid)))
hfac_t(k) = DZT(i,j,k,bid)/c2dtt(k)
end do
end do
else
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
C(i,j) = A(i,j)
A(i,j) = afac_t(k)*VDC(i,j,k,mt2,bid)
end do
end do
endif ! partial_bottom_cells
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
if (k > KMT(i,j,bid)) then
F(i,j,k) = c0
else
if (k == KMT(i,j,bid)) then
D(i,j) = hfac_t(k)+B(i,j)
else
D(i,j) = hfac_t(k)+A(i,j)+B(i,j)
endif
E(i,j,k) = A(i,j)/D(i,j)
B(i,j) = (hfac_t(k) + B(i,j))*E(i,j,k)
F(i,j,k) = C(i,j)*F(i,j,k-1)/D(i,j)
endif
end do
end do
end do ! k
!*** back substitution
do k=km-1,1,-1
do j=jb,je
do i=ib,ie
if (k < KMT(i,j,bid)) &
F(i,j,k) = F(i,j,k) + E(i,j,k)*F(i,j,k+1)
end do
end do
end do
!-----------------------------------------------------------------------
!
! The solution of the system (UZ) is DeltaTracer -- it has
! already been multiplied by dtt so we can advance the tracers
! directly
!
!-----------------------------------------------------------------------
!CDIR COLLAPSE
do k=1,km
do j=jb,je
do i=ib,ie
TNEW(i,j,k,n) = TNEW(i,j,k,n) + F(i,j,k)
end do
end do
end do
!-----------------------------------------------------------------------
!
! end tracer loop and stop timer
!
!-----------------------------------------------------------------------
end do
call timer_stop(timer_impvmixt,block_id=bid)
!-----------------------------------------------------------------------
!EOC
end subroutine impvmixt_correct
!***********************************************************************
!BOP
! !IROUTINE: impvmixu
! !INTERFACE:
subroutine impvmixu(UNEW, VNEW, this_block) 1,2
! !DESCRIPTION:
! Computes the implicit vertical mixing of momentum
! using a tridiagonal solver.
! Since the top boundary condition (given wind stress) does not
! depend on the velocity, its influence is accounted for in the
! explicit part (subroutine vdiffu). On the other hand, the bottom
! b.c. does depend on velocity ($\tau = c|{\bf\rm u}|(u,v)$).,
! but both the speed and velocity in this term are evaluated
! at the old time, so the bottom drag is also accounted for
! in the explicit part.
!
! !REVISION HISTORY:
! same as module
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km), intent(inout) :: &
UNEW, VNEW ! U,V velocities at new time level
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block info for current block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,k, &! dummy loop indices
ib,ie,jb,je, &! beg,end indices for physical domain
bid ! local block address
real (r8), dimension(nx_block,ny_block) :: &
A,B,C,D ! various temp variables
real (r8), dimension(nx_block,ny_block,km) :: &
E,F1,F2 ! various work arrays
real (r8), dimension(km) :: &
hfac_u
!-----------------------------------------------------------------------
!
! initialize
!
!-----------------------------------------------------------------------
bid = this_block%local_id
ib = this_block%ib
ie = this_block%ie
jb = this_block%jb
je = this_block%je
A=c0
B=c0
C=c0
D=c0
E=c0
F1=c0
F2=c0
call timer_start(timer_impvmixu,block_id=bid)
!*** note that this array is over-written for partial bottom cells
do k=1,km
hfac_u(k) = dz(k)/c2dtu
end do
!-----------------------------------------------------------------------
!
! solve tridiagonal system at each grid point
!
!-----------------------------------------------------------------------
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
A(i,j) = afac_u(1)*VVC(i,j,1,bid)
D(i,j) = hfac_u(1) + A(i,j)
E(i,j,1) = A(i,j)/D(i,j)
B(i,j) = hfac_u(1)*E(i,j,1)
F1(i,j,1) = hfac_u(1)*UNEW(i,j,1)/D(i,j)
F2(i,j,1) = hfac_u(1)*VNEW(i,j,1)/D(i,j)
end do
end do
do k=2,km
if (partial_bottom_cells) then
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
C(i,j) = A(i,j)
hfac_u(k) = DZU(i,j,k,bid)/c2dtu
A(i,j) = aidif*VVC(i,j,k,bid)/(p5*(DZU(i,j,k,bid) + &
DZU(i,j,k+1,bid)))
end do
end do
else
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
C(i,j) = A(i,j)
A(i,j) = afac_u(k)*VVC(i,j,k,bid)
end do
end do
endif ! partial_bottom_cells
!CDIR COLLAPSE
do j=jb,je
do i=ib,ie
if (k < KMU(i,j,bid)) then
D(i,j) = hfac_u(k) + A(i,j) + B(i,j)
E(i,j,k) = A(i,j)/D(i,j)
B(i,j) = (hfac_u(k) + B(i,j))*E(i,j,k)
F1(i,j,k) = (hfac_u(k)*UNEW(i,j,k) + &
C(i,j)*F1(i,j,k-1))/D(i,j)
F2(i,j,k) = (hfac_u(k)*VNEW(i,j,k) + &
C(i,j)*F2(i,j,k-1))/D(i,j)
else if (k == KMU(i,j,bid)) then
D(i,j) = hfac_u(k) + B(i,j)
E(i,j,k) = A(i,j)/D(i,j)
B(i,j) = (hfac_u(k) + B(i,j))*E(i,j,k)
F1(i,j,k) = (hfac_u(k)*UNEW(i,j,k) + &
C(i,j)*F1(i,j,k-1))/D(i,j)
F2(i,j,k) = (hfac_u(k)*VNEW(i,j,k) + &
C(i,j)*F2(i,j,k-1))/D(i,j)
else
F1(i,j,k) = c0
F2(i,j,k) = c0
endif
end do
end do
end do ! k
do k=km-1,1,-1
do j=jb,je
do i=ib,ie
if (k < KMU(i,j,bid)) then
F1(i,j,k) = F1(i,j,k) + E(i,j,k)*F1(i,j,k+1)
F2(i,j,k) = F2(i,j,k) + E(i,j,k)*F2(i,j,k+1)
endif
end do
end do
end do
!CDIR COLLAPSE
do k=1,km
do j=jb,je
do i=ib,ie
UNEW(i,j,k) = F1(i,j,k)
VNEW(i,j,k) = F2(i,j,k)
end do
end do
end do
!-----------------------------------------------------------------------
!
! UVEL,VVEL(newtime) now hold a modified rhs that has already been
! multiplied by dtu to advance UVEL.
!
!-----------------------------------------------------------------------
call timer_stop(timer_impvmixu,block_id=bid)
!-----------------------------------------------------------------------
!EOC
end subroutine impvmixu
!***********************************************************************
!BOP
! !IROUTINE: convad
! !INTERFACE:
subroutine convad(TNEW, RHONEW, this_block, bid) 1,7
! !DESCRIPTION:
! Convectively adjusts tracers at adjacent depth
! levels if they are gravitationally unstable.
! This routine convectively adjusts tracers at the new
! time (before they are updated at the end of step).
!
! !REVISION HISTORY:
! same as module
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(inout) :: &
TNEW ! tracers at new time level
real (r8), dimension(nx_block,ny_block,km), intent(inout) :: &
RHONEW ! density at new time level
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block info for current block
integer (int_kind), intent(in) :: &
bid ! local block index
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
nc, &! dummy index for number of adjustment passes
n, &! dummy index for tracer number
k, ks ! dummy vertical level indices
real (r8), dimension(nx_block,ny_block) :: &
RHOK, &! density of water after adiab. displacement to k+1
RHOKP, &! actual density of water at level k+1
WORK1,WORK2 ! local work space
logical (log_kind) :: &
lpec, lncnv ! logical flags for computing tavg diagnostics
!-----------------------------------------------------------------------
!
! do nothing if convective adjustment not chosen
!
!-----------------------------------------------------------------------
if (convection_itype /= convect_type_adjust) return
!-----------------------------------------------------------------------
!
! perform convective adjustment
!
!-----------------------------------------------------------------------
do nc = 1,nconvad ! repeat nconvad times
do ks = 1,2
do k = ks,km-1,2
call state(k ,k+1,TNEW(:,:,k ,1), TNEW(:,:,k ,2), &
this_block, RHOOUT=RHOK )
call state(k+1,k+1,TNEW(:,:,k+1,1), TNEW(:,:,k+1,2), &
this_block, RHOOUT=RHOKP)
if (partial_bottom_cells) then
do n = 1,nt
where ((RHOK > RHOKP) .and. (k < KMT(:,:,bid)))
TNEW(:,:,k,n) = c1/ &
(DZT(:,:,k ,bid)/dttxcel(k ) + &
DZT(:,:,k+1,bid)/dttxcel(k+1)) &
*(DZT(:,:,k ,bid)/dttxcel(k )*TNEW(:,:,k,n) &
+ DZT(:,:,k+1,bid)/dttxcel(k+1)*TNEW(:,:,k+1,n))
TNEW(:,:,k+1,n) = TNEW(:,:,k,n)
end where
enddo
else
do n = 1,nt
where ((RHOK > RHOKP) .and. (k < KMT(:,:,bid)))
TNEW(:,:,k,n) = dzwxcel(k)* &
(dztxcel(k )*TNEW(:,:,k ,n) &
+ dztxcel(k+1)*TNEW(:,:,k+1,n))
TNEW(:,:,k+1,n) = TNEW(:,:,k,n)
endwhere
enddo
endif
enddo
enddo
enddo
!-----------------------------------------------------------------------
!
! compute new density based on new tracers and compute adjustment
! diagnostic
!
!-----------------------------------------------------------------------
lpec = tavg_requested(tavg_PEC)
lncnv = tavg_requested(tavg_NCNV)
do k = 1,km
if ((lpec .or. lncnv) .and. mix_pass /= 1) then
WORK1 = RHONEW(:,:,k)
endif
call state(k,k,TNEW(:,:,k,1), TNEW(:,:,k,2), &
this_block, RHOOUT=RHONEW(:,:,k))
if (lpec .and. mix_pass /= 1) then
WORK2 = RHONEW(:,:,k) - WORK1
call accumulate_tavg_field(WORK2,tavg_PEC,bid,k)
endif
if (lncnv .and. mix_pass /=1) then
if (.not. lpec) WORK2 = RHONEW(:,:,k) - WORK1
where (abs(WORK2) > 1.e-8_r8)
WORK1 = c1
elsewhere
WORK1 = c0
end where
call accumulate_tavg_field(WORK1,tavg_NCNV,bid,k)
endif
enddo
!-----------------------------------------------------------------------
!EOC
end subroutine convad
!***********************************************************************
end module vertical_mix
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||