module ocean_nphysics_util_mod
#define COMP isc:iec,jsc:jec
!
! Stephen M. Griffies
!
!
! Tim Leslie
!
!
!
! Utilities for neutral physics, including the code to compute
! space-time dependent diffusivities.
!
!
!
! Utilities for neutral physics, including the code to compute
! space-time dependent diffusivities and many diagnostics.
!
!
!
!
!
! D.B. Chelton, R.A. deSzoeke, M.G. Schlax, K.E. Naggar, N. Siwertz
! Geographical Variability of the First Baroclinic Rossby Radius of Deformation
! Journal of Physical Oceanography (1998) vol 28 pages 433-460
!
!
!
! 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
!
!
!
! Held and Larichev
! A scaling theory for horizontally homogeneous baroclinically
! unstable flow on a beta plane
! Journal of Atmospheric Sciences (1996) vol 53 pages 946-952
!
!
!
! M. Visbeck, J.C. Marshall, T. Haine and M. Spall
! Specification of eddy transfer coefficients in coarse resolution ocean
! circulation models
! Journal of Physical Oceanography (1997) vol 27 pages 381--402
!
!
!
! D. Ferreira, J. Marshall, and P. Heimbach,
! Estimating eddy stresses by fitting dynamics to observations
! using a residual-mean ocean circulation omdel and its adjoint.
! Journal of Physical Oceanography (2005) vol 35 pages 1891-1910.
!
!
!
! K. Eden, Eddy length scales in the North Atlantic, 2007,
! Preprint.
!
!
!
! K. Eden and R. Greatbatch, 2008: Towards a mesoscale eddy closure,
! Ocean Modelling, vol. 20, pages 223-239
!
!
!
! Diffusivities can be determined in a number of manners
!
! TIME INDEPENDENT
!
! Various methods are available for specifying a time
! independent diffusivity, either globally uniform or
! with selections of spatial dependence.
!
! TIME DEPENDENT (as a function of the flow)
!
! Various methods are available for determining the
! diffusivity that changes in time according to the
! properties of the fluid. There are various means
! for specifying the length and time scales needed
! to set the diffusivity.
!
! LENGTH SCALES
!
! 1. First baroclinic Rossby radius (estimated as per Chelton etal).
! Equatorial Rossby radius is used within 5deg of the equator.
!
! 2. Width of the baroclinic zone, as done in the Hadley Centre
! model and documented in the MOM3 Manual.
!
! 3. Specified length scale set independent of the flow.
!
! TIME SCALE
!
! When using either of the above for the length scale,
! the time scale is determined by the Eady growth rate.
!
! COMBINED LENGTH/TIME SCALE
!
! Another option, used in the GFDL CM2.X coupled climate models,
! is to set the diffusivity proportional to the depth averaged
! absolute value of the horizontal density gradient.
!
!
!
!
!
!
!
!
! For printing starting and ending checksums for restarts
!
!
!
! For Time%init=.true. and wishing to ensure starting with a clean
! suite of nphysics_util fields, even if ocean_neutral.res.nc exists.
!
!
!
! For computing drhodz used in slope calculation.
! We must keep drhodz < 0 in order to maintain integrity of the
! quasi-Stokes streamfunction as well as computation of buoyancy frequency.
! Default epsln_drhodz=1e-30.
!
!
! For computing the vertical deriviative of locally referenced
! potrho as in the preferred MOM algorithm rather than the
! earlier mom4p0 approach. Default drhodz_mom4p1=.true.
!
!
! For horizontal laplacian smoothing the vertical derivative
! of density prior to its use in computing the neutral slopes.
! This smoothing helps to produce regularized slopes.
! Note that this option breaks the integrity of the triads
! and is thus NOT generally recommended.
! Default drhodz_smooth_horz=.false.
!
!
! For vertical 1-2-1 smoothing the vertical derivative of
! density prior to its use in computing the neutral slopes.
! This smoothing helps to produce regularized slopes.
! Note that this option breaks the integrity of the triads
! and is thus NOT generally recommended.
! Default drhodz_smooth_vert=.false.
!
!
!
! Velocity scale that is used for computing the MICOM Laplacian mixing
! coefficient used in smoothing of drhodzb.
! Default vel_micom_smooth=0.2.
!
!
! For number of 1-2-1 passes through to smooth drhodz or
! eady_rate in vertical. Default num_121_passes=1.
!
!
!
! Neutral diffusivity used for experiments using a constant diffusivity.
!
!
! GM-skew diffusivity used for experiments using a constant diffusivity.
!
!
! Will set aredi_array=agm_array, over-riding any other specification
! of aredi_array. Default aredi_equal_agm=.true.
!
!
!
! Value of the maximum neutral direction slope above which the neutral fluxes are
! either tapered to zero or saturated. Typical value is smax=0.01 or smaller.
!
!
! Width in slope over which use tanh with dm_taper scheme to taper fluxes in
! steep sloped regions. Typical value swidth=0.1*smax
!
!
!
! For calculation of gradients of scalars along a neutral direction, then
! when abs(slope) > smax_grad_gamma_scalar, will compute the gradient using
! only the vertical scalar gradient, since the slopes are so large they are
! effectively infinite.
! Default smax_grad_gamma_scalar=.01
!
!
! If .true., then use a horizontal diffusivity in the neutral boundary layer.
!
!
! Velocity scale that is used for computing the MICOM horizontal diffusivity
! within the neutral boundary layer.
!
!
! Constant horizontal diffusivity for the boundary layer.
! Default ah_bdy=0.0.
!
!
!
! If .true., then the GM-skew diffusivity is set according to a velocity scale
! times the grid spacing.
!
!
! Velocity scale that is used for computing the MICOM diffusivity.
!
!
!
! If true, will set agm_array as constant within two latitudinal zones.
! The idea is that one may wish to use a larger agm in the ACC than
! elsewhere.
!
!
! Boundary between agm in the south and north zones.
!
!
! Ratio between the large agm used in the southern latitudinal zone
! to that used in the north.
! agm_array(north) = agm
! agm_array(south) = agm*agm_lat_zones_ratio
!
!
!
! Set bryan_lewis_aredi=.true. when want to have aredi a function of depth
! according to the Bryan and Lewis (1979) profile. Maintained for legacy
! purposes, and not recommended for new models.
!
!
! ahs = adjustable parameter at the surface for bryan_lewis_aredi
!
!
! ahb = adjustable parameter at the bottom for bryan_lewis_aredi
!
!
!
! For those cases with agm_closure=.false. where we wish to read in
! the agm_array from restart files and keep the value from the restart.
! This approach allows us to read in a spatially dependent agm_array
! that may have been computed from another integration, but to leave
! the coefficient static in time.
! Default agm_read_restart=.false.
!
!
!
! If .true. then will compute the GM-skew diffusivity as a function of the flow.
! The length scale is determined by the Rossby radius and the time scale is
! determined by the Eady growth rate. Diffusivities are depth independent.
!
!
! Maximum GM diffusivity allowed when using agm_closure=.true.
!
!
! Minimum GM diffusivity allowed when using agm_closure=.true.
!
!
! Dimensionless tuning parameter for computing flow dependent diffusivities.
!
!
! Upper depth where start the depth integration to compute the Eady
! growth rate and/or baroclinicity.
!
!
! Deeper depth where finish the depth integration to compute the Eady
! growth rate and/or baroclinicity.
!
!
!
! For computing the agm coefficient using a 3-dimensional
! scaling by (N/Nref)^2, with N=buoyancy frequency and
! Nref the buoyancy frequency at the base of the neutral
! blayer. Default agm_closure_n2_scale=.false.
!
!
! Coefficient setting the scale for the diffusivity computed from
! agm_closure_n2_scale.
! Default agm_closure_n2_scale_coeff=1e3.
!
!
! For taking the reference buoyancy frequency as agm_closure_buoy_freq
! for the (N/Nref)^2 scaling.
! Default agm_closure_n2_scale_nref_cst=.false.
!
!
!
! For computing the agm coefficient using only the vertically
! averaged magnitude of the horizontal density gradient
! (i.e., the "baroclinicity").
!
!
! For computing the agm coefficient using only the vertically
! averaged horizontal density gradient, we need to specify a
! buoyancy frequency, which is taken to be fixed over all space-time.
!
!
!
! For setting a maximum length scale for the agm_closure calculation.
! Default agm_closure_length_cap=.false.
!
!
! Maximum length scale used for computing agm_closure.
! Default agm_closure_length_max=50e3.
!
!
!
! Use fixed length scale for computing agm_closure diffusivity
!
!
! Fixed length scale for use with agm_closure_fixed_length
!
!
!
! For computing the agm_closure length scale according to Rossby radius.
!
!
! Maximum Rossby radius used for agm_closure_length_rossby and
! the neutral_sine_taper. Default = 100e3 m.
!
!
! Minimum Rossby Radius used for agm_closure_length_rossby and
! the neutral_sine_taper. Default = 15e3 m.
!
!
!
! For computing the agm_closure length scale according to radius of baroclinic zone.
!
!
! Max number of horizontal grid points for use in computing the baroclinic zone radius.
!
!
! Critical growth rate for determining width of the baroclinic zone.
!
!
! Dimensionless scaling used to set a maximum for agm_growth.
!
!
!
! For computing the agm_closure length scale according to minimum
! of the Rhines scale and the Rossby radius, and using 3d Eady
! growth rate.
!
!
! For use in regularizing the growth rate used in the eden/greatbatch approach.
! Default agm_closure_eden_gamma=200. Setting to zero removes the regularization.
!
!
! To set the length scale for agm_closure_eden_greatbatch to constant.
! Default agm_closure_eden_length_const=.false.
!
!
! Length scale for use with agm_closure_eden_length_const=.true.
! Default agm_closure_eden_length=10e3.
!
!
!
! For vertical 1-2-1 smoothing the eady_rate
! Default agm_closure_eady_smooth=.false.
!
!
! For horizontal Laplacian smoothing of eady growth rate.
! Default agm_closure_eady_smooth_horz=.false.
!
!
! For capping the eady growth rate to avoid huge values.
! Default agm_closure_eady_cap=.false.
!
!
! To set the Eady growth rate to its average within mixed layer region.
! This is used to avoid spuriously large values which often appear just
! in the upper regions of the ocean mixed layer.
! Default agm_closure_eady_ave_mixed=.false.
!
!
!
! For smoothing the agm diffusivity in space when nonconstant diffusivity used.
! Default is agm_smooth_space=.false.
!
!
! For smoothing the agm diffusivity in time when nonconstant diffusivity used.
! Default is agm_smooth_time=.false.
!
!
! The damping time used for time smoothing agm_array.
! Default agm_damping_time=10days.
!
!
!
! For an overall scaling of the agm coefficient, according to
! the relative resolution of the grid and deformation radius.
! Default is agm_closure_grid_scaling=.false.
!
!
! Power used to scale the agm_closure diffusivity.
! Default is agm_closure_grid_scaling_power=2.0
!
!
! For an overall scaling of the aredi coefficient, according to
! the relative resolution of the grid and deformation radius.
! This option is used only when aredi_equal_agm=.false.
! Default is aredi_diffusivity_grid_scaling=.false.
!
!
!
! For drhodz used in calculation of dianeutral velocity component
! from cabbeling and thermobaricity.
! Default epsln_drhodz_diagnostics=1e-7.
!
!
! For smoothing the diagnosed dianeutral velocity component using a
! horizontal 1-2-1 smoother. Default is wdianeutral_smooth=.true.
!
!
!
! For smoothing the diagnosed contribution to steric sea level
! time tendency associated with the GM90 scheme.
! Default is smooth_eta_tend_gm90=.false.
!
!
use constants_mod, only: epsln, pi
use diag_manager_mod, only: register_diag_field, register_static_field, need_data
use fms_mod, only: FATAL, file_exist
use fms_mod, only: open_namelist_file, check_nml_error, close_file, write_version_number
use fms_io_mod, only: register_restart_field, save_restart, restore_state
use fms_io_mod, only: restart_file_type
use mpp_mod, only: input_nml_file, mpp_pe, mpp_min, mpp_error, stdout, stdlog
use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE
use mpp_domains_mod, only: mpp_update_domains, EUPDATE, NUPDATE
use time_manager_mod, only: time_type, increment_time
use ocean_domains_mod, only: get_local_indices, set_ocean_domain
use ocean_density_mod, only: calc_cabbeling_thermobaricity
use ocean_operators_mod, only: FDX_T, FDY_T, FMX, FMY, LAP_T, S2D
use ocean_parameters_mod, only: missing_value, onehalf, onefourth
use ocean_parameters_mod, only: rho0, rho0r, grav
use ocean_tracer_diag_mod, only: calc_mixed_layer_depth
use ocean_types_mod, only: ocean_grid_type, ocean_domain_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_density_type
use ocean_types_mod, only: ocean_time_type, ocean_thickness_type, ocean_time_steps_type
use ocean_types_mod, only: tracer_3d_0_nk_type, tracer_3d_1_nk_type
use ocean_util_mod, only: write_timestamp, write_chksum_3d, write_chksum_2d, write_chksum_header
use ocean_util_mod, only: diagnose_2d, diagnose_3d, diagnose_sum
use ocean_util_mod, only: write_line, write_note, write_warning, write_chksum_3d
use ocean_tracer_util_mod, only: diagnose_3d_rho
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk6
use ocean_workspace_mod, only: wrk1_v, wrk2_v, wrk3_v, wrk1_2d, wrk2_2d, wrk3_2d
implicit none
public ocean_nphysics_util_init
public ocean_nphysics_coeff_init
public ocean_nphysics_coeff_end
public tracer_derivs
public neutral_slopes
public compute_eady_rate
public compute_baroclinicity
public compute_rossby_radius
public compute_bczone_radius
public compute_diffusivity
public ocean_nphysics_util_restart
public transport_on_rho_gm
public transport_on_nrho_gm
public transport_on_theta_gm
public cabbeling_thermob_tendency
public compute_eta_tend_gm90
public watermass_diag_init
public watermass_diag
public watermass_diag_ndiffuse
public watermass_diag_sdiffuse
private calc_gradx_gamma_scalar
private calc_grady_gamma_scalar
private pressure_derivs
private cabbeling_thermob_init
private
#include
type(ocean_grid_type), pointer :: Grd => NULL()
type(ocean_domain_type), pointer :: Dom => NULL()
type(ocean_domain_type), save :: BCzone_domain
! to write some output to all processors
integer :: unit=6
! clock ids
integer :: id_clock_tracer_derivs
integer :: id_clock_neutral_slopes
integer :: id_clock_compute_eady_rate
integer :: id_clock_compute_baroclinicity
integer :: id_clock_compute_rossby_radius
integer :: id_clock_compute_bczone_radius
integer :: id_clock_compute_diffusivity
integer :: id_clock_cabbeling_thermob
integer :: id_clock_eta_tend_gm90
integer :: id_transport_on_nrho_gm
integer :: id_transport_on_rho_gm
integer :: id_transport_on_theta_gm
! for diagnostics manager
integer :: id_slope31 =-1
integer :: id_slope32 =-1
integer :: id_aredi =-1
integer :: id_agm =-1
integer :: id_aredi_3d =-1
integer :: id_agm_3d =-1
integer :: id_eady_mld =-1
integer :: id_eady_rate =-1
integer :: id_eady_rate_zave =-1
integer :: id_rossby =-1
integer :: id_rossby_radius =-1
integer :: id_rossby_equator =-1
integer :: id_rossby_nonequator =-1
integer :: id_bczone =-1
integer :: id_gw_speed =-1
integer :: id_sqrt2betaCr =-1
integer :: id_growth_rate_baroclinic =-1
integer :: id_baroclinicity =-1
integer :: id_agm_growth_rate =-1
integer :: id_agm_length =-1
integer :: id_rhines_length =-1
integer :: id_agm_qg =-1
integer :: id_N2slope =-1
integer :: id_N2_nblayer_base =-1
integer :: id_N2_for_agm =-1
integer :: id_coriolis_param =-1
integer :: id_beta_param =-1
integer :: id_grid_length =-1
integer :: id_agm_grid_scaling =-1
integer :: id_aredi_grid_scaling =-1
integer :: id_ksurf_blayer =-1
integer :: id_cabbeling_tend =-1
integer :: id_thermobaric_tend =-1
integer :: id_cabbeling_speed =-1
integer :: id_thermobaric_speed =-1
integer :: id_cabbeling_tend_intz =-1
integer :: id_thermobaric_tend_intz =-1
integer :: id_eta_tend_gm90 =-1
integer :: id_eta_tend_gm90_glob =-1
integer :: id_cabbeling_tend_intz_glob =-1
integer :: id_thermobaric_tend_intz_glob =-1
integer :: id_neut_rho_cabbeling =-1
integer :: id_neut_rho_cabbeling_on_nrho =-1
integer :: id_wdian_cabbeling =-1
integer :: id_wdian_cabbeling_on_nrho =-1
integer :: id_tform_rho_cabbel_on_nrho =-1
integer :: id_neut_rho_thermob =-1
integer :: id_neut_rho_thermob_on_nrho =-1
integer :: id_wdian_thermob =-1
integer :: id_wdian_thermob_on_nrho =-1
integer :: id_tform_rho_thermb_on_nrho =-1
! for watermass transformation diagnostics
integer :: id_neut_rho_nphysics
integer :: id_wdian_rho_nphysics
integer :: id_neut_rho_ndiff
integer :: id_wdian_rho_ndiff
integer :: id_tform_rho_ndiff
integer :: id_neut_rho_ndiff_on_nrho
integer :: id_wdian_rho_ndiff_on_nrho
integer :: id_tform_rho_ndiff_on_nrho
integer :: id_eta_tend_ndiff_tend
integer :: id_eta_tend_ndiff_tend_glob
integer :: id_neut_temp_ndiff
integer :: id_wdian_temp_ndiff
integer :: id_tform_temp_ndiff
integer :: id_neut_temp_ndiff_on_nrho
integer :: id_wdian_temp_ndiff_on_nrho
integer :: id_tform_temp_ndiff_on_nrho
integer :: id_neut_salt_ndiff
integer :: id_wdian_salt_ndiff
integer :: id_tform_salt_ndiff
integer :: id_neut_salt_ndiff_on_nrho
integer :: id_wdian_salt_ndiff_on_nrho
integer :: id_tform_salt_ndiff_on_nrho
integer :: id_neut_rho_gm
integer :: id_wdian_rho_gm
integer :: id_tform_rho_gm
integer :: id_neut_rho_gm_on_nrho
integer :: id_wdian_rho_gm_on_nrho
integer :: id_tform_rho_gm_on_nrho
integer :: id_eta_tend_gm_tend
integer :: id_eta_tend_gm_tend_glob
integer :: id_neut_temp_gm
integer :: id_wdian_temp_gm
integer :: id_tform_temp_gm
integer :: id_neut_temp_gm_on_nrho
integer :: id_wdian_temp_gm_on_nrho
integer :: id_tform_temp_gm_on_nrho
integer :: id_neut_salt_gm
integer :: id_wdian_salt_gm
integer :: id_tform_salt_gm
integer :: id_neut_salt_gm_on_nrho
integer :: id_wdian_salt_gm_on_nrho
integer :: id_tform_salt_gm_on_nrho
! for transport_on_nrho_gm
integer :: id_tx_trans_nrho_gm=-1
integer :: id_ty_trans_nrho_gm=-1
! for transport_on_rho_gm
integer :: id_tx_trans_rho_gm=-1
integer :: id_ty_trans_rho_gm=-1
! for transport_on_theta_gm
integer :: id_tx_trans_theta_gm=-1
integer :: id_ty_trans_theta_gm=-1
! for restart
type(restart_file_type), save :: Nphysics_util_restart
integer :: num_prog_tracers=0
integer :: index_temp=-1
integer :: index_salt=-1
integer :: neutralrho_nk
logical :: used
! internally set for computing watermass diagnstics
logical :: compute_watermass_diag = .false.
! 2D array of micom gm diffusivities (m^2/sec)
real, dimension(:,:), allocatable :: agm_micom
! 2D array of deformation radius scaling
real, dimension(:,:), allocatable :: agm_grid_scaling
real, dimension(:,:), allocatable :: aredi_grid_scaling
! 3D array of gm length scale (metre)
real, dimension(:,:,:), allocatable :: agm_length
! 3D array of agm_array for local purposes (m^2/sec)
real, dimension(:,:,:), allocatable :: agm_array_local
! 3D array of aredi_array for local purposes (m^2/sec)
real, dimension(:,:,:), allocatable :: aredi_array_local
! 3D array of gm growth rate scale (sec^-1)
real, dimension(:,:,:), allocatable :: agm_growth_rate
! 2D array of gm growth rate max (sec^-1)
real, dimension(:,:), allocatable :: agm_growth_rate_max
! grid length array (metre)
real, dimension(:,:), allocatable :: grid_length
! absolute value of the Coriolis parameter (sec^-1)
real, dimension(:,:), allocatable :: coriolis_param
! beta = d(Coriolis)/dy (m^-1 sec^-1)
real, dimension(:,:), allocatable :: beta_param
! for determining baroclinic zone radius (metre)
real, dimension(:,:), allocatable :: bczone_dxt
real, dimension(:,:), allocatable :: bczone_dyt
real, dimension(:,:), allocatable :: bczone_rate
! for Eady growth rate and baroclinicity calculation
real, dimension(:,:), allocatable :: count_x
real, dimension(:,:), allocatable :: count_y
! Eady growth rate (sec^-1)
real, dimension(:,:,:), allocatable :: eady_rate
real, dimension(:,:), allocatable :: eady_rate_zave
! 1st baroclinic gravity wave speed (m/s)
real, dimension(:,:), allocatable :: gravity_wave_speed
! for Eden and Greatbatch near equator (m)
real, dimension(:,:), allocatable :: sqrt2betaCr
!vertically ave abs(horiz density gradient) (kg/m^4)
real, dimension(:,:), allocatable :: baroclinicity
! (m^2/sec) for smoothing
real, dimension(:,:), allocatable :: smooth_lap
! cabbeling and thermobaricity tendencies
real, dimension(:,:,:), allocatable :: cabbeling_param
real, dimension(:,:,:), allocatable :: thermobaric_param
real, dimension(:,:,:), allocatable :: gradx_gamma_temp
real, dimension(:,:,:), allocatable :: grady_gamma_temp
real, dimension(:,:,:), allocatable :: gradx_gamma_press
real, dimension(:,:,:), allocatable :: grady_gamma_press
real, dimension(:,:,:), allocatable :: deriv_x_press
real, dimension(:,:,:), allocatable :: deriv_y_press
real, dimension(:,:,:), allocatable :: deriv_z_press
! gm contribution to steric sea level tendency
real, dimension(:,:,:), allocatable :: taper_fcn
real, dimension(:,:,:), allocatable :: slopex
real, dimension(:,:,:), allocatable :: agm_rho_slopex
real, dimension(:,:,:), allocatable :: dslopex_dx
real, dimension(:,:,:), allocatable :: dslopex_dz
real, dimension(:,:,:), allocatable :: gradx_gamma_slopex
real, dimension(:,:,:), allocatable :: slopey
real, dimension(:,:,:), allocatable :: agm_rho_slopey
real, dimension(:,:,:), allocatable :: dslopey_dy
real, dimension(:,:,:), allocatable :: dslopey_dz
real, dimension(:,:,:), allocatable :: grady_gamma_slopey
! some useful constants
real :: grav_r
real :: gravrho0r
real :: gravrho0r_buoyr
real :: gamma_damp
real :: dtime
real :: cellarea_r
! for specifying transport units
! can either be Sv or mks
character(len=32) :: transport_dims ='Sv (10^9 kg/s)'
real :: transport_convert=1.0e-9
!*************************************
! nml settings
!
! for how to compute drhodz prior to computing slopes
logical :: drhodz_mom4p1 =.true.
logical :: drhodz_smooth_horz =.false.
logical :: drhodz_smooth_vert =.false.
integer :: num_121_passes = 1
real :: vel_micom_smooth = 0.2 ! m/sec
real :: epsln_drhodz = 1e-30
! globally constant diffusivities
real :: aredi = 1.0e3 ! constant neutral diffusion tracer diffusivity (m^2/sec)
real :: agm = 1.0e3 ! constant gent-mcwilliams skew-diffusion diffusivity (m^2/sec)
! will set aredi_array=agm_array
logical :: aredi_equal_agm =.true.
! maximum neutral direction slope allowed before tapering fluxes
real :: smax=0.01
! slope width over which fluxes tapered using tanh function using dm_taper scheme
real :: swidth = 0.05*.01
! if true, apply horizontal diffusivity in neutral boundary layer for nphysicsA module
logical :: neutral_horz_mix_bdy =.false.
real :: vel_micom_bdy = 0.0 ! velocity scale (m/s) for horizontal diffusivity w/i boundary
real :: ah_bdy = 0.0 ! constant horiz diffusivity in neutral bdy
! time independent, vertically dependent neutral diffusivity according to Bryan-Lewis.
! maintained in MOM only for legacy purposes. not recommended for new models.
logical :: bryan_lewis_aredi = .false.
real :: ahs = 0.0
real :: ahb = 0.0
! for setting agm according to latitude zones
logical :: agm_lat_bands = .false.
real :: agm_lat_bands_boundary = -999. ! boundary between agm in the south and north zones
real :: agm_lat_bands_ratio = 1.0 ! ratio agm(south)/agm(north)
! set diffusivity according to local horizontal area
! of grid and a specified velocity scale (m/s)
logical :: tracer_mix_micom =.false.
real :: vel_micom = 0.0
! determining how to read in restart diffusivity
logical :: agm_read_restart=.false.
! for setting diffusivity according to flow properties
logical :: agm_closure =.false.
real :: agm_closure_scaling = 2.0 ! dimensionless tuning parameter
real :: agm_closure_max = 2.e3 ! maximum diffusivity allowed when agm_closure=.true.
real :: agm_closure_min = 2.e2 ! minimum diffusivity allowed when agm_closure=.true.
! for fixed length scale (metres) set by agm_closure_length
logical :: agm_closure_length_fixed =.false.
real :: agm_closure_length = 50.e3 ! metres
! for length scale set according to estimate of first baroclinic Rossby radius
logical :: agm_closure_length_rossby =.false.
real :: rossby_radius_max=100e3 ! metres
real :: rossby_radius_min=15e3 ! metres
! for length scale set according to radius of baroclinic zone
logical :: agm_closure_length_bczone =.false.
integer :: bczone_max_pts =10 ! max # points searched for determining baroclinic zone width
real :: agm_closure_bczone_crit_rate =1.4e-6 ! critical growth rate for determining baroclinic zone (sec^-1)
! for capping the agm_length scale
logical :: agm_closure_length_cap = .false.
real :: agm_closure_length_max = 50e3
! for scaling diffusivity by (N/Nref)^2
logical :: agm_closure_n2_scale=.false.
logical :: agm_closure_n2_scale_nref_cst=.false.
real :: agm_closure_n2_scale_coeff=1e3
! for diffusity according to Eden and Greatbatch (2008) and Eden(2007)
logical :: agm_closure_eden_greatbatch =.false.
logical :: agm_closure_eden_length_const=.false.
real :: agm_closure_eden_gamma = 200.0
real :: agm_closure_eden_length = 10e3
! for diffusivity set proportional to vertically averaged horiz density gradient
logical :: agm_closure_baroclinic = .true.
! sec^-1 buoyancy frequency for use with agm_closure_baroclinic
real :: agm_closure_buoy_freq = 0.004
! depths over which compute the Eady growth and/or horizontal density gradient
real :: agm_closure_upper_depth =100.0 ! metre
real :: agm_closure_lower_depth =2000.0 ! metre
! for smoothing and capping eady growth rate
logical :: agm_closure_eady_smooth_vert=.false.
logical :: agm_closure_eady_smooth_horz=.false.
logical :: agm_closure_eady_cap =.false.
logical :: agm_closure_eady_ave_mixed =.false.
! dimensionless number to set maximum value of agm_growth.
! agm_closure_growth_scale =0.5 yields agm_growth max=coriolis parameter.
real :: agm_closure_growth_scale =0.5
! for smoothing agm_array in space and time
logical :: agm_smooth_space=.false.
logical :: agm_smooth_time =.false.
real :: agm_damping_time=10.0
! for overall scaling of the agm coefficient
logical :: agm_closure_grid_scaling=.false.
real :: agm_closure_grid_scaling_power=2.0
! for scaling aredi_array according to size of grid and Rossby radius
logical :: aredi_diffusivity_grid_scaling=.false.
! for Tim%init=.true. and ocean_neutral.res.nc exists, but still wish
! to start from initial fields.
logical :: nphysics_util_zero_init=.false.
! for computing tendencies from thermobaricity and cabbeling
logical :: diagnose_cabbeling_thermob=.false. ! internally set
logical :: wdianeutral_smooth=.true.
real :: epsln_drhodz_diagnostics=1e-7
real :: smax_grad_gamma_scalar=.01
real :: swidthr ! inverse swidth
real :: smax_swidthr ! useful combination of terms
! for computing contribution to steric sea level evolution from GM
logical :: smooth_eta_tend_gm90 =.false.
logical :: diagnose_eta_tend_gm90 =.false. ! internally set
! for debugging
logical :: debug_this_module = .false.
!*************************************
character(len=128) :: version=&
'$Id: ocean_nphysics_util.F90,v 20.0 2013/12/14 00:14:52 fms Exp $'
character (len=128) :: tagname = &
'$Name: tikal $'
character(len=*), parameter :: FILENAME=&
__FILE__
logical :: module_is_initialized = .FALSE.
namelist /ocean_nphysics_util_nml/ debug_this_module, nphysics_util_zero_init, &
smax, swidth, epsln_drhodz, drhodz_mom4p1, &
drhodz_smooth_horz, drhodz_smooth_vert, num_121_passes, &
aredi, agm, aredi_equal_agm, tracer_mix_micom, vel_micom, &
bryan_lewis_aredi, ahs, ahb, &
neutral_horz_mix_bdy, vel_micom_bdy, ah_bdy, &
agm_lat_bands, agm_lat_bands_boundary, agm_lat_bands_ratio, &
rossby_radius_max, rossby_radius_min, agm_read_restart, &
agm_closure, agm_closure_scaling, agm_closure_max, agm_closure_min, &
agm_closure_growth_scale, agm_closure_length_fixed, agm_closure_length, &
agm_closure_length_rossby, agm_closure_length_bczone, &
bczone_max_pts, agm_closure_bczone_crit_rate, &
agm_closure_eden_greatbatch, agm_closure_eden_gamma, &
agm_closure_eden_length_const, agm_closure_eden_length, &
agm_closure_eady_smooth_vert, agm_closure_eady_smooth_horz, &
agm_closure_eady_ave_mixed, agm_closure_eady_cap, &
agm_closure_baroclinic, agm_closure_buoy_freq, &
agm_closure_upper_depth, agm_closure_lower_depth, &
agm_closure_length_cap, agm_closure_length_max, &
agm_smooth_space, vel_micom_smooth, agm_smooth_time, agm_damping_time, &
agm_closure_grid_scaling, agm_closure_grid_scaling_power, &
aredi_diffusivity_grid_scaling, &
agm_closure_n2_scale, agm_closure_n2_scale_coeff, &
agm_closure_n2_scale_nref_cst, &
smax_grad_gamma_scalar, &
epsln_drhodz_diagnostics, wdianeutral_smooth, &
smooth_eta_tend_gm90
contains
!#######################################################################
!
!
!
! Initialize the utility module for neutral physics.
!
!
subroutine ocean_nphysics_util_init(Grid, Domain, Time, Time_steps, Dens, T_prog, &
agm_closure_lower_dept, agm_closure_upper_dept, agm_closure_buoy_frq, &
smx, swidt, cmip_units, 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_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
real, intent(inout) :: agm_closure_lower_dept
real, intent(inout) :: agm_closure_upper_dept
real, intent(inout) :: agm_closure_buoy_frq
real, intent(inout) :: smx
real, intent(inout) :: swidt
logical, intent(in) :: cmip_units
logical, intent(in), optional :: debug
integer :: ioun, io_status, ierr
integer :: i,j,n
integer :: num_schemes
real :: param
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if ( module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error from ocean_nphysics_util_mod: already initialized')
endif
module_is_initialized = .TRUE.
call write_version_number(version, tagname)
Dom => Domain
Grd => Grid
! provide for namelist over-ride of default values
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_nphysics_util_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_nphysics_util_nml')
#else
ioun = open_namelist_file()
read (ioun,ocean_nphysics_util_nml,IOSTAT=io_status)
ierr = check_nml_error(io_status,'ocean_nphysics_util_nml')
call close_file (ioun)
#endif
write (stdoutunit,'(/)')
write (stdoutunit,ocean_nphysics_util_nml)
write (stdlogunit,ocean_nphysics_util_nml)
if (PRESENT(debug) .and. .not. debug_this_module) then
debug_this_module = debug
endif
if(debug_this_module) then
call write_note(FILENAME,&
'running ocean_nphysics_util_mod with debug_this_module=.true.')
endif
#ifndef MOM_STATIC_ARRAYS
call get_local_indices(Domain,isd,ied,jsd,jed,isc,iec,jsc,jec)
nk = Grid%nk
#endif
! for setting up temp and salt tracers
num_prog_tracers = size(T_prog(:))
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_nphysics_util_init')
endif
! some useful constants
dtime = Time_steps%dtime_t
grav_r = 1.0/grav
gravrho0r = grav*rho0r !for buoyancy frequency calculation
gravrho0r_buoyr = gravrho0r/agm_closure_buoy_freq !for agm_closure_baroclinic
gamma_damp = dtime/(86400.0*agm_damping_time) !for damping time dependent agm_array
swidthr = 1.0/(swidth + epsln) !for slope taper function for cabbeling/thermob diagnostic
smax_swidthr = smax_grad_gamma_scalar*swidthr !for slope taper function for cabbeling/thermob diagnostic
neutralrho_nk = size(Dens%neutralrho_ref(:))
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
if(cmip_units) then
transport_convert=1.0
transport_dims = 'kg/s'
else
transport_convert=1.0e-9
transport_dims = 'Sv (10^9 kg/s)'
endif
! for sending back to ocean_nphysics_mod
agm_closure_lower_dept = agm_closure_lower_depth
agm_closure_upper_dept = agm_closure_upper_depth
agm_closure_buoy_frq = agm_closure_buoy_freq
smx = smax
swidt = swidth
! Coriolis parameter and beta parameter
allocate (coriolis_param(isd:ied,jsd:jed))
allocate (beta_param(isd:ied,jsd:jed))
coriolis_param(:,:) = 2.0*7.292e-5*abs(sin(Grd%phit(:,:)))
beta_param(:,:) = 2.28e-11*abs(cos(Grd%phit(:,:)))
allocate (gravity_wave_speed(isd:ied,jsd:jed))
allocate (sqrt2betaCr(isd:ied,jsd:jed))
allocate (eady_rate(isd:ied,jsd:jed,nk))
allocate (eady_rate_zave(isd:ied,jsd:jed))
allocate (baroclinicity(isd:ied,jsd:jed))
sqrt2betaCr(:,:) = 0.0
gravity_wave_speed(:,:) = 0.0
eady_rate(:,:,:) = 0.0
eady_rate_zave(:,:) = 0.0
baroclinicity(:,:) = 0.0
! for computing the baroclinic zone radius using Hadley Centre search algorithm
call set_ocean_domain(BCzone_domain,Grd,xhalo=bczone_max_pts,yhalo=bczone_max_pts,name='bczone',maskmap=Dom%maskmap)
allocate (bczone_rate(isc-bczone_max_pts:iec+bczone_max_pts,jsc-bczone_max_pts:jec+bczone_max_pts))
allocate (bczone_dxt (isc-bczone_max_pts:iec+bczone_max_pts,jsc-bczone_max_pts:jec+bczone_max_pts))
allocate (bczone_dyt (isc-bczone_max_pts:iec+bczone_max_pts,jsc-bczone_max_pts:jec+bczone_max_pts))
bczone_dxt=0.0 ; bczone_dyt=0.0 ; bczone_rate=0.0
bczone_dxt(isc:iec,jsc:jec) = Grd%dxt(isc:iec,jsc:jec)
bczone_dyt(isc:iec,jsc:jec) = Grd%dyt(isc:iec,jsc:jec)
call mpp_update_domains (bczone_dxt(:,:), BCzone_domain%domain2d)
call mpp_update_domains (bczone_dyt(:,:), BCzone_domain%domain2d)
allocate (grid_length(isd:ied,jsd:jed))
grid_length(:,:) = 2.0*Grd%dxt(:,:)*Grd%dyt(:,:)/(Grd%dxt(:,:)+Grd%dyt(:,:))
allocate (agm_micom(isd:ied,jsd:jed))
if(tracer_mix_micom) then
agm_micom(:,:) = vel_micom*grid_length(:,:)
else
agm_micom(:,:) = 0.0
endif
allocate (agm_length(isd:ied,jsd:jed,nk))
agm_length(:,:,:) = agm_closure_length*Grd%tmask(:,:,:)
allocate (agm_grid_scaling(isd:ied,jsd:jed))
agm_grid_scaling(:,:) = Grd%tmask(:,:,1)
allocate (aredi_grid_scaling(isd:ied,jsd:jed))
aredi_grid_scaling(:,:) = Grd%tmask(:,:,1)
allocate (agm_growth_rate(isd:ied,jsd:jed,nk))
agm_growth_rate(:,:,:) = 0.0
allocate (agm_array_local(isd:ied,jsd:jed,nk))
agm_array_local(:,:,:) = agm*Grd%tmask(:,:,:)
allocate (aredi_array_local(isd:ied,jsd:jed,nk))
aredi_array_local(:,:,:) = aredi*Grd%tmask(:,:,:)
!max agm_growth_rate is 2.0*agm_closure_growth_scale*param
!if agm_closure_growth_scale=0.5, then agm_growth_rate is <= param
allocate (agm_growth_rate_max(isd:ied,jsd:jed))
do j=jsd,jed
do i=isd,ied
param = max(coriolis_param(i,j), sqrt2betaCr(i,j))
agm_growth_rate_max(i,j) = agm_closure_growth_scale*param
enddo
enddo
if(Time_steps%aidif /= 1.0 .and. aredi /= 0) then
call mpp_error(FATAL, &
'==>Error from ocean_nphysics_util_mod: stability requires aidif=1.0 to handle K33 implicitly.')
endif
if(agm_closure) then
num_schemes = 0
if(agm_closure_length_fixed) then
num_schemes = num_schemes + 1
call write_note(FILENAME,&
'Computing 2d flow-dependent tracer diffusivity with agm_closure_length_fixed.')
endif
if(agm_closure_length_rossby) then
num_schemes = num_schemes + 1
call write_note(FILENAME,&
'Computing 2d flow-dependent tracer diffusivity with agm_closure_length_rossby.')
endif
if(agm_closure_length_bczone) then
num_schemes = num_schemes + 1
call write_note(FILENAME,&
'Computing 2d flow-dependent tracer diffusivity with agm_closure_length_bczone.')
endif
if(agm_closure_eden_greatbatch) then
num_schemes = num_schemes + 1
call write_note(FILENAME,&
'Computing 3d flow-dependent tracer diffusivity with agm_closure_eden_greatbatch.')
endif
if(agm_closure_baroclinic) then
num_schemes = num_schemes + 1
call write_note(FILENAME,&
'Computing 2d flow-dependent tracer diffusivity with agm_closure_baroclinic.')
endif
if(agm_closure_n2_scale) then
num_schemes = num_schemes + 1
call write_note(FILENAME,&
'Computing 3d flow-dependent tracer diffusivity with agm_closure_n2_scale.')
endif
if(num_schemes == 0) then
call mpp_error(FATAL, &
'==>Error: with agm_closure=.true., must choose one of the "agm_closure" methods')
endif
if(num_schemes > 1) then
call mpp_error(FATAL, &
'==>Error: with agm_closure=.true., can choose only one of the available agm_closure methods')
endif
write(stdoutunit,'(a,e10.5)') &
' The maximum allowable diffusivity (m^2/s) is given by ',agm_closure_max
write(stdoutunit,'(a,e10.5)') &
' The minimum allowable diffusivity (m^2/s) is given by ',agm_closure_min
write(stdoutunit,'(a,e10.5,a,e10.5)') &
' Depths (m) between which compute eady growth and baroclinicity = ', &
agm_closure_upper_depth, ' ',agm_closure_lower_depth
if(agm_closure_n2_scale) then
call write_note(FILENAME,&
'Diffusivity will be computed as agm = coeff*(N/Nref)^2.')
endif
if(agm_closure_baroclinic) then
call write_note(FILENAME,&
'Length and time scales set by vertically averaged baroclinicity |grad(rho)|,')
write(stdoutunit,'(a,e10.5)') &
' as well as the constant buoyancy freq(sec^-1) = ',agm_closure_buoy_freq
write(stdoutunit,'(a,e10.5)') &
' and the constant length scale (m) = ',agm_closure_length
else
call write_note(FILENAME,&
'Eady growth rate gives inverse time scale.')
if(agm_closure_length_fixed) then
write(stdoutunit,'(a,e12.5,a)') &
' Length scale set by nml parameter agm_closure_length =',agm_closure_length,' metre'
elseif(agm_closure_length_rossby) then
call write_line('First baroclinic Rossby radius gives the length scale.')
elseif(agm_closure_length_bczone) then
call write_line('Radius of baroclinic zone gives the length scale.')
elseif(agm_closure_eden_greatbatch) then
if(agm_closure_eden_length_const) then
write(stdoutunit,'(a,e12.5,a)') &
' Take constant length scale of ',agm_closure_eden_length, 'metre'
else
call write_line('Min(rossby,rhines) gives the length scale.')
endif
endif
endif
endif ! endif for agm_closure
! for smoothing
allocate (smooth_lap(isd:ied,jsd:jed))
smooth_lap(:,:) = vel_micom_smooth*grid_length(:,:)
! for Eady growth rate and baroclinicity calculation
allocate (count_x(isd:ied,jsd:jed))
allocate (count_y(isd:ied,jsd:jed))
count_x(:,:) = 0.0
count_y(:,:) = 0.0
do i=isc,iec
do j=jsc,jec
count_x(i,j) = &
min(Grd%tmask(i+1,j,1) + Grd%tmask(i-1,j,1), 1.0/(2.0*(Grd%tmask(i+1,j,1) + Grd%tmask(i-1,j,1) + epsln)))
count_y(i,j) = &
min(Grd%tmask(i,j+1,1) + Grd%tmask(i,j-1,1), 1.0/(2.0*(Grd%tmask(i,j+1,1) + Grd%tmask(i,j-1,1) + epsln)))
enddo
enddo
! initialize clock ids
id_clock_tracer_derivs = mpp_clock_id('(Ocean neutral: tracer derivs)' ,grain=CLOCK_ROUTINE)
id_clock_neutral_slopes = mpp_clock_id('(Ocean neutral: slopes)' ,grain=CLOCK_ROUTINE)
id_clock_compute_eady_rate = mpp_clock_id('(Ocean neutral: eady rate)' ,grain=CLOCK_ROUTINE)
id_clock_compute_baroclinicity = mpp_clock_id('(Ocean neutral: baroclinic)' ,grain=CLOCK_ROUTINE)
id_clock_compute_rossby_radius = mpp_clock_id('(Ocean neutral: rossby radius)' ,grain=CLOCK_ROUTINE)
id_clock_compute_bczone_radius = mpp_clock_id('(Ocean neutral: bc zone radius)' ,grain=CLOCK_ROUTINE)
id_clock_compute_diffusivity = mpp_clock_id('(Ocean neutral: diffusivity)' ,grain=CLOCK_ROUTINE)
id_clock_cabbeling_thermob = mpp_clock_id('(Ocean neutral: cabbel/thermob)' ,grain=CLOCK_ROUTINE)
id_clock_eta_tend_gm90 = mpp_clock_id('(Ocean neutral: eta_tend_gm90)' ,grain=CLOCK_ROUTINE)
id_transport_on_nrho_gm = mpp_clock_id('(Ocean neutral: nrho-trans)' ,grain=CLOCK_ROUTINE)
id_transport_on_rho_gm = mpp_clock_id('(Ocean neutral: rho-trans)' ,grain=CLOCK_ROUTINE)
id_transport_on_theta_gm = mpp_clock_id('(Ocean neutral: theta-trans)' ,grain=CLOCK_ROUTINE)
! for diagnostics manager
id_ksurf_blayer= -1
id_ksurf_blayer= register_diag_field ('ocean_model', 'ksurf_blayer', &
Grd%tracer_axes(1:2), Time%model_time, &
'k-value at base of surf nblayer region', 'dimensionless', &
missing_value=missing_value, range=(/-1.0,1.e1/))
id_N2slope = -1
id_N2slope = register_diag_field ('ocean_model', 'N2slope', &
Grd%tracer_axes_wt(1:3), Time%model_time, &
'Squared buoyancy frequency used in neutral slope calculation', &
'1/s^2', missing_value=missing_value, range=(/-1e3,1.e10/))
id_N2_for_agm = -1
id_N2_for_agm = register_diag_field ('ocean_model', 'N2_for_agm', &
Grd%tracer_axes(1:3), Time%model_time, &
'Squared buoyancy frequency used for (N/Nref)^2 agm calculation', &
'1/s^2', missing_value=missing_value, range=(/-1e3,1.e10/))
id_N2_nblayer_base= -1
id_N2_nblayer_base= register_diag_field ('ocean_model', 'N2_nblayer_base',&
Grd%tracer_axes(1:2), Time%model_time, &
'Squared buoyancy frequency at base of nblayer', &
'1/s^2', missing_value=missing_value, range=(/-1e3,1.e10/))
id_slope31 = -1
id_slope31 = register_diag_field ('ocean_model', 'slope31', &
Grd%tracer_axes(1:3), Time%model_time, &
'neutral slope -(rho_x/rho_z)', 'dimensionless', &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_slope32 = -1
id_slope32 = register_diag_field ('ocean_model', 'slope32', &
Grd%tracer_axes(1:3), Time%model_time, &
'neutral slope -(rho_y/rho_z)', 'dimensionless', &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_eady_mld = -1
id_eady_mld = register_diag_field ('ocean_model', 'eady_mld', &
Grd%tracer_axes(1:2), Time%model_time, &
'Mixed layer depth inside of which compute ave Eady rate', &
's^-1', missing_value=missing_value, range=(/-10.0,1.e10/))
id_eady_rate_zave = -1
id_eady_rate_zave = register_diag_field ('ocean_model', 'eady_rate_zave',&
Grd%tracer_axes(1:2), Time%model_time, &
'Eady growth rate averaged over depth range for nphysics ', &
's^-1', missing_value=missing_value, range=(/-10.0,1.e10/))
id_eady_rate = -1
id_eady_rate = register_diag_field ('ocean_model', 'eady_rate',&
Grd%tracer_axes(1:3), Time%model_time, &
'Eady growth rate used in neutral physics ', 's^-1', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_rossby = -1
id_rossby = register_diag_field ('ocean_model', 'rossby', &
Grd%tracer_axes(1:2), Time%model_time, &
'Rossby radius used in neutral physics', 'm', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_rossby_radius = -1
id_rossby_radius = register_diag_field ('ocean_model', 'rossby_radius', &
Grd%tracer_axes(1:2), Time%model_time, &
'Rossby radius computed without min/max bounds', 'm', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_rossby_equator = -1
id_rossby_equator = register_diag_field ('ocean_model', 'rossby_equator', &
Grd%tracer_axes(1:2), Time%model_time, &
'Equatorial Rossby radius', 'm', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_rossby_nonequator = -1
id_rossby_nonequator = register_diag_field ('ocean_model', 'rossby_nonequator', &
Grd%tracer_axes(1:2), Time%model_time, &
'Rossby radius outside equatorial region', 'm', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_bczone = -1
id_bczone = register_diag_field ('ocean_model', 'bczone', &
Grd%tracer_axes(1:2), Time%model_time, &
'radius of baroclinic zone', 'm', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_sqrt2betaCr = -1
id_sqrt2betaCr = register_diag_field ('ocean_model', 'sqrt2betaCr', &
Grd%tracer_axes(1:2), Time%model_time, &
'sqrt(2 * 1st b/c wave speed *beta) in neutral physics ', '1/sec',&
missing_value=missing_value, range=(/-1e0,1e10/))
id_gw_speed = -1
id_gw_speed = register_diag_field ('ocean_model', 'gw_speed', &
Grd%tracer_axes(1:2), Time%model_time, &
'Gravity wave speed used in neutral physics ', 'm/s',&
missing_value=missing_value, range=(/-1e2,1e2/))
id_baroclinicity = -1
id_baroclinicity = register_diag_field ('ocean_model', 'baroclinicity', &
Grd%tracer_axes(1:2), Time%model_time, &
'vertically averaged horz density gradient', 'kg/m^4', &
missing_value=missing_value, range=(/-1e10,1.e10/))
id_growth_rate_baroclinic = -1
id_growth_rate_baroclinic = register_diag_field ('ocean_model', 'growth_rate_baroclinic', &
Grd%tracer_axes(1:3), Time%model_time, &
'growth rate using baroclinicity', 's^-1', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_agm_growth_rate = -1
id_agm_growth_rate = register_diag_field ('ocean_model', 'agm_growth_rate',&
Grd%tracer_axes(1:3), Time%model_time, &
'effective growth rate for agm', 's^-1', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_agm_length = -1
id_agm_length = register_diag_field ('ocean_model', 'agm_length', &
Grd%tracer_axes(1:3), Time%model_time, &
'effective length scale for agm', 'm', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_rhines_length = -1
id_rhines_length = register_diag_field ('ocean_model', 'rhines_length', &
Grd%tracer_axes(1:3), Time%model_time, &
'Rhines length approximated as eady_rate/beta', 'm', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_agm_qg = -1
id_agm_qg = register_diag_field ('ocean_model', 'agm_qg', &
Grd%tracer_axes(1:3), Time%model_time, &
'agm from QG theory', 'm^2/s', &
missing_value=missing_value, range=(/-10.0,1.e10/))
id_agm = -1
id_agm = register_diag_field ('ocean_model', 'agm', &
Grd%tracer_axes(1:2), Time%model_time, &
'GM diffusivity at surface', 'm^2/sec', &
missing_value=missing_value, range=(/-10.0,1.e10/),&
standard_name='ocean_tracer_bolus_laplacian_diffusivity')
id_agm_grid_scaling = register_diag_field ('ocean_model','agm_grid_scaling', &
Grd%tracer_axes(1:2), Time%model_time, &
'Scaling of AGM according to Delta(s)^2/(Delta(s)^2 + Rossby^2)',&
'dimensionless', missing_value=missing_value, range=(/-1e1,1e1/))
id_aredi_grid_scaling = register_diag_field ('ocean_model','aredi_grid_scaling', &
Grd%tracer_axes(1:2), Time%model_time, &
'Scaling of Aredi according to Delta(s)^2/(Delta(s)^2 + Rossby^2)',&
'dimensionless', missing_value=missing_value, range=(/-1e1,1e1/))
id_agm_3d = -1
id_agm_3d = register_diag_field ('ocean_model', 'agm_3d', &
Grd%tracer_axes(1:3), Time%model_time, &
'3d GM diffusivity', 'm^2/sec', &
missing_value=missing_value, range=(/-10.0,1.e10/),&
standard_name='ocean_tracer_bolus_laplacian_diffusivity')
id_aredi = -1
id_aredi = register_diag_field ('ocean_model', 'aredi', &
Grd%tracer_axes(1:2), Time%model_time, &
'neutral diffusivity at k=1', 'm^2/sec', &
missing_value=missing_value, range=(/-10.0,1.e20/),&
standard_name='ocean_tracer_epineutral_laplacian_diffusivity')
id_aredi_3d = -1
id_aredi_3d = register_diag_field ('ocean_model', 'aredi_3d', &
Grd%tracer_axes(1:3), Time%model_time, &
'3d neutral diffusivity', 'm^2/sec', &
missing_value=missing_value, range=(/-10.0,1.e20/))
id_tx_trans_nrho_gm = register_diag_field ('ocean_model','tx_trans_nrho_gm', Dens%neutralrho_axes_flux_x(1:3),&
Time%model_time, 'T-cell i-mass transport from GM on neutral rho',trim(transport_dims), &
missing_value=missing_value, range=(/-1e20,1e20/))
id_ty_trans_nrho_gm = register_diag_field ('ocean_model','ty_trans_nrho_gm', Dens%neutralrho_axes_flux_y(1:3),&
Time%model_time, 'T-cell j-mass transport from GM on neutral rho',trim(transport_dims), &
missing_value=missing_value, range=(/-1e20,1e20/))
id_tx_trans_rho_gm = register_diag_field ('ocean_model','tx_trans_rho_gm', Dens%potrho_axes_flux_x(1:3),&
Time%model_time, 'T-cell i-mass transport from GM on pot_rho',trim(transport_dims),&
missing_value=missing_value, range=(/-1e20,1e20/))
id_ty_trans_rho_gm = register_diag_field ('ocean_model','ty_trans_rho_gm', Dens%potrho_axes_flux_y(1:3),&
Time%model_time, 'T-cell j-mass transport from GM on pot_rho',trim(transport_dims),&
missing_value=missing_value, range=(/-1e20,1e20/))
id_tx_trans_theta_gm = register_diag_field ('ocean_model','tx_trans_theta_gm', Dens%theta_axes_flux_x(1:3),&
Time%model_time, 'T-cell i-mass transport from GM on theta',trim(transport_dims), &
missing_value=missing_value, range=(/-1e20,1e20/))
id_ty_trans_theta_gm = register_diag_field ('ocean_model','ty_trans_theta_gm', Dens%theta_axes_flux_y(1:3),&
Time%model_time, 'T-cell j-mass transport from GM on theta',trim(transport_dims), &
missing_value=missing_value, range=(/-1e20,1e20/))
! static fields
id_coriolis_param = register_static_field ('ocean_model', 'coriolis_param', Grd%tracer_axes(1:2), &
'Coriolis frequency on T-cell for nphysics', '1/s', &
missing_value=missing_value, range=(/-10.0,10.0/))
call diagnose_2d(Time, Grd, id_coriolis_param, coriolis_param(:,:))
id_beta_param = register_static_field ('ocean_model', 'beta_param', Grd%tracer_axes(1:2),&
'Beta=df/dy on T-cell for nphysics', '1/s', &
missing_value=missing_value, range=(/-10.0,10.0/))
call diagnose_2d(Time, Grd, id_beta_param, beta_param(:,:))
id_grid_length = register_static_field ('ocean_model', 'grid_length', Grd%tracer_axes(1:2), &
'Grid length scale used for nphysics calculations', 'm',&
missing_value=missing_value, range=(/-10.0,1e10/))
call diagnose_2d(Time, Grd, id_grid_length, grid_length(:,:))
! initialize cabbeling and thermobaricity related diagnostics
call cabbeling_thermob_init(Time, Dens)
! eta tendency calculations
id_eta_tend_gm90= -1
id_eta_tend_gm90= register_diag_field ('ocean_model','eta_tend_gm90', &
Grd%tracer_axes(1:2), Time%model_time, &
'analytic (and less accurate) form of non-Bouss steric sea level tendency from GM90', 'm/s',&
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_gm90 > 0) diagnose_eta_tend_gm90=.true.
id_eta_tend_gm90_glob= -1
id_eta_tend_gm90_glob= register_diag_field ('ocean_model', 'eta_tend_gm90_glob', &
Time%model_time, &
'global mean analytic (and less accurate) form of non-bouss steric sea level tendency from GM90',&
'm/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_gm90_glob > 0) diagnose_eta_tend_gm90=.true.
if(diagnose_eta_tend_gm90) then
call write_note(FILENAME, &
'diagnose_eta_tend_gm90==.true., diagnosing global sea level tendency from GM90.')
allocate (taper_fcn(isd:ied,jsd:jed,nk))
allocate (slopex(isd:ied,jsd:jed,nk))
allocate (agm_rho_slopex(isd:ied,jsd:jed,nk))
allocate (dslopex_dx(isd:ied,jsd:jed,nk))
allocate (dslopex_dz(isd:ied,jsd:jed,nk))
allocate (gradx_gamma_slopex(isd:ied,jsd:jed,nk))
allocate (slopey(isd:ied,jsd:jed,nk))
allocate (agm_rho_slopey(isd:ied,jsd:jed,nk))
allocate (dslopey_dy(isd:ied,jsd:jed,nk))
allocate (dslopey_dz(isd:ied,jsd:jed,nk))
allocate (grady_gamma_slopey(isd:ied,jsd:jed,nk))
taper_fcn(:,:,:) = 0.0
slopex(:,:,:) = 0.0
agm_rho_slopex(:,:,:) = 0.0
dslopex_dx(:,:,:) = 0.0
dslopex_dz(:,:,:) = 0.0
gradx_gamma_slopex(:,:,:) = 0.0
slopey(:,:,:) = 0.0
agm_rho_slopey(:,:,:) = 0.0
dslopey_dy(:,:,:) = 0.0
dslopey_dz(:,:,:) = 0.0
grady_gamma_slopey(:,:,:) = 0.0
endif
end subroutine ocean_nphysics_util_init
! NAME="ocean_nphysics_util_init"
!#######################################################################
!
!
!
! Initialize the cabbeling and thermobaricity related diagnostic
! fields and register the diagnostics to the diag manager.
!
!
subroutine cabbeling_thermob_init(Time, Dens)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type),intent(in) :: Dens
id_cabbeling_tend= -1
id_cabbeling_tend = register_diag_field ('ocean_model', 'cabbeling_tend',&
Grd%tracer_axes(1:3), Time%model_time, &
'(1/rho)*d(rho)/dt from cabbeling', 's^-1', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_cabbeling_tend > 0) diagnose_cabbeling_thermob =.true.
id_cabbeling_speed= -1
id_cabbeling_speed= register_diag_field ('ocean_model', 'cabbeling_speed',&
Grd%tracer_axes(1:3), Time%model_time, &
'(dz/rho)*d(rho)/dt from cabbeling', 'm/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_cabbeling_tend > 0) diagnose_cabbeling_thermob =.true.
id_cabbeling_tend_intz= -1
id_cabbeling_tend_intz= register_diag_field ('ocean_model', 'cabbeling_tend_intz',&
Grd%tracer_axes(1:2), Time%model_time, &
'vertical sum of (dz/rho)*d(rho)/dt from cabbeling', 'm/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_cabbeling_tend_intz > 0) diagnose_cabbeling_thermob =.true.
id_cabbeling_tend_intz_glob= -1
id_cabbeling_tend_intz_glob= register_diag_field ('ocean_model', &
'cabbeling_tend_intz_glob', Time%model_time, &
'global mean vertical sum of (dz/rho)*d(rho)/dt from cabbeling',&
'm/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_cabbeling_tend_intz_glob > 0) diagnose_cabbeling_thermob =.true.
id_thermobaric_tend= -1
id_thermobaric_tend= register_diag_field ('ocean_model', 'thermobaric_tend',&
Grd%tracer_axes(1:3), Time%model_time, &
'(1/rho)*d(rho)/dt from thermobaricity', 's^-1', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_thermobaric_tend > 0) diagnose_cabbeling_thermob =.true.
id_thermobaric_speed= -1
id_thermobaric_speed= register_diag_field ('ocean_model', 'thermobaric_speed',&
Grd%tracer_axes(1:3), Time%model_time, &
'(dz/rho)*d(rho)/dt from thermobaricity', 'm/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_thermobaric_speed > 0) diagnose_cabbeling_thermob =.true.
id_thermobaric_tend_intz= -1
id_thermobaric_tend_intz= register_diag_field ('ocean_model', 'thermobaric_tend_intz',&
Grd%tracer_axes(1:2), Time%model_time, &
'vertical sum of (dz/rho)*d(rho)/dt from thermobaricity', 'm/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_thermobaric_tend_intz > 0) diagnose_cabbeling_thermob =.true.
id_thermobaric_tend_intz_glob= -1
id_thermobaric_tend_intz_glob= register_diag_field ('ocean_model', &
'thermobaric_tend_intz_glob', Time%model_time, &
'global mean vertical sum of (dz/rho)*d(rho)/dt from thermobaricity',&
'm/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_thermobaric_tend_intz_glob > 0) diagnose_cabbeling_thermob =.true.
id_neut_rho_cabbeling= -1
id_neut_rho_cabbeling= register_diag_field ('ocean_model', 'neut_rho_cabbeling', &
Grd%tracer_axes(1:3), Time%model_time, &
'update of locally referenced potential density from cabbeling', '(kg/m^3)/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_neut_rho_cabbeling > 0) diagnose_cabbeling_thermob =.true.
id_wdian_cabbeling= -1
id_wdian_cabbeling= register_diag_field ('ocean_model', 'wdian_cabbeling',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to cabbeling', 'kg/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_wdian_cabbeling > 0) diagnose_cabbeling_thermob =.true.
id_neut_rho_thermob= -1
id_neut_rho_thermob= register_diag_field ('ocean_model', 'neut_rho_thermob', &
Grd%tracer_axes(1:3), Time%model_time, &
'update of locally referenced potential density from thermob', '(kg/m^3)/s',&
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_neut_rho_thermob > 0) diagnose_cabbeling_thermob =.true.
id_wdian_thermob= -1
id_wdian_thermob= register_diag_field ('ocean_model', 'wdian_thermob',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to thermobaricity', 'kg/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_wdian_thermob > 0) diagnose_cabbeling_thermob =.true.
id_neut_rho_cabbeling_on_nrho= -1
id_neut_rho_cabbeling_on_nrho= register_diag_field ('ocean_model', 'neut_rho_cabbeling_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'update of locally referenced potential density due to cabbeling as binned to neutral density',&
'(kg/m^3)/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_neut_rho_cabbeling_on_nrho > 0) diagnose_cabbeling_thermob =.true.
id_wdian_cabbeling_on_nrho= -1
id_wdian_cabbeling_on_nrho= register_diag_field ('ocean_model', 'wdian_cabbeling_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to cabbeling as binned to neutral density', 'kg/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_wdian_cabbeling_on_nrho > 0) diagnose_cabbeling_thermob =.true.
id_neut_rho_thermob_on_nrho= -1
id_neut_rho_thermob_on_nrho= register_diag_field ('ocean_model', 'neut_rho_thermob_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'update of locally referenced potential density due to thermobaricity as binned to neutral density',&
'(kg/m^3)/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_neut_rho_thermob_on_nrho > 0) diagnose_cabbeling_thermob =.true.
id_wdian_thermob_on_nrho= -1
id_wdian_thermob_on_nrho= register_diag_field ('ocean_model', 'wdian_thermob_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to thermobaricity as binned to neutral density', 'kg/s',&
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_wdian_thermob_on_nrho > 0) diagnose_cabbeling_thermob =.true.
id_tform_rho_cabbel_on_nrho =-1
id_tform_rho_cabbel_on_nrho = register_diag_field ('ocean_model', 'tform_rho_cabbel_on_nrho',&
Dens%neutralrho_axes(1:3), Time%model_time, &
'watermass transform due to cabbeling as binned to neutral rho layers', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_rho_cabbel_on_nrho > 0) diagnose_cabbeling_thermob =.true.
id_tform_rho_thermb_on_nrho =-1
id_tform_rho_thermb_on_nrho = register_diag_field ('ocean_model', 'tform_rho_thermb_on_nrho',&
Dens%neutralrho_axes(1:3), Time%model_time, &
'watermass transform due to thermobaricity as binned to neutral rho layers', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_rho_thermb_on_nrho > 0) diagnose_cabbeling_thermob =.true.
if(diagnose_cabbeling_thermob) then
call write_note(FILENAME, &
'diagnose_cabbeling_thermob==.true., diagnosing tendencies from cabbeling and thermobaricity.')
allocate (cabbeling_param(isd:ied,jsd:jed,nk))
allocate (thermobaric_param(isd:ied,jsd:jed,nk))
allocate (gradx_gamma_temp(isd:ied,jsd:jed,nk))
allocate (grady_gamma_temp(isd:ied,jsd:jed,nk))
allocate (gradx_gamma_press(isd:ied,jsd:jed,nk))
allocate (grady_gamma_press(isd:ied,jsd:jed,nk))
allocate (deriv_x_press(isd:ied,jsd:jed,nk))
allocate (deriv_y_press(isd:ied,jsd:jed,nk))
allocate (deriv_z_press(isd:ied,jsd:jed,nk))
cabbeling_param(:,:,:) = 0.0
thermobaric_param(:,:,:) = 0.0
gradx_gamma_temp(:,:,:) = 0.0
grady_gamma_temp(:,:,:) = 0.0
gradx_gamma_press(:,:,:) = 0.0
grady_gamma_press(:,:,:) = 0.0
deriv_x_press(:,:,:) = 0.0
deriv_y_press(:,:,:) = 0.0
deriv_z_press(:,:,:) = 0.0
endif
end subroutine cabbeling_thermob_init
! NAME="cabbeling_thermob_init"
!#######################################################################
!
!
!
! Initialize the diffusivities used in neutral physics.
! Need to initialize them after the ocean_nphysics_util_init routine,
! since need to have the domain parameters known for passing the
! array size information into the ocean_nphysics_coeff_init routine.
!
!
subroutine ocean_nphysics_coeff_init(Time, Thickness, rossby_radius, rossby_radius_raw, &
bczone_radius, agm_array, aredi_array, ah_array)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
real, dimension(isd:,jsd:), intent(inout) :: rossby_radius
real, dimension(isd:,jsd:), intent(inout) :: rossby_radius_raw
real, dimension(isd:,jsd:), intent(inout) :: bczone_radius
real, dimension(isd:,jsd:,:), intent(inout) :: agm_array
real, dimension(isd:,jsd:,:), intent(inout) :: aredi_array
real, dimension(isd:,jsd:), intent(inout) :: ah_array
integer :: i,j,k
integer :: i_delta, j_delta, k_delta
integer :: id_restart
real :: ft, deltaX, deltaY, delta, delta_iso, delta_iso0, delta_min, A_max, A_max0
character(len=64) :: file_name
integer :: stdoutunit
stdoutunit=stdout()
! initialization based on constant coefficients
agm_array(:,:,:) = agm*Grd%tmask(:,:,:)
agm_array_local(:,:,:) = agm*Grd%tmask(:,:,:)
aredi_array(:,:,:) = aredi*Grd%tmask(:,:,:)
aredi_array_local(:,:,:) = aredi*Grd%tmask(:,:,:)
! horizontal diffusivity in neutral boundary layer region
ah_array(:,:) = 0.0
if(neutral_horz_mix_bdy) then
call write_note(FILENAME,&
'Adding horizontal diffusivity in neutral boundary.')
call write_line('This method is implemented only for the case with neutral_physics_simple=.false.')
if(vel_micom_bdy > 0.0) then
ah_array(:,:) = vel_micom_bdy*grid_length(:,:)
else
ah_array(:,:) = ah_bdy
endif
do j=jsc,jec
write (stdoutunit,'(a,i4,a,e14.7,a)') &
' ah_array in neutral bdy layer at (isc,',j,',1)= ',ah_array(isc,j),' m^2/s'
enddo
endif
! Bryan-Lewis profile for aredi_array
! this profile is implemented for legacy purposes
if(bryan_lewis_aredi) then
call write_note(FILENAME,&
'Using Bryan-Lewis depth profile for the Redi diffusivity.')
call write_line('The Bryan-Lewis profile has been tested ONLY with agm==0.')
call write_line('Vertical dependence to the neutral diffusivity has not been implemented')
call write_line('according to a physical theory giving a vertical structure depending')
call write_line('on the flow field. More theoretical work is required.')
do k=1,nk
aredi_array(:,:,k) = (ahb + (ahs - ahb)*exp(-Grd%zt(k)/500.0))
enddo
do k=1,nk
write (stdoutunit,'(a,i3,e16.8)') ' k, diffusivity = ', k, aredi_array(isc,jsc,k)
enddo
endif
! Set diffusivity according to agm_lat_bands
if(agm_lat_bands) then
call write_note(FILENAME,&
'Setting agm_array according to agm_lat_bands.')
write(stdoutunit,'(a,e10.5)') &
' The ratio agm(south)/agm(north) is given by ',agm_lat_bands_boundary
write(stdoutunit,'(a,e10.5)') &
' with the latitude separating the bands given by ',agm_lat_bands_ratio
if(agm_lat_bands_boundary <= -90.0) then
write(stdoutunit,'(1x,a)') &
' Since agm_lat_bands_boundary <= -90, will default to globally constant agm value'
endif
agm_array_local(:,:,:) = agm*Grd%tmask(:,:,:)
if(agm_lat_bands_boundary > -90.) then
do j=jsc-Dom%yhalo,jec+Dom%yhalo
do i=isc-Dom%xhalo,iec+Dom%xhalo
if(Grd%yt(i,j) <= agm_lat_bands_boundary) then
agm_array_local(i,j,:) = agm*agm_lat_bands_ratio*Grd%tmask(i,j,:)
endif
enddo
enddo
endif
endif
! grid-scale dependent diffusivity suggested by that commonly used in MICOM
! vel_micom (m/s) sets velocity scale.
! space scale is set by grid size.
if(tracer_mix_micom) then
do k=1,nk
agm_array_local(:,:,k) = agm_micom(:,:)
enddo
do j=jsc,jec
write (stdoutunit,'(a,i4,a,e14.7,a)') &
' Micom agm_array at (isc,',j,',1) = ',agm_array_local(isc,j,1),' m^2/s'
enddo
if(aredi_equal_agm) then
aredi_array(:,:,:) = agm_array_local(:,:,:)
call write_note(FILENAME,&
'aredi_array = agm_array as given by Micom grid scale dependent diffusivity')
else
call write_note(FILENAME,&
'Since (agm/=aredi), aredi_array=aredi')
endif
endif
! to have the halos filled
agm_array(:,:,:) = agm_array_local(:,:,:)
! make mandatory=.false. to facilitate backward compatibility with older restart files.
file_name = 'ocean_neutral.res.nc'
id_restart = register_restart_field(Nphysics_util_restart, file_name, 'agm_array', agm_array, &
domain=Dom%domain2d, mandatory=.false.)
id_restart = register_restart_field(Nphysics_util_restart, file_name, 'aredi_array', aredi_array, &
domain=Dom%domain2d, mandatory=.false.)
id_restart = register_restart_field(Nphysics_util_restart, file_name, 'rossby_radius', rossby_radius, &
domain=Dom%domain2d, mandatory=.false.)
id_restart = register_restart_field(Nphysics_util_restart, file_name, 'rossby_radius_raw', rossby_radius_raw, &
domain=Dom%domain2d, mandatory=.false.)
id_restart = register_restart_field(Nphysics_util_restart, file_name, 'bczone_radius', bczone_radius, &
domain=Dom%domain2d, mandatory=.false.)
if(Time%init .and. nphysics_util_zero_init) then
call write_note(FILENAME,&
'Starting ocean_nphysics_util fields from raw initialization.')
elseif(.NOT. file_exist('INPUT/ocean_neutral.res.nc')) then
if (.NOT. Time%init) call mpp_error(FATAL,'Expecting file INPUT/ocean_neutral.res.nc to exist.&
&This file was not found and Time%init=.false.')
else
call restore_state(Nphysics_util_restart)
call mpp_update_domains(agm_array,Dom%domain2d)
write (stdoutunit,'(1x,a)') 'agm_array being read from restart.'
call write_chksum_3d('checksum start agm_array', agm_array(COMP,:)*Grd%tmask(COMP,:))
call mpp_update_domains(aredi_array,Dom%domain2d)
write (stdoutunit,'(1x,a)') 'aredi_array being read from restart.'
call write_chksum_3d('checksum start aredi_array', aredi_array(COMP,:)*Grd%tmask(COMP,:))
call mpp_update_domains(rossby_radius,Dom%domain2d)
write (stdoutunit,'(1x,a)') 'rossby_radius being read from restart.'
call write_chksum_2d('checksum start rossby_radius', rossby_radius(COMP)*Grd%tmask(COMP,1))
call mpp_update_domains(rossby_radius_raw,Dom%domain2d)
write (stdoutunit,'(1x,a)') 'rossby_radius_raw being read from restart.'
call write_chksum_2d('checksum start rossby_radius_raw', rossby_radius_raw(COMP)*Grd%tmask(COMP,1))
call mpp_update_domains(bczone_radius,Dom%domain2d)
write (stdoutunit,'(1x,a)') 'bczone_radius being read from restart.'
call write_chksum_2d('checksum start bczone_radius', bczone_radius(COMP)*Grd%tmask(COMP,1))
endif
! for those cases with agm_array determined by
! time invariant initialization methods.
if(.not. agm_closure) then
if(.not. agm_read_restart) then
call write_note(FILENAME,&
'agm_closure=.false. and agm_read_restart=.false.')
call write_line('=> agm_array set to static profiles.')
agm_array(:,:,:) = agm_array_local(:,:,:)
else
call write_note(FILENAME,&
'agm_closure=.false. and agm_read_restart=.true.')
call write_line('=> agm_array set to restart values.')
endif
endif
! allow aredi value to be read from namelist rather than from a restart
if(.not. agm_read_restart .and. .not. aredi_equal_agm) then
write(stdoutunit,'(1x,a)') &
'aredi_equal_agm=.false. and agm_read_restart=.false. => aredi_array set to static profiles.'
aredi_array(:,:,:) = aredi_array_local(:,:,:)
endif
if(aredi_equal_agm) then
call write_note(FILENAME,&
'aredi_equal_agm=.true. will force aredi_array to equal agm_array')
aredi_array(:,:,:) = agm_array(:,:,:)
else
call write_note(FILENAME,&
'aredi_equal_agm=.false. allows aredi_array to differ from agm_array')
endif
! compute maximum stable neutral slope available for neutral diffusion
i_delta = isc; j_delta = jsc; k_delta = 1
delta_iso = 1e30
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(Grd%tmask(i,j,k) /= 0.0) then
ft = 0.5/(aredi_array(i,j,k)*dtime + epsln)
deltaX = Grd%dxt(i,j)*Thickness%dzt(i,j,k)*ft
deltaY = Grd%dyt(i,j)*Thickness%dzt(i,j,k)*ft
if (delta_iso >= deltaX .or. delta_iso >= deltaY) then
i_delta = i; j_delta = j; k_delta = k
delta_iso = min(deltaX,deltaY)
endif
endif
enddo
enddo
enddo
delta_iso = delta_iso + 1.e-12*mpp_pe() ! to separate redundancies
delta_iso0 = delta_iso
call mpp_min (delta_iso)
! show most unstable location
if (delta_iso == delta_iso0) then
write(unit,'(a)')
write(unit,'(a)')'---Neutral direction slope check I for linear stability of neutral diffusion---'
write(unit,'(a,e14.7)')'With a neutral physics time step (secs) of ', dtime
write(unit,'(a)')'the most stringent linear stability constraint was found at the following ocean cell:'
write(unit,'(a,i4,a,i4,a,e14.7)')'long(',i_delta,',',j_delta,') = ',Grd%xt(i_delta,j_delta)
write(unit,'(a,i4,a,i4,a,e14.7)')'lat (',i_delta,',',j_delta,') = ',Grd%yt(i_delta,j_delta)
write(unit,'(a,i4,a,i4,a,i4,a,e14.7)')'thick(',i_delta,',',j_delta,',',k_delta,') = ', &
Thickness%dzt(i_delta,j_delta,k_delta)
write(unit,'(a,i4,a,i4,a,i4,a,e14.7)')'aredi(',i_delta,',',j_delta,',',k_delta,') = ', &
aredi_array(i_delta,j_delta,k_delta)
write(unit,'(a,e14.7,a)')'delta_iso = ',delta_iso,' is the maximum neutral direction slope'
write(unit,'(a)')'available for linear stability of the neutral diffusion scheme.'
write(unit,'(a)')'The namelist parameter smax should conservatively be <= delta_iso.'
if(smax >= delta_iso) then
write(unit,'(a,f10.5,a)')'==> Warning: The namelist parameter smax= ',smax, ' is >= to delta_iso.'
write(unit,'(a)')'Linear stability of the neutral diffusion scheme may be compromised.'
endif
write(unit,'(a)')
endif
! Compute maximum diffusivity available given a maximum slope of smax
i_delta = isc; j_delta = jsc; k_delta = 1
ft = 0.5/(smax*dtime + epsln)
delta_min = Thickness%dzt(i_delta,j_delta,k_delta)*Grd%dxt(i_delta,j_delta)
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(Grd%tmask(i,j,k) /= 0.0) then
delta = min(Grd%dxt(i,j),Grd%dyt(i,j))*Thickness%dzt(i,j,k)
if (delta_min > delta) then
i_delta = i; j_delta = j; k_delta = k
delta_min = delta
endif
endif
enddo
enddo
enddo
A_max = ft*delta_min + 1.e-12*mpp_pe() ! to separate redundancies
A_max0 = A_max
call mpp_min (A_max)
! show most unstable location
if (A_max == A_max0) then
write(unit,'(a)')
write(unit,'(a)')'---Neutral direction slope check II for linear stability of neutral diffusion---'
write(unit,'(a,e14.7)')'Assuming maximum Redi neutral diffusion slope of ', smax
write(unit,'(a,e14.7)')'and neutral physics time step (secs) of ', dtime
write(unit,'(a)')'the most stringent linear stability constraint was found at the following ocean cell:'
write(unit,'(a,i4,a,i4,a,e14.7)')'long(',i_delta,',',j_delta,') = ',&
Grd%xt(i_delta,j_delta)
write(unit,'(a,i4,a,i4,a,e14.7)')'lat (',i_delta,',',j_delta,') = ',&
Grd%yt(i_delta,j_delta)
write(unit,'(a,i4,a,i4,a,i4,a,e14.7)')'thick(',i_delta,',',j_delta,',',k_delta,')= ',&
Thickness%dzt(i_delta,j_delta,k_delta)
write(unit,'(a,e14.7,a)')'A_max = ',A_max,' (m^2/sec) is the maximum neutral diffusivity'
write(unit,'(a)')'available for linear stability of the neutral diffusion scheme.'
write(unit,'(a)')'Conservatively, neutral diffusivities used in the model should be less than A_max.'
write(unit,'(a)')'--------------------------------------------------------------------------------'
write(unit,'(a)')
endif
end subroutine ocean_nphysics_coeff_init
! NAME="ocean_nphysics_coeff_init"
!#######################################################################
!
!
!
! Compute the tracer derivatives.
!
! Horizontal derivatives are taken along surfaces of
! constant vertical coordinate (constant k-level)
!
! This approach ensures that when neutral physics defaults to "horizontal" physics
! next to boundaries, it will do so as horizontal, defined along surfaces of constant
! s-surfaces, and so will not generate spurious extrema.
!
! Additionally, when using generalized vertical coordinates, the neutral diffusion
! slope should be computed relative to the s-surfaces. The skew diffusion slope
! should ideally be computed with respect to z-surfaces, as z-surfaces define
! available potential energy. However, when s and z surfaces are reasonably close,
! as they are in the interior for zstar and pstar vertical coordinates, then we
! choose to to dissipate thickness as defined relative to the zstar or pstar surfaces.
! This should not be such a big deal, and it is certainly easier computationally than
! worrying about computing two separate sets of slopes. More on this detail is
! discussed in "Elements of MOM".
!
! NOTE: This approach is not appropriate for sigma-models. Indeed, many assumptions
! in the neutral physics modules need to be rethought for terrain following vertical
! coordinates.
!
! Vertical neutral density derivative for use in fz_terms
! and fz_flux, and for use in fx_flux and fy_flux.
! Note that the derivative at k=nk vanishes by definition
! since these derivatives are at the bottom of tracer cell.
! also note the use of -epsln_drhodz ensures the vertical
! derivative is always < 0. We also support the same
! approach used in the mom4p0d code for legacy purposes.
!
! Comments about smoothing drhodz:
!
! 1/ Tests in coupled 1-degree model showed extreme sensitivity
! of MOC to smoothing. GFDL users generally do NOT smooth, hence
! the default drhodz_smooth_vert=drhodz_smooth_horz=.false.
!
! 2/ Smoothing the vertical derivative of drhodzb and drhodzh helps
! is greatly needed for producing a regularized (i.e., well behaved)
! neutral slope vector.
!
! 3/ An attempt was made to smooth dTdz and dSdz rather
! than drhodz. The resulting slope was smooth, but not as
! smooth as when acting on drhodz itself.
!
!
!
subroutine tracer_derivs(Time, taum1, dtime, drhodT, drhodS, T_prog, Dens, dzwtr, &
dTdx, dTdy, dTdz, dSdx, dSdy, dSdz, drhodzb, drhodzh)
type(ocean_time_type), intent(in) :: Time
integer, intent(in) :: taum1
real, intent(in) :: dtime
real, dimension(isd:,jsd:,:), intent(in) :: drhodT
real, dimension(isd:,jsd:,:), intent(in) :: drhodS
real, dimension(isd:,jsd:,0:), intent(in) :: dzwtr
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
type(ocean_density_type), intent(in) :: Dens
type(tracer_3d_1_nk_type), intent(inout) :: dTdx(:)
type(tracer_3d_1_nk_type), intent(inout) :: dTdy(:)
type(tracer_3d_0_nk_type), intent(inout) :: dTdz(:)
type(tracer_3d_1_nk_type), intent(inout) :: dSdx
type(tracer_3d_1_nk_type), intent(inout) :: dSdy
type(tracer_3d_0_nk_type), intent(inout) :: dSdz
real, dimension(isd:,jsd:,:,0:), intent(inout) :: drhodzb
real, dimension(isd:,jsd:,:,0:), intent(inout) :: drhodzh
integer :: i, j, k, m, n
integer :: kp1, kbot
integer :: kr, kpkr, km1pkr
real :: dTdz_ijk1, dTdz_ijk2
real :: drhodzb0_prev, drhodzb1_prev
real :: drhodzh0_prev, drhodzh1_prev
real :: tmpb0, tmpb1, tmph0, tmph1
call mpp_clock_begin(id_clock_tracer_derivs)
! Horizontal derivatives taken along surfaces of
! constant vertical coordinate (constant k-level)
do n=1,num_prog_tracers
do k=1,nk
kp1 = min(k+1,nk)
wrk1(:,:,1) = T_prog(n)%field(:,:,k,taum1)
dTdx(n)%field(:,:,k) = FDX_T(wrk1(:,:,1))*FMX(Grd%tmask(:,:,k))
dTdy(n)%field(:,:,k) = FDY_T(wrk1(:,:,1))*FMY(Grd%tmask(:,:,k))
enddo
enddo
do k=1,nk
kp1 = min(k+1,nk)
wrk1(:,:,1) = Dens%rho_salinity(:,:,k,taum1)
dSdx%field(:,:,k) = FDX_T(wrk1(:,:,1))*FMX(Grd%tmask(:,:,k))
dSdy%field(:,:,k) = FDY_T(wrk1(:,:,1))*FMY(Grd%tmask(:,:,k))
enddo
! vertical derivative
do n=1,num_prog_tracers
do k=1,nk
kp1 = min(k+1,nk)
do j=jsd,jed
do i=isd,ied
dTdz(n)%field(i,j,k) = Grd%tmask(i,j,kp1)*dzwtr(i,j,k) &
*(T_prog(n)%field(i,j,k,taum1)-T_prog(n)%field(i,j,kp1,taum1))
enddo
enddo
enddo
enddo
do k=1,nk
kp1 = min(k+1,nk)
do j=jsd,jed
do i=isd,ied
dSdz%field(i,j,k) = Grd%tmask(i,j,kp1)*dzwtr(i,j,k) &
*(Dens%rho_salinity(i,j,k,taum1)-Dens%rho_salinity(i,j,kp1,taum1))
enddo
enddo
enddo
! vertical neutral density derivative
if(drhodz_mom4p1) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
dTdz_ijk1 = dTdz(index_temp)%field(i,j,k) ; dTdz_ijk2 = dSdz%field(i,j,k)
kr=0
kpkr = min(k+kr,nk)
drhodzb(i,j,k,kr) = drhodT(i,j,kpkr)*dTdz_ijk1 + drhodS(i,j,kpkr)*dTdz_ijk2
drhodzb(i,j,k,kr) = min(drhodzb(i,j,k,kr), -epsln_drhodz)
kr=1
kpkr = min(k+kr,nk)
drhodzb(i,j,k,kr) = drhodT(i,j,kpkr)*dTdz_ijk1 + drhodS(i,j,kpkr)*dTdz_ijk2
drhodzb(i,j,k,kr) = min(drhodzb(i,j,k,kr), -epsln_drhodz)
kr=0
km1pkr = max(k-1+kr,1)
drhodzh(i,j,k,kr) = drhodT(i,j,k)*dTdz(index_temp)%field(i,j,km1pkr) &
+ drhodS(i,j,k)*dSdz%field(i,j,km1pkr)
drhodzh(i,j,k,kr) = min(drhodzh(i,j,k,kr), -epsln_drhodz)
kr=1
km1pkr = max(k-1+kr,1)
drhodzh(i,j,k,kr) = drhodT(i,j,k)*dTdz(index_temp)%field(i,j,km1pkr) &
+ drhodS(i,j,k)*dSdz%field(i,j,km1pkr)
drhodzh(i,j,k,kr) = min(drhodzh(i,j,k,kr), -epsln_drhodz)
enddo
enddo
enddo
else
do k=1,nk
do j=jsd,jed
do i=isd,ied
dTdz_ijk1 = dTdz(index_temp)%field(i,j,k) ; dTdz_ijk2 = dSdz%field(i,j,k)
kr=0
kpkr = min(k+kr,nk)
drhodzb(i,j,k,kr) = drhodT(i,j,kpkr)*dTdz_ijk1 &
+ drhodS(i,j,kpkr)*dTdz_ijk2 -epsln
kr=1
kpkr = min(k+kr,nk)
drhodzb(i,j,k,kr) = drhodT(i,j,kpkr)*dTdz_ijk1 &
+ drhodS(i,j,kpkr)*dTdz_ijk2 -epsln
kr=0
km1pkr = max(k-1+kr,1)
drhodzh(i,j,k,kr) = drhodT(i,j,k)*dTdz(index_temp)%field(i,j,km1pkr) &
+ drhodS(i,j,k)*dSdz%field(i,j,km1pkr) -epsln
kr=1
km1pkr = max(k-1+kr,1)
drhodzh(i,j,k,kr) = drhodT(i,j,k)*dTdz(index_temp)%field(i,j,km1pkr) &
+ drhodS(i,j,k)*dSdz%field(i,j,km1pkr) -epsln
enddo
enddo
enddo
endif
! vertically smooth the vertical derivative of density to
! produce a smooth neutral slope vector for all flux components.
if(drhodz_smooth_vert) then
do m=1,num_121_passes
do j=jsd,jed
do i=isd,ied
drhodzb0_prev = onefourth*drhodzb(i,j,1,0)
drhodzb1_prev = onefourth*drhodzb(i,j,1,1)
drhodzh0_prev = onefourth*drhodzh(i,j,1,0)
drhodzh1_prev = onefourth*drhodzh(i,j,1,1)
kbot=Grd%kmt(i,j)
if(kbot>1) then
do k=2,kbot-2
tmpb0 = drhodzb(i,j,k,0)
drhodzb(i,j,k,0) = drhodzb0_prev + onehalf*drhodzb(i,j,k,0) &
+ onefourth*drhodzb(i,j,k+1,0)
drhodzb0_prev = onefourth*tmpb0
tmpb1 = drhodzb(i,j,k,1)
drhodzb(i,j,k,1) = drhodzb1_prev + onehalf*drhodzb(i,j,k,1) &
+ onefourth*drhodzb(i,j,k+1,1)
drhodzb1_prev = onefourth*tmpb1
tmph0 = drhodzh(i,j,k,0)
drhodzh(i,j,k,0) = drhodzh0_prev + onehalf*drhodzh(i,j,k,0) &
+ onefourth*drhodzh(i,j,k+1,0)
drhodzh0_prev = onefourth*tmph0
tmph1 = drhodzh(i,j,k,1)
drhodzh(i,j,k,1) = drhodzh1_prev + onehalf*drhodzh(i,j,k,1) &
+ onefourth*drhodzh(i,j,k+1,1)
drhodzh1_prev = onefourth*tmph1
enddo
endif
enddo
enddo
enddo
endif
! horizontally smooth the vertical derivative of density to
! produce a smooth neutral slope vector for all flux components.
! since drhodzb is needed only on the computational domain, there
! is no need to update its values in the halos.
if(drhodz_smooth_horz) then
do k=1,nk-1
drhodzb(:,:,k,0) = drhodzb(:,:,k,0) + dtime*LAP_T(drhodzb(:,:,k,0),smooth_lap(:,:))
drhodzb(:,:,k,1) = drhodzb(:,:,k,1) + dtime*LAP_T(drhodzb(:,:,k,1),smooth_lap(:,:))
drhodzh(:,:,k,0) = drhodzh(:,:,k,0) + dtime*LAP_T(drhodzh(:,:,k,0),smooth_lap(:,:))
drhodzh(:,:,k,1) = drhodzh(:,:,k,1) + dtime*LAP_T(drhodzh(:,:,k,1),smooth_lap(:,:))
enddo
call mpp_update_domains(drhodzh(:,:,:,0), Dom%domain2d, complete=.false.)
call mpp_update_domains(drhodzh(:,:,:,1), Dom%domain2d, complete=.true.)
endif
! compute squared buoyancy frequency
wrk2(:,:,:)=0.0
do k=1,nk
wrk2(:,:,k) = -onehalf*grav*rho0r*(drhodzb(:,:,k,0)+drhodzb(:,:,k,1))
enddo
call diagnose_3d(Time, Grd, id_N2slope, wrk2(:,:,:))
call mpp_clock_end(id_clock_tracer_derivs)
end subroutine tracer_derivs
! NAME="tracer_derivs"
!#######################################################################
!
!
!
! Subroutine computes the neutral slopes for the triads associated
! with the vertical flux component.
!
! Array tensor_31 initially holds the x-slope used for flux component fz.
! Array tensor_32 initially holds the y-slope used for flux component fz.
!
! In subsequent calculations, these arrays will be multipied by the
! diffusivities.
!
! No slope tapering is applied in this routine.
!
! slopes are computed over k=1,nk-1, since the slope at k=nk
! should be zero.
!
!
!
subroutine neutral_slopes(Time, dTdx, dTdy, dSdx, dSdy, drhodT, drhodS, drhodzb, tensor_31, tensor_32)
type(ocean_time_type), intent(in) :: Time
type(tracer_3d_1_nk_type), intent(in) :: dTdx(:)
type(tracer_3d_1_nk_type), intent(in) :: dTdy(:)
type(tracer_3d_1_nk_type), intent(in) :: dSdx
type(tracer_3d_1_nk_type), intent(in) :: dSdy
real, dimension(isd:,jsd:,:), intent(in) :: drhodT
real, dimension(isd:,jsd:,:), intent(in) :: drhodS
real, dimension(isd:,jsd:,:,0:), intent(in) :: drhodzb
real, dimension(isd:,jsd:,:,0:,0:), intent(inout) :: tensor_31
real, dimension(isd:,jsd:,:,0:,0:), intent(inout) :: tensor_32
integer :: i, j, k, kp1
real :: drhodT_ijk, drhodS_ijk, drhodT_ijkp1, drhodS_ijkp1
real :: tmask_ijkp1
real :: dTdx_ijk1, dTdx_ijk2, dTdx_im1jk1, dTdx_im1jk2
real :: dTdx_ijkp11, dTdx_ijkp12, dTdx_im1jkp11, dTdx_im1jkp12
real :: dTdy_ijk1, dTdy_ijk2, dTdy_ijm1k1, dTdy_ijm1k2
real :: dTdy_ijkp11, dTdy_ijkp12, dTdy_ijm1kp11, dTdy_ijm1kp12
real :: drhodzbr_ijk0, drhodzbr_ijk1
call mpp_clock_begin(id_clock_neutral_slopes)
tensor_31(:,:,:,:,:) = 0.0
tensor_32(:,:,:,:,:) = 0.0
do k=1,nk-1
kp1 = k+1
do j=jsc,jec
do i=isc,iec
tmask_ijkp1 = Grd%tmask(i,j,kp1)
drhodT_ijk = drhodT(i,j,k)
drhodS_ijk = drhodS(i,j,k)
drhodT_ijkp1 = drhodT(i,j,kp1)
drhodS_ijkp1 = drhodS(i,j,kp1)
dTdx_ijk1 = dTdx(index_temp)%field(i,j,k)
dTdx_ijk2 = dSdx%field(i,j,k)
dTdx_im1jk1 = dTdx(index_temp)%field(i-1,j,k)
dTdx_im1jk2 = dSdx%field(i-1,j,k)
dTdx_ijkp11 = dTdx(index_temp)%field(i,j,kp1)
dTdx_ijkp12 = dSdx%field(i,j,kp1)
dTdx_im1jkp11 = dTdx(index_temp)%field(i-1,j,kp1)
dTdx_im1jkp12 = dSdx%field(i-1,j,kp1)
dTdy_ijk1 = dTdy(index_temp)%field(i,j,k)
dTdy_ijk2 = dSdy%field(i,j,k)
dTdy_ijm1k1 = dTdy(index_temp)%field(i,j-1,k)
dTdy_ijm1k2 = dSdy%field(i,j-1,k)
dTdy_ijkp11 = dTdy(index_temp)%field(i,j,kp1)
dTdy_ijkp12 = dSdy%field(i,j,kp1)
dTdy_ijm1kp11 = dTdy(index_temp)%field(i,j-1,kp1)
dTdy_ijm1kp12 = dSdy%field(i,j-1,kp1)
drhodzbr_ijk0 = 1.0/drhodzb(i,j,k,0)
drhodzbr_ijk1 = 1.0/drhodzb(i,j,k,1)
! ip=jq=0
! kr=0
tensor_31(i,j,k,0,0) = -tmask_ijkp1 &
*(drhodT_ijk*dTdx_im1jk1 + drhodS_ijk*dTdx_im1jk2) &
*drhodzbr_ijk0
tensor_32(i,j,k,0,0) = -tmask_ijkp1 &
*(drhodT_ijk*dTdy_ijm1k1 + drhodS_ijk*dTdy_ijm1k2) &
*drhodzbr_ijk0
! kr=1
tensor_31(i,j,k,0,1) = -tmask_ijkp1 &
*(drhodT_ijkp1*dTdx_im1jkp11 + drhodS_ijkp1*dTdx_im1jkp12) &
*drhodzbr_ijk1
tensor_32(i,j,k,0,1) = -tmask_ijkp1 &
*(drhodT_ijkp1*dTdy_ijm1kp11 + drhodS_ijkp1*dTdy_ijm1kp12) &
*drhodzbr_ijk1
! ip=jq=1
! kr=0
tensor_31(i,j,k,1,0) = -tmask_ijkp1 &
*(drhodT_ijk*dTdx_ijk1 + drhodS_ijk*dTdx_ijk2) &
*drhodzbr_ijk0
tensor_32(i,j,k,1,0) = -tmask_ijkp1 &
*(drhodT_ijk*dTdy_ijk1 + drhodS_ijk*dTdy_ijk2) &
*drhodzbr_ijk0
! kr=1
tensor_31(i,j,k,1,1) = -tmask_ijkp1 &
*(drhodT_ijkp1*dTdx_ijkp11 + drhodS_ijkp1*dTdx_ijkp12) &
*drhodzbr_ijk1
tensor_32(i,j,k,1,1) = -tmask_ijkp1 &
*(drhodT_ijkp1*dTdy_ijkp11 + drhodS_ijkp1*dTdy_ijkp12) &
*drhodzbr_ijk1
enddo
enddo
enddo
! send to diagnostic manager
if (id_slope31 > 0) then
wrk1(isc:iec,jsc:jec,:) = onefourth* &
(tensor_31(isc:iec,jsc:jec,:,0,0) + tensor_31(isc:iec,jsc:jec,:,0,1) + &
tensor_31(isc:iec,jsc:jec,:,1,0) + tensor_31(isc:iec,jsc:jec,:,1,1))
call diagnose_3d(Time, Grd, id_slope31, wrk1(:,:,:))
endif
if (id_slope32 > 0) then
wrk1(isc:iec,jsc:jec,:) = onefourth* &
(tensor_32(isc:iec,jsc:jec,:,0,0) + tensor_32(isc:iec,jsc:jec,:,0,1) + &
tensor_32(isc:iec,jsc:jec,:,1,0) + tensor_32(isc:iec,jsc:jec,:,1,1))
call diagnose_3d(Time, Grd, id_slope32, wrk1(:,:,:))
endif
call mpp_clock_end(id_clock_neutral_slopes)
end subroutine neutral_slopes
! NAME="neutral_slopes"
!#######################################################################
!
!
!
! Finish computing eady growth rate.
!
!
subroutine compute_eady_rate(Time, Thickness, T_prog, Dens, eady_termx, eady_termy)
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(in) :: T_prog(:)
real, dimension(isd:,jsd:,:),intent(in) :: eady_termx
real, dimension(isd:,jsd:,:),intent(in) :: eady_termy
integer :: i, j, k, m, kbot, tau
real :: eady_rate_prev, tmp
real :: eady_mld(isd:ied,jsd:jed)
if(.not. agm_closure) return
call mpp_clock_begin(id_clock_compute_eady_rate)
tau = Time%tau
wrk1(:,:,:) = Grd%tmask(:,:,:)
eady_mld(:,:) = 0.0
! raw Eady growth rate, with regularization not yet applied
do k=1,nk
eady_rate(COMP,k) = Grd%tmask(COMP,k) &
*sqrt((eady_termx(COMP,k)*count_x(COMP)+eady_termy(COMP,k)*count_y(COMP))*gravrho0r)
enddo
! apply cap to the Eady growth rate
if(agm_closure_eady_cap) then
do k=1,nk
eady_rate(COMP,k) = 2.0*eady_rate(COMP,k)*agm_growth_rate_max(COMP) &
/(epsln+eady_rate(COMP,k)+agm_growth_rate_max(COMP))
enddo
endif
! vertically average Eady within surface mixed layer
if(agm_closure_eady_ave_mixed) then
call calc_mixed_layer_depth(Thickness, &
T_prog(index_salt)%field(isd:ied,jsd:jed,:,tau), &
Dens%rho_salinity(isd:ied,jsd:jed,:,tau), &
Dens%rho(isd:ied,jsd:jed,:,tau), &
Dens%pressure_at_depth(isd:ied,jsd:jed,:), &
eady_mld)
do j=jsc,jec
do i=isc,iec
eady_mld(i,j) = Grd%tmask(i,j,1)*min(eady_mld(i,j), Grd%ht(i,j))
enddo
enddo
! do not believe eady_rate at k=1, so always skip its contribution
wrk1_2d(:,:) = 0.0
wrk2_2d(:,:) = 0.0
do k=2,nk
do j=jsc,jec
do i=isc,iec
if(Thickness%depth_zt(i,j,k) < eady_mld(i,j)) then
wrk1_2d(i,j) = wrk1_2d(i,j) + Thickness%dzt(i,j,k)
wrk2_2d(i,j) = wrk2_2d(i,j) + eady_rate(i,j,k)*Thickness%dzt(i,j,k)
endif
enddo
enddo
enddo
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(Thickness%depth_zt(i,j,k) <= eady_mld(i,j)) then
eady_rate(i,j,k) = Grd%tmask(i,j,k)*wrk2_2d(i,j)/(wrk1_2d(i,j)+epsln)
endif
enddo
enddo
enddo
endif ! endif for agm_closure_eady_ave_mixed
! apply vertical 1-2-1 smoothing
if(agm_closure_eady_smooth_vert) then
do m=1,num_121_passes
do j=jsc,jec
do i=isc,iec
eady_rate_prev = onefourth*eady_rate(i,j,1)
kbot=Grd%kmt(i,j)
if(kbot>1) then
do k=2,kbot-2
tmp = eady_rate(i,j,k)
eady_rate(i,j,k) = eady_rate_prev + onehalf*eady_rate(i,j,k) &
+ onefourth*eady_rate(i,j,k+1)
eady_rate_prev = onefourth*tmp
enddo
endif
enddo
enddo
enddo
endif
! apply horizontal 1-2-1 smoothing
if(agm_closure_eady_smooth_horz) then
call mpp_update_domains (eady_rate(:,:,:), Dom%domain2d)
do k=1,nk
eady_rate(:,:,k) = S2D(eady_rate(:,:,k))
enddo
endif
! compute vertical average over specified depth
eady_rate_zave(:,:) = 0.0
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(Thickness%depth_zwt(i,j,k) >= agm_closure_upper_depth .and. &
Thickness%depth_zwt(i,j,k) <= agm_closure_lower_depth) then
eady_rate_zave(i,j) = eady_rate_zave(i,j) + eady_rate(i,j,k)*Thickness%dzwt(i,j,k)
wrk1_2d(i,j) = wrk1_2d(i,j) + Thickness%dzwt(i,j,k)
endif
enddo
enddo
enddo
eady_rate_zave(COMP) = Grd%tmask(COMP,1)*eady_rate_zave(COMP)/(wrk1_2d(COMP)+epsln)
call diagnose_3d(Time, Grd, id_eady_rate, eady_rate(:,:,:))
call diagnose_2d(Time, Grd, id_eady_rate_zave, eady_rate_zave(:,:))
call diagnose_2d(Time, Grd, id_eady_mld, eady_mld(:,:))
call mpp_clock_end(id_clock_compute_eady_rate)
end subroutine compute_eady_rate
! NAME="compute_eady_rate"
!#######################################################################
!
!
!
! Finish computing baroclinicity, which is defined to be the vertically
! averaged magnitude of the horizontal density gradient.
!
!
subroutine compute_baroclinicity(Time, baroclinic_termx, baroclinic_termy)
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:,jsd:), intent(in) :: baroclinic_termx
real, dimension(isd:,jsd:), intent(in) :: baroclinic_termy
real :: vertical_range
if(.not. agm_closure) return
call mpp_clock_begin(id_clock_compute_baroclinicity)
vertical_range = epsln + agm_closure_lower_depth - agm_closure_upper_depth
baroclinicity(COMP) = Grd%tmask(COMP,1)* &
(baroclinic_termx(COMP)*count_x(COMP)+baroclinic_termy(COMP)*count_y(COMP)) &
/vertical_range
call diagnose_2d(Time, Grd, id_baroclinicity, baroclinicity(:,:))
call mpp_clock_end(id_clock_compute_baroclinicity)
end subroutine compute_baroclinicity
! NAME="compute_baroclinicity"
!#######################################################################
!
!
!
! Subroutine computes the first baroclinic Rossby radius of deformation.
! Employ WKB approach described by Chelton et al. In particular,
! use formulae (2.2), (2.3a) and (2.3b) from their paper.
!
! Place a max and min value on the Rossby radius.
!
! Compute buoyancy frequency in terms of vertical gradient of
! locally referenced potential density. Place the reference point
! at the interface between the tracer cells, which is also where
! the vertical derivative of neutral density is located. This amounts
! to a centered difference computation similar to that used by
! Chelton et al. equation (B.4).
!
!
subroutine compute_rossby_radius(Thickness, dTdz, dSdz, Time, drhodT, drhodS, &
rossby_radius, rossby_radius_raw)
type(ocean_thickness_type), intent(in) :: Thickness
type(tracer_3d_0_nk_type), intent(in) :: dTdz(:)
type(tracer_3d_0_nk_type), intent(in) :: dSdz
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:,jsd:,:), intent(in) :: drhodT
real, dimension(isd:,jsd:,:), intent(in) :: drhodS
real, dimension(isd:,jsd:), intent(inout) :: rossby_radius
real, dimension(isd:,jsd:), intent(inout) :: rossby_radius_raw
real :: drhodzb_speed
integer :: i,j,k
call mpp_clock_begin(id_clock_compute_rossby_radius)
gravity_wave_speed(:,:) = 0.0
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
drhodzb_speed = 0.5*( (drhodT(i,j,k) + drhodT(i,j,k+1))*dTdz(index_temp)%field(i,j,k) &
+(drhodS(i,j,k) + drhodS(i,j,k+1))*dSdz%field(i,j,k) )
gravity_wave_speed(i,j) = gravity_wave_speed(i,j) &
+ Thickness%dzwt(i,j,k)*sqrt(abs(gravrho0r*drhodzb_speed))
enddo
enddo
enddo
gravity_wave_speed(:,:) = gravity_wave_speed(:,:)/pi
! for Eden near the equator
sqrt2betaCr(:,:) = sqrt(2.0*gravity_wave_speed(:,:)*beta_param(:,:))
wrk1_2d(:,:) = 0.0
wrk2_2d(:,:) = 0.0
rossby_radius_raw(:,:) = 0.0
wrk1_2d(COMP) = gravity_wave_speed(COMP)/(coriolis_param(COMP) + epsln)
wrk2_2d(COMP) = sqrt(gravity_wave_speed(COMP)/(2.0*beta_param(COMP)+epsln))
rossby_radius_raw(COMP) = min(wrk1_2d(COMP), wrk2_2d(COMP))
agm_grid_scaling(:,:) = 1.0
if(agm_closure_grid_scaling) then
! for backward bitwise compatibility
if(agm_closure_grid_scaling_power==2.0) then
agm_grid_scaling(COMP) = Grd%tmask(COMP,1) &
*grid_length(COMP)**2/(grid_length(COMP)**2 + rossby_radius_raw(COMP)**2)
else
agm_grid_scaling(COMP) = Grd%tmask(COMP,1) &
*grid_length(COMP)**agm_closure_grid_scaling_power &
/( grid_length(COMP)**agm_closure_grid_scaling_power &
+ rossby_radius_raw(COMP)**agm_closure_grid_scaling_power)
endif
endif
aredi_grid_scaling(:,:) = 1.0
if(aredi_diffusivity_grid_scaling) then
aredi_grid_scaling(COMP) = Grd%tmask(COMP,1) &
*grid_length(COMP)**agm_closure_grid_scaling_power &
/( grid_length(COMP)**agm_closure_grid_scaling_power &
+ rossby_radius_raw(COMP)**agm_closure_grid_scaling_power)
call mpp_update_domains (aredi_grid_scaling(:,:), Dom%domain2d)
endif
rossby_radius(COMP) = min(rossby_radius_max,max(rossby_radius_min,rossby_radius_raw(COMP)))
call diagnose_2d(Time, Grd, id_rossby_radius,rossby_radius_raw(:,:))
call diagnose_2d(Time, Grd, id_rossby, rossby_radius(:,:))
call diagnose_2d(Time, Grd, id_rossby_nonequator, wrk1_2d(:,:))
call diagnose_2d(Time, Grd, id_rossby_equator, wrk2_2d(:,:))
call diagnose_2d(Time, Grd, id_gw_speed, gravity_wave_speed(:,:))
call diagnose_2d(Time, Grd, id_sqrt2betaCr, sqrt2betaCr(:,:))
call diagnose_2d(Time, Grd, id_agm_grid_scaling, agm_grid_scaling(:,:))
call diagnose_2d(Time, Grd, id_aredi_grid_scaling, aredi_grid_scaling(:,:))
if(debug_this_module) then
call write_chksum_header(FILENAME, ' ', Time%model_time)
call write_chksum_2d('checksum rossby_radius', rossby_radius(COMP)*Grd%tmask(COMP,1))
call write_chksum_2d('checksum rossby_radius_raw', rossby_radius_raw(COMP)*Grd%tmask(COMP,1))
endif
call mpp_clock_end(id_clock_compute_rossby_radius)
end subroutine compute_rossby_radius
! NAME="compute_rossby_radius"
!#######################################################################
!
!
!
! Subroutine computes the radius of the baroclinic zone in a manner
! suggested by the Hadley Centre approach (Malcolm Roberts, personal
! communication).
!
! Algorithm is used in MOM3 and documented in the MOM3 Manual.
!
!
!
subroutine compute_bczone_radius(Time, bczone_radius)
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:,jsd:), intent(inout) :: bczone_radius
integer :: i,j
real :: n_zone, e_zone, s_zone, w_zone, fract, nstot, ewtot
integer :: ip, jq
if(.not. agm_closure) return
if(.not. agm_closure_length_bczone) return
call mpp_clock_begin(id_clock_compute_bczone_radius)
bczone_rate(COMP) = eady_rate_zave(COMP)
call mpp_update_domains (bczone_rate(:,:), BCzone_domain%domain2d)
do j=jsc,jec
do i=isc,iec
if (bczone_rate(i,j) > agm_closure_bczone_crit_rate) then
! search northward
n_zone = bczone_dyt(i,j)
do jq=j+1,j+bczone_max_pts
if (bczone_rate(i,jq) > agm_closure_bczone_crit_rate) then
n_zone = n_zone + bczone_dyt(i,jq)
else
exit
endif
enddo
! search southward
s_zone = bczone_dyt(i,j)
do jq=j-1,j-bczone_max_pts,-1
if (bczone_rate(i,jq) > agm_closure_bczone_crit_rate) then
s_zone = s_zone + bczone_dyt(i,jq)
else
exit
endif
enddo
! search eastward
e_zone = bczone_dxt(i,j)
do ip=i+1,i+bczone_max_pts
if (bczone_rate(ip,j) > agm_closure_bczone_crit_rate) then
e_zone = e_zone + bczone_dxt(ip,j)
else
exit
endif
enddo
! search westward
w_zone = bczone_dxt(i,j)
do ip=i-1,i-bczone_max_pts,-1
if (bczone_rate(ip,j) > agm_closure_bczone_crit_rate) then
w_zone = w_zone + bczone_dxt(ip,j)
else
exit
endif
enddo
! total radius (subtraction accounts for double-counting central point)
nstot=n_zone+s_zone-bczone_dyt(i,j)
ewtot=e_zone+w_zone-bczone_dxt(i,j)
if (nstot < ewtot) then
fract=min(n_zone,s_zone)/max(n_zone,s_zone)
bczone_radius(i,j)=fract*nstot
else
fract=min(e_zone,w_zone)/max(e_zone,w_zone)
bczone_radius(i,j)=fract*ewtot
endif
endif
enddo
enddo
call diagnose_2d(Time, Grd, id_bczone, bczone_radius(:,:))
if(debug_this_module) then
call write_chksum_header(FILENAME, ' ', Time%model_time)
call write_chksum_2d('checksum bczone_radius', bczone_radius(COMP)*Grd%tmask(COMP,1))
endif
call mpp_clock_end(id_clock_compute_bczone_radius)
end subroutine compute_bczone_radius
! NAME="compute_bczone_radius"
!#######################################################################
!
!
!
! Subroutine computes flow dependent diffusivity.
! Allow for an added dimensionless tuning factor as well as a
! minimum and maximum diffusivity.
!
!
subroutine compute_diffusivity(Time, ksurf_blayer, drhodz_zt, rossby_radius, &
bczone_radius, agm_array, aredi_array)
type(ocean_time_type), intent(in) :: Time
integer, dimension(isd:,jsd:),intent(in) :: ksurf_blayer
real, dimension(isd:,jsd:,:), intent(in) :: drhodz_zt
real, dimension(isd:,jsd:), intent(in) :: rossby_radius
real, dimension(isd:,jsd:), intent(in) :: bczone_radius
real, dimension(isd:,jsd:,:), intent(inout) :: agm_array
real, dimension(isd:,jsd:,:), intent(inout) :: aredi_array
integer :: i, j, k
real :: denom, param, active_cells, tmp
call mpp_clock_begin(id_clock_compute_diffusivity)
if(agm_closure) then
agm_growth_rate = 0.0
agm_length = 0.0
wrk1 = 0.0 !growth_rate
wrk2 = 0.0 !rhines length
wrk3 = 0.0 !agm_qg
wrk4 = 0.0 !agm_fast
wrk2_2d(:,:) = 0.0 !drhodz_ref
! set the length scale
if(agm_closure_length_fixed .or. agm_closure_baroclinic) then
agm_length(:,:,:) = Grd%tmask(:,:,:)*agm_closure_length
elseif(agm_closure_length_rossby) then
do k=1,nk
agm_length(COMP,k) = Grd%tmask(COMP,k)* &
2.0*grid_length(COMP)*rossby_radius(COMP)/(grid_length(COMP)+rossby_radius(COMP)+epsln)
enddo
elseif(agm_closure_length_bczone) then
do k=1,nk
agm_length(COMP,k) = Grd%tmask(COMP,k)* &
2.0*grid_length(COMP)*bczone_radius(COMP)/(grid_length(COMP)+bczone_radius(COMP)+epsln)
enddo
elseif(agm_closure_eden_greatbatch) then
if(agm_closure_eden_length_const) then
do k=1,nk
agm_length(COMP,k) = agm_closure_eden_length*Grd%tmask(COMP,k)
enddo
else
do k=1,nk
wrk2(COMP,k) = eady_rate(COMP,k)/(beta_param(COMP)+epsln)
agm_length(COMP,k) = min(wrk2(COMP,k),rossby_radius(COMP))
agm_length(COMP,k) = Grd%tmask(COMP,k)* &
2.0*grid_length(COMP)*agm_length(COMP,k)/(grid_length(COMP)+agm_length(COMP,k)+epsln)
enddo
endif
endif
if(agm_closure_length_cap) then
do k=1,nk
agm_length(COMP,k) = Grd%tmask(COMP,k)*min(agm_closure_length_max, agm_length(COMP,k))
enddo
endif
! set the growth rates
if(agm_closure_baroclinic) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
tmp = gravrho0r_buoyr*baroclinicity(i,j)
wrk1(i,j,k) = 2.0*tmp*agm_growth_rate_max(i,j) &
/(tmp+agm_growth_rate_max(i,j)+epsln)
enddo
enddo
enddo
elseif(agm_closure_eden_greatbatch) then
if(agm_closure_eden_gamma /= 0.0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
param = max(coriolis_param(i,j), sqrt2betaCr(i,j))
denom = sqrt(param**2 + agm_closure_eden_gamma*eady_rate(i,j,k)**2)+epsln
wrk1(i,j,k) = param*eady_rate(i,j,k)/denom
enddo
enddo
enddo
else
do k=1,nk
wrk1(COMP,k) = eady_rate(COMP,k)
enddo
endif
else
do k=1,nk
wrk1(COMP,k) = eady_rate_zave(COMP)
enddo
endif
! diffusivity computed as coeff*(N/Nref)**2
if(agm_closure_n2_scale) then
if(agm_closure_n2_scale_nref_cst) then
wrk2_2d(:,:) = rho0*grav_r*agm_closure_buoy_freq**2
else
! reference drhodz taken one level beneath blayer base,
! and no deeper than one cell from bottom.
do j=jsc,jec
do i=isc,iec
if(ksurf_blayer(i,j) > 0 .and. Grd%kmt(i,j) > 1) then
k = min(Grd%kmt(i,j)-1, ksurf_blayer(i,j)+1)
wrk2_2d(i,j) = abs(drhodz_zt(i,j,k))
endif
enddo
enddo
endif
do k=1,nk
agm_growth_rate(COMP,k) = Grd%tmask(COMP,k)*wrk1(COMP,k) ! computed just for diagnostics
wrk3(COMP,k) = agm_closure_n2_scale_coeff*Grd%tmask(COMP,k) &
*abs(drhodz_zt(COMP,k))/(epsln+wrk2_2d(COMP))
wrk4(COMP,k) = agm_grid_scaling(COMP)*(wrk3(COMP,k) + agm_micom(COMP))
wrk4(COMP,k) = Grd%tmask(COMP,k)*max(agm_closure_min, min(agm_closure_max, wrk4(COMP,k)))
enddo
! diffusivity computed as growth_rate*length_scale**2
else
do k=1,nk
agm_growth_rate(COMP,k) = Grd%tmask(COMP,k)*wrk1(COMP,k)
wrk3(COMP,k) = agm_closure_scaling*agm_growth_rate(COMP,k)*agm_length(COMP,k)**2
wrk4(COMP,k) = agm_grid_scaling(COMP)*(wrk3(COMP,k) + agm_micom(COMP))
wrk4(COMP,k) = Grd%tmask(COMP,k)*max(agm_closure_min, min(agm_closure_max, wrk4(COMP,k)))
enddo
endif
! time damping to get slowly evolving diffusivity
if(agm_smooth_time) then
do k=1,nk
agm_array(COMP,k) = &
Grd%tmask(COMP,k)*(agm_array(COMP,k) - gamma_damp*(agm_array(COMP,k)-wrk4(COMP,k)))
enddo
else
do k=1,nk
agm_array(COMP,k) = wrk4(COMP,k)
enddo
endif
! spatial smoothing
if(agm_smooth_space) then
wrk1_2d(:,:) = 0.0
call mpp_update_domains (agm_array(:,:,:), Dom%domain2d)
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(Grd%tmask(i,j,k)==1.0) then
active_cells = 4.0 +&
Grd%tmask(i-1,j,k) +&
Grd%tmask(i+1,j,k) +&
Grd%tmask(i,j-1,k) +&
Grd%tmask(i,j+1,k)
if (active_cells > 4.0) then
wrk1_2d(i,j) = &
(4.0*agm_array(i,j,k) +&
agm_array(i-1,j,k) +&
agm_array(i+1,j,k) +&
agm_array(i,j-1,k) +&
agm_array(i,j+1,k)) / active_cells
else
wrk1_2d(i,j) = agm_array(i,j,k)
endif
endif
enddo
enddo
agm_array(COMP,k) = wrk1_2d(COMP)*Grd%tmask(COMP,k)
enddo
endif
! need agm_array on full data domain
call mpp_update_domains (agm_array(:,:,:), Dom%domain2d)
call diagnose_3d(Time, Grd, id_agm_growth_rate, agm_growth_rate(:,:,:))
call diagnose_3d(Time, Grd, id_agm_length, agm_length(:,:,:))
call diagnose_3d(Time, Grd, id_rhines_length, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_agm_qg, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_growth_rate_baroclinic, wrk1(:,:,:))
if(id_N2_nblayer_base > 0) then
wrk3_2d(:,:) = grav*rho0r*wrk2_2d(:,:)
call diagnose_2d(Time, Grd, id_N2_nblayer_base, wrk2_2d(:,:))
endif
if(id_N2_for_agm > 0) then
wrk4(:,:,:) = 0.0
do k=1,nk
wrk4(:,:,k) = grav*rho0r*abs(drhodz_zt(:,:,k))
enddo
call diagnose_3d(Time, Grd, id_N2_for_agm, wrk4(:,:,:))
endif
endif ! endif of agm_closure
! set redi diffusivity
if(aredi_diffusivity_grid_scaling) then
do k=1,nk
aredi_array(:,:,k) = aredi*aredi_grid_scaling(:,:)*Grd%tmask(:,:,k)
enddo
endif
if(aredi_equal_agm) then
aredi_array(:,:,:) = agm_array(:,:,:)
endif
! diagnostics
call diagnose_2d(Time, Grd, id_agm, agm_array(:,:,1))
call diagnose_3d(Time, Grd, id_agm_3d, agm_array(:,:,:))
call diagnose_2d(Time, Grd, id_aredi, aredi_array(:,:,1))
call diagnose_3d(Time, Grd, id_aredi_3d, aredi_array(:,:,:))
if (id_ksurf_blayer > 0) then
wrk1_2d(:,:) = ksurf_blayer(:,:)
call diagnose_2d(Time, Grd, id_ksurf_blayer, wrk1_2d(:,:))
endif
if(debug_this_module) then
call write_chksum_header(FILENAME, ' ', Time%model_time)
call write_chksum_3d('checksum agm_array', agm_array(COMP,:)*Grd%tmask(COMP,:))
endif
call mpp_clock_end(id_clock_compute_diffusivity)
end subroutine compute_diffusivity
! NAME="compute_diffusivity"
!#######################################################################
!
!
!
! Classify horizontal GM mass transport according to neutral density classes.
!
! NOTE: This diagnostic works with transport integrated from bottom to
! a particular cell depth. To get transport_on_nrho_gm, a remapping is
! performed, rather than the binning done for trans_rho.
!
! Code history
! 2008: algorithm based (incorrectly) on transport_on_rho
! 2009: algorithm corrected to be consistent with remapping
! used in tracer_on_rho algorithm
!
!
subroutine transport_on_nrho_gm (Time, Dens, tx_trans_gm, ty_trans_gm)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: tx_trans_gm
real, dimension(isd:,jsd:,:), intent(in) :: ty_trans_gm
type(time_type) :: next_time
integer :: i, j, k, k_rho, neutralrho_nk
real :: work(isd:ied,jsd:jed,size(Dens%neutralrho_ref),2)
real :: W1, W2
real, dimension(jsc:jec) :: nrho_minj, nrho_maxj
real :: nrho_min, nrho_max
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error ocean_neutral_util_mod (transport_on_nrho_gm): module needs initialization ')
endif
next_time = increment_time(Time%model_time, int(dtime), 0)
if (need_data(id_tx_trans_nrho_gm,next_time) .or. need_data(id_ty_trans_nrho_gm,next_time)) then
neutralrho_nk = size(Dens%neutralrho_ref(:))
work(:,:,:,:) = 0.0
! for (i,j) points with neutralrho_ref < neutralrho(k=1), work=0
! for (i,j) points with neutralrho_ref > neutralrho(k=kmt), work=0
! these assumptions mean there is no need to specially handle the endpoints,
! since the initial value for work is 0.
! interpolate trans_gm from k-levels to neutralrho_nk-levels
do k=1 , nk-1
do j = jsc,jec
nrho_maxj(j) = maxval(Dens%neutralrho(isc:iec,j,k+1),mask=Grd%tmask(isc:iec,j,k+1)==1.)
nrho_minj(j) = minval(Dens%neutralrho(isc:iec,j,k),mask=Grd%tmask(isc:iec,j,k)==1.)
enddo
nrho_max = maxval(nrho_maxj)
nrho_min = minval(nrho_minj)
if (nrho_max == -huge(nrho_max)) exit ! only rock below this level
do k_rho = 1,neutralrho_nk
if (nrho_max < Dens%neutralrho_ref(k_rho)) cycle
if (nrho_min > Dens%neutralrho_ref(k_rho)) cycle
do j = jsc,jec
if (nrho_maxj(j) < Dens%neutralrho_ref(k_rho)) cycle
if (nrho_minj(j) > Dens%neutralrho_ref(k_rho)) cycle
do i = isc,iec
if ( Dens%neutralrho_ref(k_rho) > Dens%neutralrho(i,j,k) ) then
if (Dens%neutralrho_ref(k_rho) <= Dens%neutralrho(i,j,k+1)) then
W1= Dens%neutralrho_ref(k_rho)- Dens%neutralrho(i,j,k)
W2= Dens%neutralrho(i,j,k+1) - Dens%neutralrho_ref(k_rho)
work(i,j,k_rho,1) = (tx_trans_gm(i,j,k+1)*W1 +tx_trans_gm(i,j,k)*W2) &
/(W1 + W2 + epsln)
work(i,j,k_rho,2) = (ty_trans_gm(i,j,k+1)*W1 +ty_trans_gm(i,j,k)*W2) &
/(W1 + W2 + epsln)
endif
endif
enddo
enddo
enddo
enddo
do k_rho = 1,neutralrho_nk
work(COMP,k_rho,1) = work(COMP,k_rho,1)*Grd%tmask(COMP,1)
work(COMP,k_rho,2) = work(COMP,k_rho,2)*Grd%tmask(COMP,1)
enddo
if (id_tx_trans_nrho_gm > 0) then
call diagnose_3d(Time, Grd, id_tx_trans_nrho_gm, work(:,:,:,1), &
nk_lim=neutralrho_nk, use_mask=.false.)
endif
if (id_ty_trans_nrho_gm > 0) then
call diagnose_3d(Time, Grd, id_ty_trans_nrho_gm, work(:,:,:,2), &
nk_lim=neutralrho_nk, use_mask=.false.)
endif
endif
end subroutine transport_on_nrho_gm
! NAME="transport_on_nrho_gm"
!#######################################################################
!
!
!
! Classify horizontal GM mass transport according to potential density classes.
!
! Algorithm based on linear interpolation of function on s-surfaces to
! function on rho-surfaces.
!
! Diagnostic makes sense when potrho is monotonically increasing with
! depth, although the algorithm does not explicitly make this assumption.
!
! NOTE: This diagnostic works with transport integrated from bottom to
! a particular cell depth. To get transport_on_rho_gm, a remapping is
! performed, rather than the binning done for trans_rho.
!
! Code history
!
! 2008: algorithm based (incorrectly) on transport_on_rho
! 2009: algorithm corrected to be consistent with remapping
! used in tracer_on_rho algorithm
!
!
!
subroutine transport_on_rho_gm (Time, Dens, tx_trans_gm, ty_trans_gm)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: tx_trans_gm
real, dimension(isd:,jsd:,:), intent(in) :: ty_trans_gm
type(time_type) :: next_time
integer :: i, j, k, k_rho, potrho_nk
real :: work(isd:ied,jsd:jed,size(Dens%potrho_ref),2)
real :: W1, W2
real, dimension(jsc:jec) :: rho_minj, rho_maxj
real :: rho_min, rho_max
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error ocean_neutral_util_mod (transport_on_rho_gm): module needs initialization ')
endif
next_time = increment_time(Time%model_time, int(dtime), 0)
if (need_data(id_tx_trans_rho_gm,next_time) .or. need_data(id_ty_trans_rho_gm,next_time)) then
potrho_nk = size(Dens%potrho_ref(:))
work(:,:,:,:) = 0.0
! for (i,j) points with potrho_ref < potrho(k=1), work=0
! for (i,j) points with potrho_ref > potrho(k=kmt), work=0
! these assumptions mean there is no need to specially handle the endpoints,
! since the initial value for work is 0.
! interpolate trans_gm from k-levels to krho-levels
do k = 1,nk-1
do j = jsc,jec
rho_maxj(j) = maxval(Dens%potrho(isc:iec,j,k+1),mask=Grd%tmask(isc:iec,j,k+1)==1.)
rho_minj(j) = minval(Dens%potrho(isc:iec,j,k),mask=Grd%tmask(isc:iec,j,k)==1.)
enddo
rho_max = maxval(rho_maxj)
rho_min = minval(rho_minj)
if (rho_max == -huge(rho_max)) exit ! only rock below this level
do k_rho = 1,potrho_nk
if (rho_max < Dens%potrho_ref(k_rho)) cycle
if (rho_min > Dens%potrho_ref(k_rho)) cycle
do j = jsc,jec
if (rho_maxj(j) < Dens%potrho_ref(k_rho)) cycle
if (rho_minj(j) > Dens%potrho_ref(k_rho)) cycle
do i = isc,iec
if ( Dens%potrho_ref(k_rho) > Dens%potrho(i,j,k) ) then
if (Dens%potrho_ref(k_rho) <= Dens%potrho(i,j,k+1)) then
W1= Dens%potrho_ref(k_rho)- Dens%potrho(i,j,k)
W2= Dens%potrho(i,j,k+1) - Dens%potrho_ref(k_rho)
work(i,j,k_rho,1) = (tx_trans_gm(i,j,k+1)*W1 +tx_trans_gm(i,j,k)*W2) &
/(W1 + W2 + epsln)
work(i,j,k_rho,2) = (ty_trans_gm(i,j,k+1)*W1 +ty_trans_gm(i,j,k)*W2) &
/(W1 + W2 + epsln)
endif
endif
enddo
enddo
enddo
enddo
do k_rho = 1,potrho_nk
work(COMP,k_rho,1) = work(COMP,k_rho,1)*Grd%tmask(COMP,1)
work(COMP,k_rho,2) = work(COMP,k_rho,2)*Grd%tmask(COMP,1)
enddo
call diagnose_3d(Time, Grd, id_tx_trans_rho_gm, work(:,:,:,1), &
nk_lim=potrho_nk, use_mask=.false.)
call diagnose_3d(Time, Grd, id_ty_trans_rho_gm, work(:,:,:,2), &
nk_lim=potrho_nk, use_mask=.false.)
endif
end subroutine transport_on_rho_gm
! NAME="transport_on_rho_gm"
!#######################################################################
!
!
!
! Classify horizontal GM mass transport according to potential temp classes.
!
! Algorithm based on linear interpolation of function on s-surfaces to
! function on rho-surfaces.
!
! Diagnostic makes sense when potential temp is monotonically increasing
! with depth, although the algorithm does not explicitly make this assumption.
!
! NOTE: This diagnostic works with transport integrated from bottom to
! a particular cell depth. To get transport_on_theta_gm, a remapping is
! performed, rather than the binning done for trans_rho.
!
! Code history
!
! 2009: algorithm based on transport_on_rho_gm
!
!
!
subroutine transport_on_theta_gm (Time, Dens, Theta, tx_trans_gm, ty_trans_gm)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(in) :: Theta
real, dimension(isd:,jsd:,:), intent(in) :: tx_trans_gm
real, dimension(isd:,jsd:,:), intent(in) :: ty_trans_gm
type(time_type) :: next_time
integer :: i, j, k, k_theta, theta_nk, tau
real :: work(isd:ied,jsd:jed,size(Dens%potrho_ref),2)
real :: W1, W2
real, dimension(jsc:jec) :: theta_minj, theta_maxj
real :: theta_min, theta_max
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error ocean_neutral_util_mod (transport_on_theta_gm): module needs initialization ')
endif
tau = Time%tau
next_time = increment_time(Time%model_time, int(dtime), 0)
if (need_data(id_tx_trans_theta_gm,next_time) .or. need_data(id_ty_trans_theta_gm,next_time)) then
theta_nk = size(Dens%theta_ref(:))
work(:,:,:,:) = 0.0
! for (i,j) points with theta_ref > theta(k=1), work=0
! for (i,j) points with theta_ref < theta(k=kmt), work=0
! these assumptions mean there is no need to specially handle the endpoints,
! since the initial value for work is 0.
! interpolate trans_gm from k-levels to theta-levels
! note sign change in the if-tests relative to transport on rho and nrho
do k = 1,nk-1
do j = jsc,jec
theta_maxj(j) = maxval(Theta%field(isc:iec,j,k,tau),mask=Grd%tmask(isc:iec,j,k)==1.)
theta_minj(j) = minval(Theta%field(isc:iec,j,k+1,tau),mask=Grd%tmask(isc:iec,j,k+1)==1.)
enddo
theta_max = maxval(theta_maxj)
theta_min = minval(theta_minj)
if (theta_max == -huge(theta_max)) exit ! only rock below this level
do k_theta = 1,theta_nk
if (theta_max < Dens%theta_ref(k_theta)) cycle
if (theta_min > Dens%theta_ref(k_theta)) cycle
do j = jsc,jec
if (theta_maxj(j) < Dens%theta_ref(k_theta)) cycle
if (theta_minj(j) > Dens%theta_ref(k_theta)) cycle
do i = isc,iec
if ( Dens%theta_ref(k_theta) < Theta%field(i,j,k,tau) ) then
if (Dens%theta_ref(k_theta) >= Theta%field(i,j,k+1,tau)) then
W1= -Dens%theta_ref(k_theta) + Theta%field(i,j,k,tau)
W2= -Theta%field(i,j,k+1,tau) + Dens%theta_ref(k_theta)
work(i,j,k_theta,1) = (tx_trans_gm(i,j,k+1)*W1 +tx_trans_gm(i,j,k)*W2) &
/(W1 + W2 + epsln)
work(i,j,k_theta,2) = (ty_trans_gm(i,j,k+1)*W1 +ty_trans_gm(i,j,k)*W2) &
/(W1 + W2 + epsln)
endif
endif
enddo
enddo
enddo
enddo
do k_theta = 1,theta_nk
work(COMP,k_theta,1) = work(COMP,k_theta,1)*Grd%tmask(COMP,1)
work(COMP,k_theta,2) = work(COMP,k_theta,2)*Grd%tmask(COMP,1)
enddo
call diagnose_3d(Time, Grd, id_tx_trans_theta_gm, work(:,:,:,1), &
nk_lim=theta_nk, use_mask=.false.)
call diagnose_3d(Time, Grd, id_ty_trans_theta_gm, work(:,:,:,2), &
nk_lim=theta_nk, use_mask=.false.)
endif
end subroutine transport_on_theta_gm
! NAME="transport_on_theta_gm"
!#######################################################################
!
!
!
! Diagnose contribution to global mean sea level evolution arising
! from analytic form of Gent-McWilliams scheme.
!
! This routine computes the diagnostic based on an analytic form of the
! GM90 contribution. The raw numerical form is computed inside the respective
! nphysics modules. The raw numerical form is more accurate and thus
! recomended for purposes of sea level budgets.
!
! Compute an averaged slope using tensor_31 and tensor_32.
! Then compute the neutral divergence of this slope vector, just
! as if each component of the slope was a scalar.
!
! To avoid stencil issues with bottom cells, mask to zero
! contributions from cells next to the bottom in either of the
! three directions.
!
! Send output to diagnostic manager.
!
! Subroutine history:
! Feb2010 version 1.0: Stephen.Griffies
!
!
!
subroutine compute_eta_tend_gm90 (Time, Thickness, Dens, &
dtheta_dx, dtheta_dy, &
dsalt_dx , dsalt_dy , &
drhodzh, &
dtwedyt, dzwtr, delqc, &
tensor_31, tensor_32, agm_array)
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) :: dtheta_dx
real, dimension(isd:,jsd:,:), intent(in) :: dtheta_dy
real, dimension(isd:,jsd:,:), intent(in) :: dsalt_dx
real, dimension(isd:,jsd:,:), intent(in) :: dsalt_dy
real, dimension(isd:,jsd:,:,0:), intent(in) :: drhodzh
real, dimension(isd:,jsd:,0:), intent(in) :: dtwedyt
real, dimension(isd:,jsd:,0:), intent(in) :: dzwtr
real, dimension(isd:,jsd:,:,0:), intent(in) :: delqc
real, dimension(isd:,jsd:,:,0:,0:), intent(in) :: tensor_31
real, dimension(isd:,jsd:,:,0:,0:), intent(in) :: tensor_32
real, dimension(isd:,jsd:,:), intent(in) :: agm_array
integer :: i, j, k, kp1, tau
real :: extended_mask
real :: absslope, absslopex, absslopey
real :: drho_dz, drho_dz_rhor2, divergence
if(.not. diagnose_eta_tend_gm90) return
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error ocean_neutral_util_mod (compute_eta_tend_gm90): module needs initialization ')
endif
call mpp_clock_begin(id_clock_eta_tend_gm90)
tau = Time%tau
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
! approximate the neutral slope at the T-point as the average
! of the previously computed triad slopes contained in tensor_31 and tensor_32.
! taper to zero in regions of steep slope.
do k=1,nk
do j=jsc,jec
do i=isc,iec
slopex(i,j,k) = onefourth* &
(tensor_31(i,j,k,0,0) + tensor_31(i,j,k,0,1) + tensor_31(i,j,k,1,0) + tensor_31(i,j,k,1,1))
slopey(i,j,k) = onefourth* &
(tensor_32(i,j,k,0,0) + tensor_32(i,j,k,0,1) + tensor_32(i,j,k,1,0) + tensor_32(i,j,k,1,1))
taper_fcn(i,j,k) = Grd%tmask(i,j,k)
absslopex = abs(slopex(i,j,k))
absslopey = abs(slopey(i,j,k))
absslope = max(absslopex,absslopey)
if(absslope > smax_grad_gamma_scalar) then
taper_fcn(i,j,k) = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
slopex(i,j,k) = slopex(i,j,k)*taper_fcn(i,j,k)
slopey(i,j,k) = slopey(i,j,k)*taper_fcn(i,j,k)
endif
agm_rho_slopex(i,j,k) = slopex(i,j,k)*Dens%rho(i,j,k,tau)*agm_array(i,j,k)
agm_rho_slopey(i,j,k) = slopey(i,j,k)*Dens%rho(i,j,k,tau)*agm_array(i,j,k)
enddo
enddo
enddo
call mpp_update_domains(agm_rho_slopex(:,:,:), Dom%domain2d)
call mpp_update_domains(agm_rho_slopey(:,:,:), Dom%domain2d)
! compute derivatives of agm_rho_slope vector
do k=1,nk
kp1 = min(k+1,nk)
wrk1(:,:,1) = agm_rho_slopex(:,:,k)
dslopex_dx(:,:,k) = FDX_T(wrk1(:,:,1))*FMX(Grd%tmask(:,:,k))
wrk2(:,:,1) = agm_rho_slopey(:,:,k)
dslopey_dy(:,:,k) = FDY_T(wrk2(:,:,1))*FMY(Grd%tmask(:,:,k))
enddo
do k=1,nk
kp1 = min(k+1,nk)
dslopex_dz(:,:,k) = Grd%tmask(:,:,kp1)*dzwtr(:,:,k) &
*(agm_rho_slopex(:,:,k)-agm_rho_slopex(:,:,kp1))
dslopey_dz(:,:,k) = Grd%tmask(:,:,kp1)*dzwtr(:,:,k) &
*(agm_rho_slopey(:,:,k)-agm_rho_slopey(:,:,kp1))
enddo
! compute the i-gradient of agm_rho_slopex along a neutral direction
! units are rho*diffusivity/length
call calc_gradx_gamma_scalar(Thickness,Dens%drhodT(:,:,:),Dens%drhodS(:,:,:), &
dtheta_dx,dsalt_dx, &
drhodzh,dtwedyt,delqc, &
dslopex_dx,dslopex_dz,gradx_gamma_slopex)
! compute the j-gradient of agm_rho_slopey along a neutral direction
! units are rho*diffusivity/length
call calc_grady_gamma_scalar(Thickness,Dens%drhodT(:,:,:),Dens%drhodS(:,:,:), &
dtheta_dy,dsalt_dy, &
drhodzh,dtwedyt,delqc, &
dslopey_dy,dslopey_dz,grady_gamma_slopey)
! accumulate the vertical integral for the tendency
! unconcerned with diagnosing contributions where slope is steep,
! so taper the contributions to zero in those regions.
wrk1_2d(:,:) = 0.0
do k=1,nk
kp1=min(nk,k+1)
do j=jsc,jec
do i=isc,iec
extended_mask = Grd%tmask(i,j,k)* &
(Grd%tmask(i,j,kp1)+epsln)* &
(Grd%tmask(i+1,j,k)+epsln)* &
(Grd%tmask(i-1,j,k)+epsln)* &
(Grd%tmask(i,j+1,k)+epsln)* &
(Grd%tmask(i,j-1,k)+epsln)* &
(Grd%tmask(i+1,j,kp1)+epsln)* &
(Grd%tmask(i-1,j,kp1)+epsln)* &
(Grd%tmask(i,j+1,kp1)+epsln)* &
(Grd%tmask(i,j-1,kp1)+epsln)
drho_dz = abs(Dens%drhodz_zt(i,j,k))
drho_dz_rhor2 = extended_mask*drho_dz/(Dens%rho(i,j,k,tau)**2 + epsln)
divergence = gradx_gamma_slopex(i,j,k) + grady_gamma_slopey(i,j,k)
wrk1_2d(i,j) = wrk1_2d(i,j) + taper_fcn(i,j,k)*Thickness%dzt(i,j,k)*drho_dz_rhor2*divergence
enddo
enddo
enddo
if (id_eta_tend_gm90 > 0 .or. id_eta_tend_gm90_glob > 0) then
if(smooth_eta_tend_gm90) then
call mpp_update_domains (wrk1_2d(:,:), Dom%domain2d)
wrk1_2d(:,:) = S2D(wrk1_2d(:,:))
endif
if(id_eta_tend_gm90 > 0) then
call diagnose_2d(Time, Grd, id_eta_tend_gm90, -1.0*wrk1_2d(:,:))
endif
call diagnose_sum(Time, Grd, Dom, id_eta_tend_gm90_glob, wrk1_2d, -1.0*cellarea_r)
endif
call mpp_clock_end(id_clock_eta_tend_gm90)
end subroutine compute_eta_tend_gm90
! NAME="compute_eta_tend_gm90"
!#######################################################################
!
!
!
! Compute tendencies from cabbeling and thermobaricity.
!
! To avoid stencil issues with bottom cells, mask to zero
! contributions from cells next to the bottom in either of the
! three directions.
!
! Set cabbeling and thermobaricity to zero in regions where
! vertical stratification is gravitationally unstable. The idea is
! that in such regions, convective mixing will wash away the impact of
! along-neutral diffusion, in which case cabbeling and thermobaricity
! are not relevant anyhow.
!
! Set the vertical stratification drho_dz to a negative number
! with lower bound on its magnitude as set by epsln_drhodz_diagnostics.
!
! Send output to diagnostic manager.
!
! Subroutine history:
! Jan2010 version 1.0: initial coding by Stephen.Griffies
! Mar2010 version 2.0: tweaks on the division by drho_dz and strat_mask
!
!
!
subroutine cabbeling_thermob_tendency (Time, Thickness, T_prog, Dens, &
dtheta_dx, dtheta_dy, dtheta_dz, &
dsalt_dx, dsalt_dy, &
drhodzh, &
dxtdtsn, dtwedyt, dzwtr, delqc, &
aredi_array)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: dtheta_dx
real, dimension(isd:,jsd:,:), intent(in) :: dtheta_dy
real, dimension(isd:,jsd:,0:), intent(in) :: dtheta_dz
real, dimension(isd:,jsd:,:), intent(in) :: dsalt_dx
real, dimension(isd:,jsd:,:), intent(in) :: dsalt_dy
real, dimension(isd:,jsd:,:,0:), intent(in) :: drhodzh
real, dimension(isd:,jsd:,0:), intent(in) :: dxtdtsn
real, dimension(isd:,jsd:,0:), intent(in) :: dtwedyt
real, dimension(isd:,jsd:,0:), intent(in) :: dzwtr
real, dimension(isd:,jsd:,:,0:), intent(in) :: delqc
real, dimension(isd:,jsd:,:), intent(in) :: aredi_array
integer :: i, j, k, kp1, tau
real :: extended_mask
real :: drho_dz, strat_mask
real :: temporary
if(.not. diagnose_cabbeling_thermob) return
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error ocean_neutral_util_mod (cabbeling_thermob_tendency): module needs initialization ')
endif
call mpp_clock_begin(id_clock_cabbeling_thermob)
tau = Time%tau
! get the cabbeling and thermobaricity parameters from ocean_density.F90
call calc_cabbeling_thermobaricity(Dens%rho(:,:,:,tau), T_prog(index_salt)%field(:,:,:,tau), &
T_prog(index_temp)%field(:,:,:,tau), Dens%pressure_at_depth(:,:,:), &
Dens%drhodT(:,:,:), Dens%drhodS(:,:,:), Dens%drhodP(:,:,:), &
cabbeling_param, thermobaric_param)
! compute gradient of the pressure field
call pressure_derivs(Dens%pressure_at_depth(:,:,:), dzwtr)
! compute the i-gradient of temperature and pressure along a neutral direction
call calc_gradx_gamma_scalar(Thickness,Dens%drhodT(:,:,:),Dens%drhodS(:,:,:), &
dtheta_dx,dsalt_dx, &
drhodzh,dtwedyt,delqc, &
dtheta_dx,dtheta_dz,gradx_gamma_temp)
call calc_gradx_gamma_scalar(Thickness,Dens%drhodT(:,:,:),Dens%drhodS(:,:,:), &
dtheta_dx,dsalt_dx, &
drhodzh,dtwedyt,delqc, &
deriv_x_press,deriv_z_press,gradx_gamma_press)
! compute the j-gradient of temperature and pressure along a neutral direction
call calc_grady_gamma_scalar(Thickness,Dens%drhodT(:,:,:),Dens%drhodS(:,:,:), &
dtheta_dy,dsalt_dy, &
drhodzh,dxtdtsn,delqc, &
dtheta_dy,dtheta_dz,grady_gamma_temp)
call calc_grady_gamma_scalar(Thickness,Dens%drhodT(:,:,:),Dens%drhodS(:,:,:), &
dtheta_dy,dsalt_dy, &
drhodzh,dxtdtsn,delqc, &
deriv_y_press,deriv_z_press,grady_gamma_press)
wrk1(:,:,:) = 0.0 ! tendency from cabbeling (1/sec)
wrk2(:,:,:) = 0.0 ! tendency from thermobaricity (1/sec)
wrk3(:,:,:) = 0.0 ! speed from cabbeling (m/sec)
wrk4(:,:,:) = 0.0 ! speed from thermobaricity (m/sec)
wrk5(:,:,:) = 0.0 ! (drho/dt) from cabbeling (kg/m^3)/s
wrk6(:,:,:) = 0.0 ! (drho/dt) from thermobaricity (kg/m^3)/s
wrk1_2d(:,:) = 0.0 ! vertical sum of cabbeling speed (m/s)
wrk2_2d(:,:) = 0.0 ! vertical sum of thermobaricity speed (m/s)
wrk1_v(:,:,:,:) = 0.0 ! dianeutral mass transport (kg/s) from (cabbeling,thermobaricity)
do k=1,nk
kp1=min(nk,k+1)
do j=jsc,jec
do i=isc,iec
extended_mask = Grd%tmask(i,j,k)* &
(Grd%tmask(i,j,kp1)+epsln)* &
(Grd%tmask(i+1,j,k)+epsln)* &
(Grd%tmask(i-1,j,k)+epsln)* &
(Grd%tmask(i,j+1,k)+epsln)* &
(Grd%tmask(i,j-1,k)+epsln)* &
(Grd%tmask(i+1,j,kp1)+epsln)* &
(Grd%tmask(i-1,j,kp1)+epsln)* &
(Grd%tmask(i,j+1,kp1)+epsln)* &
(Grd%tmask(i,j-1,kp1)+epsln)
! zero those points where stratification is gravitationally unstable
strat_mask = 1.0
if(Dens%drhodz_zt(i,j,k) > 0.0) strat_mask=0.0
wrk1(i,j,k) = extended_mask*strat_mask*aredi_array(i,j,k)*cabbeling_param(i,j,k) &
*(gradx_gamma_temp(i,j,k)*gradx_gamma_temp(i,j,k)+grady_gamma_temp(i,j,k)*grady_gamma_temp(i,j,k))
wrk2(i,j,k) = extended_mask*strat_mask*aredi_array(i,j,k)*thermobaric_param(i,j,k) &
*(gradx_gamma_temp(i,j,k)*gradx_gamma_press(i,j,k)+grady_gamma_temp(i,j,k)*grady_gamma_press(i,j,k))
! provide lower bound on abs(drho_dz), and ensure drho_dz < 0
drho_dz = min(-epsln_drhodz_diagnostics,Dens%drhodz_zt(i,j,k))
temporary = Grd%tmask(i,j,k)*Grd%dat(i,j)*Dens%rho(i,j,k,tau) &
*Thickness%rho_dzt(i,j,k,tau)/(epsln+drho_dz*Thickness%dzt(i,j,k))
wrk1_v(i,j,k,1) = wrk1(i,j,k)*temporary
wrk1_v(i,j,k,2) = wrk2(i,j,k)*temporary
wrk3(i,j,k) = wrk1(i,j,k)*Thickness%dzt(i,j,k)
wrk4(i,j,k) = wrk2(i,j,k)*Thickness%dzt(i,j,k)
wrk5(i,j,k) = wrk1(i,j,k)*Dens%rho(i,j,k,tau)
wrk6(i,j,k) = wrk2(i,j,k)*Dens%rho(i,j,k,tau)
wrk1_2d(i,j) = wrk1_2d(i,j) + wrk3(i,j,k)
wrk2_2d(i,j) = wrk2_2d(i,j) + wrk4(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_cabbeling_tend, wrk1(:,:,:))
call diagnose_3d(Time, Grd, id_thermobaric_tend, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_cabbeling_speed, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_thermobaric_speed, wrk4(:,:,:))
call diagnose_2d(Time, Grd, id_cabbeling_tend_intz, wrk1_2d(:,:))
call diagnose_2d(Time, Grd, id_thermobaric_tend_intz, wrk2_2d(:,:))
call diagnose_sum(Time, Grd, Dom, id_cabbeling_tend_intz_glob, wrk1_2d, cellarea_r)
call diagnose_sum(Time, Grd, Dom, id_thermobaric_tend_intz_glob, wrk2_2d, cellarea_r)
if(wdianeutral_smooth) then
if(id_wdian_cabbeling > 0 .or. id_wdian_cabbeling_on_nrho > 0 .or. &
id_wdian_thermob > 0 .or. id_wdian_thermob_on_nrho > 0) then
call mpp_update_domains (wrk1_v(:,:,:,1), Dom%domain2d)
call mpp_update_domains (wrk1_v(:,:,:,2), Dom%domain2d)
do k=1,nk
wrk1_v(:,:,k,1) = S2D(wrk1_v(:,:,k,1))
wrk1_v(:,:,k,2) = S2D(wrk1_v(:,:,k,2))
enddo
endif
endif
call diagnose_3d(Time, Grd, id_neut_rho_cabbeling, wrk5(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_cabbeling_on_nrho, wrk5)
call diagnose_3d(Time, Grd, id_wdian_cabbeling, wrk1_v(:,:,:,1))
call diagnose_3d_rho(Time, Dens, id_wdian_cabbeling_on_nrho, wrk1_v(:,:,:,1))
call diagnose_3d(Time, Grd, id_neut_rho_thermob, wrk6(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_thermob_on_nrho, wrk6)
call diagnose_3d(Time, Grd, id_wdian_thermob, wrk1_v(:,:,:,2))
call diagnose_3d_rho(Time, Dens, id_wdian_thermob_on_nrho, wrk1_v(:,:,:,2))
if (id_tform_rho_cabbel_on_nrho > 0) then
wrk3(:,:,:) = 0.0
do k=1,nk
wrk3(COMP,k) = wrk1(COMP,k)*Dens%watermass_factor(COMP,k) &
*Dens%rho(COMP,k,tau)*Thickness%rho_dzt(COMP,k,tau)*Grd%dat(COMP)
enddo
call diagnose_3d_rho(Time, Dens, id_tform_rho_cabbel_on_nrho, wrk3)
endif
if (id_tform_rho_thermb_on_nrho > 0) then
wrk3(:,:,:) = 0.0
do k=1,nk
wrk3(COMP,k) = wrk2(COMP,k)*Dens%watermass_factor(COMP,k) &
*Dens%rho(COMP,k,tau)*Thickness%rho_dzt(COMP,k,tau)*Grd%dat(COMP)
enddo
call diagnose_3d_rho(Time, Dens, id_tform_rho_thermb_on_nrho, wrk3)
endif
call mpp_clock_end(id_clock_cabbeling_thermob)
end subroutine cabbeling_thermob_tendency
! NAME="cabbeling_thermob_tendency"
!#######################################################################
!
!
!
! Compute the pressure derivatives for use in thermobaricity diagnostic.
!
!
subroutine pressure_derivs(pressure, dzwtr)
real, dimension(isd:,jsd:,:), intent(in) :: pressure
real, dimension(isd:,jsd:,0:), intent(in) :: dzwtr
integer :: k
integer :: kp1
do k=1,nk
kp1 = min(k+1,nk)
wrk1(:,:,1) = pressure(:,:,k)
deriv_x_press(:,:,k) = FDX_T(wrk1(:,:,1))*FMX(Grd%tmask(:,:,k))
deriv_y_press(:,:,k) = FDY_T(wrk1(:,:,1))*FMY(Grd%tmask(:,:,k))
enddo
do k=1,nk
kp1 = min(k+1,nk)
deriv_z_press(:,:,k) = Grd%tmask(:,:,kp1)*dzwtr(:,:,k)*(pressure(:,:,k) - pressure(:,:,kp1))
enddo
end subroutine pressure_derivs
! NAME="pressure_derivs"
!#######################################################################
!
!
!
! Subroutine computes the i-gradient along neutral surfaces
! for a scalar field. For use in the cabbeling and thermobaricity
! diagnostic. Thus, when slope is steep, the full derivative is
! tapered, rather than just the off-diagonal term.
!
! We only compute neutral i-gradient for interior regions,
! since will be setting thermobaricity and cabbeling to zero at
! top and bottom levels; do not trust the calculation at top
! and bottom boundaries due to stencil truncations.
!
! Algorithm slightly modified from ocean_nphysC calculation of
! neutral diffusion x-flux.
!
!
!
subroutine calc_gradx_gamma_scalar(Thickness,drhodT,drhodS, &
dtheta_dx,dsalt_dx, &
drhodzh,dtwedyt,delqc, &
dscalar_dx,dscalar_dz,gradx_gamma_scalar)
type(ocean_thickness_type), intent(in) :: Thickness
real, dimension(isd:,jsd:,:), intent(in) :: drhodT
real, dimension(isd:,jsd:,:), intent(in) :: drhodS
real, dimension(isd:,jsd:,:), intent(in) :: dtheta_dx
real, dimension(isd:,jsd:,:), intent(in) :: dsalt_dx
real, dimension(isd:,jsd:,:,0:), intent(in) :: drhodzh
real, dimension(isd:,jsd:,0:), intent(in) :: dtwedyt
real, dimension(isd:,jsd:,:,0:), intent(in) :: delqc
real, dimension(isd:,jsd:,:), intent(in) :: dscalar_dx
real, dimension(isd:,jsd:,0:), intent(in) :: dscalar_dz
real, dimension(isd:,jsd:,:), intent(inout) :: gradx_gamma_scalar
integer :: i, j, k, ip
integer :: kr, kpkr
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
real :: slope, absslope
real :: drhodx, drhodz
real :: cell_mass_inv
gradx_gamma_scalar(:,:,:) = 0.0
do k=2,nk-1
! initialize arrays to zero
tensor_11(:,:,:,:) = 0.0
tensor_13(:,:,:,:) = 0.0
! scalar field 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)*dtheta_dx(i,j,k) &
+drhodS(i+ip,j,k)*dsalt_dx(i,j,k))
slope = -drhodx/drhodz
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax_grad_gamma_scalar) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
else
taperA = 1.0
endif
! fill tensor components; note tapering applied to both terms
tensor_11(i,j,ip,kr) = taperA
tensor_13(i,j,ip,kr) = taperA*slope
enddo
enddo
enddo
enddo
! scalar field dependent part of the calculation
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)*dscalar_dx(i,j,k) + tensor_13(i,j,ip,kr)*dscalar_dz(i+ip,j,k-1+kr)) &
*min(delqc(i,j,k,kr),delqc(i+1,j,k,kr))
enddo
enddo
enddo
enddo
! need to divide by mass of cell to get to a scalar gradient
do j=jsc,jec
do i=isc-1,iec
cell_mass_inv = Grd%tmask(i,j,k)*Thickness%rho_dztr(i,j,k)*Grd%dater(i,j)
gradx_gamma_scalar(i,j,k) = (sumz(i,j,0)+sumz(i,j,1))*cell_mass_inv
enddo
enddo
enddo ! enddo for k-loop
end subroutine calc_gradx_gamma_scalar
! NAME="calc_gradx_gamma_scalar"
!#######################################################################
!
!
!
! Subroutine computes the y-gradient along neutral surfaces
! for a scalar field. For use in the cabbeling and thermobaricity
! diagnostic. Thus, when slope is steep, the full derivative is
! tapered, rather than just the off-diagonal term.
!
! We only compute neutral y-gradient for interior regions,
! since will be setting thermobaricity and cabbeling to zero at
! top and bottom levels; do not trust the calculation at top
! and bottom boundaries due to stencil truncations.
!
! Algorithm slightly modified from ocean_nphysC calculation of
! neutral diffusion y-flux.
!
!
!
subroutine calc_grady_gamma_scalar(Thickness,drhodT,drhodS, &
dtheta_dy,dsalt_dy, &
drhodzh,dxtdtsn,delqc, &
dscalar_dy,dscalar_dz,grady_gamma_scalar)
type(ocean_thickness_type), intent(in) :: Thickness
real, dimension(isd:,jsd:,:), intent(in) :: drhodT
real, dimension(isd:,jsd:,:), intent(in) :: drhodS
real, dimension(isd:,jsd:,:), intent(in) :: dtheta_dy
real, dimension(isd:,jsd:,:), intent(in) :: dsalt_dy
real, dimension(isd:,jsd:,:,0:), intent(in) :: drhodzh
real, dimension(isd:,jsd:,0:), intent(in) :: dxtdtsn
real, dimension(isd:,jsd:,:,0:), intent(in) :: delqc
real, dimension(isd:,jsd:,:), intent(in) :: dscalar_dy
real, dimension(isd:,jsd:,0:), intent(in) :: dscalar_dz
real, dimension(isd:,jsd:,:), intent(inout) :: grady_gamma_scalar
integer :: i, j, k, jq
integer :: kr, kpkr
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
real :: slope, absslope
real :: drhody, drhodz
real :: cell_mass_inv
grady_gamma_scalar(:,:,:) = 0.0
do k=2,nk-1
! initialize arrays to zero
tensor_22(:,:,:,:) = 0.0
tensor_23(:,:,:,:) = 0.0
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)*dtheta_dy(i,j,k) &
+drhodS(i,j+jq,k)*dsalt_dy(i,j,k))
slope = -drhody/drhodz
absslope = abs(slope)
! taper for steep slope regions
if(absslope > smax_grad_gamma_scalar) then
taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))
else
taperA = 1.0
endif
! fill tensor components; note tapering applied to both terms
tensor_22(i,j,jq,kr) = taperA
tensor_23(i,j,jq,kr) = taperA*slope
enddo
enddo
enddo
enddo
! scalar field dependent part of the calculation
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)*dscalar_dy(i,j,k) + tensor_23(i,j,jq,kr)*dscalar_dz(i,j+jq,k-1+kr)) &
* min(delqc(i,j,k,kr),delqc(i,j+1,k,kr))
enddo
enddo
enddo
enddo
! need to divide by mass of cell to get to a scalar gradient
do j=jsc,jec
do i=isc-1,iec
cell_mass_inv = Grd%tmask(i,j,k)*Thickness%rho_dztr(i,j,k)*Grd%datnr(i,j)
grady_gamma_scalar(i,j,k) = (sumz(i,j,0)+sumz(i,j,1))*cell_mass_inv
enddo
enddo
enddo ! enddo for k-loop
end subroutine calc_grady_gamma_scalar
! NAME="calc_grady_gamma_scalar"
!#######################################################################
!
!
!
! Initialization of watermass diagnostic output files.
!
!
subroutine watermass_diag_init(Time, Dens, diagnose_gm_redi)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
logical , intent(inout) :: diagnose_gm_redi
integer :: stdoutunit
stdoutunit=stdout()
compute_watermass_diag = .false.
id_neut_rho_nphysics = -1
id_neut_rho_nphysics = register_diag_field ('ocean_model', 'neut_rho_nphysics', &
Grd%tracer_axes(1:3), Time%model_time, &
'update of locally referenced potential density from time-explicit nphysics',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_rho_nphysics > 0) compute_watermass_diag=.true.
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 mass transport due to time-explicit nphysics', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_rho_nphysics > 0) compute_watermass_diag=.true.
! full contributions to neutral diffusion
id_neut_rho_ndiff = -1
id_neut_rho_ndiff = register_diag_field ('ocean_model', 'neut_rho_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'update of locally referenced potential density time-explicit neutral diffusion',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_rho_ndiff > 0) compute_watermass_diag=.true.
id_neut_rho_ndiff_on_nrho = -1
id_neut_rho_ndiff_on_nrho = register_diag_field ('ocean_model', 'neut_rho_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'update of locally ref potrho from time-explicit neutral diffusion as binned to neutral rho layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_rho_ndiff_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_ndiff = -1
id_wdian_rho_ndiff = register_diag_field ('ocean_model', 'wdian_rho_ndiff',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to time-explicit neutral diffusion', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_rho_ndiff > 0) compute_watermass_diag=.true.
id_wdian_rho_ndiff_on_nrho = -1
id_wdian_rho_ndiff_on_nrho = register_diag_field ('ocean_model', 'wdian_rho_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to time-explicit neutral diffusion as binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_rho_ndiff_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_ndiff = -1
id_tform_rho_ndiff = register_diag_field ('ocean_model', 'tform_rho_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'watermass transform due to time-explicit neutral diffusion on levels (pre-layer binning)',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_rho_ndiff > 0) compute_watermass_diag=.true.
id_tform_rho_ndiff_on_nrho = -1
id_tform_rho_ndiff_on_nrho = register_diag_field ('ocean_model', 'tform_rho_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'watermass transform due to time-explicit neutral diffusion as binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_rho_ndiff_on_nrho > 0) compute_watermass_diag=.true.
id_eta_tend_ndiff_tend= -1
id_eta_tend_ndiff_tend= register_diag_field ('ocean_model','eta_tend_ndiff_tend', &
Grd%tracer_axes(1:2), Time%model_time, &
'numerical form of non-Bouss steric sea level tend from time explicit ndiffusion', 'm/s',&
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_ndiff_tend > 0) compute_watermass_diag=.true.
id_eta_tend_ndiff_tend_glob= -1
id_eta_tend_ndiff_tend_glob= register_diag_field ('ocean_model', 'eta_tend_ndiff_tend_glob', &
Time%model_time, &
'global mean numerical form of non-bouss steric sea level tend from time explicit ndiffusion',&
'm/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_ndiff_tend_glob > 0) compute_watermass_diag=.true.
! temp contributions to neutral diffusion
id_neut_temp_ndiff = -1
id_neut_temp_ndiff = register_diag_field ('ocean_model', 'neut_temp_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'temp related update of locally ref potrho time-explicit neutral diffusion',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_temp_ndiff > 0) compute_watermass_diag=.true.
id_neut_temp_ndiff_on_nrho = -1
id_neut_temp_ndiff_on_nrho = register_diag_field ('ocean_model', 'neut_temp_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related update of locally ref potrho from time-explicit neutral diff binned to neutral rho layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_temp_ndiff_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_temp_ndiff = -1
id_wdian_temp_ndiff = register_diag_field ('ocean_model', 'wdian_temp_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport from time-explicit neutral diffusion',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_temp_ndiff > 0) compute_watermass_diag=.true.
id_wdian_temp_ndiff_on_nrho = -1
id_wdian_temp_ndiff_on_nrho = register_diag_field ('ocean_model', 'wdian_temp_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport from time-explicit neutral diffusion binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_temp_ndiff_on_nrho > 0) compute_watermass_diag=.true.
id_tform_temp_ndiff = -1
id_tform_temp_ndiff = register_diag_field ('ocean_model', 'tform_temp_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'temp related watermass transform due to time-explicit neutral diffusion on levels',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_temp_ndiff > 0) compute_watermass_diag=.true.
id_tform_temp_ndiff_on_nrho = -1
id_tform_temp_ndiff_on_nrho = register_diag_field ('ocean_model', 'tform_temp_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related watermass transform due to time-explicit neutral diffusion binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_temp_ndiff_on_nrho > 0) compute_watermass_diag=.true.
! salinity contributions to neutral diffusion
id_neut_salt_ndiff = -1
id_neut_salt_ndiff = register_diag_field ('ocean_model', 'neut_salt_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'salt related update of locally ref potrho time-explicit neutral diffusion',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_salt_ndiff > 0) compute_watermass_diag=.true.
id_neut_salt_ndiff_on_nrho = -1
id_neut_salt_ndiff_on_nrho = register_diag_field ('ocean_model', 'neut_salt_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related update of locally ref potrho from time-explicit neutral diff binned to neutral rho layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_salt_ndiff_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_salt_ndiff = -1
id_wdian_salt_ndiff = register_diag_field ('ocean_model', 'wdian_salt_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport from time-explicit neutral diffusion',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_salt_ndiff > 0) compute_watermass_diag=.true.
id_wdian_salt_ndiff_on_nrho = -1
id_wdian_salt_ndiff_on_nrho = register_diag_field ('ocean_model', 'wdian_salt_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport from time-explicit neutral diffusion binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_salt_ndiff_on_nrho > 0) compute_watermass_diag=.true.
id_tform_salt_ndiff = -1
id_tform_salt_ndiff = register_diag_field ('ocean_model', 'tform_salt_ndiff', &
Grd%tracer_axes(1:3), Time%model_time, &
'salt related watermass transform due to time-explicit neutral diffusion on levels',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_salt_ndiff > 0) compute_watermass_diag=.true.
id_tform_salt_ndiff_on_nrho = -1
id_tform_salt_ndiff_on_nrho = register_diag_field ('ocean_model', 'tform_salt_ndiff_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related watermass transform due to time-explicit neutral diffusion binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_salt_ndiff_on_nrho > 0) compute_watermass_diag=.true.
! full contributions to GM transport
id_neut_rho_gm = -1
id_neut_rho_gm = register_diag_field ('ocean_model', 'neut_rho_gm', &
Grd%tracer_axes(1:3), Time%model_time, &
'update of locally referenced potential density from GM transport',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_rho_gm > 0) compute_watermass_diag=.true.
id_neut_rho_gm_on_nrho = -1
id_neut_rho_gm_on_nrho = register_diag_field ('ocean_model', 'neut_rho_gm_on_nrho',&
Dens%neutralrho_axes(1:3), Time%model_time, &
'update of locally ref potrho from GM transport as binned to neutral rho layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_rho_gm_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_gm = -1
id_wdian_rho_gm = register_diag_field ('ocean_model', 'wdian_rho_gm',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to GM transport', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_rho_gm > 0) compute_watermass_diag=.true.
id_wdian_rho_gm_on_nrho = -1
id_wdian_rho_gm_on_nrho = register_diag_field ('ocean_model', 'wdian_rho_gm_on_nrho',&
Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to GM transport as binned to neutral rho layers', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_rho_gm_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_gm = -1
id_tform_rho_gm = register_diag_field ('ocean_model', 'tform_rho_gm', &
Grd%tracer_axes(1:3), Time%model_time, &
'watermass transform due to GM transport on levels (pre-layer binning)',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_rho_gm > 0) compute_watermass_diag=.true.
id_tform_rho_gm_on_nrho = -1
id_tform_rho_gm_on_nrho = register_diag_field ('ocean_model', 'tform_rho_gm_on_nrho',&
Dens%neutralrho_axes(1:3), Time%model_time, &
'watermass transform due to GM transport as binned to neutral rho layers', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_rho_gm_on_nrho > 0) compute_watermass_diag=.true.
id_eta_tend_gm_tend= -1
id_eta_tend_gm_tend= register_diag_field ('ocean_model','eta_tend_gm_tend', &
Grd%tracer_axes(1:2), Time%model_time, &
'numerical form of non-Bouss steric sea level tend from GM-based tendency', 'm/s',&
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_gm_tend > 0) compute_watermass_diag=.true.
id_eta_tend_gm_tend_glob= -1
id_eta_tend_gm_tend_glob= register_diag_field ('ocean_model', 'eta_tend_gm_tend_glob', &
Time%model_time, &
'global mean numerical form of non-bouss steric sea level tend from GM-based tendency',&
'm/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_gm_tend_glob > 0) compute_watermass_diag=.true.
! temp contributions to GM transport
id_neut_temp_gm = -1
id_neut_temp_gm = register_diag_field ('ocean_model', 'neut_temp_gm', &
Grd%tracer_axes(1:3), Time%model_time, &
'temp related update of locally referenced potential density from GM transport',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_temp_gm > 0) compute_watermass_diag=.true.
id_neut_temp_gm_on_nrho = -1
id_neut_temp_gm_on_nrho = register_diag_field ('ocean_model', 'neut_temp_gm_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related update of locally ref potrho from GM transport binned to neutral rho layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_temp_gm_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_temp_gm = -1
id_wdian_temp_gm = register_diag_field ('ocean_model', 'wdian_temp_gm',&
Grd%tracer_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport due to GM transport', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_temp_gm > 0) compute_watermass_diag=.true.
id_wdian_temp_gm_on_nrho = -1
id_wdian_temp_gm_on_nrho = register_diag_field ('ocean_model', 'wdian_temp_gm_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport due to GM transport binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_temp_gm_on_nrho > 0) compute_watermass_diag=.true.
id_tform_temp_gm = -1
id_tform_temp_gm = register_diag_field ('ocean_model', 'tform_temp_gm',&
Grd%tracer_axes(1:3), Time%model_time, &
'temp related watermass transform due to GM transport on levels', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_temp_gm > 0) compute_watermass_diag=.true.
id_tform_temp_gm_on_nrho = -1
id_tform_temp_gm_on_nrho = register_diag_field ('ocean_model', 'tform_temp_gm_on_nrho',&
Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related watermass transform due to GM transport binned to neutral rho layers', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_temp_gm_on_nrho > 0) compute_watermass_diag=.true.
! salinity contributions to GM transport
id_neut_salt_gm = -1
id_neut_salt_gm = register_diag_field ('ocean_model', 'neut_salt_gm', &
Grd%tracer_axes(1:3), Time%model_time, &
'salt related update of locally referenced potential density from GM transport',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_salt_gm > 0) compute_watermass_diag=.true.
id_neut_salt_gm_on_nrho = -1
id_neut_salt_gm_on_nrho = register_diag_field ('ocean_model', 'neut_salt_gm_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related update of locally ref potrho from GM transport binned to neutral rho layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_neut_salt_gm_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_salt_gm = -1
id_wdian_salt_gm = register_diag_field ('ocean_model', 'wdian_salt_gm',&
Grd%tracer_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport due to GM transport', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_salt_gm > 0) compute_watermass_diag=.true.
id_wdian_salt_gm_on_nrho = -1
id_wdian_salt_gm_on_nrho = register_diag_field ('ocean_model', 'wdian_salt_gm_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport due to GM transport binned to neutral rho layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wdian_salt_gm_on_nrho > 0) compute_watermass_diag=.true.
id_tform_salt_gm = -1
id_tform_salt_gm = register_diag_field ('ocean_model', 'tform_salt_gm',&
Grd%tracer_axes(1:3), Time%model_time, &
'salt related watermass transform due to GM transport on levels', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_salt_gm > 0) compute_watermass_diag=.true.
id_tform_salt_gm_on_nrho = -1
id_tform_salt_gm_on_nrho = register_diag_field ('ocean_model', 'tform_salt_gm_on_nrho',&
Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related watermass transform due to GM transport binned to neutral rho layers', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tform_salt_gm_on_nrho > 0) compute_watermass_diag=.true.
if(compute_watermass_diag) then
write(stdoutunit,'(/a/)') &
'==>Note: running ocean_nphysics_util_mod w/ compute_watermass_diag=.true.'
diagnose_gm_redi=.true.
endif
end subroutine watermass_diag_init
! NAME="watermass_diag_init"
!#######################################################################
!
!
!
! Diagnose effects from nphysics on the watermass transformation.
! For use with nphysicsA and nphysicsB.
!
!
subroutine watermass_diag(Time, T_prog, Dens, tendency_redi_temp, tendency_redi_salt)
type(ocean_time_type), intent(in) :: Time
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: tendency_redi_temp
real, dimension(isd:,jsd:,:), intent(in) :: tendency_redi_salt
integer :: i,j,k,tau
real :: tmp1, tmp2
real, dimension(isd:ied,jsd:jed) :: eta_tend
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error from ocean_nphysics_util (watermass_diag): module needs initialization ')
endif
if(.not. compute_watermass_diag) return
tau = Time%tau
! impacts from full neutral physics operator
if(id_neut_rho_nphysics > 0 .or. id_wdian_rho_nphysics > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Dens%drhodT(i,j,k)*T_prog(index_temp)%wrk1(i,j,k) &
+Dens%drhodS(i,j,k)*T_prog(index_salt)%wrk1(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_rho_nphysics, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_nphysics, wrk3(:,:,:))
endif
! impacts from neutral diffusion
if(id_neut_rho_ndiff > 0 .or. id_neut_rho_ndiff_on_nrho > 0 .or. &
id_wdian_rho_ndiff > 0 .or. id_wdian_rho_ndiff_on_nrho > 0 .or. &
id_tform_rho_ndiff > 0 .or. id_tform_rho_ndiff_on_nrho > 0 .or. &
id_eta_tend_ndiff_tend > 0 .or. id_eta_tend_ndiff_tend_glob > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Dens%drhodT(i,j,k)*tendency_redi_temp(i,j,k) &
+Dens%drhodS(i,j,k)*tendency_redi_salt(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
wrk5(i,j,k) =-wrk1(i,j,k)/(epsln+Dens%rho(i,j,k,tau)**2) ! for eta_tend
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_rho_ndiff, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_ndiff, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_ndiff, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_ndiff_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_ndiff_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_ndiff_on_nrho, wrk4)
if(id_eta_tend_ndiff_tend > 0 .or. id_eta_tend_ndiff_tend_glob > 0) then
eta_tend(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
eta_tend(i,j) = eta_tend(i,j) + wrk5(i,j,k)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_ndiff_tend, eta_tend(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_ndiff_tend_glob, eta_tend, cellarea_r)
endif
endif ! endif for ndiffuse
! temperature contributions to neutral diffusion
if(id_neut_temp_ndiff > 0 .or. id_neut_temp_ndiff_on_nrho > 0 .or. &
id_wdian_temp_ndiff > 0 .or. id_wdian_temp_ndiff_on_nrho > 0 .or. &
id_tform_temp_ndiff > 0 .or. id_tform_temp_ndiff_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Dens%drhodT(i,j,k)*tendency_redi_temp(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_temp_ndiff, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_temp_ndiff, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_temp_ndiff, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_ndiff_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_temp_ndiff_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_temp_ndiff_on_nrho, wrk4)
endif ! endif for temp_ndiffuse
! salinity contributions to neutral diffusion
if(id_neut_salt_ndiff > 0 .or. id_neut_salt_ndiff_on_nrho > 0 .or. &
id_wdian_salt_ndiff > 0 .or. id_wdian_salt_ndiff_on_nrho > 0 .or. &
id_tform_salt_ndiff > 0 .or. id_tform_salt_ndiff_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Dens%drhodS(i,j,k)*tendency_redi_salt(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_salt_ndiff, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_salt_ndiff, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_salt_ndiff, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_ndiff_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_salt_ndiff_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_salt_ndiff_on_nrho, wrk4)
endif ! endif for salt_ndiffuse
! impacts from GM skewsion
if(id_neut_rho_gm > 0 .or. id_neut_rho_gm_on_nrho > 0 .or. &
id_wdian_rho_gm > 0 .or. id_wdian_rho_gm_on_nrho > 0 .or. &
id_tform_rho_gm > 0 .or. id_tform_rho_gm_on_nrho > 0 .or. &
id_eta_tend_gm_tend > 0 .or. id_eta_tend_gm_tend_glob > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
tmp1 = T_prog(index_temp)%wrk1(i,j,k) - tendency_redi_temp(i,j,k)
tmp2 = T_prog(index_salt)%wrk1(i,j,k) - tendency_redi_salt(i,j,k)
wrk1(i,j,k) = Dens%drhodT(i,j,k)*tmp1 + Dens%drhodS(i,j,k)*tmp2
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
wrk5(i,j,k) =-wrk1(i,j,k)/(epsln+Dens%rho(i,j,k,tau)**2) ! for eta_tend
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_rho_gm, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_gm, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_gm, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_gm_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_gm_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_gm_on_nrho, wrk4)
if(id_eta_tend_gm_tend > 0 .or. id_eta_tend_gm_tend_glob > 0) then
eta_tend(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
eta_tend(i,j) = eta_tend(i,j) + wrk5(i,j,k)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_gm_tend, eta_tend(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_gm_tend_glob, eta_tend, cellarea_r)
endif
endif ! endif for GM
! impacts from GM skewsion due to temperature contributions
if(id_neut_temp_gm > 0 .or. id_neut_temp_gm_on_nrho > 0 .or. &
id_wdian_temp_gm > 0 .or. id_wdian_temp_gm_on_nrho > 0 .or. &
id_tform_temp_gm > 0 .or. id_tform_temp_gm_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
tmp1 = T_prog(index_temp)%wrk1(i,j,k) - tendency_redi_temp(i,j,k)
wrk1(i,j,k) = Dens%drhodT(i,j,k)*tmp1
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_temp_gm, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_temp_gm, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_temp_gm, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_gm_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_temp_gm_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_temp_gm_on_nrho, wrk4)
endif ! endif for GM from temperature contributions
! impacts from GM skewsion due to salinity contributions
if(id_neut_salt_gm > 0 .or. id_neut_salt_gm_on_nrho > 0 .or. &
id_wdian_salt_gm > 0 .or. id_wdian_salt_gm_on_nrho > 0 .or. &
id_tform_salt_gm > 0 .or. id_tform_salt_gm_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
tmp2 = T_prog(index_salt)%wrk1(i,j,k) - tendency_redi_salt(i,j,k)
wrk1(i,j,k) = Dens%drhodS(i,j,k)*tmp2
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_salt_gm, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_salt_gm, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_salt_gm, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_gm_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_salt_gm_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_salt_gm_on_nrho, wrk4)
endif ! endif for GM from salinity contributions
end subroutine watermass_diag
! NAME="watermass_diag"
!#######################################################################
!
!
!
! Diagnose watermass transformation from neutral diffusion.
! For use with the nphysicsC scheme.
!
!
subroutine watermass_diag_ndiffuse(Time, Dens, tendency_redi_temp, tendency_redi_salt)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: tendency_redi_temp
real, dimension(isd:,jsd:,:), intent(in) :: tendency_redi_salt
integer :: i,j,k,tau
real, dimension(isd:ied,jsd:jed) :: eta_tend
if(.not. compute_watermass_diag) return
tau = Time%tau
if (id_neut_rho_ndiff > 0 .or. id_neut_rho_ndiff_on_nrho > 0 .or. &
id_wdian_rho_ndiff > 0 .or. id_wdian_rho_ndiff_on_nrho > 0 .or. &
id_tform_rho_ndiff > 0 .or. id_tform_rho_ndiff_on_nrho > 0 .or. &
id_eta_tend_ndiff_tend > 0 .or. id_eta_tend_ndiff_tend_glob > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*(Dens%drhodT(i,j,k)*tendency_redi_temp(i,j,k) &
+Dens%drhodS(i,j,k)*tendency_redi_salt(i,j,k))
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
wrk5(i,j,k) =-wrk1(i,j,k)/(epsln+Dens%rho(i,j,k,tau)**2)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_rho_ndiff,wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_ndiff, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_ndiff, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_ndiff_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_ndiff_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_ndiff_on_nrho, wrk4)
if(id_eta_tend_ndiff_tend > 0 .or. id_eta_tend_ndiff_tend_glob > 0) then
eta_tend(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
eta_tend(i,j) = eta_tend(i,j) + wrk5(i,j,k)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_ndiff_tend, eta_tend(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_ndiff_tend_glob, eta_tend, cellarea_r)
endif
endif !endif for the neut_rho, wdian, and tform diagnostics iftest
! ndiffuse scheme contributions from temperature
if (id_neut_temp_ndiff > 0 .or. id_neut_temp_ndiff_on_nrho > 0 .or. &
id_wdian_temp_ndiff > 0 .or. id_wdian_temp_ndiff_on_nrho > 0 .or. &
id_tform_temp_ndiff > 0 .or. id_tform_temp_ndiff_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*Dens%drhodT(i,j,k)*tendency_redi_temp(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_temp_ndiff,wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_temp_ndiff, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_temp_ndiff, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_ndiff_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_temp_ndiff_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_temp_ndiff_on_nrho, wrk4)
endif !endif for the neut_rho, wdian, and tform diagnostics from temperature effects
! ndiffuse scheme contributions from salinity
if (id_neut_salt_ndiff > 0 .or. id_neut_salt_ndiff_on_nrho > 0 .or. &
id_wdian_salt_ndiff > 0 .or. id_wdian_salt_ndiff_on_nrho > 0 .or. &
id_tform_salt_ndiff > 0 .or. id_tform_salt_ndiff_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*Dens%drhodS(i,j,k)*tendency_redi_salt(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_salt_ndiff,wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_salt_ndiff, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_salt_ndiff, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_ndiff_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_salt_ndiff_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_salt_ndiff_on_nrho, wrk4)
endif !endif for the neut_rho, wdian, and tform diagnostics from salinity effects
end subroutine watermass_diag_ndiffuse
! NAME="watermass_diag_ndiffuse"
!#######################################################################
!
!
!
! Diagnose watermass transformation from skew diffusion.
! For use with the neutral physics C scheme.
!
!
subroutine watermass_diag_sdiffuse(Time, Dens, tendency_gm_temp, tendency_gm_salt)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: tendency_gm_temp
real, dimension(isd:,jsd:,:), intent(in) :: tendency_gm_salt
integer :: i,j,k,tau
real, dimension(isd:ied,jsd:jed) :: eta_tend
if(.not. compute_watermass_diag) return
tau = Time%tau
if (id_neut_rho_gm > 0 .or. id_neut_rho_gm_on_nrho > 0 .or. &
id_wdian_rho_gm > 0 .or. id_wdian_rho_gm_on_nrho > 0 .or. &
id_tform_rho_gm > 0 .or. id_tform_rho_gm_on_nrho > 0 .or. &
id_eta_tend_gm_tend > 0 .or. id_eta_tend_gm_tend_glob > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*(Dens%drhodT(i,j,k)*tendency_gm_temp(i,j,k) &
+Dens%drhodS(i,j,k)*tendency_gm_salt(i,j,k))
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
wrk5(i,j,k) =-wrk1(i,j,k)/(epsln+Dens%rho(i,j,k,tau)**2) ! for eta_tend
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_rho_gm,wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_gm, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_gm, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_gm_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_gm_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_gm_on_nrho, wrk4)
if(id_eta_tend_gm_tend > 0 .or. id_eta_tend_gm_tend_glob > 0) then
eta_tend(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
eta_tend(i,j) = eta_tend(i,j) + wrk5(i,j,k)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_gm_tend, eta_tend(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_gm_tend_glob, eta_tend, cellarea_r)
endif
endif !endif for the neut_rho, wdian, and tform diagnostics iftest
! temperature contributions
if (id_neut_temp_gm > 0 .or. id_neut_temp_gm_on_nrho > 0 .or. &
id_wdian_temp_gm > 0 .or. id_wdian_temp_gm_on_nrho > 0 .or. &
id_tform_temp_gm > 0 .or. id_tform_temp_gm_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*Dens%drhodT(i,j,k)*tendency_gm_temp(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_temp_gm, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_temp_gm, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_temp_gm, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_gm_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_temp_gm_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_temp_gm_on_nrho, wrk4)
endif !endif for the neut_rho, wdian, and tform diagnostics from temperature effects
! salinity contributions
if (id_neut_salt_gm > 0 .or. id_neut_salt_gm_on_nrho > 0 .or. &
id_wdian_salt_gm > 0 .or. id_wdian_salt_gm_on_nrho > 0 .or. &
id_tform_salt_gm > 0 .or. id_tform_salt_gm_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*Dens%drhodS(i,j,k)*tendency_gm_salt(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_salt_gm, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_salt_gm, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_salt_gm, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_gm_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_salt_gm_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_salt_gm_on_nrho, wrk4)
endif !endif for the neut_rho, wdian, and tform diagnostics from salinity effects
end subroutine watermass_diag_sdiffuse
! NAME="watermass_diag_sdiffuse"
!#######################################################################
!
!
! Write out restart files registered through register_restart_file
!
subroutine ocean_nphysics_util_restart(time_stamp)
character(len=*), intent(in), optional :: time_stamp
call save_restart(Nphysics_util_restart, time_stamp)
end subroutine ocean_nphysics_util_restart
! NAME="ocean_nphysics_util_restart"
!#######################################################################
!
!
!
! Write to restart.
!
!
subroutine ocean_nphysics_coeff_end(Time, agm_array, aredi_array, &
rossby_radius, rossby_radius_raw, bczone_radius)
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:,jsd:,:), intent(in) :: agm_array
real, dimension(isd:,jsd:,:), intent(in) :: aredi_array
real, dimension(isd:,jsd:), intent(in) :: rossby_radius
real, dimension(isd:,jsd:), intent(in) :: rossby_radius_raw
real, dimension(isd:,jsd:), intent(in) :: bczone_radius
integer :: stdoutunit
stdoutunit=stdout()
write(stdoutunit,*)' '
write(stdoutunit,*) 'From ocean_nphysics_coeff_end: ending chksum'
call write_timestamp(Time%model_time)
call write_chksum_3d('ending agm_array', agm_array(COMP,:)*Grd%tmask(COMP,:))
call write_chksum_3d('ending aredi_array', aredi_array(COMP,:)*Grd%tmask(COMP,:))
call write_chksum_2d('ending rossby_radius', rossby_radius(COMP)*Grd%tmask(COMP,1))
call write_chksum_2d('ending rossby_radius_raw', rossby_radius_raw(COMP)*Grd%tmask(COMP,1))
call write_chksum_2d('ending bczone_radius', bczone_radius(COMP)*Grd%tmask(COMP,1))
end subroutine ocean_nphysics_coeff_end
! NAME="ocean_nphysics_coeff_end"
end module ocean_nphysics_util_mod