!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module hmix_gm 1,19
!BOP
! !MODULE: hmix_gm
! !DESCRIPTION:
! This module contains routines for computing horizontal mixing
! using the Gent-McWilliams eddy transport parameterization
! and isopycnal diffusion.
! !REVISION HISTORY:
! SVN:$Id: hmix_gm.F90 22835 2010-05-07 23:11:58Z njn01 $
! !USES:
use kinds_mod
use blocks
use distribution
use domain
use constants
use broadcast
use grid
use io
use vertical_mix
use vmix_kpp
use state_mod
use time_management
use tavg
use diagnostics
use exit_mod
use registry
use hmix_gm_submeso_share
#ifdef CCSMCOUPLED
use shr_sys_mod
#endif
use timers, only: timer_start, timer_stop, get_timer
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_gm, &
hdifft_gm
!EOP
!BOC
!-----------------------------------------------------------------------
!
! variables to save from one call to next
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
ktp = 1, kbt = 2 ! refer to the top and bottom halves of a
! grid cell, respectively
real (r8), dimension(:), allocatable :: &
kappa_depth ! depth dependence for KAPPA
real (r8), dimension(:,:,:), allocatable :: &
HYXW, HXYS, & ! west and south-shifted values of above
RB, & ! Rossby radius
RBR, & ! inverse of Rossby radius
BTP, & ! beta plane approximation
BL_DEPTH, & ! boundary layer depth
UIT, VIT, & ! work arrays for isopycnal mixing velocities
WTOP_ISOP, WBOT_ISOP ! vertical component of isopycnal velocities
real (r8), dimension(:,:,:,:,:,:), allocatable :: &
SF_SLX, SF_SLY ! components of the merged streamfunction
real (r8), dimension(:,:,:,:,:), allocatable :: &
SLA_SAVE ! isopycnal slopes
real (r8), dimension(:,:,:,:), allocatable :: &
FZTOP ! vertical flux
logical (log_kind), dimension(:), allocatable :: &
compute_kappa ! compute spatially varying coefficients
! this time step?
logical (log_kind) :: &
diff_tapering, & ! different tapering for two diffusivities
cancellation_occurs, & ! specified choices for the isopycnal and
! thickness diffusion coefficients result in
! cancellation of some tensor elements
read_n2_data ! if .true., use climatological N^2 data
! instead of model dependent N^2
!-----------------------------------------------------------------------
!
! KAPPA_LATERAL and KAPPA_VERTICAL are 2D and 3D arrays, respectively,
! containing the spatial variations of the isopycnal (KAPPA_ISOP)
! and thickness (KAPPA_THIC) diffusion coefficients. Except in kappa_type_eg,
! KAPPA_LATERAL has the actual diffusion coefficient values in cm^2/s,
! whereas KAPPA_VERTICAL is unitless. So, the total coefficients are
! constructed using
!
! KAPPA_ISOP or KAPPA_THIC (:,:,:,k,bid) ~ KAPPA_LATERAL(:,:,bid)
! * KAPPA_VERTICAL(:,:,k,bid)
!
! When kappa_type_eg, KAPPA_VERTICAL contains the isopycnal diffusivity
! coefficients in cm^2/s and KAPPA_LATERAL is not used!
!
!-----------------------------------------------------------------------
real (r8), dimension(:,:,:,:,:), allocatable :: &
KAPPA_ISOP, & ! 3D isopycnal diffusion coefficient
! for top and bottom half of a grid cell
KAPPA_THIC, & ! 3D thickness diffusion coefficient
! for top and bottom half of a grid cell
HOR_DIFF ! 3D horizontal diffusion coefficient
! for top and bottom half of a grid cell
real (r8), dimension(:,:,:), allocatable :: &
KAPPA_LATERAL ! horizontal variation of KAPPA in cm^2/s
real (r8), dimension(:,:,:,:), allocatable :: &
KAPPA_VERTICAL ! vertical variation of KAPPA (unitless),
! e.g. normalized buoyancy frequency dependent
! profiles at the tracer grid points
! ( = N^2 / N_ref^2 ) OR a time-independent
! user-specified function
real (r8), dimension(:,:,:,:), allocatable :: &
BUOY_FREQ_SQ, & ! N^2 defined at level interfaces
SIGMA_TOPO_MASK ! bottom topography mask used with kappa_type_eg
!-----------------------------------------------------------------------
!
! GM specific options
!
! kappa_freq = how often spatial variations of the diffusion
! coefficients are computed. Same frequency is
! used for both coefficients.
! slope_control = tanh function (Danabasoglu and McWilliams 1995) or
! DM95 with replacement function to tanh or
! slope clipping or
! method of Gerdes et al (1991)
! diag_gm_bolus = .true. for diagnostic bolus velocity computation.
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
kappa_type_const = 1, &
kappa_type_depth = 2, &
kappa_type_vmhs = 3, &
kappa_type_hdgr = 4, &
kappa_type_dradius = 5, &
kappa_type_bfreq = 6, &
kappa_type_bfreq_vmhs = 7, &
kappa_type_bfreq_hdgr = 8, &
kappa_type_bfreq_dradius = 9, &
kappa_type_eg = 10, &
slope_control_tanh = 1, &
slope_control_notanh = 2, &
slope_control_clip = 3, &
slope_control_Gerd = 4, &
kappa_freq_never = 1, &
kappa_freq_every_time_step = 2, &
kappa_freq_once_a_day = 3
integer (int_kind) :: &
kappa_isop_type, & ! choice of KAPPA_ISOP
kappa_thic_type, & ! choice of KAPPA_THIC
kappa_freq, & ! frequency of KAPPA computations
slope_control ! choice for slope control
logical (log_kind) :: &
diag_gm_bolus ! true for diagnostic bolus velocity computation
!-----------------------------------------------------------------------
!
! if use_const_ah_bkg_srfbl = .true., then the specified constant
! value of ah_bkg_srfbl is used as the "maximum" background horizontal
! diffusivity within the surface boundary layer. Otherwise,
! KAPPA_ISOP is utilized as this "maximum".
!
!-----------------------------------------------------------------------
logical (log_kind) :: &
use_const_ah_bkg_srfbl, & ! see above
transition_layer_on ! control for transition layer parameterization
real (r8) :: &
ah, & ! isopycnal diffusivity
ah_bolus, & ! thickness (GM bolus) diffusivity
ah_bkg_bottom, & ! backgroud horizontal diffusivity at k = KMT
ah_bkg_srfbl, & ! backgroud horizontal diffusivity within the
! surface boundary layer
slm_r, & ! max. slope allowed for redi diffusion
slm_b ! max. slope allowed for bolus transport
!-----------------------------------------------------------------------
!
! the following set of variables are used in Eden and Greatbatch
! (2008) KAPPA formulation. They are in the input namelist.
!
!-----------------------------------------------------------------------
real (r8) :: &
const_eg, & ! tuning parameter (unitless)
gamma_eg, & ! (> 0) effective upper limit for inverse eddy
! time scale (unitless)
kappa_min_eg, & ! minimum KAPPA (cm^2/s)
kappa_max_eg ! maximum KAPPA (cm^2/s)
!-----------------------------------------------------------------------
!
! transition layer type variables
!
!-----------------------------------------------------------------------
type tlt_info
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
DIABATIC_DEPTH, & ! depth of the diabatic region at the
! surface, i.e. mean mixed or boundary layer
! depth
THICKNESS, & ! transition layer thickness
INTERIOR_DEPTH ! depth at which the interior, adiabatic
! region starts, i.e.
! = TLT%DIABATIC_DEPTH + TLT%THICKNESS
integer (int_kind), &
dimension(nx_block,ny_block,max_blocks_clinic) :: &
K_LEVEL, & ! k level at or below which the interior,
! adiabatic region starts
ZTW ! designates if the interior region
! starts below depth zt or zw.
! ( = 1 for zt, = 2 for zw )
end type tlt_info
type (tlt_info) :: &
TLT ! transition layer thickness related fields
!-----------------------------------------------------------------------
!
! tavg ids for tavg diagnostics related to diffusivities and
! isopycnal velocities. Zonal and meridional refer here to logical
! space only.
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_UISOP, & ! zonal isopycnal velocity
tavg_VISOP, & ! meridional isopycnal velocity
tavg_WISOP, & ! vertical isopycnal velocity
tavg_KAPPA_ISOP, & ! isopycnal diffusion coefficient
tavg_KAPPA_THIC, & ! thickness diffusion coefficient
tavg_HOR_DIFF, & ! horizontal diffusion coefficient
tavg_DIA_DEPTH, & ! depth of the diabatic region at the surface
tavg_TLT, & ! transition layer thickness
tavg_INT_DEPTH, & ! depth at which the interior region starts
tavg_ADVT_ISOP, & ! vertically-integrated T eddy-induced
! advection tendency
tavg_ADVS_ISOP, & ! vertically-integrated S eddy-induced
! advection tendency
tavg_VNT_ISOP, & ! heat flux tendency in grid-y direction
! due to eddy-induced velocity
tavg_VNS_ISOP ! salt flux tendency in grid-y direction
! due to eddy-induced velocity
!-----------------------------------------------------------------------
!
! timers
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
timer_nloop ! main n loop
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_gm
! !INTERFACE:
subroutine init_gm 1,67
! !DESCRIPTION:
! Initializes various choices and allocates necessary space.
! Also computes some time-independent arrays.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
type (block) :: &
this_block ! block information for current block
character (char_len) :: &
buoyancy_freq_filename, &! file name for the time-independent
! buoyancy frequency squared
buoyancy_freq_fmt, &! format (bin or netcdf) of buoyancy file
message ! string to hold error message
type (datafile) :: &
buoyancy_freq_data_file ! data file descriptor for buoyancy freq data
type (io_field_desc) :: &
buoyancy_freq_data_in ! io field descriptor for input buoyancy freq data
type (io_dim) :: &
i_dim, j_dim, & ! dimension descriptors for horiz dims
k_dim ! dimension descriptor for depth
!-----------------------------------------------------------------------
!
! input namelist for setting various GM options
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
nml_error, & ! error flag for namelist
k, & ! vertical level index
iblock, & ! block index
i, j ! lateral indices
real (r8) :: &
kappa_depth_1, & ! parameters for variation of KAPPA
kappa_depth_2, & ! with depth used with the kappa_type_depth
kappa_depth_scale ! option
character (char_len) :: &
kappa_choice, & ! (generic) choice for KAPPA
kappa_freq_choice, & ! frequency choice for KAPPA computation
kappa_isop_choice, & ! choice for KAPPA_ISOP
kappa_thic_choice, & ! choice for KAPPA_THIC
slope_control_choice ! choice for slope control
namelist /hmix_gm_nml/ kappa_isop_choice, &
kappa_thic_choice, &
kappa_freq_choice, &
slope_control_choice, &
kappa_depth_1, kappa_depth_2, &
kappa_depth_scale, ah, ah_bolus, &
use_const_ah_bkg_srfbl, ah_bkg_srfbl, &
ah_bkg_bottom, &
slm_r, slm_b, diag_gm_bolus, &
transition_layer_on, read_n2_data, &
buoyancy_freq_filename, &
buoyancy_freq_fmt, &
const_eg, gamma_eg, &
kappa_min_eg, kappa_max_eg
!-----------------------------------------------------------------------
!
! read input namelist for additional GM options
!
! DEFAULT SETUP:
! KAPPA_ISOP : constant (= ah)
! KAPPA_THIC : constant (= ah_bolus)
! kappa_freq : never
! slope control : method by DM95 with replacing tanh by polynomial
! (= slope_control_notanh)
!
! variation of kappa with depth used with kappa_type_depth is
! kappa_depth_1 + kappa_depth_2*exp(-z/kappa_depth_scale)
!
! with kappa_depth_1 = 1, kappa_depth_2 = 0, and
! kappa_depth_scale = 150000 cm, i.e. no depth variation
!
! ah_bolus : thickness diffusion coefficient (= 0.8e07 cm^2/s)
! ah_bkg_bottom : background horizontal diffusion at k=KMT (= 0)
! use_const_ah_bkg_srfbl : use ah_bkg_srfbl value
! ah_bkg_srfbl : background horizontal diffusion within the
! surface boundary layer (= 0.8e07 cm^2/s)
! slm_r : max. slope allowed for isopycnal (redi) diffusion (= 0.3)
! slm_b : max. slope allowed for thickness (bolus) diffusion (= 0.3)
! diag_gm_bolus : .true.
! transition_layer_on: .false.
! read_n2_data : .false.
! const_eg : 1.
! gamma_eg : 300.
! kappa_min_eg : 0.35e07 cm^2/s
! kappa_max_eg : 5.00e07 cm^2/s
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! register init_gm
!
!-----------------------------------------------------------------------
call register_string('init_gm')
kappa_isop_type = kappa_type_const
kappa_thic_type = kappa_type_const
kappa_freq = kappa_freq_never
slope_control = slope_control_notanh
kappa_depth_1 = c1
kappa_depth_2 = c0
kappa_depth_scale = 150000.0_r8
ah = 0.8e7_r8
ah_bolus = 0.8e7_r8
ah_bkg_bottom = c0
ah_bkg_srfbl = 0.8e7_r8
slm_r = 0.3_r8
slm_b = 0.3_r8
diag_gm_bolus = .true.
use_const_ah_bkg_srfbl = .true.
transition_layer_on = .false.
read_n2_data = .false.
buoyancy_freq_filename = 'unknown-buoyancy'
buoyancy_freq_fmt = 'nc'
const_eg = c1
gamma_eg = 300.0_r8
kappa_min_eg = 0.35e7_r8
kappa_max_eg = 5.0e7_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=hmix_gm_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 hmix_gm_nml')
endif
if (my_task == master_task) then
write(stdout,*) ' '
write(stdout,*) ' Document Namelist Parameters:'
write(stdout,*) ' ============================ '
write(stdout,*) ' '
write(stdout, hmix_gm_nml)
write(stdout,*) ' '
write(stdout,*) ' Gent-McWilliams options:'
write(stdout,*) ' kappa_isop choice is ', &
trim(kappa_isop_choice)
write(stdout,*) ' kappa_thic choice is ', &
trim(kappa_thic_choice)
write(stdout,*) ' kappa computation frequency choice is ', &
trim(kappa_freq_choice)
write(stdout,*) ' slope control choice is ', &
trim(slope_control_choice)
write(stdout,'(a28,1pe13.6)') ' isopycnal diffusion = ', &
ah
write(stdout,'(a28,1pe13.6)') ' thickness diffusion = ', &
ah_bolus
write(stdout,'(a28,1pe13.6)') ' backgroud diff. (bottom) = ', &
ah_bkg_bottom
write(stdout,'(a28,1pe13.6)') ' backgroud diff. (srfbl) = ', &
ah_bkg_srfbl
write(stdout,'(a28,1pe13.6)') ' max. slope for redi = ', &
slm_r
write(stdout,'(a28,1pe13.6)') ' max. slope for bolus = ', &
slm_b
if ( diag_gm_bolus ) then
write(stdout,'(a47)') &
' diagnostic bolus velocity computation is on'
endif
if ( use_const_ah_bkg_srfbl ) then
write(stdout,1001)
else
write(stdout,1002)
endif
1001 format(/,' Maximum horizontal background diffusion ', &
'coefficient', /, &
' within the boundary layer is constant at ah_bkg_srfbl.')
1002 format(/,' Maximum horizontal background diffusion ', &
'coefficient', /, &
' within the boundary layer depends on KAPPA_ISOP.')
if ( transition_layer_on ) then
write(stdout,'(a33)') ' transition layer scheme is on. '
endif
select case (kappa_isop_choice(1:4))
case ('cons')
kappa_isop_type = kappa_type_const
case ('dept')
kappa_isop_type = kappa_type_depth
case ('vmhs')
kappa_isop_type = kappa_type_vmhs
case ('hdgr')
kappa_isop_type = kappa_type_hdgr
case ('drad')
kappa_isop_type = kappa_type_dradius
case ('bfre')
kappa_isop_type = kappa_type_bfreq
case ('bfvm')
kappa_isop_type = kappa_type_bfreq_vmhs
case ('bfhd')
kappa_isop_type = kappa_type_bfreq_hdgr
case ('bfdr')
kappa_isop_type = kappa_type_bfreq_dradius
case ('edgr')
kappa_isop_type = kappa_type_eg
case default
kappa_isop_type = -1000
end select
select case (kappa_thic_choice(1:4))
case ('cons')
kappa_thic_type = kappa_type_const
case ('dept')
kappa_thic_type = kappa_type_depth
case ('vmhs')
kappa_thic_type = kappa_type_vmhs
case ('hdgr')
kappa_thic_type = kappa_type_hdgr
case ('drad')
kappa_thic_type = kappa_type_dradius
case ('bfre')
kappa_thic_type = kappa_type_bfreq
case ('bfvm')
kappa_thic_type = kappa_type_bfreq_vmhs
case ('bfhd')
kappa_thic_type = kappa_type_bfreq_hdgr
case ('bfdr')
kappa_thic_type = kappa_type_bfreq_dradius
case ('edgr')
kappa_thic_type = kappa_type_eg
case default
kappa_thic_type = -1000
end select
select case (kappa_freq_choice(1:3))
case ('nev')
kappa_freq = kappa_freq_never
case ('eve')
kappa_freq = kappa_freq_every_time_step
case ('onc')
kappa_freq = kappa_freq_once_a_day
case default
kappa_freq = -1000
end select
select case (slope_control_choice(1:4))
case ('tanh')
slope_control = slope_control_tanh
case ('nota')
slope_control = slope_control_notanh
case ('clip')
slope_control = slope_control_clip
case ('Gerd')
slope_control = slope_control_Gerd
case default
slope_control = -1000
end select
if ( kappa_isop_type == kappa_type_bfreq .or. &
kappa_isop_type == kappa_type_bfreq_vmhs .or. &
kappa_isop_type == kappa_type_bfreq_hdgr .or. &
kappa_isop_type == kappa_type_bfreq_dradius .or. &
kappa_thic_type == kappa_type_bfreq .or. &
kappa_thic_type == kappa_type_bfreq_vmhs .or. &
kappa_thic_type == kappa_type_bfreq_hdgr .or. &
kappa_thic_type == kappa_type_bfreq_dradius ) then
if ( read_n2_data ) then
write(stdout,'(a32)') ' using climatological N^2 data. '
else
write(stdout,'(a34)') ' using model data to compute N^2. '
endif
endif
if ( kappa_isop_type == kappa_type_depth .or. &
kappa_thic_type == kappa_type_depth ) then
if ( kappa_isop_type == kappa_type_depth ) &
write(stdout,'(a30)') ' KAPPA_ISOP varies with depth.'
if ( kappa_thic_type == kappa_type_depth ) &
write(stdout,'(a30)') ' KAPPA_THIC varies with depth.'
write(stdout,'(a20,1pe13.6)') ' kappa_depth_1 = ', &
kappa_depth_1
write(stdout,'(a20,1pe13.6)') ' kappa_depth_2 = ', &
kappa_depth_2
write(stdout,'(a24,1pe13.6)') ' kappa_depth_scale = ', &
kappa_depth_scale
endif
if ( kappa_isop_type == kappa_type_eg .or. &
kappa_thic_type == kappa_type_eg ) then
write(stdout,'(a16,1pe13.6)') ' const_eg = ', &
const_eg
write(stdout,'(a16,1pe13.6)') ' gamma_eg = ', &
gamma_eg
write(stdout,'(a16,1pe13.6)') ' kappa_min_eg = ', &
kappa_min_eg
write(stdout,'(a16,1pe13.6)') ' kappa_max_eg = ', &
kappa_max_eg
endif
endif
call broadcast_scalar(kappa_isop_type, master_task)
call broadcast_scalar(kappa_thic_type, master_task)
call broadcast_scalar(kappa_freq, master_task)
call broadcast_scalar(slope_control, master_task)
call broadcast_scalar(kappa_depth_1, master_task)
call broadcast_scalar(kappa_depth_2, master_task)
call broadcast_scalar(kappa_depth_scale, master_task)
call broadcast_scalar(ah, master_task)
call broadcast_scalar(ah_bolus, master_task)
call broadcast_scalar(ah_bkg_bottom, master_task)
call broadcast_scalar(ah_bkg_srfbl, master_task)
call broadcast_scalar(slm_r, master_task)
call broadcast_scalar(slm_b, master_task)
call broadcast_scalar(diag_gm_bolus, master_task)
call broadcast_scalar(use_const_ah_bkg_srfbl, master_task)
call broadcast_scalar(transition_layer_on, master_task)
call broadcast_scalar(read_n2_data, master_task)
call broadcast_scalar(buoyancy_freq_filename, master_task)
call broadcast_scalar(buoyancy_freq_fmt, master_task)
call broadcast_scalar(const_eg, master_task)
call broadcast_scalar(gamma_eg, master_task)
call broadcast_scalar(kappa_min_eg, master_task)
call broadcast_scalar(kappa_max_eg, master_task)
!-----------------------------------------------------------------------
!
! error checking
!
!-----------------------------------------------------------------------
if ( kappa_isop_type == -1000 ) then
call exit_POP(sigAbort, &
'unknown type for KAPPA_ISOP in GM setup')
endif
if ( kappa_thic_type == -1000 ) then
call exit_POP(sigAbort, &
'unknown type for KAPPA_THIC in GM setup')
endif
if ( kappa_freq == -1000 ) then
call exit_POP(sigAbort, &
'unknown type for KAPPA computation frequency in GM setup')
endif
if ( slope_control == -1000 ) then
call exit_POP(sigAbort, &
'unknown slope control method in GM setup')
endif
if ( kappa_isop_type == kappa_type_depth .and. &
( kappa_thic_type /= kappa_type_const .and. &
kappa_thic_type /= kappa_type_depth ) ) then
message = 'when kappa_isop_type is kappa_type_depth, kappa_thic_type '/&
&/ 'should be either kappa_type_const or kappa_type_depth.'
call exit_POP(sigAbort, message)
endif
if ( kappa_thic_type == kappa_type_depth .and. &
( kappa_isop_type /= kappa_type_const .and. &
kappa_isop_type /= kappa_type_depth ) ) then
message = 'when kappa_thic_type is kappa_type_depth, kappa_isop_type ' /&
&/ 'should be either kappa_type_const or kappa_type_depth.'
call exit_POP(sigAbort, message)
endif
if ( ( ( kappa_isop_type /= kappa_type_const .and. &
kappa_isop_type /= kappa_type_depth ) .and. &
( kappa_thic_type /= kappa_type_const .and. &
kappa_thic_type /= kappa_type_depth ) ) .and. &
( kappa_isop_type /= kappa_thic_type ) ) then
message = 'kappa_isop_type and kappa_thic_type should be the same ' /&
&/ 'if they are BOTH model fields dependent.'
call exit_POP(sigAbort, message)
endif
if ( ( ( kappa_isop_type == kappa_type_vmhs .and. &
kappa_thic_type == kappa_type_vmhs ) .or. &
( kappa_isop_type == kappa_type_bfreq_vmhs .and. &
kappa_thic_type == kappa_type_bfreq_vmhs ) .or. &
( kappa_isop_type == kappa_type_dradius .and. &
kappa_thic_type == kappa_type_dradius ) .or. &
( kappa_isop_type == kappa_type_bfreq_dradius .and. &
kappa_thic_type == kappa_type_bfreq_dradius ) ) .and. &
ah /= ah_bolus ) then
message = 'ah and ah_bolus should be equal when vmhs or dradius ' /&
&/ ' related dependencies are used.'
call exit_POP(sigAbort, message)
endif
if ( ( kappa_isop_type == kappa_type_depth .or. &
kappa_thic_type == kappa_type_depth ) .and. &
kappa_depth_2 == c0 ) then
message = 'kappa_depth_2 should be non-zero if a time-independent ' /&
&/ 'depth variation is requested.'
call exit_POP(sigAbort, message)
endif
if ( ( kappa_isop_type == kappa_type_vmhs .or. &
kappa_thic_type == kappa_type_vmhs .or. &
kappa_isop_type == kappa_type_bfreq_vmhs .or. &
kappa_thic_type == kappa_type_bfreq_vmhs .or. &
kappa_isop_type == kappa_type_hdgr .or. &
kappa_thic_type == kappa_type_hdgr .or. &
kappa_isop_type == kappa_type_bfreq_hdgr .or. &
kappa_thic_type == kappa_type_bfreq_hdgr ) .and. &
zt(km) <= 2.0e5_r8 ) then
message = 'the max bottom depth should be > 2000.0 m' /&
&/ ' when vmhs or hdgr dependency is chosen.'
call exit_POP(sigAbort, message)
endif
if ( ( kappa_isop_type == kappa_type_eg .or. &
kappa_thic_type == kappa_type_eg ) .and. &
use_const_ah_bkg_srfbl ) then
message = 'use_const_ah_bkg_srfbl should be false when ' /&
&/ 'kappa_type_eg is used.'
call exit_POP(sigAbort, message)
endif
if ( ( kappa_isop_type == kappa_type_vmhs .or. &
kappa_thic_type == kappa_type_vmhs .or. &
kappa_isop_type == kappa_type_hdgr .or. &
kappa_thic_type == kappa_type_hdgr .or. &
kappa_isop_type == kappa_type_dradius .or. &
kappa_thic_type == kappa_type_dradius .or. &
kappa_isop_type == kappa_type_bfreq_vmhs .or. &
kappa_thic_type == kappa_type_bfreq_vmhs .or. &
kappa_isop_type == kappa_type_bfreq_hdgr .or. &
kappa_thic_type == kappa_type_bfreq_hdgr .or. &
kappa_isop_type == kappa_type_bfreq_dradius .or. &
kappa_thic_type == kappa_type_bfreq_dradius .or. &
kappa_isop_type == kappa_type_eg .or. &
kappa_thic_type == kappa_type_eg .or. &
( ( kappa_isop_type == kappa_type_bfreq .or. &
kappa_thic_type == kappa_type_bfreq ) .and. &
.not. read_n2_data ) ) .and. &
kappa_freq == kappa_freq_never ) then
message = 'kappa_freq should not be set to never when model ' /&
&/ 'fields dependent kappa types are chosen.'
call exit_POP(sigAbort, message)
endif
if ( partial_bottom_cells ) then
call exit_POP(sigAbort, &
'hmix_gm currently incompatible with partial bottom cells')
endif
if (read_n2_data .and. trim(buoyancy_freq_filename) == 'unknown-buoyancy' ) then
call exit_POP(sigAbort, &
'Must define buoyancy_freq_filename if read_n2_data is .true.')
endif
!-----------------------------------------------------------------------
!
! allocate GM arrays
!
!-----------------------------------------------------------------------
allocate (HYXW(nx_block,ny_block,nblocks_clinic), &
HXYS(nx_block,ny_block,nblocks_clinic), &
RBR (nx_block,ny_block,nblocks_clinic), &
BTP (nx_block,ny_block,nblocks_clinic), &
BL_DEPTH(nx_block,ny_block,nblocks_clinic))
allocate (SF_SLX(nx_block,ny_block,2,2,km,nblocks_clinic), &
SF_SLY(nx_block,ny_block,2,2,km,nblocks_clinic))
allocate (FZTOP(nx_block,ny_block,nt,nblocks_clinic))
allocate (kappa_depth(km))
allocate (KAPPA_ISOP(nx_block,ny_block,2,km,nblocks_clinic), &
KAPPA_THIC(nx_block,ny_block,2,km,nblocks_clinic), &
HOR_DIFF (nx_block,ny_block,2,km,nblocks_clinic))
allocate (KAPPA_LATERAL (nx_block,ny_block,nblocks_clinic), &
KAPPA_VERTICAL(nx_block,ny_block,km,nblocks_clinic))
allocate (BUOY_FREQ_SQ(nx_block,ny_block,km,nblocks_clinic))
allocate (VDC_GM(nx_block,ny_block,km,nblocks_clinic))
allocate (compute_kappa(nblocks_clinic))
HYXW = c0
HXYS = c0
RBR = c0
BTP = c0
BL_DEPTH = c0
SF_SLX = c0
SF_SLY = c0
FZTOP = c0
VDC_GM = c0
if ( transition_layer_on ) then
allocate (SLA_SAVE(nx_block,ny_block,2,km,nblocks_clinic))
allocate (RB(nx_block,ny_block,nblocks_clinic))
SLA_SAVE = c0
RB = c0
endif
!-----------------------------------------------------------------------
!
! initialize various time-independent arrays
!
!-----------------------------------------------------------------------
do k=1,km
kappa_depth(k) = kappa_depth_1 &
+ kappa_depth_2 &
*exp(-zt(k)/kappa_depth_scale)
enddo
do iblock = 1,nblocks_clinic
this_block = get_block(blocks_clinic(iblock),iblock)
KAPPA_LATERAL(:,:,iblock) = ah
KAPPA_VERTICAL(:,:,:,iblock) = c1
KAPPA_ISOP(:,:,:,:,iblock) = ah
KAPPA_THIC(:,:,:,:,iblock) = ah_bolus
HOR_DIFF(:,:,:,:,iblock) = ah
if ( kappa_isop_type == kappa_type_depth .or. &
kappa_thic_type == kappa_type_depth ) then
do k=1,km
KAPPA_VERTICAL(:,:,k,iblock) = kappa_depth(k)
enddo
endif
HYXW(:,:,iblock) = eoshift(HYX(:,:,iblock), dim=1, shift=-1)
HXYS(:,:,iblock) = eoshift(HXY(:,:,iblock), dim=2, shift=-1)
!-----------------------------------------------------------------------
! compute the Rossby radius which will be used to
! contol KAPPA [Large et al (1997), JPO, 27,
! pp 2418-2447]. Rossby radius = c/(2*omg*sin(latitude))
! where c=200cm/s is the first baroclinic wave speed.
! 15km < Rossby radius < 100km
!-----------------------------------------------------------------------
!*** Inverse of Rossby radius
RBR(:,:,iblock) = abs(FCORT(:,:,iblock)) &
/ 200.0_r8 ! |f|/Cg, Cg = 2 m/s = 200 cm/s
RBR(:,:,iblock) = min(RBR(:,:,iblock), &
c1/1.5e+6_r8) ! Cg/|f| .ge. 15 km = 1.5e+6 cm
RBR(:,:,iblock) = max(RBR(:,:,iblock), &
1.e-7_r8) ! Cg/|f| .le. 100 km = 1.e+7 cm
if ( transition_layer_on ) then
RB(:,:,iblock) = c1 / RBR(:,:,iblock)
endif
!*** beta at t-points
call ugrid_to_tgrid(BTP(:,:,iblock),ULAT(:,:,iblock),iblock)
BTP(:,:,iblock) = c2*omega*cos(BTP(:,:,iblock))/radius
enddo
!-----------------------------------------------------------------------
! HYXW, HXYS only needed in physical domain and should
! be defined correctly there. BTP invalid on south
! and westernmost ghost cells, but not needed there
! as long as number of ghost cells is >1
!-----------------------------------------------------------------------
! call update_ghost_cells(HYXW, bndy_clinic, field_loc_t, &
! field_type_scalar)
! call update_ghost_cells(HXYS, bndy_clinic, field_loc_t, &
! field_type_scalar)
! call update_ghost_cells(BTP , bndy_clinic, field_loc_t, &
! field_type_scalar)
BUOY_FREQ_SQ = c0
if ( read_n2_data ) then
buoyancy_freq_filename = '/home/bluesky/gokhan/buoyancy_freq.nc'
buoyancy_freq_data_file = construct_file ( 'nc', &
full_name=trim(buoyancy_freq_filename) )
call data_set(buoyancy_freq_data_file, 'open_read')
i_dim = construct_io_dim ( 'i', nx_global )
j_dim = construct_io_dim ( 'j', ny_global )
k_dim = construct_io_dim ( 'k', km )
buoyancy_freq_data_in = construct_io_field &
('BUOYANCY_FREQUENCY', &
dim1=i_dim, dim2=j_dim, dim3=k_dim, &
field_loc = field_loc_center, &
field_type = field_type_scalar, &
d3d_array = BUOY_FREQ_SQ)
call data_set (buoyancy_freq_data_file, 'define', &
buoyancy_freq_data_in)
call data_set (buoyancy_freq_data_file, 'read', &
buoyancy_freq_data_in)
call data_set (buoyancy_freq_data_file, 'close')
call destroy_io_field (buoyancy_freq_data_in)
call destroy_file (buoyancy_freq_data_file)
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,'(a43,a)') &
' Buoyancy frequency climatology file read: ', &
trim(buoyancy_freq_filename)
endif
endif
compute_kappa = .false.
if (slm_r /= slm_b) then
diff_tapering = .true.
else
diff_tapering = .false.
endif
cancellation_occurs = .true.
if ( ( kappa_isop_type /= kappa_thic_type ) .or. &
( diff_tapering ) .or. &
( ( kappa_isop_type == kappa_type_const .and. &
kappa_thic_type == kappa_type_const ) .and. &
ah /= ah_bolus ) .or. &
( ( kappa_isop_type == kappa_type_depth .and. &
kappa_thic_type == kappa_type_depth ) .and. &
ah /= ah_bolus ) .or. &
( ( kappa_isop_type == kappa_type_bfreq .and. &
kappa_thic_type == kappa_type_bfreq ) .and. &
ah /= ah_bolus ) ) &
cancellation_occurs = .false.
!*** for transition layer cases, the following will always be true!!
if ( transition_layer_on ) cancellation_occurs = .false.
!-----------------------------------------------------------------------
!
! initialize topography mask used with kappa_type_eg. This mask eliminates
! excessively large values of KAPPA near sloping topography.
!
!-----------------------------------------------------------------------
if ( kappa_isop_type == kappa_type_eg .or. &
kappa_thic_type == kappa_type_eg ) then
allocate (SIGMA_TOPO_MASK(nx_block,ny_block,km,nblocks_clinic))
do iblock=1,nblocks_clinic
do k=1,km
where ( k < KMT(:,:,iblock) )
SIGMA_TOPO_MASK(:,:,k,iblock) = c1
elsewhere
SIGMA_TOPO_MASK(:,:,k,iblock) = c0
endwhere
enddo
do k=1,km-1
do j=2,ny_block-1
do i=2,nx_block-1
if ( k < KMT(i,j,iblock) ) then
if ( k == KMT(i-1,j+1,iblock) .or. &
k == KMT(i ,j+1,iblock) .or. &
k == KMT(i+1,j+1,iblock) .or. &
k == KMT(i-1,j ,iblock) .or. &
k == KMT(i+1,j ,iblock) .or. &
k == KMT(i-1,j-1,iblock) .or. &
k == KMT(i ,j-1,iblock) .or. &
k == KMT(i+1,j-1,iblock) ) &
SIGMA_TOPO_MASK(i,j,k,iblock) = c0
endif
enddo
enddo
enddo
enddo
endif
!-----------------------------------------------------------------------
!
! define tavg fields related to diffusivities
!
!-----------------------------------------------------------------------
call define_tavg_field (tavg_KAPPA_ISOP, 'KAPPA_ISOP', 3, &
long_name='Isopycnal diffusion coefficient', &
units='cm^2/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field (tavg_KAPPA_THIC, 'KAPPA_THIC', 3, &
long_name='Thickness diffusion coefficient', &
units='cm^2/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field (tavg_HOR_DIFF, 'HOR_DIFF', 3, &
long_name='Horizontal diffusion coefficient', &
units='cm^2/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
!-----------------------------------------------------------------------
!
! define tavg fields related to the near-surface eddy flux
! parameterization (i.e. transition layer is on).
!
!-----------------------------------------------------------------------
if ( transition_layer_on ) then
call define_tavg_field (tavg_DIA_DEPTH, 'DIA_DEPTH', 2, &
long_name='Depth of the Diabatic Region at the Surface', &
units='cm', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field (tavg_TLT, 'TLT', 2, &
long_name='Transition Layer Thickness', &
units='cm', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field (tavg_INT_DEPTH, 'INT_DEPTH', 2, &
long_name='Depth at which the Interior Region Starts', &
units='cm', grid_loc='2110', &
coordinates='TLONG TLAT time')
endif
!-----------------------------------------------------------------------
!
! allocate and initialize bolus velocity arrays if requested
! define tavg fields related to bolus velocity
!
!-----------------------------------------------------------------------
if ( diag_gm_bolus ) then
call register_string ('diag_gm_bolus')
allocate(WTOP_ISOP(nx_block,ny_block,nblocks_clinic), &
WBOT_ISOP(nx_block,ny_block,nblocks_clinic), &
UIT(nx_block,ny_block,nblocks_clinic), &
VIT(nx_block,ny_block,nblocks_clinic))
WTOP_ISOP = c0
WBOT_ISOP = c0
UIT = c0
VIT = c0
call define_tavg_field (tavg_UISOP, 'UISOP', 3, &
long_name='Bolus Velocity in grid-x direction (diagnostic)', &
units='cm/s', grid_loc='3211', &
coordinates='ULONG TLAT z_t time')
call define_tavg_field (tavg_VISOP, 'VISOP', 3, &
long_name='Bolus Velocity in grid-y direction (diagnostic)', &
units='cm/s', grid_loc='3121', &
coordinates='TLONG ULAT z_t time')
call define_tavg_field (tavg_WISOP, 'WISOP', 3, &
long_name='Vertical Bolus Velocity (diagnostic)', &
units='cm/s', grid_loc='3112', &
coordinates='TLONG TLAT z_w time')
call define_tavg_field (tavg_ADVT_ISOP, 'ADVT_ISOP', 2, &
long_name='Vertically-Integrated T Eddy-Induced Advection Tendency (diagnostic)', &
units='cm degC/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field (tavg_ADVS_ISOP, 'ADVS_ISOP', 2, &
long_name='Vertically-Integrated S Eddy-Induced Advection Tendency (diagnostic)', &
scale_factor=1000.0_rtavg, &
units='cm gram/kilogram/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field (tavg_VNT_ISOP, 'VNT_ISOP', 3, &
long_name='Heat Flux Tendency in grid-y Dir due to Eddy-Induced Vel (diagnostic)', &
units='degC/s', grid_loc='3121', &
coordinates='TLONG ULAT z_t time')
call define_tavg_field (tavg_VNS_ISOP, 'VNS_ISOP', 3, &
long_name='Salt Flux Tendency in grid-y Dir due to Eddy-Induced Vel (diagnostic)', &
scale_factor=1000.0_rtavg, &
units='gram/kilogram/s', grid_loc='3121', &
coordinates='TLONG ULAT z_t time')
endif
call get_timer(timer_nloop,'HMIX_TRACER_GM_NLOOP', &
nblocks_clinic, distrb_clinic%nprocs)
!-----------------------------------------------------------------------
!EOC
end subroutine init_gm
!***********************************************************************
!BOP
! !IROUTINE: hdifft_gm
! !INTERFACE:
subroutine hdifft_gm (k, GTK, TMIX, UMIX, VMIX, tavg_HDIFE_TRACER, & 1,63
tavg_HDIFN_TRACER, tavg_HDIFB_TRACER, this_block)
! !DESCRIPTION:
! Gent-McWilliams eddy transport parameterization
! and isopycnal diffusion.
!
! This routine must be called successively with k = 1,2,3,...
!
! !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 at all vertical levels
! at mixing time level
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UMIX, VMIX ! U,V at all vertical levels
! at mixing time level
integer (int_kind), dimension(nt), intent(in) :: &
tavg_HDIFE_TRACER, &! tavg id for east face diffusive flux of tracer
tavg_HDIFN_TRACER, &! tavg id for north face diffusive flux of tracer
tavg_HDIFB_TRACER ! tavg id for bottom face diffusive flux of tracer
type (block), intent(in) :: &
this_block ! block info for this sub block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(out) :: &
GTK ! diffusion+bolus advection for nth tracer at level k
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
ieast = 1, iwest = 2, &
jnorth = 1, jsouth = 2
integer (int_kind) :: &
i,j,n,kk, &! dummy loop counters
kid,ktmp, &! array indices
kk_sub, kp1, &
bid, &! local block address for this sub block
tavg_stream ! temporary stream id
real (r8) :: &
fz, dz_bottom, factor
real (r8), dimension(nx_block,ny_block) :: &
CX, CY, &
RZ, &! Dz(rho)
SLA, &! absolute value of slope
WORK1, WORK2, &! local work space
WORK3, WORK4, &! local work space
KMASK, &! ocean mask
TAPER1, TAPER2, TAPER3, &! tapering factors
UIB, VIB, &! work arrays for isopycnal mixing velocities
U_ISOP, V_ISOP ! horizontal components of isopycnal velocities
real (r8), dimension(nx_block,ny_block,nt) :: &
FX, FY ! fluxes across east, north faces
real (r8), dimension(2) :: &
reference_depth
!-----------------------------------------------------------------------
!
! initialize various quantities
!
!-----------------------------------------------------------------------
bid = this_block%local_id
U_ISOP = c0
V_ISOP = c0
WORK1 = c0
WORK2 = c0
WORK3 = c0
WORK4 = c0
if ( .not. implicit_vertical_mix ) &
call exit_POP (sigAbort, &
'implicit vertical mixing must be used with GM horiz mixing')
if ( k == 1 ) then
if ( diag_gm_bolus ) then
UIB = c0
VIB = c0
UIT(:,:,bid) = c0
VIT(:,:,bid) = c0
WBOT_ISOP(:,:,bid) = c0
endif
HOR_DIFF(:,:,ktp,k,bid) = ah_bkg_srfbl
BL_DEPTH(:,:,bid) = zw(k)
if ( vmix_itype == vmix_type_kpp ) &
BL_DEPTH(:,:,bid) = KPP_HBLT(:,:,bid)
endif
if (diag_gm_bolus) WTOP_ISOP(:,:,bid) = WBOT_ISOP(:,:,bid)
CX = merge(HYX(:,:,bid)*p25, c0, (k <= KMT (:,:,bid)) &
.and. (k <= KMTE(:,:,bid)))
CY = merge(HXY(:,:,bid)*p25, c0, (k <= KMT (:,:,bid)) &
.and. (k <= KMTN(:,:,bid)))
if ( k == 1 ) then
if ( transition_layer_on ) then
if ( vmix_itype == vmix_type_kpp ) then
call smooth_hblt ( .false., .true., bid, &
SMOOTH_OUT=TLT%DIABATIC_DEPTH(:,:,bid) )
else
TLT%DIABATIC_DEPTH(:,:,bid) = zw(k)
endif
do kk=1,km
do kk_sub = ktp,kbt
kid = kk + kk_sub - 2
SLA_SAVE(:,:,kk_sub,kk,bid) = dzw(kid)*sqrt(p5*( &
(SLX(:,:,1,kk_sub,kk,bid)**2 &
+ SLX(:,:,2,kk_sub,kk,bid)**2)/DXT(:,:,bid)**2 &
+ (SLY(:,:,1,kk_sub,kk,bid)**2 &
+ SLY(:,:,2,kk_sub,kk,bid)**2)/DYT(:,:,bid)**2)) &
+ eps
enddo
enddo
call transition_layer ( this_block )
endif
!-----------------------------------------------------------------------
!
! compute isopycnal and thickness diffusion coefficients if
! they depend on the model fields
!
!-----------------------------------------------------------------------
if ( ( kappa_isop_type == kappa_type_vmhs .or. &
kappa_thic_type == kappa_type_vmhs .or. &
kappa_isop_type == kappa_type_hdgr .or. &
kappa_thic_type == kappa_type_hdgr .or. &
kappa_isop_type == kappa_type_dradius .or. &
kappa_thic_type == kappa_type_dradius .or. &
kappa_isop_type == kappa_type_bfreq .or. &
kappa_thic_type == kappa_type_bfreq .or. &
kappa_isop_type == kappa_type_bfreq_vmhs .or. &
kappa_thic_type == kappa_type_bfreq_vmhs .or. &
kappa_isop_type == kappa_type_bfreq_hdgr .or. &
kappa_thic_type == kappa_type_bfreq_hdgr .or. &
kappa_isop_type == kappa_type_bfreq_dradius .or. &
kappa_thic_type == kappa_type_bfreq_dradius .or. &
kappa_isop_type == kappa_type_eg .or. &
kappa_thic_type == kappa_type_eg ) .and. &
( ( kappa_freq == kappa_freq_every_time_step ) &
.or. ( kappa_freq == kappa_freq_once_a_day .and. eod_last ) &
.or. ( nsteps_total == 1 ) ) ) compute_kappa(bid) = .true.
if ( compute_kappa(bid) ) then
if ( kappa_isop_type == kappa_type_vmhs .or. &
kappa_thic_type == kappa_type_vmhs .or. &
kappa_isop_type == kappa_type_bfreq_vmhs .or. &
kappa_thic_type == kappa_type_bfreq_vmhs ) then
if ( nsteps_total == 1 ) then
KAPPA_LATERAL(:,:,bid) = ah
if ( kappa_isop_type == kappa_type_const ) &
KAPPA_LATERAL(:,:,bid) = ah_bolus
else
call kappa_lon_lat_vmhs (TMIX, UMIX, VMIX, this_block)
endif
endif
if ( kappa_isop_type == kappa_type_hdgr .or. &
kappa_thic_type == kappa_type_hdgr .or. &
kappa_isop_type == kappa_type_bfreq_hdgr .or. &
kappa_thic_type == kappa_type_bfreq_hdgr ) &
call kappa_lon_lat_hdgr (TMIX, this_block)
if ( kappa_isop_type == kappa_type_dradius .or. &
kappa_thic_type == kappa_type_dradius .or. &
kappa_isop_type == kappa_type_bfreq_dradius .or. &
kappa_thic_type == kappa_type_bfreq_dradius ) &
call kappa_lon_lat_dradius (this_block)
if ( kappa_isop_type == kappa_type_bfreq .or. &
kappa_thic_type == kappa_type_bfreq .or. &
kappa_isop_type == kappa_type_bfreq_vmhs .or. &
kappa_thic_type == kappa_type_bfreq_vmhs .or. &
kappa_isop_type == kappa_type_bfreq_hdgr .or. &
kappa_thic_type == kappa_type_bfreq_hdgr .or. &
kappa_isop_type == kappa_type_bfreq_dradius .or. &
kappa_thic_type == kappa_type_bfreq_dradius ) &
call buoyancy_frequency_dependent_profile (TMIX, this_block)
if ( kappa_isop_type == kappa_type_eg .or. &
kappa_thic_type == kappa_type_eg ) &
call kappa_eg (TMIX, UMIX, VMIX, this_block)
compute_kappa(bid) = .false.
endif ! end of ( compute_kappa ) if statement
!-----------------------------------------------------------------------
!
! reinitialize the diffusivity coefficients
!
!-----------------------------------------------------------------------
if ( kappa_isop_type == kappa_type_const ) then
KAPPA_ISOP(:,:,:,:,bid) = ah
elseif ( kappa_isop_type == kappa_type_eg ) then
do kk_sub=ktp,kbt
do kk=1,km
KAPPA_ISOP(:,:,kk_sub,kk,bid) = KAPPA_VERTICAL(:,:,kk,bid)
enddo
enddo
else
do kk_sub=ktp,kbt
do kk=1,km
KAPPA_ISOP(:,:,kk_sub,kk,bid) = KAPPA_LATERAL(:,:,bid) &
* KAPPA_VERTICAL(:,:,kk,bid)
enddo
enddo
endif
if ( .not. use_const_ah_bkg_srfbl ) &
HOR_DIFF(:,:,ktp,k,bid) = KAPPA_ISOP(:,:,ktp,k,bid)
if ( kappa_thic_type == kappa_type_const ) then
KAPPA_THIC(:,:,:,:,bid) = ah_bolus
else if ( kappa_thic_type == kappa_type_depth .or. &
kappa_thic_type == kappa_type_bfreq ) then
do kk_sub=ktp,kbt
do kk=1,km
KAPPA_THIC(:,:,kk_sub,kk,bid) = ah_bolus &
* KAPPA_VERTICAL(:,:,kk,bid)
enddo
enddo
else if ( kappa_thic_type == kappa_type_eg ) then
KAPPA_THIC(:,:,:,:,bid) = KAPPA_ISOP(:,:,:,:,bid)
else
do kk_sub=ktp,kbt
do kk=1,km
KAPPA_THIC(:,:,kk_sub,kk,bid) = KAPPA_LATERAL(:,:,bid) &
* KAPPA_VERTICAL(:,:,kk,bid)
enddo
enddo
endif
!-----------------------------------------------------------------------
!
! control slope of isopycnal surfaces or KAPPA
!
!-----------------------------------------------------------------------
do kk=1,km
kp1 = min(kk+1,km)
reference_depth(ktp) = zt(kp1)
reference_depth(kbt) = zw(kp1)
if ( kk == km ) reference_depth(ktp) = zw(kp1)
do kk_sub = ktp,kbt
kid = kk + kk_sub - 2
!-----------------------------------------------------------------------
!
! control KAPPA to reduce the isopycnal mixing near the
! ocean surface Large et al (1997), JPO, 27, pp 2418-2447.
! WORK1 = ratio between the depth of water parcel and
! the vertical displacement of isopycnal surfaces
! where the vertical displacement =
! Rossby radius * slope of isopycnal surfaces
!
!-----------------------------------------------------------------------
if ( transition_layer_on ) then
SLA = SLA_SAVE(:,:,kk_sub,kk,bid)
else
SLA = dzw(kid)*sqrt(p5*( &
(SLX(:,:,1,kk_sub,kk,bid)**2 &
+ SLX(:,:,2,kk_sub,kk,bid)**2)/DXT(:,:,bid)**2 &
+ (SLY(:,:,1,kk_sub,kk,bid)**2 &
+ SLY(:,:,2,kk_sub,kk,bid)**2)/DYT(:,:,bid)**2)) &
+ eps
endif
TAPER1 = c1
if ( .not. transition_layer_on ) then
if ( kk == 1 ) then
dz_bottom = c0
else
dz_bottom = zt(kk-1)
endif
if (slope_control == slope_control_tanh) then
WORK1 = min(c1,zt(kk)*RBR(:,:,bid)/SLA)
TAPER1 = p5*(c1+sin(pi*(WORK1-p5)))
! use the Rossby deformation radius tapering
! only within the boundary layer
TAPER1 = merge(TAPER1, c1, &
dz_bottom <= BL_DEPTH(:,:,bid))
else
! sine function is replaced by
! function = 4.*x*(1.-abs(x)) for |x|<0.5
WORK1 = min(c1,zt(kk)*RBR(:,:,bid)/SLA)
TAPER1 = (p5+c2*(WORK1-p5)*(c1-abs(WORK1-p5)))
TAPER1 = merge(TAPER1, c1, &
dz_bottom <= BL_DEPTH(:,:,bid))
endif
endif
!-----------------------------------------------------------------------
!
! control KAPPA for numerical stability
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! methods to control slope
!
!-----------------------------------------------------------------------
TAPER2 = c1
TAPER3 = c1
select case (slope_control)
case (slope_control_tanh)
! method by Danabasoglu & Mcwilliams (1995)
TAPER2 = merge(p5* &
(c1-tanh(c10*SLA/slm_r-c4)), c0, SLA < slm_r)
if ( diff_tapering ) then
TAPER3 = merge(p5* &
(c1-tanh(c10*SLA/slm_b-c4)), c0, SLA < slm_b)
else
TAPER3 = TAPER2
endif
case (slope_control_notanh)
! similar to DM95 except replacing tanh by
! function = x*(1.-0.25*abs(x)) for |x|<2
! = sign(x) for |x|>2
! (faster than DM95)
do j=1,ny_block
do i=1,nx_block
if (SLA(i,j) > 0.2_r8*slm_r .and. &
SLA(i,j) < 0.6_r8*slm_r) then
TAPER2(i,j) = &
p5*(c1-(2.5_r8*SLA(i,j)/slm_r-c1)* &
(c4-abs(c10*SLA(i,j)/slm_r-c4)))
else if (SLA(i,j) >= 0.6_r8*slm_r) then
TAPER2(i,j) = c0
endif
enddo
enddo
if ( diff_tapering ) then
do j=1,ny_block
do i=1,nx_block
if (SLA(i,j) > 0.2_r8*slm_b .and. &
SLA(i,j) < 0.6_r8*slm_b) then
TAPER3(i,j) = &
p5*(c1-(2.5_r8*SLA(i,j)/slm_b-c1)* &
(c4-abs(c10*SLA(i,j)/slm_b-c4)))
else if (SLA(i,j) >= 0.6_r8*slm_b) then
TAPER3(i,j) = c0
endif
end do
end do
else
TAPER3 = TAPER2
endif
case (slope_control_clip)
! slope clipping
do n=1,2
do j=1,ny_block
do i=1,nx_block
if (abs(SLX(i,j,n,kk_sub,kk,bid) &
* dzw(kid) / HUS(i,j,bid)) > slm_r) then
SLX(i,j,n,kk_sub,kk,bid) = &
sign(slm_r * HUS(i,j,bid) &
* dzwr(kid), &
SLX(i,j,n,kk_sub,kk,bid))
endif
enddo
enddo
enddo
do n=1,2
do j=1,ny_block
do i=1,nx_block
if (abs(SLY(i,j,n,kk_sub,kk,bid) &
* dzw(kid) / HUW(i,j,bid)) > slm_r) then
SLY(i,j,n,kk_sub,kk,bid) = &
sign(slm_r * HUW(i,j,bid) &
* dzwr(kid), &
SLY(i,j,n,kk_sub,kk,bid))
endif
enddo
enddo
enddo
case (slope_control_Gerd)
! method by Gerdes et al (1991)
do j=1,ny_block
do i=1,nx_block
if (SLA(i,j) > slm_r) &
TAPER2(i,j) = (slm_r/SLA(i,j))**2
enddo
enddo
if (diff_tapering) then
do j=1,ny_block
do i=1,nx_block
if (SLA(i,j) > slm_b) &
TAPER3(i,j) = (slm_b/SLA(i,j))**2
enddo
enddo
else
TAPER3 = TAPER2
endif
end select
if ( transition_layer_on ) then
TAPER2 = merge(c1, TAPER2, reference_depth(kk_sub) &
<= TLT%DIABATIC_DEPTH(:,:,bid))
TAPER3 = merge(c1, TAPER3, reference_depth(kk_sub) &
<= TLT%DIABATIC_DEPTH(:,:,bid))
endif
if ( transition_layer_on .and. use_const_ah_bkg_srfbl ) then
HOR_DIFF(:,:,kk_sub,kk,bid) = ah_bkg_srfbl
else if ( transition_layer_on .and. &
( .not. use_const_ah_bkg_srfbl .or. &
kappa_isop_type == kappa_type_eg .or. &
kappa_thic_type == kappa_type_eg ) ) then
HOR_DIFF(:,:,kk_sub,kk,bid) = KAPPA_ISOP(:,:,kk_sub,kk,bid)
else
if ( .not. ( kk == 1 .and. kk_sub == ktp ) ) then
if ( use_const_ah_bkg_srfbl ) then
HOR_DIFF(:,:,kk_sub,kk,bid) = &
merge( ah_bkg_srfbl * (c1 - TAPER1 * TAPER2) &
* KAPPA_VERTICAL(:,:,kk,bid), &
c0, dz_bottom <= BL_DEPTH(:,:,bid) )
else
HOR_DIFF(:,:,kk_sub,kk,bid) = &
merge( KAPPA_ISOP(:,:,kk_sub,kk,bid) &
* (c1 - TAPER1 * TAPER2), &
c0, dz_bottom <= BL_DEPTH(:,:,bid) )
endif
endif
endif
KAPPA_ISOP(:,:,kk_sub,kk,bid) = &
TAPER1 * TAPER2 * KAPPA_ISOP(:,:,kk_sub,kk,bid)
KAPPA_THIC(:,:,kk_sub,kk,bid) = &
TAPER1 * TAPER3 * KAPPA_THIC(:,:,kk_sub,kk,bid)
end do ! end of kk_sub loop
!-----------------------------------------------------------------------
!
! impose the boundary conditions by setting KAPPA=0
! in the quarter cells adjacent to rigid boundaries.
! bottom B.C.s are considered within the kk-loop.
!
!-----------------------------------------------------------------------
! B.C. at the bottom
where (kk == KMT(:,:,bid))
KAPPA_ISOP(:,:,kbt,kk,bid) = c0
KAPPA_THIC(:,:,kbt,kk,bid) = c0
end where
enddo ! end of kk-loop
! B.C. at the top
KAPPA_ISOP(:,:,ktp,1,bid) = c0
KAPPA_THIC(:,:,ktp,1,bid) = c0
FZTOP(:,:,:,bid) = c0 ! zero flux B.C. at the surface
if ( transition_layer_on ) then
call merged_streamfunction ( this_block )
call apply_vertical_profile_to_isop_hor_diff ( this_block )
else
TLT%DIABATIC_DEPTH(:,:,bid) = c0
TLT%THICKNESS(:,:,bid) = c0
TLT%INTERIOR_DEPTH(:,:,bid) = c0
do kk=1,km
do kk_sub=ktp,kbt
where ( kk <= KMT(:,:,bid) )
SF_SLX(:,:,1,kk_sub,kk,bid) = &
KAPPA_THIC(:,:,kk_sub,kk,bid) &
* SLX(:,:,1,kk_sub,kk,bid) * dz(kk)
SF_SLX(:,:,2,kk_sub,kk,bid) = &
KAPPA_THIC(:,:,kk_sub,kk,bid) &
* SLX(:,:,2,kk_sub,kk,bid) * dz(kk)
SF_SLY(:,:,1,kk_sub,kk,bid) = &
KAPPA_THIC(:,:,kk_sub,kk,bid) &
* SLY(:,:,1,kk_sub,kk,bid) * dz(kk)
SF_SLY(:,:,2,kk_sub,kk,bid) = &
KAPPA_THIC(:,:,kk_sub,kk,bid) &
* SLY(:,:,2,kk_sub,kk,bid) * dz(kk)
endwhere
enddo ! end of kk_sub-loop
enddo ! end of kk-loop
endif
endif ! end of k==1 if statement
KMASK = merge(c1, c0, k < KMT(:,:,bid))
!-----------------------------------------------------------------------
!
! calculate effective vertical diffusion coefficient
! NOTE: it is assumed that VDC has been set before this
! in vmix_coeffs or something similar.
!
! Dz(VDC * Dz(T)) where D is derivative rather than difference
! VDC = (Az(dz*Ax(KAPPA*HYX*SLX**2)) + Az(dz*Ay(KAPPA*HXY*SLY**2)))*
! dzw/TAREA
!
!-----------------------------------------------------------------------
if ( k < km ) then
WORK1 = dzw(k)*KMASK*TAREA_R(:,:,bid)* &
(dz(k )*p25*KAPPA_ISOP(:,:,kbt,k, bid)* &
(HYX (:,:,bid)*SLX(:,:,ieast, kbt,k, bid)**2 &
+ HYXW(:,:,bid)*SLX(:,:,iwest, kbt,k, bid)**2 &
+ HXY (:,:,bid)*SLY(:,:,jnorth,kbt,k, bid)**2 &
+ HXYS(:,:,bid)*SLY(:,:,jsouth,kbt,k, bid)**2) &
+dz(k+1)*p25*KAPPA_ISOP(:,:,ktp,k+1,bid)* &
(HYX (:,:,bid)*SLX(:,:,ieast, ktp,k+1,bid)**2 &
+ HYXW(:,:,bid)*SLX(:,:,iwest, ktp,k+1,bid)**2 &
+ HXY (:,:,bid)*SLY(:,:,jnorth,ktp,k+1,bid)**2 &
+ HXYS(:,:,bid)*SLY(:,:,jsouth,ktp,k+1,bid)**2))
do n=1,size(VDC,DIM=4)
VDC_GM(:,:,k,bid) = WORK1
VDC(:,:,k,n,bid) = VDC(:,:,k,n,bid) + WORK1
end do
end if
!-----------------------------------------------------------------------
!
! check if some horizontal diffusion needs to be added to the
! bottom half of the bottom cell
!
!-----------------------------------------------------------------------
if ( ah_bkg_bottom /= c0 ) then
where ( k == KMT(:,:,bid) )
HOR_DIFF(:,:,kbt,k,bid) = ah_bkg_bottom
endwhere
endif
!-----------------------------------------------------------------------
!
! combine isopycnal and horizontal diffusion coefficients
!
!-----------------------------------------------------------------------
do j=1,ny_block
do i=1,nx_block-1
WORK3(i,j) = KAPPA_ISOP(i, j,ktp,k,bid) &
+ HOR_DIFF (i, j,ktp,k,bid) &
+ KAPPA_ISOP(i, j,kbt,k,bid) &
+ HOR_DIFF (i, j,kbt,k,bid) &
+ KAPPA_ISOP(i+1,j,ktp,k,bid) &
+ HOR_DIFF (i+1,j,ktp,k,bid) &
+ KAPPA_ISOP(i+1,j,kbt,k,bid) &
+ HOR_DIFF (i+1,j,kbt,k,bid)
enddo
enddo
do j=1,ny_block-1
do i=1,nx_block
WORK4(i,j) = KAPPA_ISOP(i,j, ktp,k,bid) &
+ HOR_DIFF (i,j, ktp,k,bid) &
+ KAPPA_ISOP(i,j, kbt,k,bid) &
+ HOR_DIFF (i,j, kbt,k,bid) &
+ KAPPA_ISOP(i,j+1,ktp,k,bid) &
+ HOR_DIFF (i,j+1,ktp,k,bid) &
+ KAPPA_ISOP(i,j+1,kbt,k,bid) &
+ HOR_DIFF (i,j+1,kbt,k,bid)
enddo
enddo
!-----------------------------------------------------------------------
!
! start loop over tracers
!
!-----------------------------------------------------------------------
kp1 = k + 1
if ( k == km ) kp1 = k
if ( k < km ) then
dz_bottom = dz(kp1)
factor = c1
else
dz_bottom = c0
factor = c0
endif
call timer_start(timer_nloop, block_id=this_block%local_id)
do n = 1,nt
!-----------------------------------------------------------------------
!
! calculate horizontal fluxes thru vertical faces of T-cell
! FX = dz*HYX*Ax(Az(KAPPA))*Dx(T) : flux in x-direction
! FY = dz*HXY*Ay(Az(KAPPA))*Dy(T) : flux in y-direction
!
!-----------------------------------------------------------------------
FX(:,:,n) = dz(k) * CX * TX(:,:,k,n,bid) * WORK3
FY(:,:,n) = dz(k) * CY * TY(:,:,k,n,bid) * WORK4
end do
if ( .not. cancellation_occurs ) then
do j=1,ny_block
do i=1,nx_block-1
WORK1(i,j) = KAPPA_ISOP(i,j,ktp,k,bid) &
* SLX(i,j,ieast,ktp,k,bid) * dz(k) &
- SF_SLX(i,j,ieast,ktp,k,bid)
WORK2(i,j) = KAPPA_ISOP(i,j,kbt,k,bid) &
* SLX(i,j,ieast,kbt,k,bid) * dz(k) &
- SF_SLX(i,j,ieast,kbt,k,bid)
WORK3(i,j) = KAPPA_ISOP(i+1,j,ktp,k,bid) &
* SLX(i+1,j,iwest,ktp,k,bid) * dz(k) &
- SF_SLX(i+1,j,iwest,ktp,k,bid)
WORK4(i,j) = KAPPA_ISOP(i+1,j,kbt,k,bid) &
* SLX(i+1,j,iwest,kbt,k,bid) * dz(k) &
- SF_SLX(i+1,j,iwest,kbt,k,bid)
enddo
enddo
do n = 1,nt
if (n > 2 .and. k < km) &
TZ(:,:,k+1,n,bid) = TMIX(:,:,k ,n) - TMIX(:,:,k+1,n)
do j=1,ny_block
do i=1,nx_block-1
FX(i,j,n) = FX(i,j,n) - CX(i,j) &
* ( WORK1(i,j) * TZ(i,j,k,n,bid) &
+ WORK2(i,j) * TZ(i,j,kp1,n,bid) &
+ WORK3(i,j) * TZ(i+1,j,k,n,bid) &
+ WORK4(i,j) * TZ(i+1,j,kp1,n,bid) )
enddo
enddo
end do
do j=1,ny_block-1
do i=1,nx_block
WORK1(i,j) = KAPPA_ISOP(i,j,ktp,k,bid) &
* SLY(i,j,jnorth,ktp,k,bid) * dz(k) &
- SF_SLY(i,j,jnorth,ktp,k,bid)
WORK2(i,j) = KAPPA_ISOP(i,j,kbt,k,bid) &
* SLY(i,j,jnorth,kbt,k,bid) * dz(k) &
- SF_SLY(i,j,jnorth,kbt,k,bid)
WORK3(i,j) = KAPPA_ISOP(i,j+1,ktp,k,bid) &
* SLY(i,j+1,jsouth,ktp,k,bid) * dz(k) &
- SF_SLY(i,j+1,jsouth,ktp,k,bid)
WORK4(i,j) = KAPPA_ISOP(i,j+1,kbt,k,bid) &
* SLY(i,j+1,jsouth,kbt,k,bid) * dz(k) &
- SF_SLY(i,j+1,jsouth,kbt,k,bid)
enddo
enddo
do n = 1,nt
do j=1,ny_block-1
do i=1,nx_block
FY(i,j,n) = FY(i,j,n) - CY(i,j) &
* ( WORK1(i,j) * TZ(i,j,k,n,bid) &
+ WORK2(i,j) * TZ(i,j,kp1,n,bid) &
+ WORK3(i,j) * TZ(i,j+1,k,n,bid) &
+ WORK4(i,j) * TZ(i,j+1,kp1,n,bid) )
enddo
enddo
end do
endif
do n = 1,nt
!-----------------------------------------------------------------------
!
! calculate vertical fluxes thru horizontal faces of T-cell
! - Az(dz*Ax(HYX*KAPPA*SLX*TX)) - Az(dz*Ay(HXY*KAPPA*SLY*TY))
! calculate isopycnal diffusion from flux differences
! DTK = (Dx(FX)+Dy(FY)+Dz(FZ)) / volume
!
!-----------------------------------------------------------------------
GTK(:,:,n) = c0
if ( k < km ) then
WORK3 = c0
if ( .not. cancellation_occurs ) then
!pw loop split to improve performance -- 2
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK3(i,j) = WORK3(i,j) &
+ ( dz(k) * KAPPA_ISOP(i,j,kbt,k,bid) &
* ( SLX(i,j,ieast ,kbt,k ,bid) &
* HYX(i ,j ,bid) * TX(i ,j ,k ,n,bid) &
+ SLY(i,j,jnorth,kbt,k ,bid) &
* HXY(i ,j ,bid) * TY(i ,j ,k ,n,bid) &
+ SLX(i,j,iwest ,kbt,k ,bid) &
* HYX(i-1,j ,bid) * TX(i-1,j ,k ,n,bid) &
+ SLY(i,j,jsouth,kbt,k ,bid) &
* HXY(i ,j-1,bid) * TY(i ,j-1,k ,n,bid) ) )
enddo
enddo
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK3(i,j) = WORK3(i,j) &
+ ( SF_SLX(i ,j ,ieast ,kbt,k ,bid) &
* HYX(i ,j ,bid) * TX(i ,j ,k ,n,bid) &
+ SF_SLY(i ,j ,jnorth,kbt,k ,bid) &
* HXY(i ,j ,bid) * TY(i ,j ,k ,n,bid) &
+ SF_SLX(i ,j ,iwest ,kbt,k ,bid) &
* HYX(i-1,j ,bid) * TX(i-1,j ,k ,n,bid) &
+ SF_SLY(i ,j ,jsouth,kbt,k ,bid) &
* HXY(i ,j-1,bid) * TY(i ,j-1,k ,n,bid) )
enddo
enddo
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK3(i,j) = WORK3(i,j) &
+ ( dz_bottom * KAPPA_ISOP(i,j,ktp,kp1,bid) &
* ( SLX(i ,j ,ieast ,ktp,kp1,bid) &
* HYX(i ,j ,bid) * TX(i ,j ,kp1,n,bid) &
+ SLY(i ,j ,jnorth,ktp,kp1,bid) &
* HXY(i ,j ,bid) * TY(i ,j ,kp1,n,bid) &
+ SLX(i ,j ,iwest ,ktp,kp1,bid) &
* HYX(i-1,j ,bid) * TX(i-1,j ,kp1,n,bid) &
+ SLY(i ,j ,jsouth,ktp,kp1,bid) &
* HXY(i ,j-1,bid) * TY(i ,j-1,kp1,n,bid) ) )
enddo
enddo
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK3(i,j) = WORK3(i,j) &
+ ( factor &
* ( SF_SLX(i ,j ,ieast ,ktp,kp1,bid) &
* HYX(i ,j ,bid) * TX(i ,j ,kp1,n,bid) &
+ SF_SLY(i ,j ,jnorth,ktp,kp1,bid) &
* HXY(i ,j ,bid) * TY(i ,j ,kp1,n,bid) &
+ SF_SLX(i ,j ,iwest ,ktp,kp1,bid) &
* HYX(i-1,j ,bid) * TX(i-1,j ,kp1,n,bid) &
+ SF_SLY(i ,j ,jsouth,ktp,kp1,bid) &
* HXY(i ,j-1,bid) * TY(i ,j-1,kp1,n,bid) ) )
enddo
enddo
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
fz = -KMASK(i,j) * p25 * WORK3(i,j)
GTK(i,j,n) = ( FX(i,j,n) - FX(i-1,j,n) &
+ FY(i,j,n) - FY(i,j-1,n) &
+ FZTOP(i,j,n,bid) - fz )*dzr(k)*TAREA_R(i,j,bid)
FZTOP(i,j,n,bid) = fz
enddo
enddo
else
!pw loop split to improve performance
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK3(i,j) = &
( dz(k) * KAPPA_ISOP(i,j,kbt,k,bid) &
* ( SLX(i,j,ieast ,kbt,k ,bid) &
* HYX(i ,j ,bid) * TX(i ,j ,k ,n,bid) &
+ SLY(i,j,jnorth,kbt,k ,bid) &
* HXY(i ,j ,bid) * TY(i ,j ,k ,n,bid) &
+ SLX(i,j,iwest ,kbt,k ,bid) &
* HYX(i-1,j ,bid) * TX(i-1,j ,k ,n,bid) &
+ SLY(i,j,jsouth,kbt,k ,bid) &
* HXY(i ,j-1,bid) * TY(i ,j-1,k ,n,bid) ) )
enddo
enddo
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK3(i,j) = WORK3(i,j) &
+ ( dz_bottom * KAPPA_ISOP(i,j,ktp,kp1,bid) &
* ( SLX(i ,j ,ieast ,ktp,kp1,bid) &
* HYX(i ,j ,bid) * TX(i ,j ,kp1,n,bid) &
+ SLY(i ,j ,jnorth,ktp,kp1,bid) &
* HXY(i ,j ,bid) * TY(i ,j ,kp1,n,bid) &
+ SLX(i ,j ,iwest ,ktp,kp1,bid) &
* HYX(i-1,j ,bid) * TX(i-1,j ,kp1,n,bid) &
+ SLY(i ,j ,jsouth,ktp,kp1,bid) &
* HXY(i ,j-1,bid) * TY(i ,j-1,kp1,n,bid) ) )
enddo
enddo
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
fz = -KMASK(i,j) * p5 * WORK3(i,j)
GTK(i,j,n) = ( FX(i,j,n) - FX(i-1,j,n) &
+ FY(i,j,n) - FY(i,j-1,n) &
+ FZTOP(i,j,n,bid) - fz )*dzr(k)*TAREA_R(i,j,bid)
FZTOP(i,j,n,bid) = fz
enddo
enddo
endif
else ! k = km
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
GTK(i,j,n) = ( FX(i,j,n) - FX(i-1,j,n) &
+ FY(i,j,n) - FY(i,j-1,n) &
+ FZTOP(i,j,n,bid) )*dzr(k)*TAREA_R(i,j,bid)
FZTOP(i,j,n,bid) = c0
enddo
enddo
endif
!-----------------------------------------------------------------------
!
! end of tracer loop
!
!-----------------------------------------------------------------------
enddo
call timer_stop(timer_nloop, block_id=this_block%local_id)
!-----------------------------------------------------------------------
!
! diagnostic computation of the bolus velocities
!
!-----------------------------------------------------------------------
if ( diag_gm_bolus ) then
do j=1,ny_block-1
do i=1,nx_block-1
WORK1(i,j) = ( SF_SLX(i ,j,ieast,kbt,k, bid) &
+ factor * SF_SLX(i ,j,ieast,ktp,kp1,bid) &
+ SF_SLX(i+1,j,iwest,kbt,k, bid) &
+ factor * SF_SLX(i+1,j,iwest,ktp,kp1,bid)) &
* p25 * HYX(i,j,bid)
WORK2(i,j) = ( SF_SLY(i,j ,jnorth,kbt,k, bid) &
+ factor * SF_SLY(i,j ,jnorth,ktp,kp1,bid) &
+ SF_SLY(i,j+1,jsouth,kbt,k, bid) &
+ factor * SF_SLY(i,j+1,jsouth,ktp,kp1,bid)) &
* p25 * HXY(i,j,bid)
enddo
enddo
UIB = merge( WORK1, c0, k < KMT(:,:,bid) &
.and. k < KMTE(:,:,bid) )
VIB = merge( WORK2, c0, k < KMT(:,:,bid) &
.and. k < KMTN(:,:,bid) )
WORK1 = merge( UIT(:,:,bid) - UIB, c0, k <= KMT(:,:,bid) &
.and. k <= KMTE(:,:,bid) )
WORK2 = merge( VIT(:,:,bid) - VIB, c0, k <= KMT(:,:,bid) &
.and. k <= KMTN(:,:,bid) )
U_ISOP = WORK1 * dzr(k) / HTE(:,:,bid)
V_ISOP = WORK2 * dzr(k) / HTN(:,:,bid)
if ( linertial .and. k == 2 ) then
BOLUS_SP(:,:,bid) = 50.0_r8 * sqrt(U_ISOP**c2 + V_ISOP**c2)
endif
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
if ( k < KMT(i,j,bid) ) then
WBOT_ISOP(i,j,bid) = WTOP_ISOP(i,j,bid) &
+ TAREA_R(i,j,bid) &
* ( WORK1(i,j) - WORK1(i-1,j) &
+ WORK2(i,j) - WORK2(i,j-1) )
else
WBOT_ISOP(i,j,bid) = c0
endif
enddo
enddo
endif
!-----------------------------------------------------------------------
!
! update remaining bottom-face fields to top-face fields for next
! pass
!
!-----------------------------------------------------------------------
if ( diag_gm_bolus ) then
UIT(:,:,bid) = UIB
VIT(:,:,bid) = VIB
endif
!-----------------------------------------------------------------------
!
! compute isopycnal diffusion cfl diagnostics if required
!
!-----------------------------------------------------------------------
if (ldiag_cfl) then
WORK2 = p5 * (KAPPA_ISOP(:,:,ktp,k,bid) &
+ KAPPA_ISOP(:,:,kbt,k,bid))
WORK1 = merge(c4*WORK2*(DXTR(:,:,bid)**2 + DYTR(:,:,bid)**2), &
c0, KMT(:,:,bid) > k)
WORK2 = abs(WORK1)
call cfl_hdiff (k,bid,WORK2,1,this_block)
endif
!-----------------------------------------------------------------------
!
! accumulate time average if necessary
!
!-----------------------------------------------------------------------
if ( mix_pass /= 1 ) then
tavg_stream = tavg_in_which_stream(tavg_KAPPA_ISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_KAPPA_ISOP)) then
call accumulate_tavg_field &
(p5 * (KAPPA_ISOP(:,:,ktp,k,bid) &
+ KAPPA_ISOP(:,:,kbt,k,bid)), &
tavg_KAPPA_ISOP, bid, k)
endif
tavg_stream = tavg_in_which_stream(tavg_KAPPA_THIC)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_KAPPA_THIC)) then
call accumulate_tavg_field &
(p5 * (KAPPA_THIC(:,:,ktp,k,bid) &
+ KAPPA_THIC(:,:,kbt,k,bid)), &
tavg_KAPPA_THIC, bid, k)
endif
tavg_stream = tavg_in_which_stream(tavg_HOR_DIFF)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_HOR_DIFF)) then
call accumulate_tavg_field &
(p5 * (HOR_DIFF(:,:,ktp,k,bid) &
+ HOR_DIFF(:,:,kbt,k,bid)), &
tavg_HOR_DIFF, bid, k)
endif
if ( transition_layer_on .and. k == 1 ) then
tavg_stream = tavg_in_which_stream(tavg_DIA_DEPTH)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_DIA_DEPTH)) then
call accumulate_tavg_field (TLT%DIABATIC_DEPTH(:,:,bid), &
tavg_DIA_DEPTH, bid, 1)
endif
tavg_stream = tavg_in_which_stream(tavg_TLT)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_TLT)) then
call accumulate_tavg_field (TLT%THICKNESS(:,:,bid), &
tavg_TLT, bid, 1)
endif
tavg_stream = tavg_in_which_stream(tavg_INT_DEPTH)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_INT_DEPTH)) then
call accumulate_tavg_field (TLT%INTERIOR_DEPTH(:,:,bid), &
tavg_INT_DEPTH, bid, 1)
endif
endif
if ( diag_gm_bolus ) then
tavg_stream = tavg_in_which_stream(tavg_UISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_UISOP)) then
call accumulate_tavg_field (U_ISOP, tavg_UISOP, bid, k)
endif
tavg_stream = tavg_in_which_stream(tavg_VISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_VISOP)) then
call accumulate_tavg_field (V_ISOP, tavg_VISOP, bid, k)
endif
tavg_stream = tavg_in_which_stream(tavg_WISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_WISOP)) then
call accumulate_tavg_field (WTOP_ISOP(:,:,bid), tavg_WISOP,&
bid, k)
endif
tavg_stream = tavg_in_which_stream(tavg_ADVT_ISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_ADVT_ISOP)) then
WORK1 = p5 * HTE(:,:,bid) * U_ISOP * ( TMIX(:,:,k,1) &
+ eoshift(TMIX(:,:,k,1), dim=1, shift=1) )
WORK2 = eoshift(WORK1, dim=1, shift=-1)
WORK3 = WORK1 - WORK2
WORK1 = p5 * HTN(:,:,bid) * V_ISOP * ( TMIX(:,:,k,1) &
+ eoshift(TMIX(:,:,k,1), dim=2, shift=1) )
WORK2 = eoshift(WORK1, dim=2, shift=-1)
WORK3 = WORK3 + WORK1 - WORK2
WORK1 = c0
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
if ( k <= KMT(i,j,bid) ) then
WORK1(i,j) = - dz(k) * TAREA_R(i,j,bid) * WORK3(i,j)
endif
enddo
enddo
call accumulate_tavg_field (WORK1, tavg_ADVT_ISOP, bid, k)
endif
tavg_stream = tavg_in_which_stream(tavg_ADVS_ISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_ADVS_ISOP)) then
WORK1 = p5 * HTE(:,:,bid) * U_ISOP * ( TMIX(:,:,k,2) &
+ eoshift(TMIX(:,:,k,2), dim=1, shift=1) )
WORK2 = eoshift(WORK1, dim=1, shift=-1)
WORK3 = WORK1 - WORK2
WORK1 = p5 * HTN(:,:,bid) * V_ISOP * ( TMIX(:,:,k,2) &
+ eoshift(TMIX(:,:,k,2), dim=2, shift=1) )
WORK2 = eoshift(WORK1, dim=2, shift=-1)
WORK3 = WORK3 + WORK1 - WORK2
WORK1 = c0
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
if ( k <= KMT(i,j,bid) ) then
WORK1(i,j) = - dz(k) * TAREA_R(i,j,bid) * WORK3(i,j)
endif
enddo
enddo
call accumulate_tavg_field (WORK1, tavg_ADVS_ISOP, bid, k)
endif
if ( tavg_requested(tavg_VNT_ISOP) .or. &
tavg_requested(tavg_VNS_ISOP) ) then
WORK1 = p5 * V_ISOP * HTN(:,:,bid) * TAREA_R(:,:,bid)
tavg_stream = tavg_in_which_stream(tavg_VNT_ISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_VNT_ISOP)) then
WORK2 = WORK1 * ( TMIX(:,:,k,1) &
+ eoshift(TMIX(:,:,k,1), dim=2, shift=1) )
call accumulate_tavg_field (WORK2, tavg_VNT_ISOP, bid, k)
endif
tavg_stream = tavg_in_which_stream(tavg_VNS_ISOP)
if (ltavg_on(tavg_stream) .and. tavg_requested(tavg_VNS_ISOP)) then
WORK2 = WORK1 * ( TMIX(:,:,k,2) &
+ eoshift(TMIX(:,:,k,2), dim=2, shift=1) )
call accumulate_tavg_field (WORK2, tavg_VNS_ISOP, bid, k)
endif
endif
endif ! bolus velocity option on
do n = 1,nt
if (tavg_requested(tavg_HDIFE_TRACER(n))) then
if (ltavg_on(tavg_in_which_stream(tavg_HDIFE_TRACER(n)))) then
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK1(i,j) = FX(i,j,n)*dzr(k)*TAREA_R(i,j,bid)
enddo
enddo
call accumulate_tavg_field(WORK1,tavg_HDIFE_TRACER(n),bid,k)
endif
endif
if (tavg_requested(tavg_HDIFN_TRACER(n))) then
if (ltavg_on(tavg_in_which_stream(tavg_HDIFN_TRACER(n)))) then
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK1(i,j) = FY(i,j,n)*dzr(k)*TAREA_R(i,j,bid)
enddo
enddo
call accumulate_tavg_field(WORK1,tavg_HDIFN_TRACER(n),bid,k)
endif
endif
if (tavg_requested(tavg_HDIFB_TRACER(n))) then
if (ltavg_on(tavg_in_which_stream(tavg_HDIFB_TRACER(n)))) then
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK1(i,j) = FZTOP(i,j,n,bid)*dzr(k)*TAREA_R(i,j,bid)
enddo
enddo
call accumulate_tavg_field(WORK1,tavg_HDIFB_TRACER(n),bid,k)
endif
endif
enddo
endif ! mix_pass ne 1
!-----------------------------------------------------------------------
!EOC
end subroutine hdifft_gm
!***********************************************************************
!BOP
! !IROUTINE: kappa_lon_lat_vmhs
! !INTERFACE:
subroutine kappa_lon_lat_vmhs (TMIX, UMIX, VMIX, this_block) 1,5
! !DESCRIPTION:
! Variable kappa parameterization by Visbeck et al (1997):
! \begin{equation}
! KAPPA_LATERAL = C {{l^2}\over{T}},
! \end{equation}
! where $C$ is a constant, $T$ is the (Eady) time scale
! $f/\surd\overline{Ri}$ and $l$ is the Rossby radius.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TMIX ! tracers at all vertical levels
! at mixing time level
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UMIX, VMIX ! U,V at all vertical levels
! at mixing time level
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
kp1,i,j,kk,k, &
k1,k2, &
bid ! local block address for this sub block
real (r8) :: &
const, zmin1, zmin2
real (r8), dimension(nx_block,ny_block) :: &
UKT, VKT, &
UKPT, VKPT, RNUM, &
WORK1, WORK2, WORK3, WORK4, &
TK, TKP, &
RHOT, RHOS, &
GRATE, LSC
!-----------------------------------------------------------------------
! initialization
!-----------------------------------------------------------------------
bid = this_block%local_id
GRATE = c0
LSC = c0
k1 = 0
k2 = 0
vert_loop: do k=1,km
! -2000m < z < -100m
if ( zt(k) >= 1.0e4_r8 .and. zt(k) <= 2.0e5_r8 ) then
!-----------------------------------------------------------------------
!
! start computing the Richardson number and the baroclinic wave speed
! average velocities at T points of level k
!
!-----------------------------------------------------------------------
if ( k1 == 0 ) then
k1 = k ! k1 is the upper limit for integration
call ugrid_to_tgrid(UKT,UMIX(:,:,k),bid)
call ugrid_to_tgrid(VKT,VMIX(:,:,k),bid)
TK = max(-c2, TMIX(:,:,k,1))
endif
!-----------------------------------------------------------------------
!
! compute RZ=Dz(rho) at k=k
!
!-----------------------------------------------------------------------
if ( k < km ) then
kp1 = k + 1
else
kp1 = km
endif
call state (k,kp1,TMIX(:,:,k,1),TMIX(:,:,k,2), &
this_block, DRHODT=RHOT, DRHODS=RHOS)
TKP = max(-c2, TMIX(:,:,kp1,1))
WORK1 = TK - TKP
WORK2 = TMIX(:,:,k,2) - TMIX(:,:,kp1,2)
WORK3 = RHOT * WORK1 + RHOS * WORK2
WORK3 = min(WORK3,-eps2)
!-----------------------------------------------------------------------
!
! average velocities at T points of level k+1
!
!-----------------------------------------------------------------------
call ugrid_to_tgrid(UKPT,UMIX(:,:,kp1),bid)
call ugrid_to_tgrid(VKPT,VMIX(:,:,kp1),bid)
!-----------------------------------------------------------------------
!
! local Richardson number at the bottom of T box, level k
! = -g*Dz(RHO)/RHO_0/(Dz(u)**2+Dz(v)**2)
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
where (k < KMT(:,:,bid))
! RHO_0 = 1 in denominator (approximately) in cgs unit.
WORK4 = p5*(dz(k) + DZT(:,:,k+1,bid)) ! dzw equivalent
RNUM = -WORK4/((UKT - UKPT)**2 + (VKT - VKPT)**2 + eps)
GRATE = GRATE + grav*RNUM*WORK4*WORK3
LSC = LSC - grav*WORK3
end where
else ! no partial bottom cells
where (k < KMT(:,:,bid))
! RHO_0 = 1 in denominator (approximately) in cgs unit.
RNUM = -dzw(k)/((UKT - UKPT)**2 + (VKT - VKPT)**2 + eps)
GRATE = GRATE + grav*RNUM*dzw(k)*WORK3
LSC = LSC - grav*WORK3
end where
endif ! partial bottom cells
!-----------------------------------------------------------------------
!
! bottom values become top values for next pass
!
!-----------------------------------------------------------------------
UKT = UKPT
VKT = VKPT
TK = TKP
else
if ( k1 /= 0 .and. k2 == 0 ) then
k2 = k ! k2 is the lower limit for integration
exit vert_loop
endif
endif ! if(zt(k).ge.1.0e4.and.zt(k).lt.2.0e5)
end do vert_loop
do j=1,ny_block
do i=1,nx_block
kk = KMT(i,j,bid)
if ( kk > 0 ) then
zmin1 = min(zt(k1),zt(kk))
zmin2 = min(zt(k2),zt(kk))
GRATE(i,j) = GRATE(i,j)/(zmin2-zmin1+eps) ! Ri
LSC(i,j) = LSC(i,j)*(zmin2-zmin1) ! c_g^2=N^2*H^2
else
GRATE(i,j) = c0
LSC(i,j) = c0
endif
end do
end do
! equatorial inverse of time scale and square of length scale
WORK1 = sqrt(c2*sqrt(LSC)*BTP(:,:,bid)) ! sqrt(c_g*2*beta)
WORK2 = sqrt(LSC)/(c2*BTP(:,:,bid)) ! c_g/(2 beta)
! inverse of time scale
WORK1 = max(abs(FCORT(:,:,bid)),WORK1)
GRATE = WORK1/sqrt(GRATE+eps) ! 1/T = f/sqrt(Ri)
! square of length scale
LSC = LSC/(FCORT(:,:,bid)+eps)**2 ! L^2 = c_g^2/f^2
LSC = min(LSC,WORK2)
! grid size = lower limit of length scale
WORK1 = min(DXT(:,:,bid)**2,DYT(:,:,bid)**2)
LSC = max(LSC,WORK1)
!-----------------------------------------------------------------------
!
! compute KAPPA_LATERAL
!
!-----------------------------------------------------------------------
! const = 0.015_r8 ! constant taken from Visbeck et al (1997)
const = 0.13_r8
WORK1 = const*GRATE*LSC ! const*L**2/T
! KAPPA_LATERAL is bounded by 3.0e6 < KAPPA_LATERAL < 4.0e7
KAPPA_LATERAL(:,:,bid) = min(4.0e7_r8,WORK1)
KAPPA_LATERAL(:,:,bid) = max(3.0e6_r8,KAPPA_LATERAL(:,:,bid))
where (KMT(:,:,bid) <= k1) ! KAPPA_LATERAL=3.0e6 when depth < 100m
KAPPA_LATERAL(:,:,bid) = 3.0e6_r8
end where
!-----------------------------------------------------------------------
!EOC
end subroutine kappa_lon_lat_vmhs
!***********************************************************************
!BOP
! !IROUTINE: kappa_eg
! !INTERFACE:
subroutine kappa_eg (TMIX, UMIX, VMIX, this_block) 1,2
! !DESCRIPTION:
! Variable kappa parameterization by Eden and Greatbatch
! (2008, Ocean Modelling, v. 20, 223-239):
! \begin{equation}
! KAPPA_ISOP = KAPPA_THIC = c {{L^2}*{\sigma}},
! \end{equation}
! where $c$ is a tuning parameter (=const_eg), $\sigma$ is the inverse eddy
! time scale, and $L$ is the minimum of the Rossby deformation radius and
! Rhines scale.
!
! This subroutine returns KAPPA_ISOP in KAPPA_VERTICAL. Here, KAPPA_VERTICAL
! serves as a temporary array because this subroutine may be called less
! frequently than every time step and KAPPA_ISOP is modified in each time step.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TMIX ! tracers at all vertical levels
! at mixing time level
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UMIX, VMIX ! U,V at all vertical levels
! at mixing time level
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, kp1, & ! vertical level index
bid ! local block address for this sub block
real (r8), dimension(nx_block,ny_block) :: &
TK, TKP, & ! level k and k+1 TMIX
WORK1, WORK2, WORK3, & ! work arrays
C_ROSSBY, & ! first baroclinic Rossby wave speed
L_ROSSBY, & ! Rossby deformation radius
RHOT, RHOS ! dRHOdT and dRHOdS
real (r8), dimension(nx_block,ny_block,km) :: &
SIGMA, & ! inverse eddy time scale
RI_NO ! local Richardson number
!-----------------------------------------------------------------------
! initialization
!-----------------------------------------------------------------------
bid = this_block%local_id
SIGMA = c0
RI_NO = c0
BUOY_FREQ_SQ(:,:,:,bid) = c0
!-----------------------------------------------------------------------
!
! compute buoyancy frequency and Richardson number at the bottom of
! T box:
! Ri = -g*Dz(RHO)/RHO_0/(Dz(u)**2+Dz(v)**2)
!
! RHO_0 ~ 1 in cgs units.
!
!-----------------------------------------------------------------------
do k=1,km-1
if ( k == 1 ) then
TK = max(-c2, TMIX(:,:,k,1))
endif
kp1 = k+1
call state (k, kp1, TMIX(:,:,k,1), TMIX(:,:,k,2), &
this_block, DRHODT=RHOT, DRHODS=RHOS)
TKP = max(-c2, TMIX(:,:,kp1,1))
WORK1 = TK - TKP
WORK2 = TMIX(:,:,k,2) - TMIX(:,:,kp1,2)
WORK3 = RHOT * WORK1 + RHOS * WORK2
WORK3 = min(WORK3,-eps2)
WORK1 = (UMIX(:,:,k) - UMIX(:,:,kp1))**2 &
+ (VMIX(:,:,k) - VMIX(:,:,kp1))**2
call ugrid_to_tgrid (WORK2, WORK1, bid)
where (k < KMT(:,:,bid))
BUOY_FREQ_SQ(:,:,k,bid) = - grav * WORK3 * dzwr(k)
WORK1 = dzw(k)**2 / ( WORK2 + eps2 )
RI_NO(:,:,k) = WORK1 * BUOY_FREQ_SQ(:,:,k,bid)
end where
TK = TKP
end do
!-----------------------------------------------------------------------
!
! compute the first baroclinic gravity-wave phase speed.
! Computation of Rossby deformation radius follows Chelton et al.(1998)
!
!-----------------------------------------------------------------------
C_ROSSBY = c0
k = 1
where ( k < KMT(:,:,bid) )
C_ROSSBY = C_ROSSBY + sqrt(BUOY_FREQ_SQ(:,:,k,bid)) * dzw(k-1)
endwhere
do k=1,km
where ( k < KMT(:,:,bid) )
C_ROSSBY = C_ROSSBY + sqrt(BUOY_FREQ_SQ(:,:,k,bid)) * dzw(k)
endwhere
where ( k > 1 .and. k == KMT(:,:,bid) )
C_ROSSBY = C_ROSSBY + sqrt(BUOY_FREQ_SQ(:,:,k-1,bid)) * dzw(k)
endwhere
enddo
C_ROSSBY = C_ROSSBY / pi
L_ROSSBY = min( C_ROSSBY / (abs(FCORT(:,:,bid))+eps), &
sqrt( C_ROSSBY / (c2*BTP(:,:,bid)) ) )
!-----------------------------------------------------------------------
!
! compute the inverse time scale
!
!-----------------------------------------------------------------------
do k=1,km-1
WORK1 = max( abs( FCORT(:,:,bid) ), sqrt(C_ROSSBY * c2 * BTP(:,:,bid)) )
where (k < KMT(:,:,bid))
SIGMA(:,:,k) = SIGMA_TOPO_MASK(:,:,k,bid) * WORK1 &
/ sqrt( RI_NO(:,:,k) + gamma_eg )
end where
enddo
!----------------------------------------------------------------------
!
! compute KAPPA_ISOP = c_kappa * sigma * min( Rossby_radius, Rhines scales)^2.
! note that KAPPA_ISOP is stored in KAPPA_VERTICAL. KAPPA_VERTICAL is
! located at the T-grid points. in the following, KAPPA_ISOP computed at
! the lower interface of a T-grid box is assigned to the T-grid point
! just above this interface for simplicity.
!
!-----------------------------------------------------------------------
do k=1,km
WORK1 = min(L_ROSSBY, SIGMA(:,:,k)/BTP(:,:,bid))
KAPPA_VERTICAL(:,:,k,bid) = const_eg * SIGMA(:,:,k) * WORK1**2
enddo
!----------------------------------------------------------------------
!
! use below-diabatic-layer values of KAPPA_ISOP within the surface
! diabatic layer
!
!----------------------------------------------------------------------
WORK1 = BL_DEPTH(:,:,bid)
if ( transition_layer_on ) &
WORK1 = TLT%DIABATIC_DEPTH(:,:,bid)
do k=km-1,1,-1
where ( zw(k) <= WORK1 )
KAPPA_VERTICAL(:,:,k,bid) = KAPPA_VERTICAL(:,:,k+1,bid)
endwhere
enddo
!----------------------------------------------------------------------
!
! impose lower and upper limits on KAPPA
!
!----------------------------------------------------------------------
KAPPA_VERTICAL(:,:,:,bid) = max(kappa_min_eg, KAPPA_VERTICAL(:,:,:,bid))
KAPPA_VERTICAL(:,:,:,bid) = min(kappa_max_eg, KAPPA_VERTICAL(:,:,:,bid))
end subroutine kappa_eg
!***********************************************************************
!BOP
! !IROUTINE: kappa_lon_lat_hdgr
! !INTERFACE:
subroutine kappa_lon_lat_hdgr (TMIX, this_block) 1,1
! !DESCRIPTION:
! Variable kappa parameterization based on the vertically averaged
! horizontal density gradients as a measure of baroclinicity. This
! approach is from GFDL.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TMIX ! tracers at all vertical levels
! at mixing time level
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k_min, k_max, & ! min and max vertical indices for the vertical average
k, & ! vertical loop index
i, j, & ! horizontal loop indices
kk, & ! temporary value of KMT
bid ! local block address for this sub block
real (r8), dimension(nx_block,ny_block) :: &
KMASKE, KMASKN, & ! ocean masks for the east and north faces of a T-cell
WORK, & ! work array
RHOT, & ! dRHO/dT
RHOS, & ! dRHO/dS
RHO_X, & ! grid x-direction density difference
RHO_Y, & ! grid y-direction density difference
RHO_XY_AVG ! vertically averaged horizontal density gradient
real (r8), dimension(nx_block,ny_block,2) :: &
TRACER_X, & ! grid x-direction T and S differences
TRACER_Y ! grid y-direction T and S differences
real (r8) :: &
depth_min, depth_max, & ! actual min and max depths (reflecting model
! discretization) for computing vertical
! averages at every horizontal grid point
depth_w_min, & ! protection against k_min = k_max = 0 in
depth_w_max, & ! array zw
scaling_coefficient ! multiplicative factor
!-----------------------------------------------------------------------
!
! parameter values for the formulation
!
!-----------------------------------------------------------------------
real (r8), parameter :: &
tuning_const = 2.0_r8, & ! dimensionless tuning constant
length_scale_ref = 50.0e5_r8, & ! constant length scale in cm
buoyancy_freq_ref = 4.0e-3_r8, & ! constant buoyancy frequency in 1/s
vert_average_min = 1.0e4_r8, & ! min and max depth range for the
vert_average_max = 2.0e5_r8, & ! vertical average in cm
kappa_min = 3.0e6_r8, & ! min KAPPA_LATERAL value in cm^2/s
kappa_max = 4.0e7_r8 ! max KAPPA_LATERAL value in cm^2/s
!-----------------------------------------------------------------------
!
! initialization
!
!-----------------------------------------------------------------------
bid = this_block%local_id
k_min = -1
k_max = 0
TRACER_X = c0
TRACER_Y = c0
RHO_X = c0
RHO_Y = c0
RHO_XY_AVG = c0
scaling_coefficient = tuning_const * (length_scale_ref**2) * grav &
/ ( rho_sw * buoyancy_freq_ref )
vert_average_loop: do k=1,km
if ( zt(k) >= vert_average_min .and. &
zt(k) <= vert_average_max ) then
!-----------------------------------------------------------------------
!
! start computing the vertical average of the horizontal density
! gradients
!
!-----------------------------------------------------------------------
if ( k_min == -1 ) k_min = k-1 ! k_min is the upper limit for averaging
KMASKE = merge(c1, c0, k <= KMT(:,:,bid) .and. &
k <= KMTE(:,:,bid))
KMASKN = merge(c1, c0, k <= KMT(:,:,bid) .and. &
k <= KMTN(:,:,bid))
WORK = max(-c2, TMIX(:,:,k,1))
!-----------------------------------------------------------------------
!
! compute tracer horizontal differences
!
!-----------------------------------------------------------------------
do j=1,ny_block
do i=1,nx_block-1
TRACER_X(i,j,1) = KMASKE(i,j) * ( WORK(i+1,j) &
- WORK(i,j) )
TRACER_X(i,j,2) = KMASKE(i,j) &
* ( TMIX(i+1,j,k,2) - TMIX(i,j,k,2) )
enddo
enddo
do j=1,ny_block-1
do i=1,nx_block
TRACER_Y(i,j,1) = KMASKN(i,j) * ( WORK(i,j+1) &
- WORK(i,j) )
TRACER_Y(i,j,2) = KMASKN(i,j) &
* ( TMIX(i,j+1,k,2) - TMIX(i,j,k,2) )
enddo
enddo
!-----------------------------------------------------------------------
!
! now compute horizontal density differences
!
!-----------------------------------------------------------------------
call state( k, k, TMIX(:,:,k,1), TMIX(:,:,k,2), &
this_block, DRHODT=RHOT, DRHODS=RHOS )
RHO_X = RHOT * TRACER_X(:,:,1) + RHOS * TRACER_X(:,:,2)
RHO_Y = RHOT * TRACER_Y(:,:,1) + RHOS * TRACER_Y(:,:,2)
!-----------------------------------------------------------------------
!
! compute horizontal density gradients and perform their vertical
! integral
!
!-----------------------------------------------------------------------
do j=2,ny_block-1
do i=2,nx_block-1
RHO_XY_AVG(i,j) = RHO_XY_AVG(i,j) + dz(k) * sqrt( p5 &
* ( ( RHO_X(i,j)**2 + RHO_X(i-1,j)**2 ) &
* DXTR(i,j,bid)**2 ) &
+ ( ( RHO_Y(i,j)**2 + RHO_Y(i,j-1)**2 ) &
* DYTR(i,j,bid)**2 ) )
enddo
enddo
else
if ( k_min >= 0 .and. k_max == 0 ) then
k_max = k-1 ! k_max is the lower limit for averaging
exit vert_average_loop
endif
endif
enddo vert_average_loop
!-----------------------------------------------------------------------
!
! compute the vertical average of horizontal density gradients for
! the specified vertical range
!
!-----------------------------------------------------------------------
depth_w_min = c0
depth_w_max = c0
if ( k_min > 0 ) depth_w_min = zw(k_min)
if ( k_max /= 0 ) depth_w_max = zw(k_max)
do j=2,ny_block-1
do i=2,nx_block-1
kk = KMT(i,j,bid)
if ( kk > 0 ) then
depth_min = min( depth_w_min, zw(kk) )
depth_max = min( depth_w_max, zw(kk) )
RHO_XY_AVG(i,j) = RHO_XY_AVG(i,j) &
/ ( depth_max - depth_min + eps )
else
RHO_XY_AVG(i,j) = c0
endif
enddo
enddo
!-----------------------------------------------------------------------
!
! set KAPPA_LATERAL and enforce the limits
!
!-----------------------------------------------------------------------
KAPPA_LATERAL(:,:,bid) = scaling_coefficient * RHO_XY_AVG
KAPPA_LATERAL(:,:,bid) = min( kappa_max, KAPPA_LATERAL(:,:,bid) )
KAPPA_LATERAL(:,:,bid) = max( kappa_min, KAPPA_LATERAL(:,:,bid) )
where ( KMT(:,:,bid) <= k_min ) ! KAPPA_LATERAL = kappa_min
KAPPA_LATERAL(:,:,bid) = kappa_min ! when depth < vert_average_min
endwhere
!-----------------------------------------------------------------------
!EOC
end subroutine kappa_lon_lat_hdgr
!***********************************************************************
!BOP
! !IROUTINE: kappa_lon_lat_dradius
! !INTERFACE:
subroutine kappa_lon_lat_dradius (this_block) 1
! !DESCRIPTION:
! Variable kappa parameterization based on
! \begin{equation}
! KAPPA_LATERAL = KAPPA_REF {LENGTH_SCALE \over LENGTH_SCALE_REF}
! \end{equation}
! where $LENGTH_SCALE=$ min (Rossby deformation radius, sqrt(TAREA)).
! Computation of Rossby deformation radius follows Chelton et al.(1998)
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, & ! vertical loop index
bid ! local block address for this sub block
real (r8), dimension(nx_block,ny_block) :: &
WORK1, & ! work array
LENGTH_SCALE ! mostly used as a temporary array,
! but the final content is a length scale
real (r8), parameter :: &
length_scale_ref = 15.0e+5_r8, & ! reference length scale in cm
kappa_lateral_min = 0.1e+7_r8, & ! minimum kappa in cm^2/s
kappa_lateral_max = 8.0e+7_r8 ! maximum kappa in cm^2/s
!-----------------------------------------------------------------------
!
! compute the first baroclinic gravity-wave phase speed, considering
! the full-depth integral of N (consider pi later)
!
!-----------------------------------------------------------------------
bid = this_block%local_id
LENGTH_SCALE = c0
k = 1
where ( k < KMT(:,:,bid) )
LENGTH_SCALE = LENGTH_SCALE &
+ sqrt(BUOY_FREQ_SQ(:,:,k,bid)) * dzw(k-1)
endwhere
do k=1,km
where ( k < KMT(:,:,bid) )
LENGTH_SCALE = LENGTH_SCALE &
+ sqrt(BUOY_FREQ_SQ(:,:,k,bid)) * dzw(k)
endwhere
where ( k == KMT(:,:,bid) )
LENGTH_SCALE = LENGTH_SCALE &
+ sqrt(BUOY_FREQ_SQ(:,:,k,bid)) * p5 * dz(k)
endwhere
enddo
!-----------------------------------------------------------------------
!
! now compute the Rossby radius of deformation
!
!-----------------------------------------------------------------------
where ( abs(TLAT(:,:,bid)) > c5 / radian )
LENGTH_SCALE = LENGTH_SCALE / ( pi * abs(FCORT(:,:,bid)) )
elsewhere
LENGTH_SCALE = sqrt( LENGTH_SCALE / ( c2 * pi * BTP(:,:,bid) ) )
endwhere
!-----------------------------------------------------------------------
!
! consider grid size as a possible lower limit of length scale.
! note that if the vertical integral of BUOY_FREQ_SQ is zero, then
! LENGTH_SCALE = 0. we choose to enforce limits on KAPPA later, rather
! then on LENGTH_SCALE below.
!
!-----------------------------------------------------------------------
WORK1 = min( DXT(:,:,bid), DYT(:,:,bid) )
LENGTH_SCALE = min( LENGTH_SCALE, WORK1 )
!-----------------------------------------------------------------------
!
! now compute KAPPA_LATERAL
!
!-----------------------------------------------------------------------
KAPPA_LATERAL(:,:,bid) = ah * LENGTH_SCALE / length_scale_ref
if ( kappa_isop_type == kappa_type_const ) &
KAPPA_LATERAL(:,:,bid) = ah_bolus * LENGTH_SCALE &
/ length_scale_ref
KAPPA_LATERAL(:,:,bid) = min(KAPPA_LATERAL(:,:,bid), &
kappa_lateral_max)
KAPPA_LATERAL(:,:,bid) = max(KAPPA_LATERAL(:,:,bid), &
kappa_lateral_min)
!-----------------------------------------------------------------------
!EOC
end subroutine kappa_lon_lat_dradius
!***********************************************************************
!BOP
! !IROUTINE: buoyancy_frequency_dependent_profile
! !INTERFACE:
subroutine buoyancy_frequency_dependent_profile (TMIX, this_block) 1,1
! !DESCRIPTION:
! Computes normalized buoyancy frequency ($N$) dependent profiles that
! are used to represent the vertical variation of the isopycnal and/or
! thickness diffusion coefficients. We use
! \begin{equation}
! KAPPA_VERTICAL = {{N^2}\over {N_ref^2}},
! \end{equation}
! where $N_ref$ is the value of $N$ "just" below the surface diabatic
! layer (SDL) provided that $N^2 > 0$ there. Otherwise, $N_ref$ is the
! first stable $N^2$ below SDL. KAPPA_VERTICAL is set to unity for shallower
! depths. Also, 0.1 <= KAPPA_VERTICAL <= 1.0 is enforced. Note that
! this implies KAPPA_VERTICAL = 0.1 if the local $N^2 \le 0$.
! KAPPA_VERTICAL is located at the model tracer points.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TMIX ! tracers at all vertical levels
! at mixing time level
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, & ! vertical loop index
bid ! local block address for this sub block
integer (int_kind), dimension(nx_block,ny_block) :: &
K_MIN ! k index below SDL
real (r8), dimension(nx_block,ny_block) :: &
TEMP_K, & ! temperature at level k
TEMP_KP1, & ! temperature at level k+1
RHOT, & ! dRHO/dT
RHOS, & ! dRHO/dS
BUOY_FREQ_SQ_REF, & ! reference (normalization) value
! of N^2
SDL ! surface diabatic layer (see below)
real (r8), dimension(nx_block,ny_block,km) :: &
BUOY_FREQ_SQ_NORM ! normalized N^2 defined at level interfaces
!-----------------------------------------------------------------------
!
! initialize variables
!
! note that "km+1" value in K_MIN is used as a flag, i.e. it indicates
! that a minimum value has not been found yet or that it does not
! exist.
!
!-----------------------------------------------------------------------
bid = this_block%local_id
BUOY_FREQ_SQ_NORM = c0
BUOY_FREQ_SQ_REF = c0
KAPPA_VERTICAL(:,:,:,bid) = c1
K_MIN = merge( km+1, 0, KMT(:,:,bid) /= 0 )
SDL = zw(1)
if ( vmix_itype == vmix_type_kpp ) SDL = KPP_HBLT(:,:,bid)
if ( transition_layer_on ) SDL = TLT%INTERIOR_DEPTH(:,:,bid)
if ( .not. read_n2_data ) then
!-----------------------------------------------------------------------
!
! compute N^2
!
!-----------------------------------------------------------------------
do k=1,km-1
if ( k == 1 ) &
TEMP_K = max( -c2, TMIX(:,:,k,1) )
TEMP_KP1 = max( -c2, TMIX(:,:,k+1,1) )
call state( k, k+1, TMIX(:,:,k,1), TMIX(:,:,k,2), &
this_block, DRHODT=RHOT, DRHODS=RHOS )
where ( k < KMT(:,:,bid) )
BUOY_FREQ_SQ(:,:,k,bid) = max( c0, - grav * dzwr(k) &
* ( RHOT * ( TEMP_K - TEMP_KP1 ) &
+ RHOS * ( TMIX(:,:,k, 2) - TMIX(:,:,k+1,2) ) ) )
endwhere
TEMP_K = TEMP_KP1
enddo
endif
!-----------------------------------------------------------------------
!
! determine the reference buoyancy frequency and the associated
! level index (most likely just below SDL)
!
!-----------------------------------------------------------------------
do k=1,km-1
where ( ( K_MIN == km+1 ) .and. ( zw(k) > SDL ) .and. &
( k <= KMT(:,:,bid) ) .and. &
( BUOY_FREQ_SQ(:,:,k,bid) > c0 ) )
BUOY_FREQ_SQ_REF = BUOY_FREQ_SQ(:,:,k,bid)
K_MIN = k
endwhere
enddo
!-----------------------------------------------------------------------
!
! now compute the normalized profiles at the interfaces
!
!-----------------------------------------------------------------------
do k=1,km-1
where ( ( k >= K_MIN ) .and. ( k < KMT(:,:,bid) ) .and. &
( BUOY_FREQ_SQ_REF /= c0 ) )
BUOY_FREQ_SQ_NORM(:,:,k) = &
max( BUOY_FREQ_SQ(:,:,k,bid) / BUOY_FREQ_SQ_REF, 0.1_r8 )
BUOY_FREQ_SQ_NORM(:,:,k) = &
min( BUOY_FREQ_SQ_NORM(:,:,k), c1 )
elsewhere
BUOY_FREQ_SQ_NORM(:,:,k) = c1
endwhere
enddo
do k=1,km-1
where ( k == KMT(:,:,bid)-1 )
BUOY_FREQ_SQ_NORM(:,:,k+1) = BUOY_FREQ_SQ_NORM(:,:,k)
endwhere
enddo
!-----------------------------------------------------------------------
!
! finally, do not average interface values of BUOY_FREQ_SQ_NORM to
! the tracer points, but instead copy from above to preserve the
! extrema.
!
!-----------------------------------------------------------------------
do k=2,km
where ( ( k > K_MIN ) .and. ( k <= KMT(:,:,bid) ) )
KAPPA_VERTICAL(:,:,k,bid) = BUOY_FREQ_SQ_NORM(:,:,k-1)
endwhere
enddo
!-----------------------------------------------------------------------
!EOC
end subroutine buoyancy_frequency_dependent_profile
!***********************************************************************
!BOP
! !IROUTINE: transition_layer
! !INTERFACE:
subroutine transition_layer ( this_block ) 1,3
! !DESCRIPTION:
! Compute transition layer related fields. the related algorithms
! should work even with zero transition layer depth.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, kk, & ! loop indices
bid ! local block address for this sub block
integer (int_kind), dimension(nx_block,ny_block) :: &
K_START, & ! work arrays for TLT%K_LEVEL and
K_SUB ! TLT%ZTW, respectively
logical (log_kind), dimension(nx_block,ny_block) :: &
COMPUTE_TLT ! flag
real (r8), dimension(nx_block,ny_block) :: &
WORK ! work space for TLT%THICKNESS
real (r8), dimension(2) :: &
reference_depth ! zt or zw
!-----------------------------------------------------------------------
!
! initialize various quantities
!
!-----------------------------------------------------------------------
bid = this_block%local_id
K_START = 0
K_SUB = 0
TLT%THICKNESS(:,:,bid) = c0
TLT%INTERIOR_DEPTH(:,:,bid) = c0
TLT%K_LEVEL(:,:,bid) = 0
TLT%ZTW(:,:,bid) = 0
COMPUTE_TLT = merge(.true., .false., KMT(:,:,bid) /= 0)
!-----------------------------------------------------------------------
!
! initial pass to determine the minimum transition layer thickness.
! the vertical level at or below which the interior region starts
! is also computed.
!
!-----------------------------------------------------------------------
do k=1,km
where ( COMPUTE_TLT .and. &
TLT%DIABATIC_DEPTH(:,:,bid) < zw(k) )
K_START = k+1
K_SUB = ktp
TLT%THICKNESS(:,:,bid) = zw(k) - TLT%DIABATIC_DEPTH(:,:,bid)
TLT%K_LEVEL(:,:,bid) = k
TLT%ZTW(:,:,bid) = 2
COMPUTE_TLT = .false.
endwhere
where ( k /= 1 .and. K_START == k+1 .and. &
TLT%DIABATIC_DEPTH(:,:,bid) < zt(k) )
K_START = k
K_SUB = kbt
TLT%THICKNESS(:,:,bid) = zt(k) - TLT%DIABATIC_DEPTH(:,:,bid)
TLT%K_LEVEL(:,:,bid) = k
TLT%ZTW(:,:,bid) = 1
endwhere
enddo
#ifdef CCSMCOUPLED
#ifndef _HIRES
if ( any(COMPUTE_TLT) ) then
call shr_sys_abort ('Incorrect DIABATIC_DEPTH value in TLT' &
// ' computation')
endif
#endif
#endif
!-----------------------------------------------------------------------
!
! using R * |S| as the vertical displacement length scale associated
! with a horizontal displacement equivalent to the Rossby deformation
! radius, R, determine if the transition layer thickness needs to be
! modified. |S| is the absolute value of the slope of the isopycnal
! surfaces. First, start with columns where K_SUB = kbt.
!
!-----------------------------------------------------------------------
where ( KMT(:,:,bid) == 0 .or. K_START > KMT(:,:,bid) .or. &
( K_START == KMT(:,:,bid) .and. K_SUB == kbt ) )
COMPUTE_TLT = .false.
elsewhere
COMPUTE_TLT = .true.
endwhere
do k=1,km-1
WORK = c0
where ( COMPUTE_TLT .and. K_SUB == kbt .and. &
K_START < KMT(:,:,bid) .and. K_START == k )
WORK = max(SLA_SAVE(:,:,kbt,k,bid), &
SLA_SAVE(:,:,ktp,k+1,bid)) * RB(:,:,bid)
endwhere
where ( WORK /= c0 .and. &
TLT%DIABATIC_DEPTH(:,:,bid) < (zw(k) - WORK) )
COMPUTE_TLT = .false.
endwhere
where ( WORK /= c0 .and. &
TLT%DIABATIC_DEPTH(:,:,bid) >= (zw(k) - WORK) )
K_START = K_START + 1
K_SUB = ktp
TLT%THICKNESS(:,:,bid) = zw(k) - TLT%DIABATIC_DEPTH(:,:,bid)
TLT%K_LEVEL(:,:,bid) = k
TLT%ZTW(:,:,bid) = 2
endwhere
enddo
!-----------------------------------------------------------------------
!
! now consider the deeper levels
!
!-----------------------------------------------------------------------
do k=2,km
reference_depth(ktp) = zt(k)
reference_depth(kbt) = zw(k)
do kk=ktp,kbt
WORK = c0
if (kk == ktp) then
where ( COMPUTE_TLT .and. K_START <= KMT(:,:,bid) .and. &
K_START == k )
WORK = max(SLA_SAVE(:,:,ktp,k,bid), &
SLA_SAVE(:,:,kbt,k,bid)) * RB(:,:,bid)
endwhere
else
where ( COMPUTE_TLT .and. K_START < KMT(:,:,bid) .and. &
K_START == k )
WORK = max(SLA_SAVE(:,:,kbt,k,bid), &
SLA_SAVE(:,:,ktp,k+1,bid)) * RB(:,:,bid)
endwhere
where ( COMPUTE_TLT .and. K_START == KMT(:,:,bid) .and. &
K_START == k )
WORK = SLA_SAVE(:,:,kbt,k,bid) * RB(:,:,bid)
endwhere
endif
where ( WORK /= c0 .and. &
TLT%DIABATIC_DEPTH(:,:,bid) < (reference_depth(kk) - WORK) )
COMPUTE_TLT = .false.
endwhere
where ( WORK /= c0 .and. &
TLT%DIABATIC_DEPTH(:,:,bid) >= (reference_depth(kk) - WORK) )
TLT%THICKNESS(:,:,bid) = reference_depth(kk) &
- TLT%DIABATIC_DEPTH(:,:,bid)
TLT%K_LEVEL(:,:,bid) = k
TLT%ZTW(:,:,bid) = kk
endwhere
enddo
where ( COMPUTE_TLT .and. K_START == k )
K_START = K_START + 1
endwhere
enddo
#ifdef CCSMCOUPLED
#ifndef _HIRES
if ( any(COMPUTE_TLT) ) then
call shr_sys_abort ('Incorrect TLT computations')
endif
#endif
#endif
!-----------------------------------------------------------------------
!
! compute the depth at which the interior, adiabatic region starts
!
!-----------------------------------------------------------------------
do k=1,km
where ( TLT%K_LEVEL(:,:,bid) == k .and. &
TLT%ZTW(:,:,bid) == 1 )
TLT%INTERIOR_DEPTH(:,:,bid) = zt(k)
endwhere
where ( TLT%K_LEVEL(:,:,bid) == k .and. &
TLT%ZTW(:,:,bid) == 2 )
TLT%INTERIOR_DEPTH(:,:,bid) = zw(k)
endwhere
enddo
COMPUTE_TLT = .false.
where ( KMT(:,:,bid) /= 0 .and. &
TLT%INTERIOR_DEPTH(:,:,bid) == c0 ) &
COMPUTE_TLT = .true.
where ( KMT(:,:,bid) == 0 .and. &
TLT%INTERIOR_DEPTH(:,:,bid) /= c0 ) &
COMPUTE_TLT = .true.
#ifdef CCSMCOUPLED
#ifndef _HIRES
if ( any(COMPUTE_TLT) ) then
call shr_sys_abort ('Incorrect TLT%INTERIOR_DEPTH computation')
endif
#endif
#endif
!-----------------------------------------------------------------------
!EOC
end subroutine transition_layer
!***********************************************************************
!BOP
! !IROUTINE: merged_streamfunction
! !INTERFACE:
subroutine merged_streamfunction ( this_block ) 1
! !DESCRIPTION:
! Construct a merged streamfunction that has the appropriate
! behavior in the surface diabatic region, transition layer, and
! adiabatic interior
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, kk, & ! loop indices
bid ! local block address for this sub block
real (r8), dimension(nx_block,ny_block,2) :: &
WORK1, WORK2, WORK3, WORK4 ! work arrays
real (r8), dimension(nx_block,ny_block) :: &
WORK2_NEXT, WORK4_NEXT ! WORK2 or WORK4 at next level
real (r8), dimension(nx_block,ny_block) :: &
WORK5, WORK6, WORK7 ! more work arrays
logical (log_kind), dimension(nx_block,ny_block) :: &
LMASK ! flag
real (r8), dimension(2) :: &
reference_depth ! zt or zw
!-----------------------------------------------------------------------
!
! initialize various quantities
!
!-----------------------------------------------------------------------
bid = this_block%local_id
SF_SLX(:,:,:,:,:,bid) = c0
SF_SLY(:,:,:,:,:,bid) = c0
WORK1 = c0
WORK2 = c0
WORK3 = c0
WORK4 = c0
WORK5 = c0
WORK6 = c0
WORK7 = c0
WORK2_NEXT = c0
WORK4_NEXT = c0
!-----------------------------------------------------------------------
!
! compute the interior streamfunction and its first derivative at the
! INTERIOR_DEPTH level. WORK1 and WORK2 contain the streamfunction
! and its first derivative, respectively, for the zonal component
! for the east and west sides of a grid cell. WORK3 and WORK4 are
! the corresponding fields for the meridional component for the
! north and south sides of a grid cell. Note that these definitions
! include a "dz". Also, the first derivative computations assume
! that the streamfunctions are located in the middle of the top or
! bottom half of a grid cell, hence a factor of two in WORK2 and
! WORK4 calculations.
!
!-----------------------------------------------------------------------
do k=1,km-1
do kk=1,2
LMASK = TLT%K_LEVEL(:,:,bid) == k .and. &
TLT%K_LEVEL(:,:,bid) < KMT(:,:,bid) .and. &
TLT%ZTW(:,:,bid) == 1
where ( LMASK )
WORK1(:,:,kk) = KAPPA_THIC(:,:,kbt,k,bid) &
* SLX(:,:,kk,kbt,k,bid) * dz(k)
WORK2(:,:,kk) = c2 * dzwr(k) * ( WORK1(:,:,kk) &
- KAPPA_THIC(:,:,ktp,k+1,bid) * SLX(:,:,kk,ktp,k+1,bid) &
* dz(k+1) )
WORK2_NEXT = c2 * ( &
KAPPA_THIC(:,:,ktp,k+1,bid) * SLX(:,:,kk,ktp,k+1,bid) - &
KAPPA_THIC(:,:,kbt,k+1,bid) * SLX(:,:,kk,kbt,k+1,bid) )
WORK3(:,:,kk) = KAPPA_THIC(:,:,kbt,k,bid) &
* SLY(:,:,kk,kbt,k,bid) * dz(k)
WORK4(:,:,kk) = c2 * dzwr(k) * ( WORK3(:,:,kk) &
- KAPPA_THIC(:,:,ktp,k+1,bid) * SLY(:,:,kk,ktp,k+1,bid) &
* dz(k+1) )
WORK4_NEXT = c2 * ( &
KAPPA_THIC(:,:,ktp,k+1,bid) * SLY(:,:,kk,ktp,k+1,bid) - &
KAPPA_THIC(:,:,kbt,k+1,bid) * SLY(:,:,kk,kbt,k+1,bid) )
endwhere
where ( LMASK .and. abs(WORK2_NEXT) < abs(WORK2(:,:,kk)) ) &
WORK2(:,:,kk) = WORK2_NEXT
where ( LMASK .and. abs(WORK4_NEXT) < abs(WORK4(:,:,kk)) ) &
WORK4(:,:,kk) = WORK4_NEXT
LMASK = TLT%K_LEVEL(:,:,bid) == k .and. &
TLT%K_LEVEL(:,:,bid) < KMT(:,:,bid) .and. &
TLT%ZTW(:,:,bid) == 2
where ( LMASK )
WORK1(:,:,kk) = KAPPA_THIC(:,:,ktp,k+1,bid) &
* SLX(:,:,kk,ktp,k+1,bid)
WORK2(:,:,kk) = c2 * ( WORK1(:,:,kk) &
- ( KAPPA_THIC(:,:,kbt,k+1,bid) &
* SLX(:,:,kk,kbt,k+1,bid) ) )
WORK1(:,:,kk) = WORK1(:,:,kk) * dz(k+1)
WORK3(:,:,kk) = KAPPA_THIC(:,:,ktp,k+1,bid) &
* SLY(:,:,kk,ktp,k+1,bid)
WORK4(:,:,kk) = c2 * ( WORK3(:,:,kk) &
- ( KAPPA_THIC(:,:,kbt,k+1,bid) &
* SLY(:,:,kk,kbt,k+1,bid) ) )
WORK3(:,:,kk) = WORK3(:,:,kk) * dz(k+1)
endwhere
LMASK = LMASK .and. TLT%K_LEVEL(:,:,bid) + 1 < KMT(:,:,bid)
where ( LMASK )
WORK2_NEXT = c2 * dzwr(k+1) * ( &
KAPPA_THIC(:,:,kbt,k+1,bid) * SLX(:,:,kk,kbt,k+1,bid) * dz(k+1) - &
KAPPA_THIC(:,:,ktp,k+2,bid) * SLX(:,:,kk,ktp,k+2,bid) * dz(k+2))
WORK4_NEXT = c2 * dzwr(k+1) * ( &
KAPPA_THIC(:,:,kbt,k+1,bid) * SLY(:,:,kk,kbt,k+1,bid) * dz(k+1) - &
KAPPA_THIC(:,:,ktp,k+2,bid) * SLY(:,:,kk,ktp,k+2,bid) * dz(k+2))
endwhere
where ( LMASK .and. abs(WORK2_NEXT) < abs(WORK2(:,:,kk)) ) &
WORK2(:,:,kk) = WORK2_NEXT
where ( LMASK .and. abs(WORK4_NEXT) < abs(WORK4(:,:,kk)) ) &
WORK4(:,:,kk) = WORK4_NEXT
enddo
enddo
!-----------------------------------------------------------------------
!
! compute the depth independent interpolation factors used in the
! linear and quadratic interpolations within the diabatic and
! transition regions, respectively.
!
!-----------------------------------------------------------------------
WORK5 = c0
where (KMT(:,:,bid) /= 0)
WORK5(:,:) = c1 / ( c2 * TLT%DIABATIC_DEPTH(:,:,bid) &
+ TLT%THICKNESS(:,:,bid) )
endwhere
WORK6 = c0
where ((KMT(:,:,bid) /= 0) .AND. (TLT%THICKNESS(:,:,bid) > eps))
WORK6(:,:) = WORK5(:,:) / TLT%THICKNESS(:,:,bid)
endwhere
!-----------------------------------------------------------------------
!
! start of interpolation to construct the merged streamfunction
!
!-----------------------------------------------------------------------
do k=1,km
reference_depth(ktp) = zt(k) - p25 * dz(k)
reference_depth(kbt) = zt(k) + p25 * dz(k)
do kk=ktp,kbt
!-----------------------------------------------------------------------
!
! diabatic region: use linear interpolation (in streamfunction)
!
!-----------------------------------------------------------------------
where ( reference_depth(kk) <= TLT%DIABATIC_DEPTH(:,:,bid) &
.and. k <= KMT(:,:,bid) )
SF_SLX(:,:,1,kk,k,bid) = reference_depth(kk) * WORK5 &
* ( c2 * WORK1(:,:,1) + TLT%THICKNESS(:,:,bid) &
* WORK2(:,:,1) )
SF_SLX(:,:,2,kk,k,bid) = reference_depth(kk) * WORK5 &
* ( c2 * WORK1(:,:,2) + TLT%THICKNESS(:,:,bid) &
* WORK2(:,:,2) )
SF_SLY(:,:,1,kk,k,bid) = reference_depth(kk) * WORK5 &
* ( c2 * WORK3(:,:,1) + TLT%THICKNESS(:,:,bid) &
* WORK4(:,:,1) )
SF_SLY(:,:,2,kk,k,bid) = reference_depth(kk) * WORK5 &
* ( c2 * WORK3(:,:,2) + TLT%THICKNESS(:,:,bid) &
* WORK4(:,:,2) )
endwhere
!-----------------------------------------------------------------------
!
! transition layer: use quadratic interpolation (in streamfunction)
!
!-----------------------------------------------------------------------
where ( reference_depth(kk) > TLT%DIABATIC_DEPTH(:,:,bid) &
.and. reference_depth(kk) <= TLT%INTERIOR_DEPTH(:,:,bid) &
.and. k <= KMT(:,:,bid) )
WORK7 = (TLT%DIABATIC_DEPTH(:,:,bid) &
- reference_depth(kk))**2
SF_SLX(:,:,1,kk,k,bid) = - WORK7 * WORK6 &
* ( WORK1(:,:,1) + TLT%INTERIOR_DEPTH(:,:,bid) &
* WORK2(:,:,1) ) &
+ reference_depth(kk) * WORK5 &
* ( c2 * WORK1(:,:,1) + TLT%THICKNESS(:,:,bid) &
* WORK2(:,:,1) )
SF_SLX(:,:,2,kk,k,bid) = - WORK7 * WORK6 &
* ( WORK1(:,:,2) + TLT%INTERIOR_DEPTH(:,:,bid) &
* WORK2(:,:,2) ) &
+ reference_depth(kk) * WORK5 &
* ( c2 * WORK1(:,:,2) + TLT%THICKNESS(:,:,bid) &
* WORK2(:,:,2) )
SF_SLY(:,:,1,kk,k,bid) = - WORK7 * WORK6 &
* ( WORK3(:,:,1) + TLT%INTERIOR_DEPTH(:,:,bid) &
* WORK4(:,:,1) ) &
+ reference_depth(kk) * WORK5 &
* ( c2 * WORK3(:,:,1) + TLT%THICKNESS(:,:,bid) &
* WORK4(:,:,1) )
SF_SLY(:,:,2,kk,k,bid) = - WORK7 * WORK6 &
* ( WORK3(:,:,2) + TLT%INTERIOR_DEPTH(:,:,bid) &
* WORK4(:,:,2) ) &
+ reference_depth(kk) * WORK5 &
* ( c2 * WORK3(:,:,2) + TLT%THICKNESS(:,:,bid) &
* WORK4(:,:,2) )
endwhere
!-----------------------------------------------------------------------
!
! interior, adiabatic region: no interpolation is needed. note that
! "dzw" is introduced here, too, for consistency.
!
!-----------------------------------------------------------------------
where ( reference_depth(kk) > TLT%INTERIOR_DEPTH(:,:,bid) &
.and. k <= KMT(:,:,bid) )
SF_SLX(:,:,1,kk,k,bid) = KAPPA_THIC(:,:,kk,k,bid) &
* SLX(:,:,1,kk,k,bid) * dz(k)
SF_SLX(:,:,2,kk,k,bid) = KAPPA_THIC(:,:,kk,k,bid) &
* SLX(:,:,2,kk,k,bid) * dz(k)
SF_SLY(:,:,1,kk,k,bid) = KAPPA_THIC(:,:,kk,k,bid) &
* SLY(:,:,1,kk,k,bid) * dz(k)
SF_SLY(:,:,2,kk,k,bid) = KAPPA_THIC(:,:,kk,k,bid) &
* SLY(:,:,2,kk,k,bid) * dz(k)
endwhere
enddo ! end of kk-loop
enddo ! end of k-loop
!-----------------------------------------------------------------------
!EOC
end subroutine merged_streamfunction
!***********************************************************************
!BOP
! !IROUTINE: apply_vertical_profile_to_isop_hor_diff
! !INTERFACE:
subroutine apply_vertical_profile_to_isop_hor_diff ( this_block ) 1
! !DESCRIPTION:
! Apply vertical tapers to KAPPA_ISOP and HOR_DIFF based on their
! vertical location with respect to the diabatic, transition, and
! adiabatic regions.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block info for this sub block
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, kk, & ! loop indices
bid ! local block address for this sub block
real (r8), dimension(2) :: &
reference_depth
bid = this_block%local_id
!-----------------------------------------------------------------------
!
! start of tapering
!
!-----------------------------------------------------------------------
do k=1,km
reference_depth(ktp) = zt(k) - p25 * dz(k)
reference_depth(kbt) = zt(k) + p25 * dz(k)
do kk=ktp,kbt
!-----------------------------------------------------------------------
!
! diabatic region: no isopycnal diffusion
!
!-----------------------------------------------------------------------
where ( reference_depth(kk) <= TLT%DIABATIC_DEPTH(:,:,bid) &
.and. k <= KMT(:,:,bid) )
KAPPA_ISOP(:,:,kk,k,bid) = c0
endwhere
!-----------------------------------------------------------------------
!
! transition layer: a linear combination of isopcynal and horizontal
! diffusion coefficients
!
!-----------------------------------------------------------------------
where ( reference_depth(kk) > TLT%DIABATIC_DEPTH(:,:,bid) &
.and. reference_depth(kk) <= TLT%INTERIOR_DEPTH(:,:,bid) &
.and. k <= KMT(:,:,bid) .and. &
TLT%THICKNESS(:,:,bid) > eps )
HOR_DIFF(:,:,kk,k,bid) = ( TLT%INTERIOR_DEPTH(:,:,bid) &
- reference_depth(kk) ) * HOR_DIFF(:,:,kk,k,bid) &
/ TLT%THICKNESS(:,:,bid)
KAPPA_ISOP(:,:,kk,k,bid) = ( reference_depth(kk) &
- TLT%DIABATIC_DEPTH(:,:,bid) ) &
* KAPPA_ISOP(:,:,kk,k,bid) / TLT%THICKNESS(:,:,bid)
endwhere
!-----------------------------------------------------------------------
!
! interior region: no horizontal diffusion
!
!-----------------------------------------------------------------------
where ( reference_depth(kk) > TLT%INTERIOR_DEPTH(:,:,bid) &
.and. k <= KMT(:,:,bid) )
HOR_DIFF(:,:,kk,k,bid) = c0
endwhere
enddo ! end of kk-loop
enddo ! end of k-loop
!-----------------------------------------------------------------------
!EOC
end subroutine apply_vertical_profile_to_isop_hor_diff
!***********************************************************************
end module hmix_gm
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||