!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module ocean_nphysicsA_mod
!
! Stephen M. Griffies
!
!
! Russell Fiedler
!
!
!
! Thickness weighted and density weighted time tendency for tracer
! from Laplacian neutral diffusion + Laplacian GM skew-diffusion.
!
!
!
! This module computes the cell thickness weighted and density
! weighted tracer tendency from small angle Laplacian neutral diffusion
! plus Laplacian GM skew-diffusion. The algorithms are based on
! mom4p0d methods. The fundamental differences from the ocean_nphysicsB
! methods relate to the handling of fluxes near the domain boundaries.
!
!
!
!
!
! S.M. Griffies, A. Gnanadesikan, R.C. Pacanowski, V. Larichev,
! J.K. Dukowicz, and R.D. Smith
! Isoneutral diffusion in a z-coordinate ocean model
! Journal of Physical Oceanography (1998) vol 28 pages 805-830
!
!
!
! S.M. Griffies
! The Gent-McWilliams Skew-flux
! Journal of Physical Oceanography (1998) vol 28 pages 831-841
!
!
!
! S.M. Griffies
! Fundamentals of Ocean Climate Models (2004)
! Princeton University Press
!
!
!
! S.M. Griffies
! Elements of mom4p1 (2008)
!
!
!
! G. Danabasoglu and J. C. McWilliams
! Sensitivity of the global ocean circulation to
! parameterizations of mesoscale tracer transports
! Journal of Climate (1995) vol 8 pages 2967--2987
!
!
!
! Gerdes, Koberle, and Willebrand
! The influence of numerical advection schemes on the results of ocean
! general circulation models, Climate Dynamics (1991), vol. 5,
! pages 211--226.
!
!
!
! Numerical implementation of the flux components follows the triad
! approach documented in the references and implemented in MOM2 and MOM3.
! The MOM4 algorithm accounts for partial bottom cells and generalized
! orthogonal horizontal coordinates.
!
!
!
! neutral_physics_simple=.true. requires aredi_equal_agm=.true.
! neutral_physics_simple=.true. results in down-gradient
! horizontal flux components. This setting reduces the overall cost
! of the neutral physics scheme, but it is not used at GFDL
! anymore, since we favor methods whereby treatment of GM and Redi
! in the boundary layers are distinct.
!
!
!
! In steep slope regions, neutral diffusive fluxes are tapered to
! zero with the tanh taper of Danabasoglu and McWilliams (1995) or the
! quadratic scheme of Gerdes, Koberle, and Willebrand. However, if
! neutral_physics_simple=.false., the GM skew-diffusive fluxes
! can remain nonzero if have neutral_linear_gm_taper=.true.
!
!
!
!
!
!
!
! Must be true to use this module. Default is false.
!
!
! For printing starting and ending checksums for restarts
!
!
!
! Must be true to use GM skewsion. Set to false if wish to
! incorporate the "GM-effect" through form drag, as in
! ocean_form_drag module. Default use_gm_skew=.true.
!
!
!
! To compute all contributions from neutral diffusion explicitly in time, including
! the K33 diagonal piece. This approach is available only when have small time
! steps and/or running with just a single tracer. It is for testing purposes.
!
!
!
! If .true. then must have aredi_equal_agm=.true.. The horizontal fluxes are then
! computed as horizontal downgradient diffusive fluxes regardless the neutral slope.
! This approach precluds one from being able to have the GM-skew fluxes remain active
! in the steep sloped regions, thus shutting off their effects to reduce the slopes
! of isopycnals in convective and mixed layer regimes. It is for this reason that
! neutral_physics_simple=.false. is the recommended default in MOM4.
!
!
!
! When tracer falls outside a specified range, revert to horizontal
! diffusive fluxes at this cell. This is an ad hoc and incomplete attempt
! to maintain monotonicity with the neutral physics scheme.
! Default neutral_physics_limit=.true.
!
!
! If .true. then this logical reduces the neutral fluxes to
! horizontal/vertical diffusion next to boundaries.
! This approach has been found to reduce spurious
! extrema resulting from truncation of triads used to compute
! a neutral flux component. Default tmask_neutral_on=.false.
!
!
!
! Set to true to use the tanh tapering scheme of Danabasoglu and McWilliams.
! Default is true.
!
!
! Set to true to use the quadradic tapering scheme of Gerdes, Koberle, and Willebrand.
! Default is false.
!
!
!
! If .true. then with neutral_physics_simple=.false., will linearly taper GM
! skew fluxes towards the surface within regions of steep neutral slopes.
! This approach leads to a constant horizontal eddy-induced velocity in
! the steeply sloping regions and is recommended for realistic simulations.
!
!
! If .true. then with neutral_physics_simple=.false., will apply a sine-taper
! to GM and neutral diffusive fluxes in regions where the penetration depth
! of eddies is deeper than the grid point. This method is essential for
! fine vertical resolution grids.
!
!
!
! Minimum depth of a surface turbulent boundary layer
! used in the transition of the neutral physics fluxes
! to the surface. Note that in mom4p0,
! turb_blayer_min was always set to zero.
!
!
!
! Diagnose properties of the neutral physics boundary layer, whether have
! neutral_linear_gm_taper or neutral_sine_taper true or not.
!
!
!
! For cases with neutral_physics_simple=.false., then neutral_taper_diagonal=.true.
! will taper the diagonal pieces of the horizontal flux components when neutral slopes
! are steep. With neutral_taper_diagonal=.false., then the horizontal flux components will
! remain enabled for all slopes, thus producing horizontal downgradient diffusion in
! regions of vertical neutral directions.
!
!
!
! The units for writing out the transport. Either in
! Sv (10^9 kg/s) or mks (kg/s). Note the mks unit is requested
! for CMIP5 purposes.
! Default transport_units = 'Sv'.
!
!
!
use constants_mod, only: epsln, pi, grav
use diag_manager_mod, only: register_diag_field, register_static_field, send_data, need_data
use fms_mod, only: FATAL, WARNING, NOTE
use fms_mod, only: open_namelist_file, check_nml_error, close_file, write_version_number
use mpp_domains_mod, only: mpp_update_domains
use mpp_domains_mod, only: CGRID_NE, WUPDATE, SUPDATE, EUPDATE, NUPDATE
use mpp_domains_mod, only: cyclic_global_domain, global_data_domain
use mpp_mod, only: mpp_error, mpp_chksum, stdout, stdlog
use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE
use time_manager_mod, only: set_time, time_type, increment_time, operator ( + )
use ocean_domains_mod, only: get_local_indices, set_ocean_domain
use ocean_nphysics_util_mod, only: ocean_nphysics_coeff_init, ocean_nphysics_coeff_end
use ocean_nphysics_util_mod, only: neutral_slopes, tracer_derivs
use ocean_nphysics_util_mod, only: compute_eady_rate, compute_baroclinicity
use ocean_nphysics_util_mod, only: compute_rossby_radius, compute_bczone_radius
use ocean_nphysics_util_mod, only: compute_diffusivity, ocean_nphysics_util_restart
use ocean_nphysics_util_mod, only: transport_on_nrho_gm, transport_on_rho_gm, transport_on_theta_gm
use ocean_operators_mod, only: FAX, FAY, FMX, FMY, BDX_ET, BDY_NT
use ocean_parameters_mod, only: missing_value, onehalf, onefourth, oneeigth, DEPTH_BASED
use ocean_parameters_mod, only: rho0r, rho0
use ocean_sigma_transport_mod, only: tmask_sigma_on, tmask_sigma
use ocean_types_mod, only: ocean_grid_type, ocean_domain_type, ocean_density_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_thickness_type
use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type
use ocean_types_mod, only: tracer_2d_type, tracer_3d_0_nk_type, tracer_3d_1_nk_type
use ocean_util_mod, only: write_timestamp
use ocean_workspace_mod, only: wrk1, wrk3, wrk1_v, wrk1_2d
implicit none
public ocean_nphysicsA_init
public ocean_nphysicsA_end
public nphysicsA
public ocean_nphysicsA_restart
private fx_flux
private fy_flux
private fz_flux
private fz_terms
private neutral_blayer
private slope_function_gm
private gm_velocity
private
type(ocean_grid_type), pointer :: Grd => NULL()
type(ocean_domain_type), pointer :: Dom => NULL()
type(ocean_domain_type), save :: Dom_flux
integer :: num_prog_tracers = 0
! clock ids
integer :: id_clock_neutral_blayer
integer :: id_clock_neutral_blayer_1
integer :: id_clock_neutral_blayer_2
integer :: id_clock_neutral_blayer_3
integer :: id_clock_fz_terms
integer :: id_clock_fx_flux
integer :: id_clock_fy_flux
integer :: id_clock_fz_flux
! diagnostic manager ids
logical :: used
integer :: id_k33_explicit =-1
integer :: id_ah_bdy =-1
integer :: id_ustar =-1
integer :: id_vstar =-1
integer :: id_wstar =-1
integer :: id_tx_trans_gm =-1
integer :: id_ty_trans_gm =-1
integer :: id_depth_blayer_base =-1
integer :: id_eddy_depth =-1
integer :: id_steep_depth =-1
integer :: id_slope_blayer_base =-1
integer :: id_grav_agm_dz_sx_drhodx =-1
integer :: id_grav_agm_dz_sy_drhody =-1
integer :: id_gm_eddy_ke_source =-1
integer :: id_slopex_drhodx =-1
integer :: id_slopey_drhody =-1
integer :: id_neutral_rho_nphysics
integer :: id_wdian_rho_nphysics
integer, dimension(:), allocatable :: id_neutral_physics
integer, dimension(:), allocatable :: id_k33_implicit
integer, dimension(:), allocatable :: id_flux_x ! for i-directed heat flux from neutral physics
integer, dimension(:), allocatable :: id_flux_y ! for j-directed heat flux from neutral physics
integer, dimension(:), allocatable :: id_flux_x_int_z ! for vertically integrated i-directed tracer flux
integer, dimension(:), allocatable :: id_flux_y_int_z ! for vertically integrated j-directed tracer flux
#include
#ifdef MOM4_STATIC_ARRAYS
real, dimension(isd:ied,jsd:jed,nk,0:1) :: delqc !density weighted (kg/m^3) quarter cell thickness(m)
real, dimension(isd:ied,jsd:jed,0:nk) :: dzwtr !(1/dzwt)(m^-1)
real, dimension(isd:ied,jsd:jed,0:1) :: dtew !grid distances from T-point to cell faces (m)
real, dimension(isd:ied,jsd:jed,0:1) :: dtns !grid distances from T-point to cell faces (m)
real, dimension(isd:ied,jsd:jed,0:1) :: dtwedyt !horizontal areas (m^2) of quarter cell
real, dimension(isd:ied,jsd:jed,0:1) :: dxtdtsn !horizontal areas (m^2) of quarter cell
real, dimension(isd:ied,jsd:jed,nk) :: slopex_drhodx !3D array of slopex * drhodx for diagnostics
real, dimension(isd:ied,jsd:jed,nk) :: slopey_drhody !3D array of slopey * drhody for diagnostics
real, dimension(isd:ied,jsd:jed,nk) :: grav_agm_dz_sx_drhodx !3D array of grav*agm*dz*slopex*drhodx for diagnostics
real, dimension(isd:ied,jsd:jed,nk) :: grav_agm_dz_sy_drhody !3D array of grav*agm*dz*slopey*drhody for diagnostics
real, dimension(isd:ied,jsd:jed,nk) :: aredi_array !3D array of redi diffusivities (m^2/sec)
real, dimension(isd:ied,jsd:jed,nk) :: agm_array !3D array of gm diffusivities (m^2/sec)
real, dimension(isd:ied,jsd:jed) :: ah_array !2D array of micom horizontal diffusivities (m^2/sec)
real, dimension(isd:ied,jsd:jed) :: bczone_radius !for bczone calculation (m)
real, dimension(isd:ied,jsd:jed) :: rossby_radius !first baroclinic Rossby radius (m)
real, dimension(isd:ied,jsd:jed) :: rossby_radius_raw !first baroclinic Rossby radius (m) without min/max settings
real, dimension(isd:ied,jsd:jed,nk) :: eady_termx !rho_z*(S_x)^2 for computing Eady growth rate
real, dimension(isd:ied,jsd:jed,nk) :: eady_termy !rho_z*(S_y)^2 for computing Eady growth rate
real, dimension(isd:ied,jsd:jed) :: baroclinic_termx !intermediate term for computing vert ave baroclinicity
real, dimension(isd:ied,jsd:jed) :: baroclinic_termy !intermediate term for computing vert ave baroclinicity
real, dimension(isd:ied,jsd:jed) :: grid_length !grid length scale (m)
real, dimension(isd:ied,jsd:jed,nk) :: drhodT !drho/dtheta (kg/m^3/C)
real, dimension(isd:ied,jsd:jed,nk) :: drhodS !drho/dsalinity (kg/m^3/psu)
real, dimension(isd:ied,jsd:jed,nk,0:1) :: drhodzb !vertical neutral density gradient (kg/m^3/m)
real, dimension(isd:ied,jsd:jed,nk,0:1) :: drhodzh !vertical neutral density gradient (kg/m^3/m)
real, dimension(isd:ied,jsd:jed,nk) :: K33_implicit !density weighted (kg/m^3) implicit in time
!diagonal term in redi tensor (m^2/sec)
real, dimension(isd:ied,jsd:jed,nk) :: K33_explicit !density weighted (kg/m^3) explicit in time
!diagonal term in redi tensor (m^2/sec)
real, dimension(isd:ied,jsd:jed,nk,0:1,0:1) :: tensor_31 !tracer independent portion of mixing tensor
real, dimension(isd:ied,jsd:jed,nk,0:1,0:1) :: tensor_32 !tracer independent portion of mixing tensor
real, dimension(isd:ied,jsd:jed,nk) :: tx_trans_gm !for diagnosing i-transport due to GM (Sv)
real, dimension(isd:ied,jsd:jed,nk) :: ty_trans_gm !for diagnosing j-transport due to GM (Sv)
real, dimension(isd:ied,jsd:jed) :: depth_blayer_base ! depth (m) of boundary layer base for neutral physics
real, dimension(isd:ied,jsd:jed) :: slope_blayer_base ! abs(slope) at base of neutral boundary layer
real, dimension(isd:ied,jsd:jed) :: eddy_depth ! max of depth(m) mesoscale eddies penetrate & kpp bldepth
real, dimension(isd:ied,jsd:jed) :: turb_blayer_depth ! depth (m) of surface turbulent boundary layer
real, dimension(isd:ied,jsd:jed) :: N2_blayer_base ! squared buoyancy freq (1/s^2) at base of nblayer
integer, dimension(isd:ied,jsd:jed) :: ksurf_blayer ! k-value at base of surface nblayer
#else
real, dimension(:,:,:,:), allocatable :: delqc !density weighted (kg/m^3) quarter cell thickness(m)
real, dimension(:,:,:), allocatable :: dzwtr !(1/dzwt)(m^-1)
real, dimension(:,:,:), allocatable :: dtew !grid distances from T-point to cell faces (m)
real, dimension(:,:,:), allocatable :: dtns !grid distances from T-point to cell faces (m)
real, dimension(:,:,:), allocatable :: dtwedyt !horizontal areas (m^2) of quarter cell
real, dimension(:,:,:), allocatable :: dxtdtsn !horizontal areas (m^2) of quarter cell
real, dimension(:,:,:), allocatable :: slopex_drhodx !3D array of slopex * drhodx for diagnostics
real, dimension(:,:,:), allocatable :: slopey_drhody !3D array of slopey * drhody for diagnostics
real, dimension(:,:,:), allocatable :: grav_agm_dz_sx_drhodx !3D array of grav*agm*dz*slopex*drhodx for diagnostics
real, dimension(:,:,:), allocatable :: grav_agm_dz_sy_drhody !3D array of grav*agm*dz*slopey*drhody for diagnostics
real, dimension(:,:,:), allocatable :: aredi_array !3D array of redi diffusivities (m^2/sec)
real, dimension(:,:,:), allocatable :: agm_array !3D array of gm diffusivities (m^2/sec)
real, dimension(:,:), allocatable :: ah_array !2D array of micom horizontal diffusivities (m^2/sec)
real, dimension(:,:), allocatable :: bczone_radius !for bzcone calculation (m)
real, dimension(:,:), allocatable :: rossby_radius !first baroclinic Rossby radius (m)
real, dimension(:,:), allocatable :: rossby_radius_raw !first baroclinic Rossby radius (m) without min/max settings
real, dimension(:,:,:), allocatable :: eady_termx !rho_z*(S_x)^2 for computing Eady growth rate
real, dimension(:,:,:), allocatable :: eady_termy !rho_z*(S_y)^2 for computing Eady growth rate
real, dimension(:,:), allocatable :: baroclinic_termx !intermediate term for computing vert ave baroclinicity
real, dimension(:,:), allocatable :: baroclinic_termy !intermediate term for computing vert ave baroclinicity
real, dimension(:,:), allocatable :: grid_length !grid length scale (m)
real, dimension(:,:,:), allocatable :: drhodT !drho/dtheta (kg/m^3/C)
real, dimension(:,:,:), allocatable :: drhodS !drho/dsalinity (kg/m^3/psu)
real, dimension(:,:,:,:), allocatable :: drhodzb !vertical neutral density gradient (kg/m^3/m)
real, dimension(:,:,:,:), allocatable :: drhodzh !vertical neutral density gradient (kg/m^3/m)
real, dimension(:,:,:), allocatable :: K33_implicit !density weighted (kg/m^3) implicit in time
!diagonal term in redi tensor (m^2/sec)
real, dimension(:,:,:), allocatable :: K33_explicit !density weighted (kg/m^3) explicit in time
real, dimension(:,:,:,:,:), allocatable :: tensor_31 !tracer independent portion of mixing tensor
real, dimension(:,:,:,:,:), allocatable :: tensor_32 !tracer independent portion of mixing tensor
real, dimension(:,:,:), allocatable :: tx_trans_gm !for diagnosing i-transport due to GM (Sv)
real, dimension(:,:,:), allocatable :: ty_trans_gm !for diagnosing j-transport due to GM (Sv)
real, dimension(:,:), allocatable :: depth_blayer_base ! depth (m) of boundary layer base for neutral physics
real, dimension(:,:), allocatable :: slope_blayer_base ! abs(slope) at base of neutral boundary layer
real, dimension(:,:), allocatable :: eddy_depth ! max of depth (m) mesoscale eddies penetrate & kpp bldepth
real, dimension(:,:), allocatable :: turb_blayer_depth ! depth (m) of surface turbulent boundary layer
real, dimension(:,:), allocatable :: N2_blayer_base ! squared buoyancy freq (1/s^2) at base of nblayer
integer, dimension(:,:), allocatable :: ksurf_blayer ! k-value at base of surface nblayer
#endif
! introduce following derived types so that do not need to know num_prog_tracers at compile time
type(tracer_3d_1_nk_type), dimension(:), allocatable :: dTdx ! tracer partial derivative (tracer/m)
type(tracer_3d_1_nk_type), dimension(:), allocatable :: dTdy ! tracer partial derivative (tracer/m)
type(tracer_3d_0_nk_type), dimension(:), allocatable :: dTdz ! tracer partial derivative (tracer/m)
type(tracer_2d_type), dimension(:), allocatable :: fz1 ! z-flux component for tracers at particular k
type(tracer_2d_type), dimension(:), allocatable :: fz2 ! z-flux component for tracers at particular k
type(tracer_3d_1_nk_type), dimension(:), allocatable :: flux_x ! i-flux component for tracers for all k
type(tracer_3d_1_nk_type), dimension(:), allocatable :: flux_y ! j-flux component for tracers for all k
integer :: index_temp
integer :: index_salt
character(len=128) :: version=&
'$Id: ocean_nphysicsA.F90,v 1.1.2.20.22.5.18.1.18.1.32.1.4.1 2009/10/21 20:07:27 smg Exp $'
character (len=128) :: tagname = &
'$Name: mom4p1_pubrel_dec2009_nnz $'
logical :: module_is_initialized = .FALSE.
! time step settings
real :: dtime
real :: two_dtime_inv
! vertical coordinate
integer :: vert_coordinate_class
! lower and upper depth for vertically averaging ocean properties.
! read into ocean_nphysics_util_nml
real :: agm_closure_upper_depth
real :: agm_closure_lower_depth
!**************nml settings**************
! for using form drag rather than GM-skewsion
! gm_skew=1.0 if use_gm_skew=.true.
! gm_skew=0.0 if use_gm_skew=.false.
logical :: use_gm_skew=.true.
real :: gm_skew=1.0
! for setting the slope tapering methods
real :: smax ! set in ocean_neutral_util_nml
real :: swidth ! set in ocean_neutral_util_nml
logical :: dm_taper = .true. ! tanh tapering scheme of Danabasoglu and McWilliams
real :: swidthr ! inverse swidth
real :: smax_swidthr ! useful combination of terms
real :: dm_taper_const = 1.0 ! internally set to unity when dm_taper=.true.
logical :: gkw_taper = .false. ! quadratic tapering of Gerdes, Koberle, and Willebrand
real :: gkw_taper_const = 0.0 ! internally set to unity when gkw_taper=.true.
! neutral_physics_simple=.false. is commonly used at GFDL with mom4
! neutral_physics_simple=.true. was the setting more common in MOM3.
logical :: neutral_physics_simple=.false.
! for neutral blayer
real :: turb_blayer_min = 0.0 ! metres (was always set to zero in mom4p0)
! for diagnosing properties at the base of the "neutral physics boundary layer"
logical :: neutral_blayer_diagnose=.false.
! for maintaining horizontal GM velocity constant with depth within neutral boundary layer
logical :: neutral_linear_gm_taper=.true.
! for tapering neutral physics over penetration depth of
! eddies as determined by slope and Rossby radius
logical :: neutral_sine_taper=.true.
! to remove diagonal elements to neutral diffusion within the neutral boundar layer
logical :: neutral_taper_diagonal=.false.
real :: taper_diagonal=0.0 ! taper_diagonal=1 if neutral_taper_diagonal=.true. ; =0 otherwise
! to reduce neutral fluxes to horz/vert diffusion next to model boundaries
logical :: tmask_neutral_on=.false.
! to compute K33 explicitly in time.
! diffusion_all_explicit=.false. for realistic simulations.
logical :: diffusion_all_explicit=.false.
! revert to horizontal diffusion when tracer falls outside specified range
logical :: neutral_physics_limit=.true.
! for specifying transport units
! can either be Sv or mks
character(len=32) :: transport_units='Sv'
character(len=32) :: transport_dims ='Sv (10^9 kg/s)'
real :: transport_convert=1.0e-9
! for the module as a whole
logical :: use_this_module = .false.
logical :: debug_this_module = .false.
!**************end of nml settings**************
namelist /ocean_nphysicsA_nml/ use_this_module, debug_this_module, &
neutral_physics_limit, use_gm_skew, &
dm_taper, gkw_taper, tmask_neutral_on, &
neutral_physics_simple, neutral_taper_diagonal, &
neutral_linear_gm_taper, neutral_sine_taper, &
diffusion_all_explicit, neutral_blayer_diagnose, turb_blayer_min, &
transport_units
contains
!#######################################################################
!
!
!
! Initialize the neutral physics module by registering fields for
! diagnostic output and performing some numerical checks to see
! that namelist settings are appropriate.
!
!
subroutine ocean_nphysicsA_init(Grid, Domain, Time, Time_steps, Thickness, T_prog, &
ver_coordinate_class, agm_closure_lower_dept, agm_closure_upper_dept, &
smx, swidt, debug)
type(ocean_grid_type), intent(in), target :: Grid
type(ocean_domain_type), intent(in), target :: Domain
type(ocean_time_type), intent(in) :: Time
type(ocean_time_steps_type), intent(in) :: Time_steps
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
integer, intent(in) :: ver_coordinate_class
real, intent(in) :: agm_closure_lower_dept
real, intent(in) :: agm_closure_upper_dept
real, intent(in) :: smx
real, intent(in) :: swidt
logical, intent(in), optional :: debug
integer :: ioun, io_status, ierr
integer :: i, j, n
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if ( module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error from ocean_nphysicsA_mod (ocean_nphysicsA_init):already initialized')
endif
module_is_initialized = .TRUE.
call write_version_number( version, tagname )
num_prog_tracers = size(T_prog(:))
dtime = Time_steps%dtime_t
Dom => Domain
Grd => Grid
vert_coordinate_class = ver_coordinate_class
agm_closure_lower_depth = agm_closure_lower_dept
agm_closure_upper_depth = agm_closure_upper_dept
smax = smx
swidth = swidt
! provide for namelist over-ride of default values
ioun = open_namelist_file()
read (ioun,ocean_nphysicsA_nml,IOSTAT=io_status)
write (stdoutunit,'(/)')
write (stdoutunit,ocean_nphysicsA_nml)
write (stdlogunit,ocean_nphysicsA_nml)
ierr = check_nml_error(io_status,'ocean_nphysicsA_nml')
call close_file (ioun)
#ifndef MOM4_STATIC_ARRAYS
call get_local_indices(Domain,isd,ied,jsd,jed,isc,iec,jsc,jec)
nk = Grid%nk
#endif
if(use_this_module) then
call mpp_error(NOTE, &
'==> from ocean_nphysicsA_mod: USING ocean_nphysicsA_mod.')
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_nphysicsA_mod: USING ocean_nphysicsA.'
else
call mpp_error(NOTE, &
'==> from ocean_nphysicsA_mod: NOT using ocean_nphysicsA_mod.')
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_nphysicsA_mod: NOT using ocean_nphysicsA.'
return
endif
if (PRESENT(debug) .and. .not. debug_this_module) then
debug_this_module = debug
endif
if(debug_this_module) then
write(stdoutunit,'(a)') '==>Note: running ocean_nphysicsA_mod with debug_this_module=.true.'
endif
write(stdoutunit,'(/1x,a,f10.2)') &
'==> Note from ocean_nphysicsA_mod: using forward time step of (secs)', dtime
if(transport_units=='Sv') then
transport_convert=1.0e-9
transport_dims = 'Sv (10^9 kg/s)'
else
transport_convert=1.0
transport_dims = 'kg/s'
endif
! Useful constants
two_dtime_inv = 0.5/dtime !for explicit piece of K33
swidthr = 1.0/(swidth + epsln) !for slope taper function when dm_taper used
smax_swidthr = smax*swidthr !for slope taper function when dm_taper used
if(use_gm_skew) then
gm_skew=1.0
write(stdoutunit,'(1x,a)') &
' ==> Note from ocean_nphysicsA_init: use_gm_skew=.true. so will use GM-skewsion.'
else
gm_skew=0.0
write(stdoutunit,'(1x,a)') &
' ==> Note from ocean_nphysicsA_init: use_gm_skew=.false. so will NOT use GM-skewsion.'
if(neutral_physics_simple) then
call mpp_error(FATAL, &
'==>Error from ocean_nphysicsA_mod: use_gm_skew=.false. incompatible with neutral_physics_simple=.true.')
endif
endif
if(neutral_physics_simple) then
write(stdoutunit,'(1x,a)') &
' ==> Note from ocean_nphysics_closure_mod: neutral_physics_simple=.true.'
write(stdoutunit,'(7x,a)') &
'means the horizontal SGS tracer fluxes are downgradient for all neutral slopes.'
write(stdoutunit,'(7x,a)') &
'Assumes aredi_array = agm_array.'
write(stdoutunit,'(7x,a)') &
'Analogous to MOM3 approach.'
endif
if(neutral_physics_limit) then
write(stdoutunit,'(1x,a)') &
' ==>Note from ocean_nphysicsA_mod: neutral_physics_limit=.true.'
write(stdoutunit,'(7x,a)') &
'Will revert to horizontal diffusion for points where tracer is outside specified range.'
endif
if(dm_taper .and. .not. gkw_taper) then
write(stdoutunit,'(1x,a)') &
' ==>Note from ocean_nphysicsA_mod: dm_taper=.true. Will use the tanh scheme '
write(stdoutunit,'(7x,a)') &
'of Danabasoglu and McWilliams to taper neutral physics in steep sloped regions'
dm_taper_const =1.0
gkw_taper_const=0.0
endif
if(gkw_taper .and. .not. dm_taper) then
write(stdoutunit,'(1x,a)') &
' ==>Note from ocean_nphysicsA_mod: gkw_taper=.true. Will use the quadratic scheme'
write(stdoutunit,'(7x,a)') &
'of Gerdes, Koberle, and Willebrand to taper neutral physics in steep sloped regions'
dm_taper_const =0.0
gkw_taper_const=1.0
endif
if(gkw_taper .and. dm_taper) then
dm_taper_const =0.0
gkw_taper_const=0.0
call mpp_error(FATAL, &
'==>Error from ocean_nphysicsA_mod: gkw_taper and dm_taper cannot both be set true--choose only one.')
endif
if(.not. neutral_physics_simple) then
if(neutral_linear_gm_taper) then
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_nphysicsA_mod: neutral_linear_gm_taper=.true., so will linearly'
write(stdoutunit,'(7x,a)') &
'taper GM towards the surface when reaching steep neutral slopes in surface bdy.'
else
call mpp_error(NOTE, &
' ==> Note from ocean_nphysicsA_mod: GM exponentially tapered to zero in neutral bdy layer region.')
endif
if(neutral_sine_taper) then
write(stdoutunit,'(1x,a)') &
' ==> Note: Running with neutral_sine_taper=.true., and so will use sine-taper'
write(stdoutunit,'(7x,a)') &
'on fluxes where eddy penetration depth and/or KPP hblt exceeds grid depth.'
endif
endif
if(neutral_blayer_diagnose) then
write(stdoutunit,'(1x,a)') &
' ==> Note: Running with neutral_blayer_diagnose=.true., so will diagnose properties'
write(stdoutunit,'(7x,a)')'of the neutral physics boundary layer.'
endif
if(neutral_linear_gm_taper) then
write(stdoutunit,'(1x,a)')' ==>Running with a nontrivial GM transport in steep neutral slope regions.'
endif
if(diffusion_all_explicit) then
write(stdoutunit,'(1x,a)') &
' ==> Warning: Running w/ diffusion_all_explicit=.true., which means compute K33 contribution'
write(stdoutunit,'(7x,a)') &
'to neutral diffusion explicitly in time. This method is stable ONLY if taking'
write(stdoutunit,'(7x,a)') &
'very small time steps and/or running with just a single tracer.'
endif
if(neutral_taper_diagonal) then
taper_diagonal = 1.0
write(stdoutunit,'(1x,a,f5.2)') &
' ==>Note: Running w/ neutral_taper_diagonal=.true. and so taper_diagonal = ',taper_diagonal
write(stdoutunit,'(7x,a)') &
'Will taper all pieces of horiz neutral flux components in steep neutral slope regions.'
write(stdoutunit,'(7x,a)') &
'neutral_taper_diagonal=.false. will alternatively keep the diagonal pieces untapered.'
endif
if(.not. neutral_taper_diagonal) then
taper_diagonal = 0.0
write(stdoutunit,'(1x,a,f5.2)') &
' ==> Note: Running w/ neutral_taper_diagonal=.false. so taper_diagonal = ',taper_diagonal
write(stdoutunit,'(7x,a)') &
'Diagonal pieces of horizontal flux components are untapered regardless the neutral slope.'
endif
do n=1,num_prog_tracers
T_prog(n)%neutral_physics_limit = neutral_physics_limit
enddo
allocate( dTdx(num_prog_tracers) )
allocate( dTdy(num_prog_tracers) )
allocate( dTdz(num_prog_tracers) )
allocate( fz1(num_prog_tracers) )
allocate( fz2(num_prog_tracers) )
allocate( flux_x(num_prog_tracers) )
allocate( flux_y(num_prog_tracers) )
call set_ocean_domain(Dom_flux, Grid, xhalo=Dom%xhalo, yhalo = Dom%yhalo,name='flux dom neutral',maskmap=Dom%maskmap)
#ifndef MOM4_STATIC_ARRAYS
allocate (dtew(isd:ied,jsd:jed,0:1))
allocate (dtns(isd:ied,jsd:jed,0:1))
allocate (dtwedyt(isd:ied,jsd:jed,0:1))
allocate (dxtdtsn(isd:ied,jsd:jed,0:1))
allocate (grid_length(isd:ied,jsd:jed))
allocate (delqc(isd:ied,jsd:jed,nk,0:1))
allocate (dzwtr(isd:ied,jsd:jed,0:nk))
allocate (aredi_array(isd:ied,jsd:jed,nk))
allocate (agm_array(isd:ied,jsd:jed,nk))
allocate (ah_array(isd:ied,jsd:jed))
allocate (bczone_radius(isd:ied,jsd:jed))
allocate (rossby_radius(isd:ied,jsd:jed))
allocate (rossby_radius_raw(isd:ied,jsd:jed))
allocate (eady_termx(isd:ied,jsd:jed,nk))
allocate (eady_termy(isd:ied,jsd:jed,nk))
allocate (baroclinic_termx(isd:ied,jsd:jed))
allocate (baroclinic_termy(isd:ied,jsd:jed))
allocate (tx_trans_gm(isd:ied,jsd:jed,nk))
allocate (ty_trans_gm(isd:ied,jsd:jed,nk))
allocate (grav_agm_dz_sx_drhodx(isd:ied,jsd:jed,nk))
allocate (grav_agm_dz_sy_drhody(isd:ied,jsd:jed,nk))
allocate (slopex_drhodx(isd:ied,jsd:jed,nk))
allocate (slopey_drhody(isd:ied,jsd:jed,nk))
allocate (drhodT(isd:ied,jsd:jed,nk))
allocate (drhodS(isd:ied,jsd:jed,nk))
allocate (drhodzb(isd:ied,jsd:jed,nk,0:1))
allocate (drhodzh(isd:ied,jsd:jed,nk,0:1))
allocate (tensor_31(isd:ied,jsd:jed,nk,0:1,0:1))
allocate (tensor_32(isd:ied,jsd:jed,nk,0:1,0:1))
allocate (K33_implicit(isd:ied,jsd:jed,nk))
allocate (K33_explicit(isd:ied,jsd:jed,nk))
do n=1,num_prog_tracers
allocate ( dTdx(n)%field(isd:ied,jsd:jed,nk) )
allocate ( dTdy(n)%field(isd:ied,jsd:jed,nk) )
allocate ( dTdz(n)%field(isd:ied,jsd:jed,0:nk) )
allocate ( fz1(n)%field(isd:ied,jsd:jed) )
allocate ( fz2(n)%field(isd:ied,jsd:jed) )
allocate ( flux_x(n)%field(isd:ied,jsd:jed,nk) )
allocate ( flux_y(n)%field(isd:ied,jsd:jed,nk) )
enddo
allocate(depth_blayer_base(isd:ied,jsd:jed))
allocate(slope_blayer_base(isd:ied,jsd:jed))
allocate(eddy_depth(isd:ied,jsd:jed))
allocate(turb_blayer_depth(isd:ied,jsd:jed))
allocate(N2_blayer_base(isd:ied,jsd:jed))
allocate(ksurf_blayer(isd:ied,jsd:jed))
#endif
do n=1,num_prog_tracers
dTdx(n)%field(:,:,:) = 0.0
dTdy(n)%field(:,:,:) = 0.0
dTdz(n)%field(:,:,:) = 0.0
fz1(n)%field(:,:) = 0.0
fz2(n)%field(:,:) = 0.0
flux_x(n)%field(:,:,:) = 0.0
flux_y(n)%field(:,:,:) = 0.0
enddo
depth_blayer_base = 0.0
slope_blayer_base = 0.0
eddy_depth = 0.0
turb_blayer_depth = 0.0
N2_blayer_base = 0.0
ksurf_blayer = 1
tx_trans_gm = 0.0
ty_trans_gm = 0.0
dtew(:,:,0) = Grd%dtw(:,:)
dtew(:,:,1) = Grd%dte(:,:)
dtns(:,:,0) = Grd%dts(:,:)
dtns(:,:,1) = Grd%dtn(:,:)
dtwedyt(:,:,:) = 0.0
dtwedyt(:,:,0) = Grd%dte(:,:)*Grd%dyt(:,:)
do i=isc-1,iec
dtwedyt(i,:,1) = Grd%dtw(i+1,:)*Grd%dyt(i+1,:)
enddo
dxtdtsn(:,:,:) = 0.0
dxtdtsn(:,:,0) = Grd%dxt(:,:)*Grd%dtn(:,:)
do j=jsc-1,jec
dxtdtsn(:,j,1) = Grd%dxt(:,j+1)*Grd%dts(:,j+1)
enddo
grid_length(:,:) = 2.0*Grd%dxt(:,:)*Grd%dyt(:,:)/(Grd%dxt(:,:)+Grd%dyt(:,:))
eady_termx(:,:,:) = 0.0
eady_termy(:,:,:) = 0.0
baroclinic_termx(:,:) = 0.0
baroclinic_termy(:,:) = 0.0
grav_agm_dz_sx_drhodx(:,:,:) = 0.0
grav_agm_dz_sy_drhody(:,:,:) = 0.0
slopex_drhodx(:,:,:) = 0.0
slopey_drhody(:,:,:) = 0.0
rossby_radius(:,:) = 0.0
rossby_radius_raw(:,:) = 0.0
bczone_radius(:,:) = 0.0
agm_array(:,:,:) = 0.0
aredi_array(:,:,:) = 0.0
ah_array(:,:) = 0.0
call ocean_nphysics_coeff_init(Time, Thickness, rossby_radius, rossby_radius_raw, &
bczone_radius, agm_array, aredi_array, ah_array)
drhodT(:,:,:) = 0.0
drhodS(:,:,:) = 0.0
drhodzb(:,:,:,:) = 0.0
drhodzh(:,:,:,:) = 0.0
tensor_31(:,:,:,:,:) = 0.0
tensor_32(:,:,:,:,:) = 0.0
K33_implicit(:,:,:) = 0.0
K33_explicit(:,:,:) = 0.0
index_temp=-1;index_salt=-1
do n=1,num_prog_tracers
if (T_prog(n)%name == 'temp') index_temp = n
if (T_prog(n)%name == 'salt') index_salt = n
enddo
if (index_temp == -1 .or. index_salt == -1) then
call mpp_error(FATAL, &
'==>Error: temp and/or salt not identified in call to ocean_nphysicsA_init')
endif
! initialize clock ids
id_clock_neutral_blayer = mpp_clock_id('(Ocean neutral: blayer)' ,grain=CLOCK_ROUTINE)
id_clock_neutral_blayer_1 = mpp_clock_id('(Ocean neutral: blayer1)' ,grain=CLOCK_ROUTINE)
id_clock_neutral_blayer_2 = mpp_clock_id('(Ocean neutral: blayer2)' ,grain=CLOCK_ROUTINE)
id_clock_neutral_blayer_3 = mpp_clock_id('(Ocean neutral: blayer3)' ,grain=CLOCK_ROUTINE)
id_clock_fz_terms = mpp_clock_id('(Ocean neutral: fz-terms)',grain=CLOCK_ROUTINE)
id_clock_fx_flux = mpp_clock_id('(Ocean neutral: fx-flux)' ,grain=CLOCK_ROUTINE)
id_clock_fy_flux = mpp_clock_id('(Ocean neutral: fy-flux)' ,grain=CLOCK_ROUTINE)
id_clock_fz_flux = mpp_clock_id('(Ocean neutral: fz-flux)' ,grain=CLOCK_ROUTINE)
! register fields for diagnostic output
id_ah_bdy = -1
id_ah_bdy = register_static_field ('ocean_model', 'ah_bdy', &
Grd%tracer_axes(1:2), 'Horz diffusivity in surface neutral bdy', &
'm^2/sec', missing_value=missing_value, range=(/-10.0,1.e10/))
if (id_ah_bdy > 0) then
used = send_data(id_ah_bdy, ah_array(:,:), Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
id_k33_explicit = -1
id_k33_explicit = register_diag_field ('ocean_model', 'k33_explicit', &
Grd%tracer_axes_wt(1:3), Time%model_time, &
'K33_explicit tensor element', 'm^2/sec', &
missing_value=missing_value, range=(/-10.0,1.e20/))
id_tx_trans_gm = -1
id_tx_trans_gm = register_diag_field ('ocean_model', 'tx_trans_gm', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, &
'T-cell mass i-transport from GM',trim(transport_dims), &
missing_value=missing_value, range=(/-1e8,1.e8/))
id_ty_trans_gm = -1
id_ty_trans_gm = register_diag_field ('ocean_model', 'ty_trans_gm', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, &
'T-cell mass j-transport from GM',trim(transport_dims),&
missing_value=missing_value, range=(/-1e8,1.e8/))
id_depth_blayer_base = -1
id_depth_blayer_base = register_diag_field ('ocean_model','depth_blayer_base', &
Grd%tracer_axes(1:2), Time%model_time, &
'Depth of neutral physics surface bdy-layer', 'm', &
missing_value=missing_value, range=(/-1.0,1.e8/))
id_slope_blayer_base = -1
id_slope_blayer_base = register_diag_field ('ocean_model','slope_blayer_base', &
Grd%tracer_axes(1:2), Time%model_time, &
'abs(slope) at neutral physics bdy-layer base', 'dimensionless', &
missing_value=missing_value, range=(/-1.0,1.e8/))
id_eddy_depth = -1
id_eddy_depth = register_diag_field ('ocean_model','eddy_depth', &
Grd%tracer_axes(1:2), Time%model_time, &
'eddy penetration depth', 'm', &
missing_value=missing_value, range=(/-1.0,1.e8/))
id_steep_depth = -1
id_steep_depth = register_diag_field ('ocean_model','steep_depth', &
Grd%tracer_axes(1:2), Time%model_time, &
'Depth at base of steep slope surface bdy layer', 'm', &
missing_value=missing_value, range=(/-1.0,1.e8/))
id_grav_agm_dz_sx_drhodx = -1
id_grav_agm_dz_sx_drhodx = register_diag_field ('ocean_model', 'grav_agm_dz_sx_drhodx', &
Grd%tracer_axes(1:3), Time%model_time, &
'grav*dz*agm*neutral x-slope times drho_dx', 'W/m^2', &
missing_value=missing_value, range=(/-1.e15,1.e15/))
id_grav_agm_dz_sy_drhody = -1
id_grav_agm_dz_sy_drhody = register_diag_field ('ocean_model', 'grav_agm_dz_sy_drhody', &
Grd%tracer_axes(1:3), Time%model_time, &
'grav*dz*agm*neutral y-slope times drho_dy', 'W/m^2', &
missing_value=missing_value, range=(/-1.e15,1.e15/))
id_gm_eddy_ke_source = -1
id_gm_eddy_ke_source = register_diag_field ('ocean_model', 'gm_eddy_ke_source',&
Grd%tracer_axes(1:3), Time%model_time, &
'rho0*dz*agm*(slope*N)^2 = eddy ke source from GM', 'W/m^2', &
missing_value=missing_value, range=(/-1.e1,1.e15/), &
standard_name='tendency_of_ocean_eddy_kinetic_energy_content_due_to_bolus_transport')
id_slopex_drhodx = -1
id_slopex_drhodx = register_diag_field ('ocean_model', 'slopex_drhodx', &
Grd%tracer_axes(1:3), Time%model_time, &
'neutral x-slope times drho_dx', 'kg/m^4', &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_slopey_drhody = -1
id_slopey_drhody = register_diag_field ('ocean_model', 'slopey_drhody', &
Grd%tracer_axes(1:3), Time%model_time, &
'neutral y-slope times drho_dy', 'kg/m^4', &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_ustar = -1
id_ustar = register_diag_field ('ocean_model', 'ugm', &
Grd%vel_axes_uv(1:3), Time%model_time, &
'GM zonal velocity', 'm/sec', &
missing_value=missing_value, range=(/-10.0,10.0/))
id_vstar = -1
id_vstar = register_diag_field ('ocean_model', 'vgm', &
Grd%vel_axes_uv(1:3), Time%model_time, &
'GM merid velocity', 'm/sec', &
missing_value=missing_value, range=(/-10.0,10.0/))
id_wstar = -1
id_wstar = register_diag_field ('ocean_model', 'wgm', &
Grd%vel_axes_wt(1:3), Time%model_time, &
'GM vert velocity (T-cell bottom)', 'm/sec', &
missing_value=missing_value, range=(/-10.0,10.0/))
id_neutral_rho_nphysics = -1
id_neutral_rho_nphysics = register_diag_field ('ocean_model', 'neutral_rho_nphysics', &
Grd%tracer_axes(1:3), Time%model_time, &
'update of neutral density from explicit in time nphysics', &
'rho*rho_dz/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
id_wdian_rho_nphysics = -1
id_wdian_rho_nphysics = register_diag_field ('ocean_model', 'wdian_rho_nphysics', &
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral velocity component due to explicit in time nphysics',&
'm/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
allocate (id_k33_implicit(num_prog_tracers))
allocate (id_neutral_physics(num_prog_tracers))
id_k33_implicit = -1
id_neutral_physics = -1
do n=1,num_prog_tracers
id_k33_implicit(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_k33_implicit', &
Grd%tracer_axes_wt(1:3), Time%model_time, &
'K33_implicit tensor element for '//trim(T_prog(n)%name), &
'm^2/sec', missing_value=missing_value, range=(/-10.0,1.e10/))
if (T_prog(n)%name == 'temp') then
id_neutral_physics(n) = register_diag_field ('ocean_model', 'neutral_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), Time%model_time, &
'rho*dzt*cp*explicit neutral tendency (heating)', &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e10,1.e10/))
else
id_neutral_physics(n) = register_diag_field ('ocean_model', 'neutral_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), Time%model_time, &
'rho*dzt*explicit neutral tendency for '//trim(T_prog(n)%name), &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e10,1.e10/))
endif
enddo
allocate (id_flux_x(num_prog_tracers))
allocate (id_flux_y(num_prog_tracers))
allocate (id_flux_x_int_z(num_prog_tracers))
allocate (id_flux_y_int_z(num_prog_tracers))
id_flux_x = -1
id_flux_y = -1
id_flux_x_int_z = -1
id_flux_y_int_z = -1
do n=1,num_prog_tracers
if(n == index_temp) then
id_flux_x(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_xflux_neutral', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, 'cp*neutral_xflux*dyt*rho_dzt*temp', &
'Watt', missing_value=missing_value, range=(/-1.e18,1.e18/))
id_flux_y(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_yflux_neutral', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, 'cp*neutral_yflux*dxt*rho_dzt*temp', &
'Watt', missing_value=missing_value, range=(/-1.e18,1.e18/))
id_flux_x_int_z(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_xflux_neutral_int_z', &
Grd%tracer_axes_flux_x(1:2), Time%model_time, 'z-integral cp*neutral_xflux*dyt*rho_dzt*temp', &
'Watt', missing_value=missing_value, range=(/-1.e18,1.e18/))
id_flux_y_int_z(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_yflux_neutral_int_z', &
Grd%tracer_axes_flux_y(1:2), Time%model_time, 'z-integral cp*neutral_yflux*dxt*rho_dzt*temp', &
'Watt', missing_value=missing_value, range=(/-1.e18,1.e18/))
else
id_flux_x(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_xflux_neutral', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, &
'neutral_xflux*dyt*rho_dzt*tracer for'//trim(T_prog(n)%name), &
trim(T_prog(n)%units)//' kg/sec', missing_value=missing_value, range=(/-1.e18,1.e18/))
id_flux_y(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_yflux_neutral', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, &
'neutral_yflux*dxt*rho_dzt*tracer for'//trim(T_prog(n)%name), &
trim(T_prog(n)%units)//' kg/sec', missing_value=missing_value, range=(/-1.e18,1.e18/))
id_flux_x_int_z(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_xflux_neutral_int_z', &
Grd%tracer_axes_flux_x(1:2), Time%model_time, &
'z-integral neutral_xflux*dyt*rho_dzt*tracer for'//trim(T_prog(n)%name), &
trim(T_prog(n)%units)//' kg/sec', missing_value=missing_value, range=(/-1.e18,1.e18/))
id_flux_y_int_z(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_yflux_neutral_int_z', &
Grd%tracer_axes_flux_y(1:2), Time%model_time, &
'z-integral neutral_yflux*dxt*rho_dzt*tracer for'//trim(T_prog(n)%name), &
trim(T_prog(n)%units)//' kg/sec', missing_value=missing_value, range=(/-1.e18,1.e18/))
endif
enddo
end subroutine ocean_nphysicsA_init
! NAME="ocean_nphysicsA_init"
!#######################################################################
!
!
!
! This function computes the thickness weighted and density weighted
! time tendency for tracer from neutral physics. Full discussion
! and details are provided by Griffies (2004).
!
! Here is a brief summary.
!
!---How the neutral diffusive flux components are computed:
!
! The vertical flux component is split into diagonal (3,3) and
! off-diagonal (3,1) and (3,2) terms. The off-diagonal (3,1) and (3,2)
! terms are included explicitly in time. The main contribution from the
! (3,3) term to the time tendency is included implicitly in time
! along with the usual contribution from diapycnal processes
! (vertical mixing schemes). This is the K33_implicit term.
! This approach is necessary with high vertical resolution, as
! noted by Cox (1987). However, splitting the vertical flux into
! an implicit and explicit piece compromises the
! integrity of the vertical flux component (see Griffies et al. 1998).
! So to minimize the disparity engendered by this split, the portion of
! K33 that can be stably included explicitly in time is computed along
! with the (3,1) and (3,2) terms.
!
! All other terms in the mixing tensor are included explicitly in time
! using a forward time step as required for temporal stability of
! numerical diffusive processes.
!
! The off-diagonal terms in the horizontal flux components, and all terms
! in the vertical flux component, are tapered in regions of steep neutral
! slope according to the requirements of linear stability. MOM4 allows for
! choice of two tapering schemes:
! (a) the tanh taper of Danabasoglu and McWilliams (1995)
! (b) the quadratic scheme of Gerdes, Koberle, and Willebrand (1991)
! Linear stability is far less stringent on the diagonal (1,1) and (2,2)
! part of the horizontal flux. Indeed, these terms in practice need
! not be tapered in steep sloped regions. The namelist
! neutral_taper_diagonal=.false. keeps the diagnonal terms maintained
! for all neutral slopes. This approach assists in reducing numerical
! noise in regions where the physical system experiences a lot of
! diapycnal mixing anyhow.
!
!---How the skew diffusive flux components are computed:
!
! The GM skew flux components are purely off-diagonal.
! They are generally tapered when neutral slope
! is large (neutral_physics_simple=.false).
! Doing so maintains a nontrivial GM slumping effect even when the
! neutral slopes are vertical. The alternative neutral_physics_simple=.true.
! is the approach used in MOM3, whereby GM effects are removed
! in steep sloped regions. neutral_physics_simple=.false. is
! less efficient, but has been seen to yield superior simulations.
!
!
!
subroutine nphysicsA (Time, Thickness, Dens, rho, T_prog, &
surf_blthick, gm_diffusivity, &
baroclinic_rossby_radius)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: rho
real, dimension(isd:,jsd:), intent(in) :: surf_blthick
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
real, dimension(isd:,jsd:,:), intent(inout) :: gm_diffusivity
real, dimension(isd:,jsd:), intent(inout) :: baroclinic_rossby_radius
real, dimension(isd:ied,jsd:jed) :: tchg
real, dimension(isd:ied,jsd:jed) :: tmp_flux
integer :: i, j, k, n, nn
integer :: tau, taum1
real :: temporary
if (.not. use_this_module) return
if ( .not. module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error from ocean_nphysicsA (neutral_physics): needs initialization')
endif
if (size(T_prog(:)) /= num_prog_tracers) then
call mpp_error(FATAL, &
'==>Error from ocean_nphysicsA (neutral_physics): inconsistent size of T_prog')
endif
taum1 = Time%taum1
tau = Time%tau
! time dependent delqc geometric factor
do k=1,nk
delqc(:,:,k,0) = Grd%fracdz(k,0)*Thickness%rho_dzt(:,:,k,tau)
delqc(:,:,k,1) = Grd%fracdz(k,1)*Thickness%rho_dzt(:,:,k,tau)
enddo
! time dependent inverse dzwt
do k=0,nk
dzwtr(:,:,k) = 1.0/Thickness%dzwt(:,:,k)
enddo
! compute density derivatives
do k=1,nk
do j=jsd,jed
do i=isd,ied
drhodT(i,j,k) = Dens%drhodT(i,j,k)
drhodS(i,j,k) = Dens%drhodS(i,j,k)
enddo
enddo
enddo
call tracer_derivs(Time, taum1, dtime, drhodT, drhodS, T_prog, dzwtr, &
dTdx, dTdy, dTdz, drhodzb, drhodzh)
call neutral_slopes(Time, dTdx, dTdy, drhodT, drhodS, drhodzb, tensor_31, tensor_32)
call neutral_blayer(Time, surf_blthick)
call fz_terms(Time, Thickness, Dens, T_prog, rho)
do n=1,num_prog_tracers
fz1(n)%field(:,:) = 0.0
fz2(n)%field(:,:) = 0.0
flux_x(n)%field(:,:,:) = 0.0
flux_y(n)%field(:,:,:) = 0.0
T_prog(n)%wrk1(:,:,:) = 0.0
enddo
if(neutral_physics_simple) then
! horizontal flux components are downgradient diffusive
do k=1,nk
do nn=1,num_prog_tracers
flux_x(nn)%field(:,:,k) = &
agm_array(:,:,k)*dTdx(nn)%field(:,:,k)*FMX(Thickness%rho_dzt(:,:,k,tau)*Grd%tmask(:,:,k))
flux_y(nn)%field(:,:,k) = &
agm_array(:,:,k)*dTdy(nn)%field(:,:,k)*FMY(Thickness%rho_dzt(:,:,k,tau)*Grd%tmask(:,:,k))
enddo
enddo
if (Grd%tripolar) then
do nn=1,num_prog_tracers
call mpp_update_domains(flux_x(nn)%field(:,:,:), flux_y(nn)%field(:,:,:), Dom_flux%domain2d, &
gridtype=CGRID_NE, complete=T_prog(nn)%complete)
enddo
endif
do k=1,nk
call mpp_clock_begin(id_clock_fz_flux)
call fz_flux(T_prog,k)
call mpp_clock_end(id_clock_fz_flux)
do nn=1,num_prog_tracers
tchg(:,:) = (BDX_ET(flux_x(nn)%field(:,:,k)) + BDY_NT(flux_y(nn)%field(:,:,k)))
do j=jsc,jec
do i=isc,iec
T_prog(nn)%wrk1(i,j,k) = Grd%tmask(i,j,k)*(tchg(i,j) + (fz1(nn)%field(i,j)-fz2(nn)%field(i,j)))
T_prog(nn)%th_tendency(i,j,k) = T_prog(nn)%th_tendency(i,j,k) + T_prog(nn)%wrk1(i,j,k)
fz1(nn)%field(i,j) = fz2(nn)%field(i,j)
flux_x(nn)%field(i,j,k) = Grd%dyte(i,j)*flux_x(nn)%field(i,j,k)
flux_y(nn)%field(i,j,k) = Grd%dxtn(i,j)*flux_y(nn)%field(i,j,k)
enddo
enddo
enddo
enddo
else
do k=1,nk
call fx_flux(Time, Thickness, T_prog, k)
call fy_flux(Time, Thickness, T_prog, k)
enddo
if (Grd%tripolar) then
do nn=1,num_prog_tracers
call mpp_update_domains(flux_x(nn)%field(:,:,:), flux_y(nn)%field(:,:,:), Dom_flux%domain2d, &
gridtype=CGRID_NE, complete=T_prog(nn)%complete)
enddo
endif
do k=1,nk
call mpp_clock_begin(id_clock_fz_flux)
call fz_flux(T_prog,k)
call mpp_clock_end(id_clock_fz_flux)
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc,iec
T_prog(nn)%wrk1(i,j,k) = &
Grd%tmask(i,j,k) &
*(fz1(nn)%field(i,j)-fz2(nn)%field(i,j) &
+( flux_x(nn)%field(i,j,k)-flux_x(nn)%field(i-1,j,k) &
+flux_y(nn)%field(i,j,k)-flux_y(nn)%field(i,j-1,k) )*Grd%datr(i,j) &
)
T_prog(nn)%th_tendency(i,j,k) = T_prog(nn)%th_tendency(i,j,k) + T_prog(nn)%wrk1(i,j,k)
fz1(nn)%field(i,j) = fz2(nn)%field(i,j)
enddo
enddo
enddo
enddo
endif ! endif for neutral_physics_simple
! compute Eady growth rate and baroclinicity for use in next time step
call compute_eady_rate(Time, Thickness, T_prog, Dens, eady_termx, eady_termy)
call compute_baroclinicity(Time%model_time, baroclinic_termx, baroclinic_termy)
! compute rossby radius for use in next time step
call compute_rossby_radius(Thickness, dTdz, Time%model_time, drhodT, drhodS, &
rossby_radius, rossby_radius_raw)
baroclinic_rossby_radius(:,:) = rossby_radius_raw(:,:)
! compute width of baroclinic zone for use in next time step
call compute_bczone_radius(Time%model_time, bczone_radius)
! update closure-based diffusivity for next time step
call compute_diffusivity(Time%model_time, ksurf_blayer, Dens%drhodz_zt, rossby_radius, &
bczone_radius, agm_array, aredi_array)
do nn=1,num_prog_tracers
if(id_neutral_physics(nn) > 0) then
used = send_data (id_neutral_physics(nn), T_prog(nn)%wrk1(:,:,:)*T_prog(nn)%conversion, &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
! send fluxes to diag_manager
! minus sign is due to a MOM-convention for physics fluxes
if(id_flux_x(nn) > 0) then
used = send_data (id_flux_x(nn), -1.0*T_prog(nn)%conversion*flux_x(nn)%field(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_flux_y(nn) > 0) then
used = send_data (id_flux_y(nn), -1.0*T_prog(nn)%conversion*flux_y(nn)%field(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_flux_x_int_z(nn) > 0) then
tmp_flux = 0.0
do k=1,nk
tmp_flux(isc:iec,jsc:jec) = tmp_flux(isc:iec,jsc:jec) + flux_x(nn)%field(isc:iec,jsc:jec,k)
enddo
used = send_data (id_flux_x_int_z(nn), -1.0*T_prog(nn)%conversion*tmp_flux(:,:), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if(id_flux_y_int_z(nn) > 0) then
tmp_flux = 0.0
do k=1,nk
tmp_flux(isc:iec,jsc:jec) = tmp_flux(isc:iec,jsc:jec) + flux_y(nn)%field(isc:iec,jsc:jec,k)
enddo
used = send_data (id_flux_y_int_z(nn), -1.0*T_prog(nn)%conversion*tmp_flux(:,:), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
enddo
! more diagnostic output
if(id_neutral_rho_nphysics > 0) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)* &
(drhodT(i,j,k)*T_prog(index_temp)%wrk1(i,j,k) &
+drhodS(i,j,k)*T_prog(index_salt)%wrk1(i,j,k))
enddo
enddo
enddo
used = send_data (id_neutral_rho_nphysics, wrk1(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_wdian_rho_nphysics > 0) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
temporary = Dens%drhodz_zt(i,j,k)*Thickness%rho_dzt(i,j,k,tau)
wrk1(i,j,k) = Grd%tmask(i,j,k)* &
(drhodT(i,j,k)*T_prog(index_temp)%wrk1(i,j,k) &
+drhodS(i,j,k)*T_prog(index_salt)%wrk1(i,j,k)) &
/(epsln+temporary)
enddo
enddo
enddo
used = send_data (id_wdian_rho_nphysics, wrk1(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if (id_k33_explicit > 0) used = send_data (id_k33_explicit, rho0r*K33_explicit(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
! gm_diffusivity passed to neutral_physics for use in computing form drag
gm_diffusivity(:,:,:) = 0.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
gm_diffusivity(i,j,k) = agm_array(i,j,k)
enddo
enddo
enddo
call gm_velocity(Thickness, Time)
end subroutine nphysicsA
! NAME="nphysicsA"
!#######################################################################
!
!
!
! Subroutine computes the boundary layer as determined by
! 1. steep neutral slopes
! 2. depth within which typical mesoscale eddies are partially outcropped
! 3. depth within which vertical mixing scheme (e.g., kpp) computes a boundary layer
!
! Note: Only consider surface boundary layers here.
!
! Scheme originally coded for mom4p0 by Stephen.Griffies@noaa.gov
! with help for optimization by Russell.Fiedler@csiro.au.
!
!
!
subroutine neutral_blayer(Time, surf_blthick)
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:,jsd:), intent(in) :: surf_blthick
logical :: within_interior
logical :: slopeok(isc:iec)
integer :: i, j, k, ip, jq, kb, kr, kk, kkpkr
real :: depth, slope, absslope
real :: steep_depth(isd:ied,jsd:jed)
real :: slope_blayer_baseA(isd:ied,jsd:jed)
real :: slope_blayer_baseB(isd:ied,jsd:jed)
real :: thick_31(isd:ied,jsd:jed,0:1,0:1)
real :: thick_32(isd:ied,jsd:jed,0:1,0:1)
real :: thick_13(isd:ied,jsd:jed,0:1,0:1)
real :: thick_23(isd:ied,jsd:jed,0:1,0:1)
real :: slope_31(isd:ied,jsd:jed,0:1,0:1)
real :: slope_32(isd:ied,jsd:jed,0:1,0:1)
real :: slope_13(isd:ied,jsd:jed,0:1,0:1)
real :: slope_23(isd:ied,jsd:jed,0:1,0:1)
call mpp_clock_begin(id_clock_neutral_blayer)
thick_31 = Grd%zw(1)
thick_32 = Grd%zw(1)
thick_13 = Grd%zw(1)
thick_23 = Grd%zw(1)
slope_31 = 0.0
slope_32 = 0.0
slope_13 = 0.0
slope_23 = 0.0
slope_blayer_baseA = 0.0
slope_blayer_baseB = 0.0
steep_depth = 0.0
eddy_depth = 0.0
slope_blayer_base = 0.0
depth_blayer_base = 0.0
if(neutral_sine_taper .or. neutral_blayer_diagnose) then
! Determine depth over which mesoscale eddies feel the ocean
! surface. This depth is a function of the neutral slope
! and the Rossby radius. This depth is called "eddy_depth".
! The algorithm for computing this depth is taken from
! the appendix to Large etal, 1997 JPO vol 27, 2418-2447.
!
! In addition to considering mesoscale eddy lengths,
! include the possibility that the diabatic vertical
! mixing (e.g., kpp) produces a mixed layer depth that is
! deeper than the depth that mesoscale eddies feel the ocean
! surface. Include this surf_blthick in the considerations so
! to determine the depth of this generalized "boundary layer"
! and the neutral slope at the base of the boundary layer.
call mpp_clock_begin(id_clock_neutral_blayer_1)
! 31-triads
do ip=0,1
do kr=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop31a: do kk=1,nk-1
do i=isc,iec
absslope = abs(tensor_31(i,j,kk,ip,kr))
depth = max(rossby_radius(i,j)*absslope, max(turb_blayer_min,surf_blthick(i,j)))
if(depth > Grd%zw(kk) .and. .not. slopeok(i)) then
thick_31(i,j,ip,kr) = Grd%zw(kk)
slope_31(i,j,ip,kr) = absslope
else
slopeok(i)=.true.
endif
enddo
if(all(slopeok(:))) exit kkloop31a
enddo kkloop31a
enddo
enddo
enddo
! 13-triads
do kr=0,1
do ip=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop13a: do kk=1,nk
kkpkr = min(kk+kr,nk)
do i=isc,iec
slope = -Grd%tmask(i+ip,j,kkpkr)&
*(drhodT(i+ip,j,kk)*dTdx(index_temp)%field(i,j,kk) +&
drhodS(i+ip,j,kk)*dTdx(index_salt)%field(i,j,kk)) &
/(drhodT(i+ip,j,kk)*dTdz(index_temp)%field(i+ip,j,kk-1+kr) +&
drhodS(i+ip,j,kk)*dTdz(index_salt)%field(i+ip,j,kk-1+kr) - epsln)
absslope = abs(slope)
depth = max(rossby_radius(i,j)*absslope, max(turb_blayer_min,surf_blthick(i,j)))
if(depth > Grd%zw(kk) .and. .not.slopeok(i)) then
thick_13(i,j,ip,kr) = Grd%zw(kk)
slope_13(i,j,ip,kr) = absslope
else
slopeok(i)=.true.
endif
enddo
if(all(slopeok(:))) exit kkloop13a
enddo kkloop13a
enddo
enddo
enddo
! 32-triads
do jq=0,1
do kr=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop32a: do kk=1,nk-1
do i=isc,iec
absslope = abs(tensor_32(i,j,kk,jq,kr))
depth = max(rossby_radius(i,j)*absslope, max(turb_blayer_min,surf_blthick(i,j)))
if(depth > Grd%zw(kk) .and. .not. slopeok(i) ) then
thick_32(i,j,jq,kr) = Grd%zw(kk)
slope_32(i,j,jq,kr) = absslope
else
slopeok(i) = .true.
endif
enddo
if(all(slopeok(:))) exit kkloop32a
enddo kkloop32a
enddo
enddo
enddo
! 23-triads
do kr=0,1
do jq=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop23a: do kk=1,nk
do i=isc,iec
kkpkr = min(kk+kr,nk)
slope = -Grd%tmask(i,j+jq,kkpkr)&
*(drhodT(i,j+jq,kk)*dTdy(index_temp)%field(i,j,kk)+ &
drhodS(i,j+jq,kk)*dTdy(index_salt)%field(i,j,kk)) &
/(drhodT(i,j+jq,kk)*dTdz(index_temp)%field(i,j+jq,kk-1+kr) +&
drhodS(i,j+jq,kk)*dTdz(index_salt)%field(i,j+jq,kk-1+kr) - epsln)
absslope = abs(slope)
depth = max(rossby_radius(i,j)*absslope, max(turb_blayer_min,surf_blthick(i,j)))
if(depth > Grd%zw(kk) .and. .not.slopeok(i)) then
thick_23(i,j,jq,kr) = Grd%zw(kk)
slope_23(i,j,jq,kr) = absslope
else
slopeok(i)=.true.
endif
enddo
if(all(slopeok(:))) exit kkloop23a
enddo kkloop23a
enddo
enddo
enddo
call mpp_clock_end(id_clock_neutral_blayer_1)
call mpp_clock_begin(id_clock_neutral_blayer_2)
! max of triad depth defines eddy_depth. average slope
! defines slope at boundary layer base. This method
! maximizes eddy_depth, and regulates slope at base.
do j=jsc,jec
do i=isc,iec
eddy_depth(i,j) = max(Grd%zw(1), &
thick_31(i,j,0,0), thick_31(i,j,0,1), &
thick_31(i,j,1,0), thick_31(i,j,1,1), &
thick_32(i,j,0,0), thick_32(i,j,0,1), &
thick_32(i,j,1,0), thick_32(i,j,1,1), &
thick_13(i,j,0,0), thick_13(i,j,0,1), &
thick_13(i,j,1,0), thick_13(i,j,1,1), &
thick_23(i,j,0,0), thick_23(i,j,0,1), &
thick_23(i,j,1,0), thick_23(i,j,1,1))
! make sure that none of the slopes are larger than smax
do kr=0,1
do ip=0,1
slope_31(i,j,ip,kr) = min(smax, slope_31(i,j,ip,kr))
slope_13(i,j,ip,kr) = min(smax, slope_13(i,j,ip,kr))
slope_32(i,j,ip,kr) = min(smax, slope_32(i,j,ip,kr))
slope_23(i,j,ip,kr) = min(smax, slope_23(i,j,ip,kr))
enddo
enddo
! average of the 16 slopes defines slope at boundary layer base
slope_blayer_baseA(i,j) = &
(slope_31(i,j,0,0) + slope_13(i,j,0,0) + &
slope_31(i,j,1,0) + slope_13(i,j,1,0) + &
slope_31(i,j,0,1) + slope_13(i,j,0,1) + &
slope_31(i,j,1,1) + slope_13(i,j,1,1) + &
slope_32(i,j,0,0) + slope_23(i,j,0,0) + &
slope_32(i,j,0,1) + slope_23(i,j,0,1) + &
slope_32(i,j,1,0) + slope_23(i,j,1,0) + &
slope_32(i,j,1,1) + slope_23(i,j,1,1)) /16.0
enddo
enddo
call mpp_clock_end(id_clock_neutral_blayer_2)
endif ! endif for neutral_sine_taper .or. neutral_blayer_diagnose
if(neutral_linear_gm_taper .or. neutral_blayer_diagnose) then
! Determine depth of surface boundary layer as defined
! by depth where neutral slopes are larger than smax.
thick_31 = Grd%zw(1)
thick_32 = Grd%zw(1)
thick_13 = Grd%zw(1)
thick_23 = Grd%zw(1)
! 31-triads
do ip=0,1
do kr=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop31b: do kk=1,nk-1
do i=isc,iec
absslope = abs(tensor_31(i,j,kk,ip,kr))
if(absslope > smax .and. .not. slopeok(i)) then
thick_31(i,j,ip,kr) = Grd%zw(kk)
slope_31(i,j,ip,kr) = absslope
else
slopeok(i)=.true.
endif
enddo
if(all(slopeok(:))) exit kkloop31b
enddo kkloop31b
enddo
enddo
enddo
! 13-triads
do kr=0,1
do ip=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop13b: do kk=1,nk
kkpkr = min(kk+kr,nk)
do i=isc,iec
slope = -Grd%tmask(i+ip,j,kkpkr)&
*(drhodT(i+ip,j,kk)*dTdx(index_temp)%field(i,j,kk) +&
drhodS(i+ip,j,kk)*dTdx(index_salt)%field(i,j,kk)) &
/(drhodT(i+ip,j,kk)*dTdz(index_temp)%field(i+ip,j,kk-1+kr) +&
drhodS(i+ip,j,kk)*dTdz(index_salt)%field(i+ip,j,kk-1+kr) - epsln)
absslope = abs(slope)
if(absslope > smax .and. .not.slopeok(i)) then
thick_13(i,j,ip,kr) = Grd%zw(kk)
slope_13(i,j,ip,kr) = absslope
else
slopeok(i)=.true.
endif
enddo
if(all(slopeok(:))) exit kkloop13b
enddo kkloop13b
enddo
enddo
enddo
! 32-triads
do jq=0,1
do kr=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop32b: do kk=1,nk-1
do i=isc,iec
absslope = abs(tensor_32(i,j,kk,jq,kr))
if(absslope > smax .and. .not. slopeok(i)) then
thick_32(i,j,jq,kr) = Grd%zw(kk)
slope_32(i,j,jq,kr) = absslope
else
slopeok(i)=.true.
endif
enddo
if(all(slopeok(:))) exit kkloop32b
enddo kkloop32b
enddo
enddo
enddo
! 23-triads
do kr=0,1
do jq=0,1
do j=jsc,jec
slopeok(:)=.false.
kkloop23b: do kk=1,nk
do i=isc,iec
kkpkr = min(kk+kr,nk)
slope = -Grd%tmask(i,j+jq,kkpkr)&
*(drhodT(i,j+jq,kk)*dTdy(index_temp)%field(i,j,kk)+ &
drhodS(i,j+jq,kk)*dTdy(index_salt)%field(i,j,kk)) &
/(drhodT(i,j+jq,kk)*dTdz(index_temp)%field(i,j+jq,kk-1+kr) +&
drhodS(i,j+jq,kk)*dTdz(index_salt)%field(i,j+jq,kk-1+kr) - epsln)
absslope = abs(slope)
if(absslope > smax .and. .not.slopeok(i)) then
thick_23(i,j,jq,kr) = Grd%zw(kk)
slope_23(i,j,jq,kr) = absslope
else
slopeok(i)=.true.
endif
enddo
if(all(slopeok(:))) exit kkloop23b
enddo kkloop23b
enddo
enddo
enddo
call mpp_clock_begin(id_clock_neutral_blayer_3)
! max of triad depth defines steep_depth. average slope
! defines slope at boundary layer base. This method
! maximizes steep_depth, and regulates slope at base.
do j=jsc,jec
do i=isc,iec
steep_depth(i,j) = max(Grd%zw(1), &
thick_31(i,j,0,0), thick_31(i,j,0,1), &
thick_31(i,j,1,0), thick_31(i,j,1,1), &
thick_32(i,j,0,0), thick_32(i,j,0,1), &
thick_32(i,j,1,0), thick_32(i,j,1,1), &
thick_13(i,j,0,0), thick_13(i,j,0,1), &
thick_13(i,j,1,0), thick_13(i,j,1,1), &
thick_23(i,j,0,0), thick_23(i,j,0,1), &
thick_23(i,j,1,0), thick_23(i,j,1,1))
! make sure that none of the slopes are larger than smax
do kr=0,1
do ip=0,1
slope_31(i,j,ip,kr) = min(smax, slope_31(i,j,ip,kr))
slope_13(i,j,ip,kr) = min(smax, slope_13(i,j,ip,kr))
slope_32(i,j,ip,kr) = min(smax, slope_32(i,j,ip,kr))
slope_23(i,j,ip,kr) = min(smax, slope_23(i,j,ip,kr))
enddo
enddo
! average of the 16 slopes defines slope at boundary layer base
slope_blayer_baseB(i,j) = 0.0
do kr=0,1
do ip=0,1
slope_blayer_baseB(i,j) = slope_blayer_baseB(i,j) + &
slope_31(i,j,ip,kr) + slope_13(i,j,ip,kr) + &
slope_32(i,j,ip,kr) + slope_23(i,j,ip,kr)
enddo
enddo
slope_blayer_baseB(i,j) = slope_blayer_baseB(i,j)/16.0
enddo
enddo
call mpp_clock_end(id_clock_neutral_blayer_3)
! The maximum of steep_depth and eddy_depth define
! neutral_blayer_depth. For depths shallower than
! neutral_blayer_depth, horizontal gm velocity is depth
! independent and vertical gm velocity linearly decreases
! to zero as the surface is reached.
do j=jsc,jec
do i=isc,iec
if(steep_depth(i,j) > eddy_depth(i,j)) then
slope_blayer_base(i,j) = slope_blayer_baseB(i,j)
depth_blayer_base(i,j) = steep_depth(i,j)
else
slope_blayer_base(i,j) = slope_blayer_baseA(i,j)
depth_blayer_base(i,j) = eddy_depth(i,j)
endif
enddo
enddo
! save the k-level surface transition layer.
! initialize ksurf_blayer to 1 rather than 0, to avoid
! accessing the zeroth element of an array.
ksurf_blayer(:,:) = 1
wrk1_2d(:,:) = 0.0
do j=jsd,jed
do i=isd,ied
kb=Grd%kmt(i,j)
within_interior=.false.
ksurf_blayer_loop: do k=1,kb
if(Grd%zw(k) > depth_blayer_base(i,j)) then
ksurf_blayer(i,j) = k
wrk1_2d(i,j) = float(k)
exit ksurf_blayer_loop
endif
enddo ksurf_blayer_loop
enddo
enddo
endif ! endif for neutral_linear_gm_taper .or. neutral_blayer_diagnose
if (id_depth_blayer_base > 0) then
used = send_data (id_depth_blayer_base, depth_blayer_base(:,:), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_slope_blayer_base > 0) then
used = send_data (id_slope_blayer_base, slope_blayer_base(:,:), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_eddy_depth > 0) then
used = send_data (id_eddy_depth, eddy_depth(:,:), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_steep_depth > 0) then
used = send_data (id_steep_depth, steep_depth(:,:), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
! update domains or zero the depths, depending on the chosen method
if(neutral_linear_gm_taper) then
call mpp_update_domains(depth_blayer_base, Dom%domain2d)
call mpp_update_domains(slope_blayer_base, Dom%domain2d)
else
depth_blayer_base = 0.0
endif
if(neutral_sine_taper) then
call mpp_update_domains(eddy_depth, Dom%domain2d)
else
eddy_depth = 0.0
endif
call mpp_clock_end(id_clock_neutral_blayer)
end subroutine neutral_blayer
! NAME="neutral_blayer"
!#######################################################################
!
!
!
! Subroutine computes the tracer independent pieces of the vertical
! flux component. As a result of this routine,
! Array tensor_31 = x-diffusivity*slope (m^2/sec) for fz
! Array tensor_32 = y-diffusivity*slope (m^2/sec) for fz
!
! K33 is the (3,3) term in small angle Redi diffusion tensor.
! It is broken into an explicit in time piece and implicit
! in time piece. It is weighted by density for non-Boussinesq
! and rho0 for Boussinesq.
!
! K33 has units (kg/m^3)*m^2/sec.
!
! Also will compute the squared Eady growth rate, with the maximum
! slope contributing to this growth rate set by smax.
!
!
subroutine fz_terms(Time, Thickness, Dens, T_prog, rho)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(inout) :: T_prog(num_prog_tracers)
real, dimension(isd:,jsd:,:), intent(in) :: rho
integer :: i, j, k, n, ip, jq, kr, tau
real :: baroclinic_triad, K33, K33_crit
real :: sumx(0:1), sumy(0:1)
real :: slope, absslope, depth_ratio
real :: taperA, taperB
real :: taper_slope, taper_slope2, gm_taper_slope
real :: aredi_scalar, agm_scalar, aredi_plus_agm
real :: absdrhodzb(0:1)
real :: mindelqc(0:1,0:1), delqc_ijk_1, delqc_ijkp1_0
call mpp_clock_begin(id_clock_fz_terms)
tau = Time%tau
eady_termx(:,:,:) = 0.0
eady_termy(:,:,:) = 0.0
baroclinic_termx(:,:) = 0.0
baroclinic_termy(:,:) = 0.0
! Main loop. Short ip,jq,kr loops are explicitly unrolled
! in order to expose independence and allow vectorisation
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
aredi_scalar = aredi_array(i,j,k)
agm_scalar = gm_skew*agm_array(i,j,k)
aredi_plus_agm = aredi_scalar + agm_scalar
absdrhodzb(0) = abs(drhodzb(i,j,k,0))
absdrhodzb(1) = abs(drhodzb(i,j,k,1))
delqc_ijk_1 = delqc(i,j,k,1)
delqc_ijkp1_0 = delqc(i,j,k+1,0)
!mindelqc(ip,kr) = min(delqc(i-1+ip,j,k+kr,1-kr),delqc(i+ip,j,k+kr,1-kr))
ip=0 ; kr=0
mindelqc(ip,kr) = min(delqc(i-1,j,k,1),delqc_ijk_1)
ip=0 ; kr=1
mindelqc(ip,kr) = min(delqc(i-1,j,k+1,0),delqc_ijkp1_0)
ip=1 ; kr=0
mindelqc(ip,kr) = min(delqc_ijk_1,delqc(i+1,j,k,1))
ip=1 ; kr=1
mindelqc(ip,kr) = min(delqc_ijkp1_0,delqc(i+1,j,k+1,0))
tx_trans_gm(i,j,k) = 0.0
baroclinic_triad = 0.0
ip=0
sumx(ip) = 0.0
kr=0
slope = tensor_31(i,j,k,ip,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
tx_trans_gm(i,j,k) = tx_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_31(i,j,k,ip,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumx(ip) = sumx(ip) + mindelqc(ip,kr)*taper_slope2
eady_termx(i,j,k) = eady_termx(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
kr=1
slope = tensor_31(i,j,k,ip,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
tx_trans_gm(i,j,k) = tx_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_31(i,j,k,ip,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumx(ip) = sumx(ip) + mindelqc(ip,kr)*taper_slope2
eady_termx(i,j,k) = eady_termx(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
sumx(ip) = sumx(ip)*dtew(i,j,ip)
ip=1
sumx(ip) = 0.0
kr=0
slope = tensor_31(i,j,k,ip,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
tx_trans_gm(i,j,k) = tx_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_31(i,j,k,ip,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumx(ip) = sumx(ip) + mindelqc(ip,kr)*taper_slope2
eady_termx(i,j,k) = eady_termx(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
kr=1
slope = tensor_31(i,j,k,ip,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
tx_trans_gm(i,j,k) = tx_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_31(i,j,k,ip,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumx(ip) = sumx(ip) + mindelqc(ip,kr)*taper_slope2
eady_termx(i,j,k) = eady_termx(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
sumx(ip) = sumx(ip)*dtew(i,j,ip)
if(Thickness%depth_zt(i,j,k) >= agm_closure_upper_depth .and. &
Thickness%depth_zt(i,j,k) <= agm_closure_lower_depth) then
baroclinic_termx(i,j) = baroclinic_termx(i,j) + baroclinic_triad*Thickness%dzwt(i,j,k)
endif
tx_trans_gm(i,j,k) = 0.25*tx_trans_gm(i,j,k)*Grd%dyu(i,j)
!mindelqc(jq,kr) = min(delqc(i,j-1+jq,k+kr,1-kr),delqc(i,j+jq,k+kr,1-kr))
jq=0 ; kr=0
mindelqc(jq,kr) = min(delqc(i,j-1,k,1),delqc_ijk_1)
jq=0 ; kr=1
mindelqc(jq,kr) = min(delqc(i,j-1,k+1,0),delqc_ijkp1_0)
jq=1 ; kr=0
mindelqc(jq,kr) = min(delqc_ijk_1,delqc(i,j+1,k,1))
jq=1 ; kr=1
mindelqc(jq,kr) = min(delqc_ijkp1_0,delqc(i,j+1,k+1,0))
ty_trans_gm(i,j,k) = 0.0
baroclinic_triad = 0.0
jq=0
sumy(jq) = 0.0
kr=0
slope = tensor_32(i,j,k,jq,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
ty_trans_gm(i,j,k) = ty_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_32(i,j,k,jq,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumy(jq) = sumy(jq) + mindelqc(jq,kr)*taper_slope2
eady_termy(i,j,k) = eady_termy(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
kr=1
slope = tensor_32(i,j,k,jq,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
ty_trans_gm(i,j,k) = ty_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_32(i,j,k,jq,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumy(jq) = sumy(jq) + mindelqc(jq,kr)*taper_slope2
eady_termy(i,j,k) = eady_termy(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
sumy(jq) = sumy(jq)*dtns(i,j,jq)
jq=1
sumy(jq) = 0.0
kr=0
slope = tensor_32(i,j,k,jq,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
ty_trans_gm(i,j,k) = ty_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_32(i,j,k,jq,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumy(jq) = sumy(jq) + mindelqc(jq,kr)*taper_slope2
eady_termy(i,j,k) = eady_termy(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
kr=1
slope = tensor_32(i,j,k,jq,kr)
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
taper_slope2 = slope*taper_slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
ty_trans_gm(i,j,k) = ty_trans_gm(i,j,k) + agm_scalar*gm_taper_slope
tensor_32(i,j,k,jq,kr) = aredi_scalar*taper_slope + agm_scalar*gm_taper_slope
sumy(jq) = sumy(jq) + mindelqc(jq,kr)*taper_slope2
eady_termy(i,j,k) = eady_termy(i,j,k) + absdrhodzb(kr)*absslope*absslope
baroclinic_triad = baroclinic_triad + absdrhodzb(kr)*abs(taper_slope)
sumy(jq) = sumy(jq)*dtns(i,j,jq)
if(Thickness%depth_zt(i,j,k) >= agm_closure_upper_depth .and. &
Thickness%depth_zt(i,j,k) <= agm_closure_lower_depth) then
baroclinic_termy(i,j) = baroclinic_termy(i,j) + baroclinic_triad*Thickness%dzwt(i,j,k)
endif
ty_trans_gm(i,j,k) = 0.25*ty_trans_gm(i,j,k)*Grd%dxu(i,j)
K33 = aredi_scalar*Grd%tmask(i,j,k+1)*(Grd%dxtr(i,j)*(sumx(0)+sumx(1)) &
+ Grd%dytr(i,j)*(sumy(0)+sumy(1)))*dzwtr(i,j,k)
! determine part of K33 included explicitly in time and that part implicitly.
! explicit calculation is more accurate than implicit, so aim to compute as much
! as stably possible via the explicit method.
K33_explicit(i,j,k) = K33
K33_implicit(i,j,k) = 0.0
if(.not. diffusion_all_explicit) then
K33_crit = two_dtime_inv*Thickness%rho_dzt(i,j,k,tau)*Thickness%dzt(i,j,k)
if(K33 >= K33_crit) then
K33_explicit(i,j,k) = K33_crit
K33_implicit(i,j,k) = K33-K33_crit
endif
endif
enddo
enddo
enddo
if(tmask_sigma_on) then
do j=jsc,jec
do i=isc,iec
if(tmask_sigma(i,j) > 0.0) then
k = Grd%kmt(i,j)-1
K33_implicit(i,j,k) = K33_implicit(i,j,k)*(1.0-tmask_sigma(i,j))
K33_explicit(i,j,k) = K33_explicit(i,j,k)*(1.0-tmask_sigma(i,j))
endif
enddo
enddo
endif
if(tmask_neutral_on) then
do j=jsc,jec
do i=isc,iec
K33_implicit(i,j,1) = 0.0
K33_explicit(i,j,1) = 0.0
tx_trans_gm(i,j,1) = 0.0
ty_trans_gm(i,j,1) = 0.0
if(Grd%kmt(i,j) > 1) then
k = Grd%kmt(i,j)-1
K33_implicit(i,j,k) = 0.0
K33_explicit(i,j,k) = 0.0
tx_trans_gm(i,j,k) = 0.0
ty_trans_gm(i,j,k) = 0.0
endif
enddo
enddo
endif
do n=1,num_prog_tracers
T_prog(n)%K33_implicit(:,:,:) = K33_implicit(:,:,:)
enddo
if(neutral_physics_limit) then
do n=1,num_prog_tracers
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(T_prog(n)%tmask_limit(i,j,k)==1.0) T_prog(n)%K33_implicit(i,j,k) = 0.0
enddo
enddo
enddo
enddo
endif
do n=1,num_prog_tracers
if (id_k33_implicit(n) > 0) then
used = send_data (id_k33_implicit(n), rho0r*T_prog(n)%K33_implicit(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
enddo
! rho factor to get mass transport.
if(vert_coordinate_class==DEPTH_BASED) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
tx_trans_gm(i,j,k) = tx_trans_gm(i,j,k)*rho0
ty_trans_gm(i,j,k) = ty_trans_gm(i,j,k)*rho0
enddo
enddo
enddo
else
do k=1,nk
do j=jsc,jec
do i=isc,iec
tx_trans_gm(i,j,k) = tx_trans_gm(i,j,k)*rho(i,j,k)
ty_trans_gm(i,j,k) = ty_trans_gm(i,j,k)*rho(i,j,k)
enddo
enddo
enddo
endif
call transport_on_rho_gm (Time, Dens, &
transport_convert*tx_trans_gm(:,:,:), transport_convert*ty_trans_gm(:,:,:))
call transport_on_nrho_gm(Time, Dens, &
transport_convert*tx_trans_gm(:,:,:), transport_convert*ty_trans_gm(:,:,:))
call transport_on_theta_gm(Time, Dens, T_prog(index_temp), &
transport_convert*tx_trans_gm(:,:,:), transport_convert*ty_trans_gm(:,:,:))
! transport_convert converts kg/s to chosen dimensions (either kg/s or Sv (10^9 kg/s))
if(id_tx_trans_gm > 0) then
call mpp_update_domains(tx_trans_gm, Dom%domain2d,flags=EUPDATE)
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
tx_trans_gm(i,j,k) = transport_convert*Grd%dxuer(i,j)* &
(tx_trans_gm(i+1,j,k)*Grd%due(i,j) + &
tx_trans_gm(i,j,k) *Grd%duw(i+1,j))
enddo
enddo
enddo
used = send_data (id_tx_trans_gm, tx_trans_gm(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
! transport_convert converts kg/s to chosen dimensions (either kg/s or Sv (10^9 kg/s))
if(id_ty_trans_gm > 0) then
call mpp_update_domains(ty_trans_gm, Dom%domain2d,flags=NUPDATE)
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
ty_trans_gm(i,j,k) = transport_convert*Grd%dyunr(i,j)* &
(ty_trans_gm(i,j+1,k)*Grd%dun(i,j) + &
ty_trans_gm(i,j,k) *Grd%dus(i,j+1))
enddo
enddo
enddo
used = send_data (id_ty_trans_gm, ty_trans_gm(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
call mpp_clock_end(id_clock_fz_terms)
end subroutine fz_terms
! NAME="fz_terms"
!#######################################################################
!
!
!
! Subroutine computes the zonal neutral physics tracer flux component.
! Compute this component for all tracers at level k.
!
! fx has physical dimensions (area*diffusivity*density*tracer gradient)
!
!
!
subroutine fx_flux(Time, Thickness, T_prog, k)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_prog_tracer_type), intent(in) :: T_prog(num_prog_tracers)
integer, intent(in) :: k
integer :: nn, i, j, ip, tau
integer :: kr, kpkr
real :: triad_weight(isd:ied,jsd:jed)
real :: tensor_11(isd:ied,jsd:jed,0:1,0:1)
real :: tensor_13(isd:ied,jsd:jed,0:1,0:1)
real :: sumz(isd:ied,jsd:jed,0:1)
real :: taperA, taperB
real :: slope, absslope, taper_slope, gm_taper_slope
real :: depth_ratio
real :: drhodx, drhodz
call mpp_clock_begin(id_clock_fx_flux)
tau = Time%tau
! initialize arrays to zero
tensor_11(:,:,:,:) = 0.0
tensor_13(:,:,:,:) = 0.0
grav_agm_dz_sx_drhodx(:,:,k) = 0.0
slopex_drhodx(:,:,k) = 0.0
triad_weight(:,:) = 0.0
! tracer-independent part of the calculation
do kr=0,1
kpkr = min(k+kr,nk)
do ip=0,1
do j=jsc,jec
do i=isc-1,iec
drhodz = drhodzh(i+ip,j,k,kr)
drhodx = Grd%tmask(i+ip,j,kpkr) &
*( drhodT(i+ip,j,k)*dTdx(index_temp)%field(i,j,k) &
+drhodS(i+ip,j,k)*dTdx(index_salt)%field(i,j,k))
slope = -drhodx/drhodz
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
! for diagnosing slope dot grad rho
triad_weight(i,j) = triad_weight(i,j) + Grd%tmask(i+ip,j,kpkr)
slopex_drhodx(i,j,k) = slopex_drhodx(i,j,k) + gm_taper_slope*drhodx
! fill tensor components
tensor_11(i,j,ip,kr) = aredi_array(i,j,k)*( (1.0-taper_diagonal) + taperA*taperB*taper_diagonal) &
+ ah_array(i,j)*(1.0-taperB)
tensor_13(i,j,ip,kr) = aredi_array(i,j,k)*taper_slope-gm_skew*agm_array(i,j,k)*gm_taper_slope
enddo
enddo
enddo
enddo
! tracer-dependent part of the calculation
do nn=1,num_prog_tracers
sumz(:,:,:) = 0.0
do kr=0,1
do ip=0,1
do j=jsc,jec
do i=isc-1,iec
sumz(i,j,kr) = sumz(i,j,kr) + dtwedyt(i,j,ip)*Grd%tmask(i,j,k)*Grd%tmask(i+1,j,k) &
*(tensor_11(i,j,ip,kr)*dTdx(nn)%field(i,j,k) + tensor_13(i,j,ip,kr)*dTdz(nn)%field(i+ip,j,k-1+kr)) &
* min(delqc(i,j,k,kr),delqc(i+1,j,k,kr))
enddo
enddo
enddo
enddo
do j=jsc,jec
do i=isc-1,iec
flux_x(nn)%field(i,j,k) = Grd%dxter(i,j)*(sumz(i,j,0)+sumz(i,j,1))
enddo
enddo
enddo
! apply some masks
if(tmask_neutral_on) then
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc-1,iec
if(k==Grd%kmt(i,j)) then
flux_x(nn)%field(i,j,k) = aredi_array(i,j,k)*dTdx(nn)%field(i,j,k)*Grd%dyte(i,j)* &
min(Thickness%rho_dzt(i,j,k,tau),Thickness%rho_dzt(i+1,j,k,tau))
endif
enddo
enddo
enddo
if(k==1) then
do nn=1,num_prog_tracers
flux_x(nn)%field(:,:,k) = aredi_array(:,:,k)*dTdx(nn)%field(:,:,k)*Grd%dyte(:,:) &
*FMX(Thickness%rho_dzt(:,:,k,tau))
enddo
endif
endif
if(tmask_sigma_on) then
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc-1,iec
if(k==Grd%kmt(i,j)) flux_x(nn)%field(i,j,k) = flux_x(nn)%field(i,j,k)*(1.0-tmask_sigma(i,j))
enddo
enddo
enddo
endif
if(neutral_physics_limit) then
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc-1,iec
if(T_prog(nn)%tmask_limit(i,j,k)==1.0) then
flux_x(nn)%field(i,j,k) = aredi_array(i,j,k)*dTdx(nn)%field(i,j,k)*Grd%dyte(i,j) &
*min(Thickness%rho_dzt(i,j,k,tau),Thickness%rho_dzt(i+1,j,k,tau))
endif
enddo
enddo
enddo
endif
! for diagnosing slopex_drhodx and grav_agm_dz_sx_drhodx
do j=jsc,jec
do i=isc,iec
if(triad_weight(i,j) > 0.0) then
triad_weight(i,j) = 1.0/(epsln + triad_weight(i,j))
slopex_drhodx(i,j,k) = triad_weight(i,j)*slopex_drhodx(i,j,k)
grav_agm_dz_sx_drhodx(i,j,k) = slopex_drhodx(i,j,k) &
*grav*agm_array(i,j,k)*Thickness%dzt(i,j,k)
else
slopex_drhodx(i,j,k) = 0.0
grav_agm_dz_sx_drhodx(i,j,k) = 0.0
endif
enddo
enddo
if(k==nk) then
if(id_grav_agm_dz_sx_drhodx > 0) then
used = send_data (id_grav_agm_dz_sx_drhodx, grav_agm_dz_sx_drhodx(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_slopex_drhodx > 0) then
used = send_data (id_slopex_drhodx, slopex_drhodx(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
endif
call mpp_clock_end(id_clock_fx_flux)
end subroutine fx_flux
! NAME="fx_flux"
!#######################################################################
!
!
!
! Subroutine computes the meridional neutral physics tracer flux component.
! Compute this component for all tracers at level k.
!
! fy has physical dimensions (area*diffusivity*density*tracer gradient)
!
!
!
subroutine fy_flux(Time, Thickness, T_prog, k)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_prog_tracer_type), intent(in) :: T_prog(num_prog_tracers)
integer, intent(in) :: k
integer :: nn, i, j, jq, tau
integer :: kr, kpkr
real :: triad_weight(isd:ied,jsd:jed)
real :: tensor_22(isd:ied,jsd:jed,0:1,0:1)
real :: tensor_23(isd:ied,jsd:jed,0:1,0:1)
real :: sumz(isd:ied,jsd:jed,0:1)
real :: taperA, taperB
real :: slope, absslope, taper_slope, gm_taper_slope
real :: depth_ratio
real :: drhody, drhodz
call mpp_clock_begin(id_clock_fy_flux)
tau = Time%tau
! initialize arrays to zero
tensor_22(:,:,:,:) = 0.0
tensor_23(:,:,:,:) = 0.0
grav_agm_dz_sy_drhody(:,:,k) = 0.0
slopey_drhody(:,:,k) = 0.0
triad_weight(:,:) = 0.0
! tracer-independent part of the calculation
do kr=0,1
kpkr = min(k+kr,nk)
do jq=0,1
do j=jsc-1,jec
do i=isc,iec
drhodz = drhodzh(i,j+jq,k,kr)
drhody = Grd%tmask(i,j+jq,kpkr) &
*( drhodT(i,j+jq,k)*dTdy(index_temp)%field(i,j,k) &
+drhodS(i,j+jq,k)*dTdy(index_salt)%field(i,j,k))
slope = -drhody/drhodz
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
else
taperA = 1.0
endif
! taper for grid depths shallower than eddy_depth
if(eddy_depth(i,j) >= Grd%zt(k)) then
depth_ratio = Grd%zt(k)/(epsln+eddy_depth(i,j))
taperB = 0.5*(1.0 + sin(pi*(depth_ratio-0.5)))
else
taperB = 1.0
endif
! taper times slope for use with neutral diffusion
taper_slope = taperA*taperB*slope
! taper times slope for use with GM skewsion
if(depth_blayer_base(i,j) >= Grd%zt(k)) then
gm_taper_slope = (Grd%zt(k)/depth_blayer_base(i,j))*slope_blayer_base(i,j)*sign(1.0,slope)
else
gm_taper_slope = taper_slope
endif
! for diagnosing slope dot grad rho
triad_weight(i,j) = triad_weight(i,j) + Grd%tmask(i,j+jq,kpkr)
slopey_drhody(i,j,k) = slopey_drhody(i,j,k) + gm_taper_slope*drhody
! fill tensor components
tensor_22(i,j,jq,kr) = aredi_array(i,j,k)*( (1.0-taper_diagonal) + taperA*taperB*taper_diagonal) &
+ ah_array(i,j)*(1.0-taperB)
tensor_23(i,j,jq,kr) = aredi_array(i,j,k)*taper_slope-gm_skew*agm_array(i,j,k)*gm_taper_slope
enddo
enddo
enddo
enddo
! tracer-dependent part of the calculation
do nn=1,num_prog_tracers
sumz(:,:,:) = 0.0
do kr=0,1
do jq=0,1
do j=jsc-1,jec
do i=isc,iec
sumz(i,j,kr) = sumz(i,j,kr) + dxtdtsn(i,j,jq)*Grd%tmask(i,j,k)*Grd%tmask(i,j+1,k) &
*(tensor_22(i,j,jq,kr)*dTdy(nn)%field(i,j,k) + tensor_23(i,j,jq,kr)*dTdz(nn)%field(i,j+jq,k-1+kr)) &
* min(delqc(i,j,k,kr),delqc(i,j+1,k,kr))
enddo
enddo
enddo
enddo
do j=jsc-1,jec
do i=isc,iec
flux_y(nn)%field(i,j,k) = Grd%dytnr(i,j)*(sumz(i,j,0)+sumz(i,j,1))
enddo
enddo
enddo
! apply some masks
if(tmask_neutral_on) then
do nn=1,num_prog_tracers
do j=jsc-1,jec
do i=isc,iec
if(k==Grd%kmt(i,j)) then
flux_y(nn)%field(i,j,k) = aredi_array(i,j,k)*dTdy(nn)%field(i,j,k)*Grd%dxtn(i,j) &
*min(Thickness%rho_dzt(i,j,k,tau),Thickness%rho_dzt(i,j+1,k,tau))
endif
enddo
enddo
enddo
if(k==1) then
do nn=1,num_prog_tracers
flux_y(nn)%field(:,:,k) = &
aredi_array(:,:,k)*dTdy(nn)%field(:,:,k)*Grd%dxtn(:,:)*FMY(Thickness%rho_dzt(:,:,k,tau))
enddo
endif
endif
if(tmask_sigma_on) then
do nn=1,num_prog_tracers
do j=jsc-1,jec
do i=isc,iec
if(k==Grd%kmt(i,j)) flux_y(nn)%field(i,j,k) = flux_y(nn)%field(i,j,k)*(1.0-tmask_sigma(i,j))
enddo
enddo
enddo
endif
if(neutral_physics_limit) then
do nn=1,num_prog_tracers
do j=jsc-1,jec
do i=isc,iec
if(T_prog(nn)%tmask_limit(i,j,k)==1.0) then
flux_y(nn)%field(i,j,k) = aredi_array(i,j,k)*dTdy(nn)%field(i,j,k)*Grd%dxtn(i,j) &
*min(Thickness%rho_dzt(i,j,k,tau),Thickness%rho_dzt(i,j+1,k,tau))
endif
enddo
enddo
enddo
endif
! for diagnosing slopey_drhody and grav_agm_dz_sy_drhody
do j=jsc,jec
do i=isc,iec
if(triad_weight(i,j) > 0.0) then
triad_weight(i,j) = 1.0/(epsln + triad_weight(i,j))
slopey_drhody(i,j,k) = triad_weight(i,j)*slopey_drhody(i,j,k)
grav_agm_dz_sy_drhody(i,j,k) = slopey_drhody(i,j,k) &
*grav*agm_array(i,j,k)*Thickness%dzt(i,j,k)
else
slopey_drhody(i,j,k) = 0.0
grav_agm_dz_sy_drhody(i,j,k) = 0.0
endif
enddo
enddo
if(k==nk) then
if(id_grav_agm_dz_sy_drhody > 0) then
used = send_data (id_grav_agm_dz_sy_drhody, grav_agm_dz_sy_drhody(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_slopey_drhody > 0) then
used = send_data (id_slopey_drhody, slopey_drhody(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_gm_eddy_ke_source > 0) then
used = send_data (id_gm_eddy_ke_source, grav_agm_dz_sx_drhodx(:,:,:) &
+grav_agm_dz_sy_drhody(:,:,:),&
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
endif
call mpp_clock_end(id_clock_fy_flux)
end subroutine fy_flux
! NAME="fy_flux"
!#######################################################################
!
!
!
! Subroutine computes the vertical neutral physics tracer flux component.
! Compute this component for all tracers at level k.
! Surface and bottom boundary condition fz(k=0)=fz(k=kmt(i,j))=0
!
! fz has physical dimensions (density*diffusivity*tracer gradient)
!
!
!
subroutine fz_flux(T_prog, k)
type(ocean_prog_tracer_type), intent(in) :: T_prog(num_prog_tracers)
integer, intent(in) :: k
integer :: nn, i, j, ip, jq, kr
real :: sumx_0, sumx_1, sumy_0, sumy_1
real :: temparray31(isc:iec,jsc:jec,0:1,0:1)
real :: temparray32(isc:iec,jsc:jec,0:1,0:1)
if(tmask_neutral_on .and. k==1) then
do nn = 1, num_prog_tracers
fz2(nn)%field(:,:) = 0.0
enddo
return
endif
if(k==nk) then
do nn = 1, num_prog_tracers
do j=jsc,jec
do i=isc,iec
fz2(nn)%field(i,j) = 0.0
enddo
enddo
enddo
return
elseif(k < nk) then
! tracer-independent part of the calculation
do kr=0,1
do ip=0,1
do j=jsc,jec
do i=isc,iec
temparray31(i,j,ip,kr) = tensor_31(i,j,k,ip,kr)*dtew(i,j,ip) &
*min(delqc(i-1+ip,j,k+kr,1-kr),delqc(i+ip,j,k+kr,1-kr))
enddo
enddo
enddo
do jq=0,1
do j=jsc,jec
do i=isc,iec
temparray32(i,j,jq,kr) = tensor_32(i,j,k,jq,kr)*dtns(i,j,jq) &
*min(delqc(i,j-1+jq,k+kr,1-kr),delqc(i,j+jq,k+kr,1-kr))
enddo
enddo
enddo
enddo
! tracer-dependent part of the calculation
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc,iec
sumx_0 = temparray31(i,j,0,0)*dTdx(nn)%field(i-1,j,k) &
+ temparray31(i,j,0,1)*dTdx(nn)%field(i-1,j,k+1)
sumx_1 = temparray31(i,j,1,0)*dTdx(nn)%field(i,j,k) &
+ temparray31(i,j,1,1)*dTdx(nn)%field(i,j,k+1)
sumy_0 = temparray32(i,j,0,0)*dTdy(nn)%field(i,j-1,k) &
+ temparray32(i,j,0,1)*dTdy(nn)%field(i,j-1,k+1)
sumy_1 = temparray32(i,j,1,0)*dTdy(nn)%field(i,j,k) &
+ temparray32(i,j,1,1)*dTdy(nn)%field(i,j,k+1)
! compute time explicit portion of the vertical flux
fz2(nn)%field(i,j) = Grd%tmask(i,j,k+1) &
*( Grd%dxtr(i,j)*(sumx_0+sumx_1) &
+Grd%dytr(i,j)*(sumy_0+sumy_1)) &
*dzwtr(i,j,k) &
+ K33_explicit(i,j,k)*dTdz(nn)%field(i,j,k)
enddo
enddo
enddo
if(tmask_neutral_on) then
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc,iec
if(k==(Grd%kmt(i,j)-1)) fz2(nn)%field(i,j) = 0.0
enddo
enddo
enddo
endif
if(tmask_sigma_on) then
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc,iec
if(tmask_sigma(i,j) > 0.0 .and. k==(Grd%kmt(i,j)-1) ) then
fz2(nn)%field(i,j) = fz2(nn)%field(i,j)*(1.0-tmask_sigma(i,j))
endif
enddo
enddo
enddo
endif
if(neutral_physics_limit) then
do nn=1,num_prog_tracers
do j=jsc,jec
do i=isc,iec
if(T_prog(nn)%tmask_limit(i,j,k)==1.0) fz2(nn)%field(i,j) = 0.0
enddo
enddo
enddo
endif
endif !if-test for k-level
end subroutine fz_flux
! NAME="fz_flux"
!#######################################################################
!
!
!
! Subroutine computes GM eddy-induced velocity field for diagnostics.
! Compute ustar and vstar at U-cell point, and wstar at T-cell bottom.
!
! Do a two-point average rather than more democratic four-point avg
! in order to avoid having to call mpp_update domains on tensor_31 and
! tensor_32. The 0.5 factor is due to the two-point average.
!
! Note that this algorithm is ad hoc. Researchers interested in this
! field may wish to test alternatives.
!
!
subroutine gm_velocity(Thickness, Time)
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:ied,jsd:jed) :: ugm, vgm
real, dimension(0:1) :: agmSx, agmSy
real :: normalize=0.5
integer :: i, j, k, km1, kp1
integer :: dtint
type(time_type) :: next_time
dtint = nint(dtime)
next_time = increment_time(Time%model_time, dtint, 0)
if(need_data(id_ustar, next_time) .or. &
need_data(id_vstar, next_time) .or. &
need_data(id_wstar, next_time)) then
wrk1_v(:,:,:,:) = 0.0
wrk3(:,:,:) = 0.0
ugm(:,:) = 0.0
vgm(:,:) = 0.0
do k=1,nk
km1 = max(1,k-1)
kp1 = min(nk,k+1)
do j=jsc,jec
do i=isc,iec
agmSx(0) = agm_array(i,j,k)*( &
slope_function_gm(tensor_31(i,j,km1,1,0)) + &
slope_function_gm(tensor_31(i,j,km1,1,1)) )
agmSx(1) = agm_array(i,j,k)*( &
slope_function_gm(tensor_31(i,j,k,1,0)) + &
slope_function_gm(tensor_31(i,j,k,1,1)) )
ugm(i,j) = normalize*(agmSx(0)-agmSx(1))/Thickness%dzt(i,j,k)*Grd%tmask(i,j,k)*Grd%tmask(i+1,j,k)
agmSy(0) = agm_array(i,j,k)*( &
slope_function_gm(tensor_32(i,j,km1,1,0)) + &
slope_function_gm(tensor_32(i,j,km1,1,1)) )
agmSy(1) = agm_array(i,j,k)*( &
slope_function_gm(tensor_32(i,j,k,1,0)) + &
slope_function_gm(tensor_32(i,j,k,1,1)) )
vgm(i,j) = normalize*(agmSy(0)-agmSy(1))/Thickness%dzt(i,j,k)*Grd%tmask(i,j,k)*Grd%tmask(i,j+1,k)
enddo
enddo
call mpp_update_domains (ugm(:,:), Dom%domain2d, flags=EUPDATE)
call mpp_update_domains (vgm(:,:), Dom%domain2d, flags=NUPDATE)
wrk1_v(:,:,k,1) = FAY(ugm(:,:))
wrk1_v(:,:,k,2) = FAX(vgm(:,:))
enddo
call mpp_update_domains (wrk1_v(:,:,:,1), Dom%domain2d, flags=WUPDATE)
call mpp_update_domains (wrk1_v(:,:,:,2), Dom%domain2d, flags=SUPDATE)
k=1
wrk3(:,:,k) = BDX_ET(wrk1_v(:,:,k,1)*Thickness%dzu(:,:,k)) + &
BDY_NT(wrk1_v(:,:,k,2)*Thickness%dzu(:,:,k))
do k=2,nk
wrk3(:,:,k) = wrk3(:,:,k-1) + BDX_ET(wrk1_v(:,:,k,1)*Thickness%dzu(:,:,k)) + &
BDY_NT(wrk1_v(:,:,k,2)*Thickness%dzu(:,:,k))
enddo
if(need_data(id_ustar, next_time)) then
if (id_ustar > 0) used = send_data (id_ustar, wrk1_v(:,:,:,1), &
Time%model_time, rmask=Grd%umask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(need_data(id_vstar, next_time)) then
if (id_vstar > 0) used = send_data (id_vstar, wrk1_v(:,:,:,2), &
Time%model_time, rmask=Grd%umask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(need_data(id_wstar, next_time)) then
if (id_wstar > 0) used = send_data (id_wstar, wrk3(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
else
return
endif
end subroutine gm_velocity
! NAME="gm_velocity"
!#######################################################################
!
!
!
! Function for defining effective slope in diagnostic GM velocity
! calculation. Used only for diagnostic purposes.
!
!
function slope_function_gm (slope)
real, intent(in) :: slope
real :: absslope, slope_function_gm
absslope = abs(slope)
if(absslope >= smax) then
slope_function_gm = smax*sign(1.0,slope)
else
slope_function_gm = slope
endif
end function slope_function_gm
! NAME="slope_function_gm"
!#######################################################################
!
!
! Write out restart files registered through register_restart_file
!
subroutine ocean_nphysicsA_restart(time_stamp)
character(len=*), intent(in), optional :: time_stamp
if(.not. use_this_module) return
call ocean_nphysics_util_restart(time_stamp)
end subroutine ocean_nphysicsA_restart
! NAME="ocean_nphysicsA_restart"
!#######################################################################
!
!
!
! Write to restart.
!
!
subroutine ocean_nphysicsA_end(Time)
type(ocean_time_type), intent(in) :: Time
if(.not. use_this_module) return
call ocean_nphysics_coeff_end(Time, agm_array, aredi_array, rossby_radius, rossby_radius_raw, bczone_radius)
end subroutine ocean_nphysicsA_end
! NAME="ocean_nphysicsA_end"
end module ocean_nphysicsA_mod