module ocean_tracer_mod
#define COMP isc:iec,jsc:jec
!
! Stephen Griffies
!
!
! Matt Harrison
!
!
! Richard D. Slater (initialization)
!
!
!
! This module time steps the tracer fields.
!
!
!
! This module time steps the tracer fields.
! Initialization for the tracer packages is done as well.
!
!
!
!
!
! R.C. Pacanowski and S.M. Griffies
! The MOM3 Manual (1999)
!
!
!
! S.M. Griffies, M.J. Harrison, R.C. Pacanowski, and A. Rosati
! A Technical Guide to MOM4 (2004)
!
!
!
! S.M. Griffies, Fundamentals of ocean climate models (2004)
!
!
!
!
!
!
!
! If true, then will freeze the tracer fields.
!
!
!
! To remove the T_prog%source contribution to tracer
! evolution. For debugging purposes. Default
! zero_tracer_source=.false.
!
!
!
! Limit the values of age tracer to be less than
! total run time and greater than zero.
! Default limit_age_tracer=.false.
!
!
!
! Initial maximum age tracer. This nml provides the ability to
! start an integration with an age tracer that is not initialized
! to zero, say if we took an age tracer from another spin-up.
! Default age_tracer_max_init=0.0.
!
!
!
! For remapping initial tracer distributions, generally determined
! according to depth vertical coordinates using the mom preprocessing
! schemes, onto s-coordinates. This method is of use for initializing
! terrain following coordinate simulations with mom.
!
!
!
! For computing frazil heating before the implicit vertical physics
! (which includes boundary fluxes), and before vertical convection.
! This is the order that CM2.0 and CM2.1 performed their calculations
! of frazil. It is arguable that one should NOT do frazil until the
! end of a time step, after vertical physics and after surface
! boundary fluxes.
! Default frazil_heating_before_vphysics=.false.
!
!
!
! For computing frazil heating after the implicit vertical physics
! (which includes boundary fluxes), and after vertical convection.
! This is the recommended method.
! Default frazil_heating_after_vphysics=.false.
!
!
!
! tmask_limit is derived separately for the tracers. However,
! it may be appropriate to have the mask be the same for temp
! and salinity, in which case the neutral physics fluxes are
! self-consistent. But for some cases, such as when running with
! linear eos, may not wish to have the temp and salinity coupled
! when computing the mask.
!
!
!
! For updating the tmaks_limit array. This calculation is
! recommended for the following physics and advection schemes:
! 1/ quicker advection
! 2/ neutral physics
! 3/ submesoscale closure.
! The default is compute_tmask_limit_on=.true., but if none
! of the above schemes is used, then some time savings can be
! realized by setting compute_tmask_limit_on=.false.
!
!
!
! For debugging the tracer module
!
!
! Set true to write a restart. False setting only for rare
! cases where wish to benchmark model without measuring the cost
! of writing restarts and associated chksums.
! Default is write_a_restart=.true.
!
!
!
! To linear interpolate the initial conditions for prognostic
! tracers to the partial bottom cells. Default
! interpolate_tprog_to_pbott=.true.
!
!
!
! To linear interpolate the initial conditions for diagnostic
! tracers to the partial bottom cells. Default
! interpolate_tdiag_to_pbott=.false.
!
!
!
! For adding an inflow transport from the northern boundary
! which brings in temp and salinity according to inflow data
! files. Default is inflow_nboundary=.false.
!
!
!
! For debugging ocean tracer package manager.
!
!
!
! To call a check to see that temperature and salinity
! are within their pre-selected range.
! Default use_tempsalt_check_range=.false. since this
! check may incur some cost that users should be aware of.
!
!
!
use constants_mod, only: epsln, kelvin
use diag_manager_mod, only: register_diag_field, send_data
use field_manager_mod, only: fm_string_len, fm_type_name_len
use field_manager_mod, only: fm_get_length
use field_manager_mod, only: fm_dump_list, fm_change_list, fm_loop_over_list
use fms_mod, only: read_data, file_exist, field_exist
use fms_mod, only: open_namelist_file, check_nml_error, close_file
use fms_mod, only: FATAL, WARNING, NOTE, stdout, stdlog
use fms_io_mod, only: register_restart_field, save_restart, restore_state
use fms_io_mod, only: restart_file_type, reset_field_pointer, reset_field_name
use fms_io_mod, only: field_size
use fm_util_mod, only: fm_util_get_real, fm_util_get_integer
use fm_util_mod, only: fm_util_get_logical, fm_util_get_string
use fm_util_mod, only: fm_util_get_string_array
use fm_util_mod, only: fm_util_check_for_bad_fields, fm_util_set_value
use mpp_domains_mod, only: mpp_update_domains
use mpp_domains_mod, only: mpp_global_sum, NON_BITWISE_EXACT_SUM
use mpp_io_mod, only: mpp_open, fieldtype
use mpp_io_mod, only: MPP_NETCDF, MPP_OVERWR, MPP_ASCII, MPP_RDONLY, MPP_SINGLE, MPP_MULTI
use mpp_mod, only: input_nml_file, mpp_error, mpp_pe, mpp_root_pe, mpp_broadcast, ALL_PES, mpp_max, mpp_min, mpp_chksum
use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end
use mpp_mod, only: CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, CLOCK_MODULE, CLOCK_ROUTINE
use platform_mod, only: i8_kind
use time_manager_mod, only: time_type, set_time, increment_time, operator( + )
use transport_matrix_mod, only: do_transport_matrix, transport_matrix_store_explicit
use ocean_blob_mod, only: ocean_blob_implicit, adjust_L_thickness
use ocean_convect_mod, only: convection
use ocean_domains_mod, only: get_local_indices
use ocean_density_mod, only: neutral_density, update_ocean_density_salinity
use ocean_frazil_mod, only: compute_frazil_heating
use ocean_bih_tracer_mod, only: bih_tracer
#ifdef USE_OCEAN_BGC
use ocean_generic_mod, only: ocean_generic_get_field, ocean_generic_get_field_pointer
use ocean_generic_mod, only: ocean_generic_set_pointer
#endif
use ocean_lap_tracer_mod, only: lap_tracer
use ocean_obc_mod, only: ocean_obc_tracer, ocean_obc_update_boundary, ocean_obc_tracer_init
use ocean_parameters_mod, only: ADVECT_UPWIND, ADVECT_2ND_ORDER, ADVECT_4TH_ORDER, ADVECT_6TH_ORDER
use ocean_parameters_mod, only: ADVECT_QUICKER, ADVECT_QUICKMOM3, ADVECT_MDFL_SUP_B, ADVECT_MDFL_SWEBY
use ocean_parameters_mod, only: ADVECT_PSOM, ADVECT_MDPPM, ADVECT_DST_LINEAR, ADVECT_MDMDT_TEST
use ocean_parameters_mod, only: ADVECT_MDPPM_TEST, ADVECT_MDFL_SWEBY_TEST, ADVECT_DST_LINEAR_TEST
use ocean_parameters_mod, only: missing_value, sec_in_yr_r, rho0, rho0r, cp_ocean
use ocean_parameters_mod, only: TWO_LEVEL, THREE_LEVEL
use ocean_parameters_mod, only: CONSERVATIVE_TEMP, POTENTIAL_TEMP
use ocean_parameters_mod, only: QUASI_HORIZONTAL, TERRAIN_FOLLOWING
use ocean_parameters_mod, only: GEOPOTENTIAL
use ocean_parameters_mod, only: DEPTH_BASED
use ocean_passive_mod, only: passive_tracer_init, update_tracer_passive
use ocean_shortwave_mod, only: ocean_irradiance_init
use ocean_tempsalt_mod, only: contemp_from_pottemp, pottemp_from_contemp, tempsalt_check_range
use ocean_tpm_mod, only: ocean_tpm_init
use ocean_tpm_util_mod, only: otpm_set_tracer_package
use ocean_tracer_advect_mod, only: horz_advect_tracer, vert_advect_tracer
use ocean_tracer_diag_mod, only: send_tracer_variance
use ocean_tracer_util_mod, only: rebin_onto_rho, diagnose_mass_of_layer
use ocean_tracer_util_mod, only: tracer_prog_chksum, tracer_diag_chksum, tracer_min_max
use ocean_thickness_mod, only: update_E_thickness
use ocean_types_mod, only: ocean_grid_type, ocean_domain_type, ocean_thickness_type
use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_diag_tracer_type
use ocean_types_mod, only: ocean_adv_vel_type, ocean_external_mode_type
use ocean_types_mod, only: ocean_public_type, ocean_density_type, ocean_options_type
use ocean_types_mod, only: ocean_lagrangian_type, ocean_velocity_type, blob_diag_type
use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_sum, write_chksum_3d
use ocean_tracer_util_mod, only: diagnose_3d_rho
use ocean_vert_mix_mod, only: vert_diffuse, vert_diffuse_implicit
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk6, wrk1_2d
implicit none
private
logical :: prog_module_initialized = .false.
logical :: diag_module_initialized = .false.
character(len=256) :: version='CVS $Id: ocean_tracer.F90,v 20.0 2013/12/14 00:17:20 fms Exp $'
character(len=256) :: tagname='Tag $Name: tikal $'
character(len=48), parameter :: mod_name = 'ocean_tracer_mod'
integer :: num_tracers =0
integer :: num_prog_tracers =0
integer :: num_diag_tracers =0
integer :: num_family_tracers=0
#include
type(ocean_grid_type), pointer :: Grd =>NULL()
type(ocean_domain_type), pointer :: Dom =>NULL()
! for time steps and implicit vertical mixing
integer :: tendency = 0
real :: dtime = 0.0
real :: dtimer = 0.0
real :: dtts = 0.0
real :: dtuv = 0.0
real :: cp_oceanr
integer :: index_temp =-1
integer :: index_salt =-1
integer :: index_frazil =-1
integer :: index_temp_sq =-1
integer :: index_salt_sq =-1
integer :: index_diag_temp =-1
integer :: index_added_heat =-1
integer :: index_redist_heat =-1
! for obc
logical :: have_obc=.false.
! for global normalization
real :: cellarea_r
! for vertical coordinate class
integer :: vert_coordinate_class
! for setting the temperature variable
integer :: prog_temp_variable=0
! for CMIP units (K rather than degrees C if cmip_version<=5)
real :: cmip_offset = 0.0
character(len=32) :: temp_units='degrees C'
! for possible inflow tracer settings
real, dimension(:,:,:), allocatable :: mask_inflow
real, dimension(:,:,:), allocatable :: salt_inflow
real, dimension(:,:,:), allocatable :: temp_inflow
! identification numbers for mpp clocks
integer :: id_clock_bih_tracer
integer :: id_clock_lap_tracer
integer :: id_clock_vert_diffuse
integer :: id_clock_vert_diffuse_implicit
integer :: id_clock_tracer_advect
integer :: id_clock_frazil
integer :: id_clock_convection
integer :: id_clock_blob
integer :: id_clock_blob_implicit
integer :: id_clock_adjust_L_thickness
integer :: id_clock_adjust_E_thickness
! for restart
integer :: num_tracer_restart = 0
integer, allocatable :: id_tracer_restart(:)
integer, allocatable :: id_tracer_type(:)
character(len=64), allocatable :: tracer_restart_file(:)
type(restart_file_type), allocatable :: Tracer_restart(:)
!work array on neutral density space
integer :: neutralrho_nk
real, dimension(:,:,:), allocatable :: nrho_work
real, dimension(:,:,:), allocatable :: nrho_work2
real, dimension(:,:,:), allocatable :: nrho_mass
! internally set for computing watermass diagnstics
logical :: compute_watermass_diag = .false.
! for diagnostics
logical :: used
integer, allocatable, dimension(:) :: id_eta_smooth
integer, allocatable, dimension(:) :: id_eta_smooth_on_nrho
integer, allocatable, dimension(:) :: id_pbot_smooth
integer, allocatable, dimension(:) :: id_prog
integer, allocatable, dimension(:) :: id_progT
integer, allocatable, dimension(:) :: id_prog_explicit
integer, allocatable, dimension(:) :: id_prog_rhodzt
integer, allocatable, dimension(:) :: id_prog_rhodztT
integer, allocatable, dimension(:) :: id_prog_int_rhodz
integer, allocatable, dimension(:) :: id_prog_int_rhodzT
integer, allocatable, dimension(:) :: id_prog_on_depth
integer, allocatable, dimension(:) :: id_tendency_conc
integer, allocatable, dimension(:) :: id_tendency_concL
integer, allocatable, dimension(:) :: id_tendency_concT
integer, allocatable, dimension(:) :: id_tendency
integer, allocatable, dimension(:) :: id_tendency_on_nrho
integer, allocatable, dimension(:) :: id_tendencyL
integer, allocatable, dimension(:) :: id_tendencyT
integer, allocatable, dimension(:) :: id_tendency_expl
integer, allocatable, dimension(:) :: id_tendency_explL
integer, allocatable, dimension(:) :: id_tendency_explT
integer, allocatable, dimension(:) :: id_surf_tracer
integer, allocatable, dimension(:) :: id_surf_tracerT
integer, allocatable, dimension(:) :: id_surf_tracer_sq
integer, allocatable, dimension(:) :: id_diag_surf_tracer
integer, allocatable, dimension(:) :: id_diag_surf_tracer_sq
integer, allocatable, dimension(:) :: id_bott_tracer
integer, allocatable, dimension(:) :: id_bott_tracerT
integer, allocatable, dimension(:) :: id_diag
integer, allocatable, dimension(:) :: id_diag_total
integer, allocatable, dimension(:) :: id_tmask_limit
integer :: id_neut_rho_tendency
integer :: id_pot_rho_tendency
integer :: id_neut_rho_tendency_on_nrho
integer :: id_wdian_rho_tendency
integer :: id_wdian_rho_tendency_on_nrho
integer :: id_tform_rho_tendency
integer :: id_tform_rho_tendency_on_nrho
integer :: id_neut_temp_tendency
integer :: id_neut_temp_tendency_on_nrho
integer :: id_wdian_temp_tendency
integer :: id_wdian_temp_tendency_on_nrho
integer :: id_tform_temp_tendency
integer :: id_tform_temp_tendency_on_nrho
integer :: id_neut_salt_tendency
integer :: id_neut_salt_tendency_on_nrho
integer :: id_wdian_salt_tendency
integer :: id_wdian_salt_tendency_on_nrho
integer :: id_tform_salt_tendency
integer :: id_tform_salt_tendency_on_nrho
integer :: id_neut_rho_smooth
integer :: id_pot_rho_smooth
integer :: id_neut_rho_smooth_on_nrho
integer :: id_wdian_rho_smooth
integer :: id_wdian_rho_smooth_on_nrho
integer :: id_tform_rho_smooth
integer :: id_tform_rho_smooth_on_nrho
integer :: id_eta_tend_smooth
integer :: id_eta_tend_smooth_glob
integer :: id_neut_temp_smooth
integer :: id_neut_temp_smooth_on_nrho
integer :: id_wdian_temp_smooth
integer :: id_wdian_temp_smooth_on_nrho
integer :: id_tform_temp_smooth
integer :: id_tform_temp_smooth_on_nrho
integer :: id_neut_salt_smooth
integer :: id_neut_salt_smooth_on_nrho
integer :: id_wdian_salt_smooth
integer :: id_wdian_salt_smooth_on_nrho
integer :: id_tform_salt_smooth
integer :: id_tform_salt_smooth_on_nrho
integer :: id_neut_rho_pme
integer :: id_pot_rho_pme
integer :: id_wdian_rho_pme
integer :: id_tform_rho_pme
integer :: id_neut_rho_pme_on_nrho
integer :: id_wdian_rho_pme_on_nrho
integer :: id_tform_rho_pme_on_nrho
integer :: id_neut_rho_pbl_pme_kn
integer :: id_neut_rho_pbl_pme_kn_on_nrho
integer :: id_wdian_rho_pbl_pme_kn
integer :: id_wdian_rho_pbl_pme_kn_on_nrho
integer :: id_tform_rho_pbl_pme_kn
integer :: id_tform_rho_pbl_pme_kn_on_nrho
integer :: id_neut_temp_pbl_pme_kn
integer :: id_neut_temp_pbl_pme_kn_on_nrho
integer :: idwdian_temp_pbl_pme_kn
integer :: idwdian_temp_pbl_pme_kn_on_nrho
integer :: idtform_temp_pbl_pme_kn
integer :: idtform_temp_pbl_pme_kn_on_nrho
integer :: id_neut_salt_pbl_pme_kn
integer :: id_neut_salt_pbl_pme_kn_on_nrho
integer :: idwdian_salt_pbl_pme_kn
integer :: idwdian_salt_pbl_pme_kn_on_nrho
integer :: idtform_salt_pbl_pme_kn
integer :: idtform_salt_pbl_pme_kn_on_nrho
integer :: id_neut_rho_pbl_pme_pr
integer :: id_neut_rho_pbl_pme_pr_on_nrho
integer :: id_wdian_rho_pbl_pme_pr
integer :: id_wdian_rho_pbl_pme_pr_on_nrho
integer :: id_tform_rho_pbl_pme_pr
integer :: id_tform_rho_pbl_pme_pr_on_nrho
integer :: id_neut_temp_pbl_pme_pr
integer :: id_neut_temp_pbl_pme_pr_on_nrho
integer :: idwdian_temp_pbl_pme_pr
integer :: idwdian_temp_pbl_pme_pr_on_nrho
integer :: idtform_temp_pbl_pme_pr
integer :: idtform_temp_pbl_pme_pr_on_nrho
integer :: id_neut_salt_pbl_pme_pr
integer :: id_neut_salt_pbl_pme_pr_on_nrho
integer :: idwdian_salt_pbl_pme_pr
integer :: idwdian_salt_pbl_pme_pr_on_nrho
integer :: idtform_salt_pbl_pme_pr
integer :: idtform_salt_pbl_pme_pr_on_nrho
integer :: id_frazil_on_nrho
integer :: id_neut_rho_frazil
integer :: id_pot_rho_frazil
integer :: id_wdian_rho_frazil
integer :: id_tform_rho_frazil
integer :: id_neut_rho_frazil_on_nrho
integer :: id_wdian_rho_frazil_on_nrho
integer :: id_tform_rho_frazil_on_nrho
integer :: id_eta_tend_frazil
integer :: id_eta_tend_frazil_glob
integer :: id_mass_t_on_nrho
integer :: id_mass_t_tendency_on_nrho
integer :: id_mass_nrho_layer
integer :: id_mass_nrho_tendency_layer
! for ascii output
integer :: unit=6
public update_ocean_tracer
public ocean_prog_tracer_init
public ocean_diag_tracer_init
public ocean_tracer_diagnostics_init
public ocean_tracer_end
public compute_tmask_limit
public ocean_tracer_restart
private update_advection_only
private remap_s_to_depth
private remap_depth_to_s
private inflow_nboundary_init
private watermass_diag
private send_tracer_diagnostics
!---------------nml settings---------------
logical :: zero_tendency = .false.
logical :: zero_tracer_source = .false.
logical :: debug_this_module = .false.
logical :: tmask_limit_ts_same = .true.
logical :: write_a_restart = .true.
logical :: remap_depth_to_s_init = .false.
logical :: inflow_nboundary = .false.
logical :: interpolate_tprog_to_pbott = .true.
logical :: interpolate_tdiag_to_pbott = .false.
logical :: limit_age_tracer = .false.
logical :: frazil_heating_before_vphysics = .false.
logical :: frazil_heating_after_vphysics = .false.
logical :: compute_tmask_limit_on = .true.
logical :: use_tempsalt_check_range = .false.
real :: age_tracer_max_init = 0.0
! for tracer package manager
logical :: ocean_tpm_debug = .false.
namelist /ocean_tracer_nml/ debug_this_module, zero_tendency, zero_tracer_source, write_a_restart, &
ocean_tpm_debug, tmask_limit_ts_same, remap_depth_to_s_init, &
inflow_nboundary, interpolate_tprog_to_pbott, interpolate_tdiag_to_pbott,&
limit_age_tracer, age_tracer_max_init, &
frazil_heating_before_vphysics, frazil_heating_after_vphysics, &
compute_tmask_limit_on, use_tempsalt_check_range
contains
!#######################################################################
!
!
!
! Initialization code for prognostic tracers, returning a pointer to
! the T_prog array.
!
!
function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, Time_steps, &
num_prog, vert_coordinate_type, obc, cmip_units, cmip_version, use_blobs, debug) &
result (T_prog)
type(ocean_grid_type), intent(in), target :: Grid
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_options_type), intent(inout) :: Ocean_options
type(ocean_domain_type), intent(in), target :: Domain
type(ocean_time_type), intent(in) :: Time
type(ocean_time_steps_type), intent(in) :: Time_steps
integer, intent(out) :: num_prog
integer, intent(in) :: vert_coordinate_type
logical, intent(in) :: obc
logical, intent(in) :: cmip_units
integer, intent(in) :: cmip_version
logical, intent(in) :: use_blobs
logical, intent(in), optional :: debug
! return value
type(ocean_prog_tracer_type), dimension(:), pointer :: T_prog
integer :: i, j, k, n, kb, l
integer :: ierr, num_diag
integer :: tau, taum1, taup1
integer :: ioun, io_status
integer, dimension(4) :: siz
integer :: frazil_heating_order=0
real :: fact
character(len=32) :: name
character(len=128) :: filename
logical :: initialize_as_a_passive_tracer=.false.
character(len=33) :: prog_name
character(len=138) :: prog_longname
character(len=48), parameter :: sub_name = 'ocean_prog_tracer_init'
character(len=256), parameter :: error_header = '==>Error from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
character(len=256), parameter :: warn_header = '==>Warning from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
character(len=256), parameter :: note_header = '==>Note from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
! variables for tracer package
integer :: ind
real, dimension(2) :: range_array
character(len=64) :: caller_str
character(len=fm_string_len) :: string_fm
character(len=fm_type_name_len) :: typ
character(len=fm_string_len), pointer, dimension(:) :: good_list
integer :: stdoutunit, stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if (prog_module_initialized) then
call mpp_error(FATAL, trim(error_header) // ' Prognostic tracers already initialized')
endif
nullify(T_prog)
write( stdlogunit,'(/a/)') trim(version)
have_obc = obc
dtts = Time_steps%dtts
dtuv = Time_steps%dtuv
cp_oceanr = 1.0/cp_ocean
! provide for namelist over-ride
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_tracer_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_tracer_nml')
#else
ioun = open_namelist_file()
read (ioun, ocean_tracer_nml,iostat=io_status)
ierr = check_nml_error(io_status,'ocean_tracer_nml')
call close_file (ioun)
#endif
write (stdoutunit,'(/)')
write (stdoutunit, ocean_tracer_nml)
write (stdlogunit, ocean_tracer_nml)
if (PRESENT(debug) .and. .not. debug_this_module) then
debug_this_module = debug
endif
if(debug_this_module) then
write(stdoutunit,'(a)') '==>Note: running with debug_this_module=.true. so will print lots of checksums.'
endif
if(cmip_units .and. cmip_version <= 5) then
cmip_offset = kelvin
temp_units = 'K'
else
cmip_offset = 0.0
temp_units = 'degrees C'
endif
if(zero_tendency) then
call mpp_error(NOTE, trim(note_header) // ' zero_tendency=true so will not time step tracer fields.')
Ocean_options%tracer_tendency = 'Did NOT time step prognostic tracer fields.'
else
Ocean_options%tracer_tendency = 'Time stepped the prognostic tracer fields.'
endif
if(zero_tracer_source) then
call mpp_error(NOTE, &
trim(note_header) // ' zero_tracer_source=true so remove T_prog%source from evolution.')
endif
if(.not. write_a_restart) then
write(stdoutunit,'(a)') '==>Note: running ocean_tracer with write_a_restart=.false.'
write(stdoutunit,'(a)') ' Will NOT write restart file, and so cannot restart the run.'
endif
if(frazil_heating_before_vphysics) then
frazil_heating_order=frazil_heating_order+1
write(stdoutunit,'(a)') '==>Note: frazil heating called before vertical physics and before boundary fluxes.'
write(stdoutunit,'(a)') ' This method is retained for legacy purposes: it is NOT recommended for new runs. '
write(stdoutunit,'(a)') ' '
endif
if(frazil_heating_after_vphysics) then
frazil_heating_order=frazil_heating_order+1
write(stdoutunit,'(a)') '==>Note: frazil heating called after vertical physics and after boundary fluxes.'
write(stdoutunit,'(a)') ' This is the recommended method. '
write(stdoutunit,'(a)') ' '
endif
if(frazil_heating_order>1) then
write(stdoutunit,'(a)') '==>Error from ocean_tracer_mod: choose just one temporal order for frazil heating.'
write(stdoutunit,'(a)') ' '
call mpp_error(FATAL, &
trim(note_header) // ' can only choose one temporal order for frazil heating.')
endif
if(frazil_heating_order==0) then
write(stdoutunit,'(a)') '==>Error from ocean_tracer_mod: MUST specify order for frazil heating in ocean_tracer.'
write(stdoutunit,'(a)') ' '
call mpp_error(FATAL, &
trim(note_header) // ' Need to specify an order for calling frazil heating from ocean_tracer_mod.')
endif
! some time step information
if (dtts /= dtuv) then
write (stdoutunit,'(/a)') trim(warn_header) // ' Asynchronous timesteps (dtts > dtuv) imply inaccurate transients.'
write (stdoutunit,'(a)') ' and total tracer (i.e. heat content) is not conserved.'
write (stdoutunit,'(a,f5.2,a/)') ' dtts =',dtts/dtuv,' times larger than dtuv.'
else
call mpp_error(NOTE, trim(note_header) // ' Synchronous timesteps have been specified (dtts = dtuv).')
endif
tendency = Time_steps%tendency
dtime = Time_steps%dtime_t
dtimer = 1.0/(dtime+epsln)
tau = Time%tau
taum1 = Time%taum1
taup1 = Time%taup1
! initialize clock ids
id_clock_bih_tracer = mpp_clock_id('(Ocean tracer: bih tracer) ' ,grain=CLOCK_MODULE)
id_clock_lap_tracer = mpp_clock_id('(Ocean tracer: lap tracer) ' ,grain=CLOCK_MODULE)
id_clock_vert_diffuse = mpp_clock_id('(Ocean tracer: vert diffuse) ' ,grain=CLOCK_MODULE)
id_clock_vert_diffuse_implicit = mpp_clock_id('(Ocean tracer: vert diffuse impl) ',grain=CLOCK_MODULE)
id_clock_tracer_advect = mpp_clock_id('(Ocean tracer: advection) ' ,grain=CLOCK_MODULE)
id_clock_frazil = mpp_clock_id('(Ocean tracer: frazil) ' ,grain=CLOCK_MODULE)
id_clock_convection = mpp_clock_id('(Ocean tracer: convection) ' ,grain=CLOCK_MODULE)
id_clock_blob = mpp_clock_id('(Ocean tracer: Lagrangian blobs) ' ,grain=CLOCK_MODULE)
id_clock_blob_implicit = mpp_clock_id('(Ocean tracer: implicit blobs)' ,grain=CLOCK_ROUTINE)
id_clock_adjust_L_thickness = mpp_clock_id('(Ocean tracer: adjust the L thickness)',grain=CLOCK_ROUTINE)
id_clock_adjust_E_thickness = mpp_clock_id('(Ocean tracer: adjust the E thickness)',grain=CLOCK_ROUTINE)
! perform requested debugging for ocean tracer package manager
if (ocean_tpm_debug) then
write (stdoutunit,*) ' '
write (stdoutunit,*) 'Dumping field tree at start of ocean_prog_tracer_init'
call write_timestamp(Time%model_time)
if (.not. fm_dump_list('/', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer tree')
endif
endif
! call routine to initialize the required tracer package
if (otpm_set_tracer_package('required', caller=trim(mod_name)//'('//trim(sub_name)//')') .le. 0) then
call mpp_error(FATAL, trim(error_header) // ' Could not set required packages list')
endif
! perform requested debugging for ocean tracer package manager
if (ocean_tpm_debug) then !{
write (stdoutunit,*) ' '
write (stdoutunit,*) 'Dumping /ocean_mod field tree before ocean_tpm_init'
call write_timestamp(Time%model_time)
if (.not. fm_dump_list('/ocean_mod', recursive = .true.)) then !{
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer tree')
endif !}
endif !}
! call the initialization routine for the ocean tracer packages
call ocean_tpm_init(Domain, Grid, Time, Time_steps, &
Ocean_options, debug)
! perform requested debugging for ocean tracer package manager
if (ocean_tpm_debug) then
write (stdoutunit,*) ' '
write (stdoutunit,*) 'Dumping /ocean_mod field tree after ocean_tpm_init'
call write_timestamp(Time%model_time)
if (.not. fm_dump_list('/ocean_mod', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer tree')
endif
endif
! check for any errors in the number of fields in the tracer_packages list
good_list => fm_util_get_string_array('/ocean_mod/GOOD/good_tracer_packages', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then
call fm_util_check_for_bad_fields('/ocean_mod/tracer_packages', good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
else
call mpp_error(FATAL,trim(error_header) // ' Empty "good_tracer_packages" list')
endif
! check for any errors in the number of fields in the prog_tracers list
good_list => fm_util_get_string_array('/ocean_mod/GOOD/good_prog_tracers', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then
call fm_util_check_for_bad_fields('/ocean_mod/prog_tracers', good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
else
call mpp_error(FATAL,trim(error_header) // ' Empty "good_prog_tracers" list')
endif
! check for any errors in the number of fields in the namelists list
good_list => fm_util_get_string_array('/ocean_mod/GOOD/good_namelists', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then
call fm_util_check_for_bad_fields('/ocean_mod/namelists', good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
!else
!call mpp_error(FATAL,trim(error_header) // ' Empty "good_namelists" list')
endif
! get the number of tracers
! for now, we will not try to dynamically allocate the tracer arrays here
num_prog_tracers = fm_get_length('/ocean_mod/prog_tracers')
num_diag = fm_get_length('/ocean_mod/diag_tracers')
num_prog=num_prog_tracers
! allocate arrays based on the number of prognostic tracers
allocate( T_prog (num_prog_tracers) )
allocate( id_eta_smooth (num_prog_tracers) )
allocate( id_eta_smooth_on_nrho(num_prog_tracers) )
allocate( id_pbot_smooth (num_prog_tracers) )
allocate( id_prog (num_prog_tracers) )
allocate( id_prog_explicit (num_prog_tracers) )
allocate( id_prog_rhodzt (num_prog_tracers) )
allocate( id_prog_int_rhodz (num_prog_tracers) )
allocate( id_prog_on_depth (num_prog_tracers) )
allocate( id_tendency_conc (num_prog_tracers) )
allocate( id_tendency (num_prog_tracers) )
allocate( id_tendency_on_nrho(num_prog_tracers) )
allocate( id_tendency_expl (num_prog_tracers) )
allocate( id_surf_tracer (num_prog_tracers) )
allocate( id_surf_tracer_sq (num_prog_tracers) )
allocate( id_bott_tracer (num_prog_tracers) )
if (use_blobs) then
allocate( id_progT (num_prog_tracers) )
allocate( id_prog_rhodztT (num_prog_tracers) )
allocate( id_prog_int_rhodzT(num_prog_tracers) )
allocate( id_surf_tracerT (num_prog_tracers) )
allocate( id_bott_tracerT (num_prog_tracers) )
allocate( id_tendency_concL (num_prog_tracers) )
allocate( id_tendency_concT (num_prog_tracers) )
allocate( id_tendencyL (num_prog_tracers) )
allocate( id_tendencyT (num_prog_tracers) )
allocate( id_tendency_explL (num_prog_tracers) )
allocate( id_tendency_explT (num_prog_tracers) )
endif
allocate( Tracer_restart (num_prog+num_diag) )
allocate( id_tracer_restart (num_prog+num_diag) )
allocate( tracer_restart_file (num_prog+num_diag) )
allocate( id_tracer_type (num_prog+num_diag) )
id_eta_smooth(:) = -1
id_eta_smooth_on_nrho(:)= -1
id_pbot_smooth(:) = -1
id_prog(:) = -1
id_prog_explicit(:) = -1
id_prog_rhodzt(:) = -1
id_prog_int_rhodz(:) = -1
id_prog_on_depth(:) = -1
id_tendency_conc(:) = -1
id_tendency(:) = -1
id_tendency_on_nrho(:)= -1
id_tendency_expl(:) = -1
id_surf_tracer(:) = -1
id_surf_tracer_sq(:) = -1
id_bott_tracer(:) = -1
if (use_blobs) then
id_progT(:) = -1
id_prog_rhodztT(:) = -1
id_prog_int_rhodzT(:)= -1
id_surf_tracerT(:) = -1
id_bott_tracerT(:) = -1
id_tendency_concL(:) = -1
id_tendency_concT(:) = -1
id_tendencyL(:) = -1
id_tendencyT(:) = -1
id_tendency_explL(:) = -1
id_tendency_explT(:) = -1
endif
! set logical to determine when to mpp_update tracers
do n=1,num_prog_tracers-1
T_prog(n)%complete=.false.
enddo
T_prog(num_prog_tracers)%complete=.true.
! dump the lists for the tracer packages
write (stdoutunit,*)
write (stdoutunit,*) 'Dumping tracer_packages tracer tree'
if (.not. fm_dump_list('/ocean_mod/tracer_packages', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer_packages tracer tree')
endif
write (stdoutunit,*)
write (stdoutunit,*) 'Dumping prog_tracers tracer tree'
if (.not. fm_dump_list('/ocean_mod/prog_tracers', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping prog_tracers tracer tree')
endif
write (stdoutunit,*)
write (stdoutunit,*) 'Dumping namelists tracer tree'
if (.not. fm_dump_list('/ocean_mod/namelists', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping namelists tracer tree')
endif
!### finished with initializing t_prog arrays
!### finished with call to tracer startup routine
! set local array indices
Grd => Grid
Dom => Domain
#ifndef MOM_STATIC_ARRAYS
call get_local_indices(Dom, isd, ied, jsd, jed, isc, iec, jsc, jec)
nk=Grd%nk
#endif
do n=1,num_prog_tracers
#ifndef MOM_STATIC_ARRAYS
allocate( T_prog(n)%field(isd:ied,jsd:jed,nk,3))
allocate( T_prog(n)%th_tendency(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%tendency(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%source(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%eta_smooth(isd:ied,jsd:jed))
allocate( T_prog(n)%pbot_smooth(isd:ied,jsd:jed))
allocate( T_prog(n)%wrk1(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%tmask_limit(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%K33_implicit(isd:ied,jsd:jed,nk))
if (use_blobs) then
allocate( T_prog(n)%fieldT(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%sum_blob(isd:ied,jsd:jed,nk,3))
allocate( T_prog(n)%tend_blob(isd:ied,jsd:jed,nk))
endif
#endif
T_prog(n)%field(:,:,:,:) = 0.0
T_prog(n)%th_tendency(:,:,:) = 0.0
T_prog(n)%tendency(:,:,:) = 0.0
T_prog(n)%source(:,:,:) = 0.0
T_prog(n)%eta_smooth(:,:) = 0.0
T_prog(n)%pbot_smooth(:,:) = 0.0
T_prog(n)%wrk1(:,:,:) = 0.0
T_prog(n)%tmask_limit(:,:,:) = 0.0
T_prog(n)%K33_implicit(:,:,:) = 0.0
T_prog(n)%neutral_physics_limit = .false.
if (use_blobs) then
!note that fieldT is initialised in ocean_blob_init
T_prog(n)%fieldT(:,:,:) = 0.0
T_prog(n)%sum_blob(:,:,:,:) = 0.0
T_prog(n)%tend_blob(:,:,:) = 0.0
endif
enddo
if(inflow_nboundary) then
call inflow_nboundary_init
endif
! for global surface area normalization
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
! fill the field table entries for the prognostic tracers
n = 0
do while (fm_loop_over_list('/ocean_mod/prog_tracers', name, typ, ind)) !{
if (typ .ne. 'list') then !{
call mpp_error(FATAL, trim(error_header) // ' ' // trim(name) // ' is not a list')
else !}{
n = n + 1 ! increment the array index
if (n .ne. ind) then !{
write (stdoutunit,*) trim(warn_header), ' Tracer index, ', ind, &
' does not match array index, ', n, ' for ', trim(name)
endif !}
! save the name
T_prog(n)%name = name
if (.not. fm_change_list('/ocean_mod/prog_tracers/' // trim(name))) then !{
call mpp_error(FATAL, trim(error_header) // ' Problem changing to ' // trim(name))
endif !}
caller_str = 'ocean_tracer_mod(ocean_prog_tracer_init)'
! save the units
T_prog(n)%units = fm_util_get_string('units', caller = caller_str, scalar = .true.)
! save the type
T_prog(n)%type = fm_util_get_string('type', caller = caller_str, scalar = .true.)
! save the longname
T_prog(n)%longname = fm_util_get_string('longname', caller = caller_str, scalar = .true.)
! save the conversion
T_prog(n)%conversion = fm_util_get_real('conversion', caller = caller_str, scalar = .true.)
! save the offset
T_prog(n)%offset = fm_util_get_real('offset', caller = caller_str, scalar = .true.)
! get the min and max of the tracer
T_prog(n)%min_tracer = fm_util_get_real('min_tracer', caller = caller_str, scalar = .true.)
T_prog(n)%max_tracer = fm_util_get_real('max_tracer', caller = caller_str, scalar = .true.)
! get the min and max of the range for analysis
T_prog(n)%min_range = fm_util_get_real('min_range', caller = caller_str, scalar = .true.)
T_prog(n)%max_range = fm_util_get_real('max_range', caller = caller_str, scalar = .true.)
! get the flux unit
T_prog(n)%flux_units = fm_util_get_string('flux_units', caller = caller_str, scalar = .true.)
! get the min and max of the flux range for analysis
T_prog(n)%min_flux_range = fm_util_get_real('min_flux_range', caller = caller_str, scalar = .true.)
T_prog(n)%max_flux_range = fm_util_get_real('max_flux_range', caller = caller_str, scalar = .true.)
! save the restart file
T_prog(n)%restart_file = fm_util_get_string('restart_file', caller = caller_str, scalar = .true.)
! save flag for whether the tracer must have a value in the restart file
T_prog(n)%const_init_tracer = fm_util_get_logical('const_init_tracer', caller = caller_str, scalar = .true.)
! save value to globally initialize this tracer (optional)
T_prog(n)%const_init_value = fm_util_get_real('const_init_value', caller = caller_str, scalar = .true.)
! get the horizontal-advection-scheme
string_fm = fm_util_get_string('horizontal-advection-scheme', caller = caller_str, scalar = .true.)
select case (trim(string_fm))
case ('upwind')
T_prog(n)%horz_advect_scheme = ADVECT_UPWIND
case ('2nd_order')
T_prog(n)%horz_advect_scheme = ADVECT_2ND_ORDER
case ('4th_order')
T_prog(n)%horz_advect_scheme = ADVECT_4TH_ORDER
case ('quicker')
T_prog(n)%horz_advect_scheme = ADVECT_QUICKER
case ('quickMOM3')
T_prog(n)%horz_advect_scheme = ADVECT_QUICKMOM3
case ('6th_order')
T_prog(n)%horz_advect_scheme = ADVECT_6TH_ORDER
case ('mdfl_sup_b')
T_prog(n)%horz_advect_scheme = ADVECT_MDFL_SUP_B
case ('mdfl_sweby')
T_prog(n)%horz_advect_scheme = ADVECT_MDFL_SWEBY
case ('mdfl_sweby_test')
T_prog(n)%horz_advect_scheme = ADVECT_MDFL_SWEBY_TEST
case ('dst_linear')
T_prog(n)%horz_advect_scheme = ADVECT_DST_LINEAR
case ('dst_linear_test')
T_prog(n)%horz_advect_scheme = ADVECT_DST_LINEAR_TEST
case ('mdppm')
T_prog(n)%horz_advect_scheme = ADVECT_MDPPM
T_prog(n)%ppm_hlimiter = fm_util_get_integer('ppm_hlimiter', caller = caller_str, scalar = .true.)
case ('mdppm_test')
T_prog(n)%horz_advect_scheme = ADVECT_MDPPM_TEST
T_prog(n)%ppm_hlimiter = fm_util_get_integer('ppm_hlimiter', caller = caller_str, scalar = .true.)
case ('mdmdt_test')
T_prog(n)%horz_advect_scheme = ADVECT_MDMDT_TEST
T_prog(n)%mdt_scheme = fm_util_get_integer('mdt_scheme', caller = caller_str, scalar = .true.)
case ('psom')
T_prog(n)%horz_advect_scheme = ADVECT_PSOM
! set psom_limit when using psom
T_prog(n)%psom_limit = fm_util_get_logical('psom_limit', caller = caller_str, scalar = .true.)
case default
call mpp_error(FATAL, trim(error_header) // ' Invalid horz-advect-scheme ' // trim(string_fm) // ' for ' // trim(name))
end select
! get the vertical-advection-scheme
string_fm = fm_util_get_string('vertical-advection-scheme', caller = caller_str, scalar = .true.)
select case (trim(string_fm))
case ('upwind')
T_prog(n)%vert_advect_scheme = ADVECT_UPWIND
case ('2nd_order')
T_prog(n)%vert_advect_scheme = ADVECT_2ND_ORDER
case ('4th_order')
T_prog(n)%vert_advect_scheme = ADVECT_4TH_ORDER
case ('6th_order')
T_prog(n)%vert_advect_scheme = ADVECT_6TH_ORDER
case ('quicker')
T_prog(n)%vert_advect_scheme = ADVECT_QUICKER
case ('quickMOM3')
T_prog(n)%vert_advect_scheme = ADVECT_QUICKMOM3
case ('mdfl_sup_b')
T_prog(n)%vert_advect_scheme = ADVECT_MDFL_SUP_B
case ('mdfl_sweby')
T_prog(n)%vert_advect_scheme = ADVECT_MDFL_SWEBY
case ('mdfl_sweby_test')
T_prog(n)%vert_advect_scheme = ADVECT_MDFL_SWEBY_TEST
case ('dst_linear')
T_prog(n)%vert_advect_scheme = ADVECT_DST_LINEAR
case ('dst_linear_test')
T_prog(n)%vert_advect_scheme = ADVECT_DST_LINEAR_TEST
case ('psom')
T_prog(n)%vert_advect_scheme = ADVECT_PSOM
case ('mdppm')
T_prog(n)%vert_advect_scheme = ADVECT_MDPPM
T_prog(n)%ppm_vlimiter = fm_util_get_integer('ppm_vlimiter', caller = caller_str, scalar = .true.)
case ('mdppm_test')
T_prog(n)%vert_advect_scheme = ADVECT_MDPPM_TEST
T_prog(n)%ppm_vlimiter = fm_util_get_integer('ppm_vlimiter', caller = caller_str, scalar = .true.)
case ('mdmdt_test')
T_prog(n)%vert_advect_scheme = ADVECT_MDMDT_TEST
T_prog(n)%mdt_scheme = fm_util_get_integer('mdt_scheme', caller = caller_str, scalar = .true.)
case default
call mpp_error(FATAL, trim(error_header) // ' Invalid vert-advect scheme ' // trim(string_fm) // ' for ' // trim(name))
end select
! save the max_tracer_limit
T_prog(n)%max_tracer_limit = fm_util_get_real('max_tracer_limit', caller = caller_str, scalar = .true.)
! save the min_tracer_limit
T_prog(n)%min_tracer_limit = fm_util_get_real('min_tracer_limit', caller = caller_str, scalar = .true.)
! save flag for whether to transport tracer using only advection
T_prog(n)%use_only_advection = fm_util_get_logical('use_only_advection', caller = caller_str, scalar = .true.)
if(T_prog(n)%use_only_advection) then
call mpp_error(NOTE, &
trim(note_header) // ' Will evolve '// trim(T_prog(n)%name) // ' via advection alone.')
endif
! check consistency in selection of three-dimensional advection schemes
if(T_prog(n)%horz_advect_scheme==ADVECT_MDFL_SUP_B .and. &
T_prog(n)%vert_advect_scheme /= ADVECT_MDFL_SUP_B) then
call mpp_error(FATAL,&
trim(error_header) // ' mdfl_sup_b advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%horz_advect_scheme==ADVECT_MDFL_SWEBY .and. &
T_prog(n)%vert_advect_scheme /= ADVECT_MDFL_SWEBY) then
call mpp_error(FATAL,&
trim(error_header) // ' mdfl_sweby advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%horz_advect_scheme==ADVECT_DST_LINEAR .and. &
T_prog(n)%vert_advect_scheme /= ADVECT_DST_LINEAR) then
call mpp_error(FATAL,&
trim(error_header) // ' dst_linear advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%horz_advect_scheme==ADVECT_DST_LINEAR_TEST .and. &
T_prog(n)%vert_advect_scheme /= ADVECT_DST_LINEAR_TEST) then
call mpp_error(FATAL,&
trim(error_header) // ' dst_linear_legacy advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_MDFL_SUP_B .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_MDFL_SUP_B) then
call mpp_error(FATAL,&
trim(error_header) // ' mdfl_sup_b advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_MDFL_SWEBY .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_MDFL_SWEBY) then
call mpp_error(FATAL,&
trim(error_header) // ' mdfl_sweby advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_MDFL_SWEBY_TEST .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_MDFL_SWEBY_TEST) then
call mpp_error(FATAL,&
trim(error_header) // ' mdfl_sweby_test advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_DST_LINEAR .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_DST_LINEAR) then
call mpp_error(FATAL,&
trim(error_header) // ' dst_linear advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_DST_LINEAR_TEST .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_DST_LINEAR_TEST) then
call mpp_error(FATAL,&
trim(error_header) // ' dst_linear_test advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_PSOM .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_PSOM) then
call mpp_error(FATAL,&
trim(error_header) // ' psom advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_MDPPM .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_MDPPM) then
call mpp_error(FATAL,&
trim(error_header) // ' mdppm advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_MDPPM_TEST .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_MDPPM_TEST) then
call mpp_error(FATAL,&
trim(error_header) // ' mdppm_test advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_MDMDT_TEST .and. &
T_prog(n)%horz_advect_scheme /= ADVECT_MDMDT_TEST) then
call mpp_error(FATAL,&
trim(error_header) // ' mdmdt_test advection for ' // trim(name) // 'must be for **BOTH** horz & vert')
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_2ND_ORDER .or. &
T_prog(n)%horz_advect_scheme==ADVECT_2ND_ORDER) then
if(tendency==TWO_LEVEL) then
call mpp_error(FATAL,&
trim(error_header) // ' 2nd_order advection for ' // trim(name) // 'is unstable with TWO_LEVEL time')
endif
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_4TH_ORDER .or. &
T_prog(n)%horz_advect_scheme==ADVECT_4TH_ORDER) then
if(tendency==TWO_LEVEL) then
call mpp_error(FATAL,&
trim(error_header) // ' 4th_order advection for ' // trim(name) // 'is unstable with TWO_LEVEL time')
endif
endif
if(T_prog(n)%vert_advect_scheme==ADVECT_6TH_ORDER .or. &
T_prog(n)%horz_advect_scheme==ADVECT_6TH_ORDER) then
if(tendency==TWO_LEVEL) then
call mpp_error(FATAL,&
trim(error_header) // ' 6th_order advection for ' // trim(name) // 'is unstable with TWO_LEVEL time')
endif
endif
endif !}
enddo !}
allocate(id_tmask_limit(num_prog_tracers))
do n=1,num_prog_tracers
id_tmask_limit(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_tmask_limit', &
Grd%tracer_axes(1:3), Time%model_time, &
'use upwind (not quicker) and horz diff (not neutral)', &
'none', missing_value=missing_value, range=(/-2.0,2.0/))
if (use_blobs) then
id_prog_rhodzt(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_rhodztE', &
Grd%tracer_axes(1:3), &
Time%model_time, 'E system '//trim(T_prog(n)%longname)//' * rho_dzt', &
trim(T_prog(n)%units)//'*(kg/m^3)*m', &
missing_value=missing_value, range=(/-1e6,1e12/))
id_prog_int_rhodz(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_int_rhodzE', &
Grd%tracer_axes(1:2), &
Time%model_time, 'E system vertical sum of '//trim(T_prog(n)%longname)//' * rho_dzt', &
trim(T_prog(n)%units)//'*(kg/m^3)*m', &
missing_value=missing_value, range=(/-1e20,1e20/))
else !not use_blobs
id_prog_rhodzt(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_rhodzt', &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname)//' * rho_dzt', &
trim(T_prog(n)%units)//'*(kg/m^3)*m', &
missing_value=missing_value, range=(/-1e6,1e12/))
id_prog_int_rhodz(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_int_rhodz', &
Grd%tracer_axes(1:2), &
Time%model_time, 'vertical sum of '//trim(T_prog(n)%longname)//' * rho_dzt', &
trim(T_prog(n)%units)//'*(kg/m^3)*m', &
missing_value=missing_value, range=(/-1e20,1e20/))
endif
range_array(1) = T_prog(n)%min_range
range_array(2) = T_prog(n)%max_range
if (use_blobs) then
prog_name = trim(T_prog(n)%name)//'E'
prog_longname = 'E system '//trim(T_prog(n)%longname)
else
prog_name = T_prog(n)%name
prog_longname = T_prog(n)%longname
endif
if(T_prog(n)%longname=='Potential temperature') then
id_prog(n) = register_diag_field ('ocean_model', &
trim(prog_name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(prog_longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_potential_temperature')
id_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(prog_longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_surface_temperature')
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(prog_longname), &
'squared '//trim(temp_units), &
missing_value=missing_value, range=(/0.0,1e10/), &
standard_name='square_of_sea_surface_temperature')
elseif(T_prog(n)%longname=='Conservative temperature') then
id_prog(n) = register_diag_field ('ocean_model', &
trim(prog_name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(prog_longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_conservative_temperature')
id_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(prog_longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_surface_conservative_temperature')
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(prog_longname), &
'squared '//trim(temp_units), &
missing_value=missing_value, range=(/0.0,1e10/), &
standard_name='square_of_sea_surface_conservative_temperature')
elseif(T_prog(n)%name=='salt') then
id_prog(n) = register_diag_field ('ocean_model', &
trim(prog_name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_salinity')
id_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_surface_salinity')
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(prog_longname), &
'squared '//trim(T_prog(n)%units), &
missing_value=missing_value, range=(/0.0,1e10/), &
standard_name='square_of_sea_surface_salinity')
elseif(T_prog(n)%name=='age_global') then
id_prog(n) = register_diag_field ('ocean_model', &
trim(prog_name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_age_since_surface_contact')
id_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(prog_longname), &
'squared '//trim(T_prog(n)%units), &
missing_value=missing_value, range=(/0.0,1e20/))
elseif(T_prog(n)%name=='cfc_11') then
id_prog(n) = register_diag_field ('ocean_model', &
trim(prog_name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='moles_per_unit_mass_of_cfc11_in_sea_water')
id_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(prog_longname), &
'squared '//trim(T_prog(n)%units), &
missing_value=missing_value, range=(/0.0,1e20/))
else
id_prog(n) = register_diag_field ('ocean_model', trim(prog_name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(prog_longname), trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
id_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(prog_longname), &
'squared '//trim(T_prog(n)%units), &
missing_value=missing_value, range=(/0.0,1e20/))
endif
id_prog_explicit(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_explicit', &
Grd%tracer_axes(1:3), Time%model_time, &
trim(T_prog(n)%longname)//' at taup1 updated just with explicit tendencies',&
trim(T_prog(n)%units), missing_value=missing_value, range=range_array)
id_prog_on_depth(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_on_depth', Grd%tracer_axes_depth(1:3), &
Time%model_time, trim(T_prog(n)%longname), trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
id_bott_tracer(n) = register_diag_field ('ocean_model', &
'bottom_'//trim(prog_name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(prog_longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
! register diagnostics for some tracer tendency terms
if (T_prog(n)%max_flux_range .gt. T_prog(n)%min_flux_range) then
range_array(1) = T_prog(n)%min_flux_range
range_array(2) = T_prog(n)%max_flux_range
id_tendency_conc(n) = register_diag_field ('ocean_model', &
trim(prog_name)//'_tendency_conc', Grd%tracer_axes(1:3), &
Time%model_time, 'time tendency of tracer concentration for '//trim(prog_longname),&
trim(T_prog(n)%units)//' per second', missing_value=missing_value, range=range_array)
id_tendency(n) = register_diag_field ('ocean_model', &
trim(prog_name)//'_tendency', Grd%tracer_axes(1:3), &
Time%model_time, 'time tendency for tracer '//trim(prog_longname), &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=range_array)
id_eta_smooth(n) = register_diag_field ('ocean_model', &
trim(prog_name)//'_eta_smooth', Grd%tracer_axes(1:2), &
Time%model_time, 'surface smoother for ' // trim(prog_name), &
trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=range_array)
id_pbot_smooth(n) = register_diag_field ('ocean_model', &
trim(prog_name)//'_pbot_smooth', Grd%tracer_axes(1 :2), &
Time%model_time, 'bottom smoother for ' // trim(prog_name), &
trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=range_array)
! The explicit tendency is the total system, whether blobs are used or not.
! So, we make no distinction between T and E systems.
id_tendency_expl(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_tendency_expl', Grd%tracer_axes(1:3), &
Time%model_time, 'explicit in time tendency for tracer ' // trim(T_prog(n)%longname), &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=range_array)
else
id_tendency_conc(n) = register_diag_field ('ocean_model', &
trim(prog_name)//'_tendency_conc', Grd%tracer_axes(1:3), &
Time%model_time, 'tendency of tracer concentration for '//trim(prog_longname),&
trim(T_prog(n)%units)//' per second', missing_value=missing_value, range=(/-1.e10,1e10/))
id_tendency(n) = register_diag_field ('ocean_model', &
trim(prog_name)//'_tendency', Grd%tracer_axes(1:3), &
Time%model_time, 'time tendency for tracer '//trim(prog_longname), &
trim(T_prog(n)%flux_units), missing_value=missing_value)
id_eta_smooth(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_eta_smooth', Grd%tracer_axes(1:2), &
Time%model_time, 'surface smoother for ' // trim(T_prog(n)%name), &
trim(T_prog(n)%flux_units), &
missing_value=missing_value)
id_pbot_smooth(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_pbot_smooth', Grd%tracer_axes(1:2), &
Time%model_time, 'bottom smoother for ' // trim(T_prog(n)%name), &
trim(T_prog(n)%flux_units), &
missing_value=missing_value)
! The explicit tendency is the total system, whether blobs are used or not.
! So, we make no distinction between T and E systems.
id_tendency_expl(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_tendency_expl', Grd%tracer_axes(1:3), &
Time%model_time, 'explicit in time tendency for ' // trim(T_prog(n)%longname), &
trim(T_prog(n)%flux_units), missing_value=missing_value)
endif
!Diagnostics for the Lagrangian blob framework
if (use_blobs) then
id_prog_rhodztT(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_rhodzt',&
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname)//' * rho_dzt', &
trim(T_prog(n)%units)//'*(kg/m^3)*m', &
missing_value=missing_value, range=(/-1e6,1e12/))
id_prog_int_rhodzT(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_int_rhodz',&
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname)//' * rho_dzt', &
trim(T_prog(n)%units)//'*(kg/m^3)*m', &
missing_value=missing_value, range=(/-1e20,1e20/))
if(T_prog(n)%longname=='Potential temperature') then
id_progT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_potential_temperature')
id_surf_tracerT(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_surface_temperature')
elseif(T_prog(n)%longname=='Conservative temperature') then
id_progT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_conservative_temperature')
id_surf_tracerT(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname), &
trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_surface_conservative_temperature')
elseif(T_prog(n)%name=='salt') then
id_progT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_salinity')
id_surf_tracerT(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_surface_salinity')
elseif(T_prog(n)%name=='age_global') then
id_progT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_age_since_surface_contact')
id_surf_tracerT(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
elseif(T_prog(n)%name=='cfc_11') then
id_progT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='moles_per_unit_mass_of_cfc11_in_sea_water')
id_surf_tracerT(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
else
id_progT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
id_surf_tracerT(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
endif
id_bott_tracerT(n) = register_diag_field ('ocean_model', &
'bottom_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_prog(n)%longname), &
trim(T_prog(n)%units), &
missing_value=missing_value, range=range_array)
id_tendency_concT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_tendency_conc', Grd%tracer_axes(1:3), &
Time%model_time, 'tendency of tracer concentration for '//trim(T_prog(n)%longname),&
trim(T_prog(n)%units)//' per second', missing_value=missing_value, range=(/-1.e10,1e10/))
id_tendency_concL(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_tendency_concL', Grd%tracer_axes(1:3), &
Time%model_time, 'L contribution to tendency of tracer concentration for '//trim(T_prog(n)%longname),&
trim(T_prog(n)%units)//' per second', missing_value=missing_value, range=(/-1.e10,1e10/))
id_tendencyT(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_tendency', Grd%tracer_axes(1:3), &
Time%model_time, 'time tendency for tracer '//trim(T_prog(n)%longname), &
trim(T_prog(n)%flux_units), missing_value=missing_value)
id_tendencyL(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'L_tendency', Grd%tracer_axes(1:3), &
Time%model_time, 'time tendency for tracer '//trim(T_prog(n)%longname)//'L', &
trim(T_prog(n)%flux_units), missing_value=missing_value)
endif ! endif for use_blobs
if (T_prog(n)%name == 'temp') index_temp = n
if (T_prog(n)%name == 'salt') index_salt = n
if (T_prog(n)%name == 'temp_sq') index_temp_sq = n
if (T_prog(n)%name == 'salt_sq') index_salt_sq = n
if (T_prog(n)%name == 'redist_heat') index_redist_heat = n
if (T_prog(n)%name == 'added_heat') index_added_heat = n
if (T_prog(n)%longname == 'Conservative temperature') prog_temp_variable=CONSERVATIVE_TEMP
if (T_prog(n)%longname == 'Potential temperature') prog_temp_variable=POTENTIAL_TEMP
enddo ! n-loop
if(prog_temp_variable==CONSERVATIVE_TEMP) then
write (stdoutunit,'(/1x,a)') &
' ==> Note from ocean_tracer_mod: prognostic temperature = conservative temperature.'
write (stdoutunit,'(1x,a/)') &
' diagnostic temperature = potential temperature.'
elseif(prog_temp_variable==POTENTIAL_TEMP) then
write (stdoutunit,'(/1x,a)') &
' ==> Note from ocean_tracer_mod: prognostic temperature = potential temperature.'
write (stdoutunit,'(1x,a/)') &
' diagnostic temperature = conservative temperature.'
else
call mpp_error(FATAL, &
'==>Error in ocean_tracer_mod: model temperature variable remains unspecified')
endif
if (index_temp == -1 .or. index_salt == -1 ) then
call mpp_error(FATAL, trim(error_header) // &
' temp or salt not included in table entries. MOM needs both.')
endif
#ifdef USE_OCEAN_BGC
!Point the %field4d_ptr and %tendency for "generic" tracers to the corresponding T_prog(n)%field and T_prog(n)%K33_implicit
!to utilize the already allocated memory in MOM and to avoid copying the arrays back and forth between generic tracers and MOM
do n=1,num_prog_tracers
if(T_prog(n)%type .eq. 'generic') then
call ocean_generic_set_pointer(T_prog(n)%name, 'field', T_prog(n)%field, isd, jsd)
!T_prog(n)%K33_implicit is used in vertdiff method below for calculating vertical diffusivity
call ocean_generic_set_pointer(T_prog(n)%name, 'tendency', T_prog(n)%K33_implicit, isd, jsd)
endif
enddo
#endif
! read prognostic tracer initial conditions or restarts
write (stdoutunit,*) ' '
write (stdoutunit,*) trim(note_header), &
' Reading prognostic tracer initial conditions or restarts'
do n=1,num_prog_tracers !{
write (stdoutunit,*)
write (stdoutunit,*) &
'Initializing tracer number', n,' at time level tau. This tracer is called ',trim(T_prog(n)%name)
if(.not. Time%init) then
if(tendency==TWO_LEVEL) then
write (stdoutunit,'(/a)')'Expecting only one time record from the tracer restart.'
elseif(tendency==THREE_LEVEL) then
write (stdoutunit,'(/a)')'Expecting two time records from the tracer restart.'
endif
endif
T_prog(n)%field(:,:,:,taum1) = 0.0
! initialization logic
!
! 1. Always call passive_tracer_init once.
! If Time%init=.false. then return w/o doing anything
! If not using passive tracers, then return w/o doing anything
!
! 2. If filename exists with field in it, then fill the tracer with field in filename.
!
! 3. If filename exists but there is no field in it, then fill the tracer with const_init_value
! if const_init_tracer is true, otherwise end with a fatal error.
!
! 4. If filename does not exist, then fill the tracer with const_init_value
! if const_init_tracer is true, otherwise end with a fatal error.
!
filename = T_prog(n)%restart_file
call passive_tracer_init(Time%init, T_prog(n), initialize_as_a_passive_tracer)
do l = 1, num_tracer_restart
if(trim(filename) == tracer_restart_file(l)) exit
end do
if(l>num_tracer_restart) then
num_tracer_restart = num_tracer_restart + 1
tracer_restart_file(l) = trim(filename)
end if
id_tracer_type(n) = l
if(tendency==THREE_LEVEL) then
id_tracer_restart(n) = register_restart_field(Tracer_restart(l), filename, T_prog(n)%name, &
T_prog(n)%field(:,:,:,tau), T_prog(n)%field(:,:,:,taup1), Domain%domain2d)
else
id_tracer_restart(n) = register_restart_field(Tracer_restart(l), filename, T_prog(n)%name, &
T_prog(n)%field(:,:,:,tau), Domain%domain2d)
end if
filename = 'INPUT/'//trim(T_prog(n)%restart_file)
if( Time%init .and. T_prog(n)%const_init_tracer ) then
write (stdoutunit,*) &
'Initializing the tracer ',trim(T_prog(n)%name),' to the constant ',T_prog(n)%const_init_value
T_prog(n)%field(:,:,:,1) = T_prog(n)%const_init_value*Grd%tmask(:,:,:)
T_prog(n)%field(:,:,:,2) = T_prog(n)%const_init_value*Grd%tmask(:,:,:)
T_prog(n)%field(:,:,:,3) = T_prog(n)%const_init_value*Grd%tmask(:,:,:)
elseif (field_exist(trim(filename), trim(T_prog(n)%name))) then
call field_size(filename,T_prog(n)%name, siz)
if (tendency==TWO_LEVEL .and. siz(4) > 1) then
write(stdoutunit,'(/a)') &
'==>WARMING: Attempt to read restart for tracer from a 3-level time scheme (2 time records)'
write(stdoutunit,'(a)') &
' when running MOM with 2-level timestepping (only need 1 time record in restart).'
write(stdoutunit,'(a)') &
' Reduce restart file to a single time record in order to avoid confusion.'
call mpp_error(WARNING, &
'Reading 3-time level ocean tracer (w/ 2 time records) while using 2-level (needs only 1 record)')
endif
write (stdoutunit,*) 'Reading restart for prog tracer ', trim(T_prog(n)%name), ' from file ', &
trim(T_prog(n)%restart_file)
call restore_state(Tracer_restart(l), id_tracer_restart(n))
! modify initial conditions
if (Time%init) then
if(remap_depth_to_s_init .and. vert_coordinate_type == QUASI_HORIZONTAL) then
call mpp_error(WARNING, trim(note_header) // &
'No need to remap init-conditions from depth to s-coord with quasi-horizontal vert-coord.')
endif
if(.not. remap_depth_to_s_init .and. vert_coordinate_type == TERRAIN_FOLLOWING) then
call mpp_error(WARNING, trim(note_header) // &
'May need to remap init-conditions from depth to s-coord with terrain following vert-coord.')
endif
if(remap_depth_to_s_init .and. vert_coordinate_type == TERRAIN_FOLLOWING) then
call mpp_error(NOTE, trim(note_header) // &
'Remap init-tracer from depth to s-coord for terrain following vert-coord.')
endif
! remap the depth-based initial conditions onto s-coordinates
if(remap_depth_to_s_init) then
write (stdoutunit,*) &
'After reading tracer init, remapping depth to s-coord for tracer ',trim(T_prog(n)%name)
call remap_depth_to_s(Thickness, Time, T_prog(n)%field(:,:,:,tau))
endif
! modifications due to partial bottom cell
if(vert_coordinate_type == QUASI_HORIZONTAL .and. interpolate_tprog_to_pbott) then
write (stdoutunit,*) &
'After reading ic, linearly interpolate ',trim(T_prog(n)%name),' to partial cell bottom.'
do j=jsc,jec
do i=isc,iec
kb = Grd%kmt(i,j)
if (kb .gt. 1) then
if(use_blobs) then
fact = Thickness%dzwtT(i,j,kb-1)/Grd%dzw(kb-1)
else
fact = Thickness%dzwt(i,j,kb-1)/Grd%dzw(kb-1)
endif
T_prog(n)%field(i,j,kb,tau) = &
T_prog(n)%field(i,j,kb-1,tau) - &
fact*(T_prog(n)%field(i,j,kb-1,tau) - &
T_prog(n)%field(i,j,kb,tau))
endif
! zero out tracers in land points
do k = kb+1,nk
T_prog(n)%field(i,j,k,tau) = 0.0
enddo
enddo
enddo
endif
endif
else
if(.not. initialize_as_a_passive_tracer) then
call mpp_error(FATAL, &
trim(T_prog(n)%name) // ' is not initialized as an ideal passive tracer, ' // &
' it is not initialized as constant, and it does not exist in the file ' &
// trim(filename) // &
'All tracers must have initialization specified. There is no default.')
endif
endif
write (stdoutunit,*) &
'Completed initialization of tracer ',trim(T_prog(n)%name),' at time level tau'
! fill halos for tau tracer value
call mpp_update_domains(T_prog(n)%field(:,:,:,tau), Dom%domain2d)
! read in second time level when running with threelevel scheme
if(tendency==THREE_LEVEL .and. .not.Time%init) then
write (stdoutunit,*) 'Since Time%init=.false. and using threelevel tendency, read taup1 for ',trim(T_prog(n)%name)
write (stdoutunit,*) 'Completed read of restart for ', trim(T_prog(n)%name),' at time taup1'
call mpp_update_domains(T_prog(n)%field(:,:,:,taup1), Dom%domain2d)
else
! initialize tracer at taup1 to tracer at tau
T_prog(n)%field(:,:,:,taup1) = T_prog(n)%field(:,:,:,tau)
endif
enddo !} n-loop
write (stdoutunit,*) ' '
write (stdoutunit,*) trim(note_header), ' finished reading prognostic tracer restarts.'
if(have_obc) call ocean_obc_tracer_init(Time, T_prog, num_prog_tracers, debug_this_module)
prog_module_initialized = .true.
if (debug_this_module) then
do n=1,num_prog_tracers
write(stdoutunit,'(a)') ' '
write(stdoutunit,*) 'From ocean_tracer_mod: initial ',trim(T_prog(n)%name), ' chksum (tau) ==>'
call tracer_prog_chksum(Time, T_prog(n), tau)
write(stdoutunit,*) 'From ocean_tracer_mod: initial ',trim(T_prog(n)%name), ' chksum (taup1) ==>'
call tracer_prog_chksum(Time, T_prog(n), taup1)
call tracer_min_max(Time, Thickness, T_prog(n))
enddo
endif
end function ocean_prog_tracer_init
! NAME="ocean_prog_tracer_init">
!#######################################################################
!
!
!
! Initialization code for diagnostic tracers, returning a pointer to the T_diag array
!
!
function ocean_diag_tracer_init (Time, Thickness, vert_coordinate_type, num_diag, use_blobs) &
result (T_diag) !{
type(ocean_time_type), intent(in), target :: Time
type(ocean_thickness_type), intent(in) :: Thickness
integer, intent(in) :: vert_coordinate_type
integer, intent(out) :: num_diag
logical, intent(in) :: use_blobs
!
! Return type
!
type(ocean_diag_tracer_type), dimension(:), pointer :: T_diag
integer :: kb, i, j, k, n, l
character(len=32) :: name
character(len=128) :: filename
real :: fact
character(len=48), parameter :: sub_name = 'ocean_diag_tracer_init'
character(len=256), parameter :: error_header = '==>Error from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
character(len=256), parameter :: warn_header = '==>Warning from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
character(len=256), parameter :: note_header = '==>Note from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
! variables for tracer package
integer :: ind, id_restart
real, dimension(2) :: range_array
character(len=64) :: caller_str
character(len=fm_type_name_len) :: typ
character(len=fm_string_len), pointer, dimension(:) :: good_list
integer :: stdoutunit
stdoutunit=stdout()
if (diag_module_initialized) then
call mpp_error(FATAL, trim(error_header) // ' Diagnostic tracers already initialized')
endif
nullify(T_diag)
! calls to internally generated (i.e., no field table entries) diagnostic tracers
call ocean_irradiance_init
!
! Check for any errors in the number of fields in the diag_tracers list
!
good_list => fm_util_get_string_array('/ocean_mod/GOOD/good_diag_tracers', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then !{
call fm_util_check_for_bad_fields('/ocean_mod/diag_tracers', good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
endif !}
! Get the number of tracers
! For now, we will not try to dynamically allocate the tracer arrays here
num_diag_tracers = fm_get_length('/ocean_mod/diag_tracers')
num_diag = num_diag_tracers
!
! Only do the following if there are diag tracers
!
if (num_diag_tracers .gt. 0) then !{
write (stdoutunit,*)
write (stdoutunit,*) trim(note_header), ' ', num_diag_tracers, ' diagnostic tracers requested.'
!
! Allocate the T_diag array
!
allocate( T_diag(num_diag_tracers) )
allocate( id_diag(num_diag_tracers) )
allocate( id_diag_total(num_diag_tracers) )
allocate( id_diag_surf_tracer (num_diag_tracers) )
allocate( id_diag_surf_tracer_sq(num_diag_tracers) )
id_diag(:) = -1
id_diag_total(:) = -1
id_diag_surf_tracer(:) = -1
id_diag_surf_tracer_sq(:) = -1
!
! Dump the diagnostic tracer tree
!
write (stdoutunit,*)
write (stdoutunit,*) 'Dumping ocean diag field tree after reading diag tracer tree'
if (.not. fm_dump_list('/ocean_mod/diag_tracers', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping diag_tracers tracer tree')
endif
!### finished with initializing t_diag arrays
!### finished with call to tracer startup routine
! allocate tracer type arrays and fill in tracer table information
if (num_diag_tracers .gt. 0 ) then
do n=1,num_diag_tracers
#ifndef MOM_STATIC_ARRAYS
allocate(T_diag(n)%field(isd:ied,jsd:jed,nk))
#endif
T_diag(n)%field(:,:,:) = 0.0
enddo
endif
!
! process the T_diag array
!
n = 0
do while (fm_loop_over_list('/ocean_mod/diag_tracers', name, typ, ind)) !{
if (typ .ne. 'list') then !{
call mpp_error(FATAL, trim(error_header) // ' ' // trim(name) // ' is not a list')
else !}{
n = n + 1
if (n .ne. ind) then
write (stdoutunit,*) trim(warn_header), ' Tracer index, ', ind, &
' does not match array index, ', n, ' for ', trim(name)
endif
! save the name (required)
T_diag(n)%name = name
if (.not. fm_change_list('/ocean_mod/diag_tracers/' // trim(name))) then
call mpp_error(FATAL, trim(error_header) // ' Problem changing to ' // trim(name))
endif
caller_str = 'ocean_tracer_mod(ocean_diag_tracer_init)'
! save the units (optional)
T_diag(n)%units = fm_util_get_string('units', caller = caller_str, scalar = .true.)
! save the type (optional)
T_diag(n)%type = fm_util_get_string('type', caller = caller_str, scalar = .true.)
! save the longname (optional)
T_diag(n)%longname = fm_util_get_string('longname', caller = caller_str, scalar = .true.)
! save the conversion (optional)
T_diag(n)%conversion = fm_util_get_real('conversion', caller = caller_str, scalar = .true.)
! save the offset (optional)
T_diag(n)%offset = fm_util_get_real('offset', caller = caller_str, scalar = .true.)
! get the min and max of the tracer for analysis (optional)
T_diag(n)%min_tracer = fm_util_get_real('min_tracer', caller = caller_str, scalar = .true.)
T_diag(n)%max_tracer = fm_util_get_real('max_tracer', caller = caller_str, scalar = .true.)
! get the min and max of the range for analysis (optional)
T_diag(n)%min_range = fm_util_get_real('min_range', caller = caller_str, scalar = .true.)
T_diag(n)%max_range = fm_util_get_real('max_range', caller = caller_str, scalar = .true.)
! save the input restart file (optional)
T_diag(n)%restart_file = fm_util_get_string('restart_file', caller = caller_str, scalar = .true.)
! save flag for whether to initialize this tracer (optional)
T_diag(n)%const_init_tracer = fm_util_get_logical('const_init_tracer', caller = caller_str, scalar = .true.)
! save value to globally initialize this tracer (optional)
T_diag(n)%const_init_value = fm_util_get_real('const_init_value', caller = caller_str, scalar = .true.)
endif !}
enddo !}
! register diagnostics for the tracer
do n=1,num_diag_tracers
range_array(1) = T_diag(n)%min_range
range_array(2) = T_diag(n)%max_range
if(T_diag(n)%longname=='Potential temperature') then
id_diag(n) = register_diag_field ('ocean_model', trim(T_diag(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_diag(n)%longname), trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_water_potential_temperature')
id_diag_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_diag(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_diag(n)%longname), trim(temp_units), &
missing_value=missing_value, range=range_array, &
standard_name='sea_surface_temperature')
id_diag_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_diag(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_diag(n)%longname), &
'squared '//trim(temp_units), &
missing_value=missing_value, range=(/0.0,1e10/), &
standard_name='square_of_sea_surface_temperature')
elseif(T_diag(n)%longname=='Conservative temperature') then
id_diag(n) = register_diag_field ('ocean_model', trim(T_diag(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_diag(n)%longname), trim(temp_units), &
missing_value=missing_value, range=range_array)
id_diag_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_diag(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_diag(n)%longname), trim(temp_units), &
missing_value=missing_value, range=range_array)
id_diag_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_diag(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_diag(n)%longname), &
'squared '//trim(temp_units), &
missing_value=missing_value, range=(/0.0,1e10/))
else
id_diag(n) = register_diag_field ('ocean_model', trim(T_diag(n)%name), &
Grd%tracer_axes(1:3), &
Time%model_time, trim(T_diag(n)%longname), trim(T_diag(n)%units), &
missing_value=missing_value)
id_diag_surf_tracer(n) = register_diag_field ('ocean_model', &
'surface_'//trim(T_diag(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, trim(T_diag(n)%longname), trim(T_diag(n)%units), &
missing_value=missing_value, range=range_array)
id_diag_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_diag(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_diag(n)%longname), &
'squared '//trim(T_diag(n)%units), &
missing_value=missing_value, range=(/0.0,1e20/))
endif
id_diag_total(n) = register_diag_field ('ocean_model','total_ocean'//trim(T_diag(n)%name), &
Time%model_time, 'Total ocean mass for tracer '//trim(T_diag(n)%name), 'kg/1e18', &
missing_value=missing_value, range=(/0.0,1e20/))
enddo
#ifdef USE_OCEAN_BGC
!Get the %filed for "generic" tracers as it might have already been set.
if (num_diag_tracers .gt. 0 ) then
do n=1,num_diag_tracers
if(T_diag(n)%type .eq. 'generic') then
!call ocean_generic_get_field(T_diag(n)%name,T_diag(n)%field)
call ocean_generic_set_pointer(T_diag(n)%name, 'field', T_diag(n)%field, isd, jsd)
!nnz: find a way to use their already allocated field pointer directly.
!#ifndef MOM_STATIC_ARRAYS
!! deallocate(T_diag(n)%field)
!!#endif
!! call ocean_generic_get_field_pointer(T_diag(n)%name,T_diag(n)%field)
endif
enddo
endif
#endif
! read diagnostic tracer restarts
write (stdoutunit,*) trim(note_header), ' Reading diagnostic tracer initial conditions and/or restarts'
do n=1,num_diag_tracers !{
write (stdoutunit,*)
write (stdoutunit,*) &
'Initializing tracer number', n,' at time level tau. This tracer is called ',trim(T_diag(n)%name)
if (T_diag(n)%restart_file .eq. ' ') then !{
write (stdoutunit,*) 'Skipping tracer ', trim(T_diag(n)%name)
else !}{
filename = T_diag(n)%restart_file
do l = 1, num_tracer_restart
if(trim(filename) == tracer_restart_file(l)) exit
end do
if (l>num_tracer_restart) then
num_tracer_restart = num_tracer_restart + 1
tracer_restart_file(l) = trim(filename)
end if
id_restart = register_restart_field(Tracer_restart(l), filename, T_diag(n)%name, &
T_diag(n)%field(:,:,:), Dom%domain2d)
filename = 'INPUT/'//trim(T_diag(n)%restart_file)
if (Time%init .and. T_diag(n)%const_init_tracer ) then
write (stdoutunit,*) &
'Initializing diagnostic tracer ', trim(T_diag(n)%name),' to constant ', T_diag(n)%const_init_value
T_diag(n)%field(:,:,:) = T_diag(n)%const_init_value
elseif (field_exist(filename, T_diag(n)%name)) then
write (stdoutunit,*) 'Reading restart for diag tracer ', trim(T_diag(n)%name), ' from file ', &
trim(T_diag(n)%restart_file)
call restore_state(Tracer_restart(l), id_restart)
write (stdoutunit,*) 'Completed read of ic/restart for diagnostic tracer ', trim(T_diag(n)%name)
if (Time%init) then !{
if(remap_depth_to_s_init .and. vert_coordinate_type == QUASI_HORIZONTAL) then
call mpp_error(WARNING, trim(note_header) // &
'No need to remap init-conditions from depth to s-coord with quasi-horizontal vert-coord.')
endif
if(.not. remap_depth_to_s_init .and. vert_coordinate_type == TERRAIN_FOLLOWING) then
call mpp_error(WARNING, trim(note_header) // &
'May need to remap init-conditions from depth to s-coord with terrain following vert-coord.')
endif
if(remap_depth_to_s_init .and. vert_coordinate_type == TERRAIN_FOLLOWING) then
call mpp_error(NOTE, trim(note_header) // &
'Remap init-tracer from depth to s-coord for terrain following vert-coord.')
endif
! remap the depth-based initial conditions onto s-coordinates
if(remap_depth_to_s_init) then
write (stdoutunit,*) &
'After reading tracer init, remapping depth to s-coord for tracer ',trim(T_diag(n)%name)
call remap_depth_to_s(Thickness, Time, T_diag(n)%field(:,:,:))
endif
! modifications due to partial bottom cell
if (vert_coordinate_type == QUASI_HORIZONTAL .and. interpolate_tdiag_to_pbott) then
! linearly interpolate tracers to partial bottom cells
write (stdoutunit,*) 'Linearly interpolate tracer ',trim(T_diag(n)%name),' to partial cell bottom.'
do j = jsc, jec !{
do i = isc, iec !{
kb = Grd%kmt(i,j)
if (kb .gt. 1) then !{
if(use_blobs) then
fact = Thickness%dzwtT(i,j,kb-1)/Grd%dzw(kb-1)
else
fact = Thickness%dzwt(i,j,kb-1)/Grd%dzw(kb-1)
endif
T_diag(n)%field(i,j,kb) = &
(1.-fact)*T_diag(n)%field(i,j,kb-1) + fact*T_diag(n)%field(i,j,kb)
endif !}
! zero out tracers in land points
do k = kb+1,nk !{
T_diag(n)%field(i,j,k) = 0.0
enddo !} k
enddo !} i
enddo !} j
endif
endif !}
else
call mpp_error(FATAL, trim(T_diag(n)%name) &
// ' is not initialized as constant, and it does not exist in the file ' &
// trim(filename) // &
'. All tracers must have initialization specified. There is no default.')
endif
call mpp_update_domains(T_diag(n)%field(:,:,:), Dom%domain2d)
endif !}
enddo !} n
write (stdoutunit,*) trim(note_header), ' Finished reading diagnostic tracer restarts.'
diag_module_initialized = .true.
if (debug_this_module) then
write(stdoutunit,'(a)') ' '
do n=1,num_diag_tracers
if (T_diag(n)%restart_file .ne. ' ') then !{
write(stdoutunit,*) 'From ocean_tracer_mod: initial ',trim(T_diag(n)%name), ' chksum (diag) ==>'
call tracer_diag_chksum(Time, T_diag(n))
endif !}
enddo
endif
do n=1,num_diag_tracers
if (T_diag(n)%longname == 'Conservative temperature') index_diag_temp = n
if (T_diag(n)%longname == 'Potential temperature') index_diag_temp = n
enddo
if (index_diag_temp == -1) then
call mpp_error(FATAL, trim(error_header) // ' diagnostic temperature not included. MOM needs it.')
endif
else !}{
write (stdoutunit,*)
write (stdoutunit,*) trim(note_header), 'No diagnostic tracers requested.'
endif !}
end function ocean_diag_tracer_init !}
! NAME="ocean_diag_tracer_init">
!#######################################################################
!
!
!
! Update value of tracer concentration to time taup1.
!
! Note that T_prog(n)%source is added at the very end
! of the time step, after the rho_dzt factor has been
! divided.
!
!
subroutine update_ocean_tracer (Time, Dens, Adv_vel, Thickness, pme, diff_cbt, &
T_prog, T_diag, L_system, &
Velocity, Ext_mode, EL_diag, use_blobs)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(inout) :: Dens
type(ocean_adv_vel_type), intent(in) :: Adv_vel
type(ocean_thickness_type), intent(inout) :: Thickness
real, dimension(isd:,jsd:), intent(in) :: pme
real, dimension(isd:,jsd:,:,:), intent(in) :: diff_cbt
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_diag_tracer_type), intent(inout) :: T_diag(:)
type(ocean_lagrangian_type), intent(inout) :: L_system
type(ocean_velocity_type), intent(in) :: Velocity
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(blob_diag_type), intent(inout) :: EL_diag(0:)
logical, intent(in) :: use_blobs
type(time_type) :: next_time
type(time_type) :: time_step
real :: age_tracer_max
integer :: i, j, k, kbot, n
integer :: taum1, tau, taup1
integer :: stdoutunit
stdoutunit=stdout()
taum1 = Time%taum1
tau = Time%tau
taup1 = Time%taup1
if (size(T_prog(:)) < num_prog_tracers) then
call mpp_error(FATAL,&
'==>Error from ocean_tracer_mod (update_ocean_tracer): T_prog array size too small')
endif
if (zero_tendency) then
! hold tracer fields constant in time
do n=1,num_prog_tracers
T_prog(n)%field(:,:,:,taup1) = T_prog(n)%field(:,:,:,tau)
enddo
else
time_step = set_time(int(dtts),0)
next_time = Time%model_time + time_step
! add contributions to T_prog%th_tendency terms arising
! from time explicit advection and diffusion
do n=1,num_prog_tracers
if(n==index_temp_sq .or. n==index_salt_sq) cycle
call mpp_clock_begin(id_clock_tracer_advect)
call horz_advect_tracer(Time, Adv_vel, Thickness, Dens, &
T_prog(1:num_prog_tracers), T_prog(n), n, dtime)
call vert_advect_tracer(Time, Adv_vel, Dens, Thickness, &
T_prog(1:num_prog_tracers), T_prog(n), n, dtime)
call mpp_clock_end(id_clock_tracer_advect)
call mpp_clock_begin(id_clock_lap_tracer)
call lap_tracer(Time, Thickness, T_prog(n), n)
call mpp_clock_end(id_clock_lap_tracer)
call mpp_clock_begin(id_clock_bih_tracer)
call bih_tracer(Time, Thickness, T_prog(n), n)
call mpp_clock_end(id_clock_bih_tracer)
call mpp_clock_begin(id_clock_vert_diffuse)
call vert_diffuse(Time, Thickness, n, T_prog(n), diff_cbt)
call mpp_clock_end(id_clock_vert_diffuse)
! --add pme contributions to the k=1 cell.
! --river contributions are in T_prog%th_tendency from ocean_rivermix_mod
! so long as the rivermix module is enabled.
! --add eta smoother contributions due to surface height smoothing.
do j=jsc,jec
do i=isc,iec
T_prog(n)%th_tendency(i,j,1) = T_prog(n)%th_tendency(i,j,1) + &
Grd%tmask(i,j,1)*(pme(i,j)*T_prog(n)%tpme(i,j) + T_prog(n)%eta_smooth(i,j) )
enddo
enddo
! --add bottom smoother contributions due to bottom pressure smoothing
do j=jsc,jec
do i=isc,iec
kbot = Grd%kmt(i,j)
if(kbot>0) then
T_prog(n)%th_tendency(i,j,kbot) = T_prog(n)%th_tendency(i,j,kbot) + &
Grd%tmask(i,j,kbot)*T_prog(n)%pbot_smooth(i,j)
endif
enddo
enddo
! update tracer concentration to taup1 using explicit tendency update
do k=1,nk
do j=jsc,jec
do i=isc,iec
T_prog(n)%field(i,j,k,taup1) = &
(Thickness%rho_dzt(i,j,k,taum1)*T_prog(n)%field(i,j,k,taum1) &
+ dtime*T_prog(n)%th_tendency(i,j,k)) &
*Thickness%rho_dztr(i,j,k)
enddo
enddo
enddo
! Add the contribution from entraining/detraining/destroyed blobs
if (use_blobs) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
T_prog(n)%field(i,j,k,taup1) = T_prog(n)%field(i,j,k,taup1) &
+ dtime*T_prog(n)%tend_blob(i,j,k)*Thickness%rho_dztr(i,j,k)
enddo
enddo
enddo
endif
! this side boundary condition does not conserve tracer...
! it simply inserts tracer to the domain with a specified value.
if(inflow_nboundary) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(mask_inflow(i,j,k)==1.0) then
T_prog(index_temp)%field(i,j,k,taup1) = temp_inflow(i,j,k)*Grd%tmask(i,j,k)
T_prog(index_salt)%field(i,j,k,taup1) = salt_inflow(i,j,k)*Grd%tmask(i,j,k)
endif
enddo
enddo
enddo
endif
if(have_obc) then
call ocean_obc_tracer(T_prog(n)%field, Adv_vel%uhrho_et, Adv_vel%vhrho_nt, &
Thickness, pme, taum1, tau, taup1, Time, n)
endif
if (debug_this_module) then
write(stdoutunit,'(a)') ' '
write(stdoutunit,*) 'From ocean_tracer_mod: tracers after explicit tendency update (taup1)'
call tracer_prog_chksum(Time, T_prog(n), taup1)
call tracer_min_max(Time, Thickness, T_prog(n))
endif
if (id_tendency_expl(n)> 0) then
call diagnose_3d(Time, Grd, id_tendency_expl(n), T_prog(n)%th_tendency(:,:,:)*T_prog(n)%conversion)
endif
if (id_eta_smooth(n)> 0) then
call diagnose_2d(Time, Grd, id_eta_smooth(n), T_prog(n)%eta_smooth(:,:)*T_prog(n)%conversion)
endif
if (id_eta_smooth_on_nrho(n)> 0) then
wrk1(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = T_prog(n)%eta_smooth(i,j)*T_prog(n)%conversion
enddo
enddo
call diagnose_3d_rho(Time, Dens, id_eta_smooth_on_nrho(n), wrk1)
end if
if (id_pbot_smooth(n)> 0) then
call diagnose_2d(Time, Grd, id_pbot_smooth(n), T_prog(n)%pbot_smooth(:,:)*T_prog(n)%conversion)
endif
call diagnose_3d(Time, Grd, id_prog_explicit(n), T_prog(n)%field(:,:,:,taup1))
enddo ! enddo for n=1,num_prog_tracers
if (do_transport_matrix) then
! accumulate time explicit pieces for transport matrix
call transport_matrix_store_explicit(Time, T_prog, isd, ied, jsd, jed, nk, isc, iec, jsc, jec, dtts, Grd%tmask)
endif
! compute ocean heating due to frazil ice formation
! This is when frazil was computed for CM2.0 and CM2.1 simulations at GFDL
if(frazil_heating_before_vphysics) then
call mpp_clock_begin(id_clock_frazil)
call update_ocean_density_salinity(T_prog, taup1, Dens)
call compute_frazil_heating(Time, Thickness, Dens, &
T_prog(1:num_prog_tracers), T_diag(1:num_diag_tracers))
call mpp_clock_end(id_clock_frazil)
endif
! compute time implicit vertical diffusion of tracer concentration
call mpp_clock_begin(id_clock_vert_diffuse_implicit)
call vert_diffuse_implicit(diff_cbt, index_salt, Time, Thickness, Dens, T_prog(1:num_prog_tracers))
call mpp_clock_end(id_clock_vert_diffuse_implicit)
! call the time implicit part of the Lagrangian blobs
call mpp_clock_begin(id_clock_blob)
! Call the routines that form/move/destroy blobs implicitly in time
call mpp_clock_begin(id_clock_blob_implicit)
call ocean_blob_implicit(Time, Thickness, T_prog(1:num_prog_tracers), Dens, Adv_vel, Velocity, EL_diag(:))
call mpp_clock_end(id_clock_blob_implicit)
! We need to update the Lagrangian system and the partition of Thickness.
call mpp_clock_begin(id_clock_adjust_L_thickness)
call adjust_L_thickness(Time, Thickness, T_prog, L_system)
call mpp_clock_end(id_clock_adjust_L_thickness)
! We need to adjust the E system thickness to reflect the new L system thickness
call mpp_clock_begin(id_clock_adjust_E_thickness)
call update_E_thickness(Time, Grd, Thickness, Dens, Ext_mode)
call mpp_clock_end(id_clock_adjust_E_thickness)
call mpp_clock_end(id_clock_blob)
! convectively adjust water columns
call mpp_clock_begin(id_clock_convection)
call convection (Time, Thickness, T_prog(1:num_prog_tracers), Dens)
call mpp_clock_end(id_clock_convection)
! compute ocean heating due to frazil ice formation
! Computing frazil here ensures that the ocean temperature
! at the end of a time step will never be below freezing.
if(frazil_heating_after_vphysics) then
call mpp_clock_begin(id_clock_frazil)
call update_ocean_density_salinity(T_prog, taup1, Dens)
call compute_frazil_heating(Time, Thickness, Dens, &
T_prog(1:num_prog_tracers), T_diag(1:num_diag_tracers))
call mpp_clock_end(id_clock_frazil)
endif
if (debug_this_module) then
write(stdoutunit,*) 'From ocean_tracer_mod: tracers at end of update_ocean_tracer (taup1)'
call tracer_prog_chksum(Time, T_prog(index_temp), taup1)
call tracer_min_max(Time, Thickness, T_prog(index_temp))
call tracer_prog_chksum(Time, T_prog(index_salt), taup1)
call tracer_min_max(Time, Thickness, T_prog(index_salt))
call tracer_prog_chksum(Time, T_prog(index_added_heat), taup1)
call tracer_min_max(Time, Thickness, T_prog(index_added_heat))
call tracer_prog_chksum(Time, T_prog(index_redist_heat), taup1)
call tracer_min_max(Time, Thickness, T_prog(index_redist_heat))
endif
if(have_obc) then
do n=1,num_prog_tracers
call mpp_update_domains(T_prog(n)%field(:,:,:,taup1), Dom%domain2d, complete=T_prog(n)%complete)
enddo
do n=1,num_prog_tracers
call ocean_obc_update_boundary(T_prog(n)%field(:,:,:,taup1), 'T')
enddo
endif
! for numerically testing tracer advection schemes
do n=1,num_prog_tracers
if(T_prog(n)%use_only_advection) then
call update_advection_only(Time, Adv_vel, Dens, Thickness, T_prog, n)
endif
enddo
! add tracer sources, which have dimensions tracer concentration per time
if(.not. zero_tracer_source) then
do n=1,num_prog_tracers
do k=1,nk
do j=jsd,jed
do i=isd,ied
T_prog(n)%field(i,j,k,taup1) = T_prog(n)%field(i,j,k,taup1) &
+ dtts*T_prog(n)%source(i,j,k)
enddo
enddo
enddo
enddo
endif
! For age tracer, limit the lower bound at 0
! and upper bound at age_tracer_max.
age_tracer_max = age_tracer_max_init + Time%itt0*dtts*sec_in_yr_r
if(limit_age_tracer) then
do n=1,num_prog_tracers
if(T_prog(n)%name(1:3) =='age') then
do k=1,nk
do j=jsd,jed
do i=isd,ied
T_prog(n)%field(i,j,k,taup1) = max(0.0,T_prog(n)%field(i,j,k,taup1))
T_prog(n)%field(i,j,k,taup1) = min(age_tracer_max,T_prog(n)%field(i,j,k,taup1))
enddo
enddo
enddo
endif
enddo
endif
! update squared tracer
call update_tracer_passive(num_prog_tracers, T_prog, T_diag, dtts, taup1)
endif ! endif for zero_tendency
! update for Dens%rho_salinity, now that all of the
! updates have been applied to the prognostic salinity.
call update_ocean_density_salinity(T_prog, taup1, Dens)
! perform range check for temperature and salinity
if(use_tempsalt_check_range .or. debug_this_module) then
call tempsalt_check_range(T_prog(index_salt)%field(:,:,:,taup1),&
T_prog(index_temp)%field(:,:,:,taup1),&
Thickness%depth_zt(:,:,:))
endif
! send some tracer diagnostics at time tau
call send_tracer_diagnostics(Time, T_prog, T_diag, Thickness, Dens, use_blobs)
! compute watermass diagnostics
call watermass_diag(Time, T_prog, T_diag, Dens, Thickness, pme)
! update the diagnostic temperature variable to taup1. This must be done after
! diagnostics are sent at time tau.
if(prog_temp_variable==CONSERVATIVE_TEMP) then
T_diag(index_diag_temp)%field(:,:,:) = pottemp_from_contemp(Dens%rho_salinity(:,:,:,taup1), &
T_prog(index_temp)%field(:,:,:,taup1)) &
*Grd%tmask(:,:,:)
else
T_diag(index_diag_temp)%field(:,:,:) = contemp_from_pottemp(Dens%rho_salinity(:,:,:,taup1), &
T_prog(index_temp)%field(:,:,:,taup1)) &
*Grd%tmask(:,:,:)
endif
! Update the combined tracer concentration to be at time taup1.
! Do so here so that diagnostics for fieldT are at tau.
if (use_blobs) then
do n=1,num_prog_tracers
do k=1,nk
do j=jsd,jed
do i=isd,ied
T_prog(n)%fieldT(i,j,k) = T_prog(n)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
+ T_prog(n)%sum_blob(i,j,k,taup1)*Grd%datr(i,j)
T_prog(n)%fieldT(i,j,k) = Grd%tmask(i,j,k)*T_prog(n)%fieldT(i,j,k)/(epsln+Thickness%rho_dztT(i,j,k,taup1))
enddo
enddo
enddo
call mpp_update_domains(T_prog(n)%fieldT(:,:,:), Dom%domain2d, complete=T_prog(n)%complete)
enddo
if (debug_this_module) then
write(stdoutunit, '(/,a)') 'Checksums for tracer fields (E and L systems)'
do n=1,num_prog_tracers
call write_chksum_3d('Combined '//trim(T_prog(n)%name), Grd%tmask(COMP,:)*T_prog(n)%fieldT(COMP,:))
call write_chksum_3d('Sum of L system '//trim(T_prog(n)%name), Grd%tmask(COMP,:)*T_prog(n)%sum_blob(COMP,:,taup1))
enddo
endif
endif
end subroutine update_ocean_tracer
! NAME="update_ocean_tracer">
!#######################################################################
!
!
!
!
! Redo tracer updates for those that use only advection--nothing else.
! This method is useful for testing advection schemes.
!
! T_prog(n)%use_only_advection==.true. ignores all boundary forcing
! and sources, so if T_prog(n)%stf or pme, rivers, sources
! are nonzero, tracer diagnostics will spuriously indicate
! non-conservation.
!
! Assume for these tests that
! (1) vertical advection is done fully explictly in time
! (2) pme, rivers, stf, btf, and other sources are zero
! (3) do not use advect_sweby_all
!
!
!
subroutine update_advection_only(Time, Adv_vel, Dens, Thickness, T_prog, n)
type(ocean_time_type), intent(in) :: Time
type(ocean_adv_vel_type), intent(in) :: Adv_vel
type(ocean_density_type), intent(in) :: Dens
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
integer, intent(in) :: n
integer :: i, j, k
integer :: taup1, taum1
taup1 = Time%taup1
taum1 = Time%taum1
T_prog(n)%th_tendency = 0.0
call horz_advect_tracer(Time, Adv_vel, Thickness, Dens, &
T_prog(1:num_prog_tracers), T_prog(n), n, dtime, store_flux=.FALSE.)
call vert_advect_tracer(Time, Adv_vel, Dens, Thickness, &
T_prog(1:num_prog_tracers), T_prog(n), n, dtime)
do k=1,nk
do j=jsc,jec
do i=isc,iec
T_prog(n)%field(i,j,k,taup1) = &
(Thickness%rho_dzt(i,j,k,taum1)*T_prog(n)%field(i,j,k,taum1) &
+ dtime*T_prog(n)%th_tendency(i,j,k)) * Thickness%rho_dztr(i,j,k)
enddo
enddo
enddo
end subroutine update_advection_only
! NAME="update_advection_only">
!#######################################################################
!
!
! Write out restart files registered through register_restart_file
!
subroutine ocean_tracer_restart(Time, T_prog, time_stamp)
type(ocean_time_type), intent(in) :: Time
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
character(len=*), intent(in), optional :: time_stamp
integer :: tau, taup1, n
tau = Time%tau
taup1 = Time%taup1
do n=1,num_prog_tracers
if(tendency==THREE_LEVEL) then
call reset_field_pointer(Tracer_restart(id_tracer_type(n)), id_tracer_restart(n), T_prog(n)%field(:,:,:,tau), &
T_prog(n)%field(:,:,:,taup1) )
else
call reset_field_pointer(Tracer_restart(id_tracer_type(n)), id_tracer_restart(n), T_prog(n)%field(:,:,:,taup1) )
endif
enddo
do n=1, num_tracer_restart
call save_restart(Tracer_restart(n), time_stamp)
enddo
end subroutine ocean_tracer_restart
! NAME="ocean_tracer_restart"
!#######################################################################
!
!
!
! Write ocean tracer restarts
!
!
subroutine ocean_tracer_end(Time, T_prog, T_diag, use_blobs)
type(ocean_time_type), intent(in) :: Time
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
type(ocean_diag_tracer_type), intent(in) :: T_diag(:)
logical, intent(in) :: use_blobs
integer :: n, tau, taup1, unit
integer(i8_kind) :: chksum
character(len=128) :: filename
real :: fieldsum
integer :: stdoutunit
stdoutunit=stdout()
tau = Time%tau
taup1 = Time%taup1
if(.not. write_a_restart) then
write(stdoutunit,'(/a)') '==>Warning from ocean_tracer_mod (ocean_tracer_end): NO restart written.'
call mpp_error(WARNING,'==>Warning from ocean_tracer_mod (ocean_tracer_end): NO restart written.')
return
endif
! open ascii file containing checksums used to check for reproducibility
filename = 'RESTART/ocean_tracer.res'
call mpp_open(unit, trim(filename),form=MPP_ASCII,&
action=MPP_OVERWR,threading=MPP_SINGLE,fileset=MPP_SINGLE,nohdrs=.true.)
call ocean_tracer_restart(Time, T_prog)
! write the restart fields to the file T_prog(n)%restart_file
! write the chksums to the ascii file RESTART/ocean_tracer.res
do n=1,num_prog_tracers
if(tendency==THREE_LEVEL) then
write(stdoutunit,*) ' '
write(stdoutunit,*) 'Ending ',trim(T_prog(n)%name), ' chksum (tau) ==>'
call tracer_prog_chksum(Time, T_prog(n), tau, chksum)
if (mpp_pe() == mpp_root_pe()) write(unit,*) trim(T_prog(n)%name)//' tau chksum =', chksum
endif
write(stdoutunit,*) ' '
write(stdoutunit,*) 'Ending ',trim(T_prog(n)%name), ' chksum (taup1) ==>'
call tracer_prog_chksum(Time, T_prog(n), taup1, chksum)
if (mpp_pe() == mpp_root_pe()) write(unit,*) trim(T_prog(n)%name)//' taup1 chksum =', chksum
if (use_blobs) then
fieldsum = mpp_chksum(Grd%tmask(isc:iec,jsc:jec,:)*T_prog(n)%fieldT(isc:iec,jsc:jec,:))
if (mpp_pe() == mpp_root_pe()) write(stdoutunit,*) trim(T_prog(n)%name)//' Total field =', fieldsum
endif
enddo
! write the restart fields to the file T_diag(n)%restart_file
! write the chksums to the ascii file RESTART/ocean_tracer.res where writing restart
do n=1,num_diag_tracers
if (T_diag(n)%restart_file .ne. ' ') then
write(stdoutunit,*) ' '
write(stdoutunit,*) 'Ending ',trim(T_diag(n)%name), ' chksum (diag) ==>'
call tracer_diag_chksum(Time, T_diag(n), chksum)
if (mpp_pe() == mpp_root_pe()) then
write(unit,*) trim(T_diag(n)%name)//' chksum =', chksum
endif
endif
enddo
return
end subroutine ocean_tracer_end
! NAME="ocean_tracer_end">
!#######################################################################
!
!
!
! Provide for possibility that quicker advection reverts to
! first order upwind when tracer is outside a specified range.
! Likewise, may wish to revert neutral physics to horizontal diffusion.
!
! For this purpose, we define a mask which is set to unity where
! fluxes revert to first order upwind advection
! (if using quicker) and horizontal diffusion (if using neutral).
!
! This method is very ad hoc. What is preferred for advection is to use
! a monotonic scheme, such as mdfl_sweby or mdppm. For neutral physics,
! no analogous monotonic scheme has been implemented. Such could be
! useful, especially for passive tracers. In the meantime, tmask_limit
! provides a very rough limiter for neutral physics to help keep tracers
! within specified bounds.
!
!
!
subroutine compute_tmask_limit(Time, T_prog)
type(ocean_time_type), intent(in) :: Time
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
integer :: tau
integer :: i, j, k, n
integer :: stdoutunit
stdoutunit=stdout()
if(.not. compute_tmask_limit_on) return
tau = Time%tau
do n=1,num_prog_tracers
do k=1,nk
do j=jsc-1,jec
do i=isc-1,iec
T_prog(n)%tmask_limit(i,j,k) = 0.0
! is tracer outside range?
if(Grd%tmask(i,j,k) == 1.0) then
if(T_prog(n)%field(i,j,k,tau) < T_prog(n)%min_tracer_limit .or. &
T_prog(n)%field(i,j,k,tau) > T_prog(n)%max_tracer_limit) then
T_prog(n)%tmask_limit(i,j,k) = Grd%tmask(i,j,k)
endif
endif
! is east tracer outside range?
if(Grd%tmask(i+1,j,k) == 1.0) then
if(T_prog(n)%field(i+1,j,k,tau) < T_prog(n)%min_tracer_limit .or. &
T_prog(n)%field(i+1,j,k,tau) > T_prog(n)%max_tracer_limit) then
T_prog(n)%tmask_limit(i,j,k) = Grd%tmask(i,j,k)
endif
endif
! is north tracer outside range?
if(Grd%tmask(i,j+1,k) == 1.0) then
if(T_prog(n)%field(i,j+1,k,tau) < T_prog(n)%min_tracer_limit .or. &
T_prog(n)%field(i,j+1,k,tau) > T_prog(n)%max_tracer_limit) then
T_prog(n)%tmask_limit(i,j,k) = Grd%tmask(i,j,k)
endif
endif
enddo
enddo
enddo
! is deeper tracer outside range?
do k=1,nk-1
do j=jsc-1,jec
do i=isc-1,iec
if(Grd%tmask(i,j,k+1) == 1.0) then
if(T_prog(n)%field(i,j,k+1,tau) < T_prog(n)%min_tracer_limit .or. &
T_prog(n)%field(i,j,k+1,tau) > T_prog(n)%max_tracer_limit) then
T_prog(n)%tmask_limit(i,j,k) = Grd%tmask(i,j,k)
endif
endif
enddo
enddo
enddo
enddo ! end of loop n=1,num_prog_tracers
! insist that if temperature or salinity tmask_limit is 1.0, then both must be 1.0
if(tmask_limit_ts_same) then
do k=1,nk
do j=jsc-1,jec
do i=isc-1,iec
if( T_prog(index_temp)%tmask_limit(i,j,k)==1.0 &
.or. T_prog(index_salt)%tmask_limit(i,j,k)==1.0) then
T_prog(index_temp)%tmask_limit(i,j,k)=1.0
T_prog(index_salt)%tmask_limit(i,j,k)=1.0
endif
enddo
enddo
enddo
endif
! ensure that redist_heat and added_heat masks are the same as temp
if ( index_redist_heat /= -1 ) T_prog(index_redist_heat)%tmask_limit = T_prog(index_temp)%tmask_limit
if ( index_added_heat /= -1 ) T_prog(index_added_heat)%tmask_limit = T_prog(index_temp)%tmask_limit
! debugging and diagnostics
do n=1,num_prog_tracers
if(debug_this_module) then
write(stdoutunit,*) ' '
call write_timestamp(Time%model_time)
call write_chksum_3d('tmask_limit', T_prog(n)%tmask_limit(COMP,:))
endif
call diagnose_3d(Time, Grd, id_tmask_limit(n), T_prog(n)%tmask_limit(:,:,:))
enddo
end subroutine compute_tmask_limit
! NAME="compute_tmask_limit"
!#######################################################################
!
!
!
! Remap in the vertical from s-coordinate to depth and then send to
! diagnostic manager.
!
! This routine is mostly of use for terrain following vertical
! coordinates, which generally deviate a lot from depth or pressure
! coordinates. The zstar and pstar coordinates are very similar
! to z or pressure, so there is no need to do the remapping for
! purposes of visualization.
!
! The routine needs to be made more general and faster.
! It also has been found to be problematic, so it is NOT
! recommended. It remains here as a template for a better algorithm.
! Remapping methods in Ferret are much better.
!
! Use rho_dzt weighting to account for nonBoussinesq.
!
! Author: Stephen.Griffies
!
!
!
subroutine remap_s_to_depth(Thickness, Time, array_in, ntracer)
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:,jsd:,:), intent(in) :: array_in
integer, intent(in) :: ntracer
integer :: i, j, k, kk, tau
tau = Time%tau
wrk1 = 0.0
wrk2 = 0.0
wrk3 = 0.0
k=1
do kk=1,nk
do j=jsd,jed
do i=isd,ied
if(Thickness%depth_zt(i,j,kk) < Grd%zw(k)) then
wrk1(i,j,k) = wrk1(i,j,k) + array_in(i,j,kk)*Thickness%rho_dzt(i,j,kk,tau)
wrk2(i,j,k) = wrk2(i,j,k) + Thickness%rho_dzt(i,j,kk,tau)
endif
enddo
enddo
enddo
do k=1,nk-1
do kk=1,nk
do j=jsd,jed
do i=isd,ied
if(Grd%zw(k) <= Thickness%depth_zt(i,j,kk) .and. Thickness%depth_zt(i,j,kk) < Grd%zw(k+1)) then
wrk1(i,j,k+1) = wrk1(i,j,k+1) + array_in(i,j,kk)*Thickness%rho_dzt(i,j,kk,tau)
wrk2(i,j,k+1) = wrk2(i,j,k+1) + Thickness%rho_dzt(i,j,kk,tau)
endif
enddo
enddo
enddo
enddo
k=nk
do kk=1,nk
do j=jsd,jed
do i=isd,ied
if(Thickness%depth_zt(i,j,kk) > Grd%zw(k) ) then
wrk1(i,j,k) = wrk1(i,j,k) + array_in(i,j,kk)*Thickness%rho_dzt(i,j,kk,tau)
wrk2(i,j,k) = wrk2(i,j,k) + Thickness%rho_dzt(i,j,kk,tau)
endif
enddo
enddo
enddo
do k=1,nk
do j=jsd,jed
do i=isd,ied
wrk3(i,j,k) = Grd%tmask_depth(i,j,k)*wrk1(i,j,k)/(wrk2(i,j,k)+epsln)
enddo
enddo
enddo
used = send_data (id_prog_on_depth(ntracer), wrk3(:,:,:), &
Time%model_time,rmask=Grd%tmask_depth(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
end subroutine remap_s_to_depth
! NAME="remap_s_to_depth"
!#######################################################################
!
!
!
! Remap in the vertical from depth to s-coordinate. This routine is
! used for initializing terrain following coordinate models given an
! initial tracer field generated assuming depth-like vertical coordinate.
!
! This routine is of use for terrain following vertical coordinates,
! which generally deviate a lot from z, zstar, pressure, or pstar
! coordinates.
!
! Algorithm is very rudimentary and can be made better.
!
! Author: Stephen.Griffies
!
!
!
subroutine remap_depth_to_s(Thickness, Time, array)
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_time_type), intent(in) :: Time
real, dimension(isd:,jsd:,:), intent(inout) :: array
integer :: i, j, k, kk
integer :: stdoutunit
stdoutunit=stdout()
wrk1 = 0.0
k=1
do kk=1,nk
do j=jsd,jed
do i=isd,ied
if(Thickness%depth_zt(i,j,kk) < Grd%zw(k)) then
wrk1(i,j,kk) = array(i,j,k)
endif
enddo
enddo
enddo
do k=1,nk-1
do kk=1,nk
do j=jsd,jed
do i=isd,ied
if(Grd%zw(k) <= Thickness%depth_zt(i,j,kk) .and. Thickness%depth_zt(i,j,kk) < Grd%zw(k+1)) then
wrk1(i,j,kk) = array(i,j,k+1)
endif
enddo
enddo
enddo
enddo
k=nk
do kk=1,nk
do j=jsd,jed
do i=isd,ied
if(Thickness%depth_zt(i,j,kk) > Grd%zw(k) ) then
wrk1(i,j,kk) = array(i,j,k)
endif
enddo
enddo
enddo
if(debug_this_module) then
write(stdoutunit,'(/a)')'==>from remap_depth_to_s'
call write_timestamp(Time%model_time)
do k=1,nk
do j=jsc,jsc
do i=isc,isc
write(stdoutunit,'(a,i3,a,i3,a,i3,a,f12.4)') &
'array_in (',i+Dom%ioff,',',j+Dom%joff,',',k,') = ',array(i,j,k)
write(stdoutunit,'(a,i3,a,i3,a,i3,a,f12.4)') &
'array_out(',i+Dom%ioff,',',j+Dom%joff,',',k,') = ',wrk1(i,j,k)
enddo
enddo
enddo
endif
array(:,:,:) = wrk1(:,:,:)
end subroutine remap_depth_to_s
! NAME="remap_depth_to_s"
!#######################################################################
!
!
!
! Initialize the mask, temp, and salt values at the northern boundary
! where we are specifying values for an inflow boundary condition.
!
!
!
subroutine inflow_nboundary_init
character(len=128) :: filename
integer :: stdoutunit
stdoutunit=stdout()
allocate(mask_inflow(isd:ied,jsd:jed,nk))
allocate(temp_inflow(isd:ied,jsd:jed,nk))
allocate(salt_inflow(isd:ied,jsd:jed,nk))
write(stdoutunit,'(a)') &
'==>Note: running with inflow_nboundary. Specify tracer at north boundary. This is nonconservative.'
mask_inflow(:,:,:) = 0.0
temp_inflow(:,:,:) = 0.0
salt_inflow(:,:,:) = 0.0
filename = 'INPUT/mask_inflow.nc'
if (file_exist(trim(filename))) then
call read_data(filename,'mask_inflow',mask_inflow,Dom%domain2d,timelevel=1)
call mpp_update_domains(mask_inflow(:,:,:), Dom%domain2d)
else
call mpp_error(FATAL,&
'==>Error from ocean_tracer_mod: inflow_nboundary_init cannot find mask_inflow file.')
endif
filename = 'INPUT/temp_inflow.nc'
if (file_exist(trim(filename))) then
call read_data(filename,'temp_inflow',temp_inflow,Dom%domain2d,timelevel=1)
call mpp_update_domains(temp_inflow(:,:,:), Dom%domain2d)
else
call mpp_error(FATAL,&
'==>Error from ocean_tracer_mod: inflow_nboundary_init cannot find temp_inflow file.')
endif
filename = 'INPUT/salt_inflow.nc'
if (file_exist(trim(filename))) then
call read_data(filename,'salt_inflow',salt_inflow,Dom%domain2d,timelevel=1)
call mpp_update_domains(salt_inflow(:,:,:), Dom%domain2d)
else
call mpp_error(FATAL,&
'==>Error from ocean_tracer_mod: inflow_nboundary_init cannot find salt_inflow file.')
endif
end subroutine inflow_nboundary_init
! NAME="inflow_nboundary_init"
!#######################################################################
!
!
!
! Do some initialization required for diagnostics.
! Density must already be initialized for this routine to be used.
! Hence, this routine is called from ocean_model.F90.
!
!
subroutine ocean_tracer_diagnostics_init(Time, Dens, T_diag, T_prog, ver_coordinate_class)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
type(ocean_diag_tracer_type), intent(in) :: T_diag(:)
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
integer, intent(in) :: ver_coordinate_class
integer :: n
integer :: stdoutunit
stdoutunit=stdout()
compute_watermass_diag = .false.
vert_coordinate_class = ver_coordinate_class
neutralrho_nk = size(Dens%neutralrho_ref(:))
allocate( nrho_work(isd:ied,jsd:jed,neutralrho_nk) )
allocate( nrho_work2(isd:ied,jsd:jed,neutralrho_nk) )
allocate( nrho_mass(isd:ied,jsd:jed,neutralrho_nk) )
nrho_work(:,:,:) = 0.0
nrho_work2(:,:,:) = 0.0
nrho_mass(:,:,:) = 0.0
do n=1,num_diag_tracers
if (T_diag(n)%name == 'frazil') index_frazil = n
enddo
do n=1,num_prog_tracers
! Neutral density binning of heat budget diagnostics
id_tendency_on_nrho(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_tendency_on_nrho', Dens%neutralrho_axes(1:3), &
Time%model_time, 'time tendency for tracer '//trim(T_prog(n)%longname)//" binned to neutral density", &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e20,1.e20/))
id_eta_smooth_on_nrho(n) = register_diag_field ('ocean_model', &
trim(T_prog(n)%name)//'_eta_smooth_on_nrho', Dens%neutralrho_axes(1:3), &
Time%model_time, 'surface smoother for ' // trim(T_prog(n)%longname)//" binned to neutral density", &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e20,1.e20/))
enddo
! time tendency for locally referenced potential density
id_neut_rho_tendency = register_diag_field ('ocean_model','neut_rho_tendency', &
Grd%tracer_axes(1:3), Time%model_time,'total time tendency of locally ref potrho',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_tendency > 0) compute_watermass_diag=.true.
! time tendency for potential density
id_pot_rho_tendency = register_diag_field ('ocean_model','pot_rho_tendency', &
Grd%tracer_axes(1:3), Time%model_time,'total time tendency of depth ref potrho',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_pot_rho_tendency > 0) compute_watermass_diag=.true.
id_neut_rho_tendency_on_nrho = register_diag_field ('ocean_model', &
'neut_rho_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'total time tendency of locally ref potrho as binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_tendency = register_diag_field ('ocean_model','wdian_rho_tendency',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to time tendency of locally ref potrho', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_tendency > 0) compute_watermass_diag=.true.
id_wdian_rho_tendency_on_nrho = register_diag_field ('ocean_model', &
'wdian_rho_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to total time tendency as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_tendency = register_diag_field ('ocean_model', &
'tform_rho_tendency', Grd%tracer_axes(1:3), Time%model_time, &
'water mass transformation from total tendency on levels (pre-layer binning)',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_tendency > 0) compute_watermass_diag=.true.
id_tform_rho_tendency_on_nrho = register_diag_field ('ocean_model', &
'tform_rho_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'water mass transformation from total tendency as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_tendency_on_nrho > 0) compute_watermass_diag=.true.
! time tendency from temperature effects
id_neut_temp_tendency = register_diag_field ('ocean_model','neut_temp_tendency', &
Grd%tracer_axes(1:3), Time%model_time,'temp related time tendency of locally ref potrho',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_tendency > 0) compute_watermass_diag=.true.
id_neut_temp_tendency_on_nrho = register_diag_field ('ocean_model', &
'neut_temp_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related time tendency of locally ref potrho binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_temp_tendency = register_diag_field ('ocean_model','wdian_temp_tendency', &
Grd%tracer_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport due to time tendency of locally ref potrho',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_temp_tendency > 0) compute_watermass_diag=.true.
id_wdian_temp_tendency_on_nrho = register_diag_field ('ocean_model', &
'wdian_temp_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport due to time tendency as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_temp_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_tform_temp_tendency = register_diag_field ('ocean_model', &
'tform_temp_tendency', Grd%tracer_axes(1:3), Time%model_time, &
'temp related water mass transformation from tendency on levels',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_temp_tendency > 0) compute_watermass_diag=.true.
id_tform_temp_tendency_on_nrho = register_diag_field ('ocean_model', &
'tform_temp_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related water mass transformation from tendency binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_temp_tendency_on_nrho > 0) compute_watermass_diag=.true.
! time tendency from salinity effects
id_neut_salt_tendency = register_diag_field ('ocean_model','neut_salt_tendency', &
Grd%tracer_axes(1:3), Time%model_time,'salt related time tendency of locally ref potrho',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_tendency > 0) compute_watermass_diag=.true.
id_neut_salt_tendency_on_nrho = register_diag_field ('ocean_model', &
'neut_salt_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related time tendency of locally ref potrho binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_salt_tendency = register_diag_field ('ocean_model','wdian_salt_tendency', &
Grd%tracer_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport due to time tendency of locally ref potrho',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_salt_tendency > 0) compute_watermass_diag=.true.
id_wdian_salt_tendency_on_nrho = register_diag_field ('ocean_model', &
'wdian_salt_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport due to time tendency as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_salt_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_tform_salt_tendency = register_diag_field ('ocean_model', &
'tform_salt_tendency', Grd%tracer_axes(1:3), Time%model_time, &
'salt related water mass transformation from tendency on levels',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_salt_tendency > 0) compute_watermass_diag=.true.
id_tform_salt_tendency_on_nrho = register_diag_field ('ocean_model', &
'tform_salt_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related water mass transformation from tendency binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_salt_tendency_on_nrho > 0) compute_watermass_diag=.true.
! smoothing operator for locally referenced potential density
id_neut_rho_smooth = register_diag_field ('ocean_model','neut_rho_smooth',&
Grd%tracer_axes(1:3), Time%model_time, &
'time tendency of local ref potrho from eta/pbot smoother', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_smooth > 0) compute_watermass_diag=.true.
! smoothing operator for potential density
id_pot_rho_smooth = register_diag_field ('ocean_model','pot_rho_smooth',&
Grd%tracer_axes(1:3), Time%model_time, &
'time tendency of depth ref potrho from eta/pbot smoother', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_pot_rho_smooth > 0) compute_watermass_diag=.true.
id_neut_rho_smooth_on_nrho = register_diag_field ('ocean_model', &
'neut_rho_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'time tendency of local ref potrho from eta/pbot smoother binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_smooth = register_diag_field ('ocean_model','wdian_rho_smooth',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to eta/pbot smoother', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_smooth > 0) compute_watermass_diag=.true.
id_wdian_rho_smooth_on_nrho = register_diag_field ('ocean_model', &
'wdian_rho_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to eta/pbot smoother binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_smooth = register_diag_field ('ocean_model', &
'tform_rho_smooth', Dens%neutralrho_axes(1:3), Time%model_time,&
'water mass transformation from smoother on levels', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_smooth_on_nrho = register_diag_field ('ocean_model', &
'tform_rho_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time,&
'water mass transformation from smoother binned to neutral density', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_eta_tend_smooth = register_diag_field ('ocean_model','eta_tend_smooth',&
Grd%tracer_axes(1:2), Time%model_time, &
'non-Bouss steric sea level tendency from eta/pbot smoother', &
'm/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_eta_tend_smooth > 0) compute_watermass_diag=.true.
id_eta_tend_smooth_glob = register_diag_field ('ocean_model', &
'eta_tend_smooth_glob', Time%model_time, &
'global mean non-Bouss steric sea level tendency from eta/pbot smoother',&
'm/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_eta_tend_smooth_glob > 0) compute_watermass_diag=.true.
! smoothing operator from temperature effects
id_neut_temp_smooth = register_diag_field ('ocean_model','neut_temp_smooth',&
Grd%tracer_axes(1:3), Time%model_time, &
'temp related time tendency of local ref potrho from eta/pbot smoother', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_smooth > 0) compute_watermass_diag=.true.
id_neut_temp_smooth_on_nrho = register_diag_field ('ocean_model', &
'neut_temp_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related time tendency of local ref potrho from eta/pbot smoother binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_temp_smooth = register_diag_field ('ocean_model','wdian_temp_smooth',&
Grd%tracer_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport due to eta/pbot smoother', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_temp_smooth > 0) compute_watermass_diag=.true.
id_wdian_temp_smooth_on_nrho = register_diag_field ('ocean_model', &
'wdian_temp_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport due to eta/pbot smoother binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_temp_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_tform_temp_smooth = register_diag_field ('ocean_model', &
'tform_temp_smooth', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related water mass transformation from smoother on levels',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_temp_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_tform_temp_smooth_on_nrho = register_diag_field ('ocean_model', &
'tform_temp_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related water mass transformation from smoother binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_temp_smooth_on_nrho > 0) compute_watermass_diag=.true.
! smoothing operator from salinity effects
id_neut_salt_smooth = register_diag_field ('ocean_model','neut_salt_smooth',&
Grd%tracer_axes(1:3), Time%model_time, &
'salt related time tendency of local ref potrho from eta/pbot smoother', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_smooth > 0) compute_watermass_diag=.true.
id_neut_salt_smooth_on_nrho = register_diag_field ('ocean_model', &
'neut_salt_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related time tendency of local ref potrho from eta/pbot smoother binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_salt_smooth = register_diag_field ('ocean_model','wdian_salt_smooth',&
Grd%tracer_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport due to eta/pbot smoother', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_salt_smooth > 0) compute_watermass_diag=.true.
id_wdian_salt_smooth_on_nrho = register_diag_field ('ocean_model', &
'wdian_salt_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport due to eta/pbot smoother binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_salt_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_tform_salt_smooth = register_diag_field ('ocean_model', &
'tform_salt_smooth', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related water mass transformation from smoother on levels',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_salt_smooth_on_nrho > 0) compute_watermass_diag=.true.
id_tform_salt_smooth_on_nrho = register_diag_field ('ocean_model', &
'tform_salt_smooth_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related water mass transformation from smoother binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_salt_smooth_on_nrho > 0) compute_watermass_diag=.true.
! contribution of frazil to heat budget binned on neutral density
id_frazil_on_nrho = register_diag_field ('ocean_model', 'frazil_on_nrho', Dens%neutralrho_axes(1:3), &
Time%model_time, 'ocn frazil heat flux over time step binned to neutral density', 'W/m^2',&
missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_frazil_on_nrho > 0) compute_watermass_diag=.true.
! contributions from frazil formation to local referenced potential density
id_neut_rho_frazil = register_diag_field ('ocean_model', 'neut_rho_frazil',&
Grd%tracer_axes(1:3), Time%model_time, &
'update of local ref potrho from frazil formation', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_frazil > 0) compute_watermass_diag=.true.
! contributions from frazil formation to potential density
id_pot_rho_frazil = register_diag_field ('ocean_model', 'pot_rho_frazil',&
Grd%tracer_axes(1:3), Time%model_time, &
'update of depth ref potrho from frazil formation', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_pot_rho_frazil > 0) compute_watermass_diag=.true.
id_neut_rho_frazil_on_nrho = register_diag_field ('ocean_model', &
'neut_rho_frazil_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'update of local ref potrho from frazil formation as binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_frazil_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_frazil = register_diag_field ('ocean_model', 'wdian_rho_frazil',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to frazil formation', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_frazil > 0) compute_watermass_diag=.true.
id_wdian_rho_frazil_on_nrho = register_diag_field ('ocean_model', &
'wdian_rho_frazil_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to frazil formation as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_frazil_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_frazil_on_nrho = register_diag_field ('ocean_model', &
'tform_rho_frazil_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time,&
'water mass transformation from frazil as binned to neutral density', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_frazil_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_frazil = register_diag_field ('ocean_model', &
'tform_rho_frazil', Dens%neutralrho_axes(1:3), Time%model_time, &
'water mass transformation from frazil on levels (pre-layer binning)',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_frazil > 0) compute_watermass_diag=.true.
id_eta_tend_frazil = register_diag_field ('ocean_model','eta_tend_frazil',&
Grd%tracer_axes(1:2), Time%model_time, &
'non-Bouss steric sea level tendency from frazil', &
'm/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_eta_tend_frazil > 0) compute_watermass_diag=.true.
id_eta_tend_frazil_glob = register_diag_field ('ocean_model', &
'eta_tend_frazil_glob', Time%model_time, &
'global mean non-Bouss steric sea level tendency from frazil',&
'm/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_eta_tend_frazil_glob > 0) compute_watermass_diag=.true.
! contributions from pme in flux-form
id_neut_rho_pme = register_diag_field ('ocean_model','neut_rho_pme', &
Grd%tracer_axes(1:3), Time%model_time,'time tendency of local ref potrho from pme',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_pme > 0) compute_watermass_diag=.true.
id_pot_rho_pme = register_diag_field ('ocean_model','pot_rho_pme', &
Grd%tracer_axes(1:3), Time%model_time,'time tendency of depth ref potrho from pme',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_pot_rho_pme > 0) compute_watermass_diag=.true.
id_neut_rho_pme_on_nrho = register_diag_field ('ocean_model', &
'neut_rho_pme_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'time tendency of local ref potrho from pme as binned to neutral density',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_pme_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_pme = register_diag_field ('ocean_model','wdian_rho_pme',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to pme', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_pme > 0) compute_watermass_diag=.true.
id_wdian_rho_pme_on_nrho = register_diag_field ('ocean_model', &
'wdian_rho_pme_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport due to pme as binned to neutral density', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_pme_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_pme = register_diag_field ('ocean_model', &
'tform_rho_pme', Dens%neutralrho_axes(1:3), Time%model_time, &
'water mass transformation from pme on levels (pre-layer binning)',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_pme > 0) compute_watermass_diag=.true.
id_tform_rho_pme_on_nrho = register_diag_field ('ocean_model', &
'tform_rho_pme_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time,&
'water mass transformation from pme as binned to neutral density', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_pme_on_nrho > 0) compute_watermass_diag=.true.
! contributions from pme in advective form using kinematic based method
id_neut_rho_pbl_pme_kn = register_diag_field ('ocean_model', &
'neut_rho_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'kinematic advective-form material time derivative from pme', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_pbl_pme_kn > 0) compute_watermass_diag=.true.
id_neut_rho_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'neut_rho_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'kinematic advective-form material time derivative from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_pbl_pme_kn = register_diag_field ('ocean_model', &
'wdian_rho_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'kinematic advective-form dianeutral transport from pme', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_pbl_pme_kn > 0) compute_watermass_diag=.true.
id_wdian_rho_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'wdian_rho_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'kinematic advective-form dianeutral transport from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_pbl_pme_kn = register_diag_field ('ocean_model', &
'tform_rho_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'kinematic advective-form transform from pme pre-binning to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_pbl_pme_kn > 0) compute_watermass_diag=.true.
id_tform_rho_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'tform_rho_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'kinematic advective-form transform from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
! contributions from pme temp pieces in advective form using kinematic based method
id_neut_temp_pbl_pme_kn = register_diag_field ('ocean_model', &
'neut_temp_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'temp related kinematic advective-form material time derivative from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_pbl_pme_kn > 0) compute_watermass_diag=.true.
id_neut_temp_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'neut_temp_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related kinematic advective-form material time derivative from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
idwdian_temp_pbl_pme_kn = register_diag_field ('ocean_model', &
'wdian_temp_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'temp related kinematic advective-form dianeutral transport from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_temp_pbl_pme_kn > 0) compute_watermass_diag=.true.
idwdian_temp_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'wdian_temp_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related kinematic advective-form dianeutral transport from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_temp_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
idtform_temp_pbl_pme_kn = register_diag_field ('ocean_model', &
'tform_temp_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'temp related kinematic advective-form transform from pme pre-binning to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_temp_pbl_pme_kn > 0) compute_watermass_diag=.true.
idtform_temp_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'tform_temp_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related kinematic advective-form transform from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_temp_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
! contributions from pme salt pieces in advective form using kinematic based method
id_neut_salt_pbl_pme_kn = register_diag_field ('ocean_model', &
'neut_salt_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'salt related kinematic advective-form material time derivative from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_pbl_pme_kn > 0) compute_watermass_diag=.true.
id_neut_salt_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'neut_salt_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related kinematic advective-form material time derivative from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
idwdian_salt_pbl_pme_kn = register_diag_field ('ocean_model', &
'wdian_salt_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'salt related kinematic advective-form dianeutral transport from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_salt_pbl_pme_kn > 0) compute_watermass_diag=.true.
idwdian_salt_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'wdian_salt_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related kinematic advective-form dianeutral transport from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_salt_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
idtform_salt_pbl_pme_kn = register_diag_field ('ocean_model', &
'tform_salt_pbl_pme_kn', Grd%tracer_axes(1:3), Time%model_time, &
'salt related kinematic advective-form transform from pme pre-binning to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_salt_pbl_pme_kn > 0) compute_watermass_diag=.true.
idtform_salt_pbl_pme_kn_on_nrho = register_diag_field ('ocean_model', &
'tform_salt_pbl_pme_kn_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related kinematic advective-form transform from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_salt_pbl_pme_kn_on_nrho > 0) compute_watermass_diag=.true.
! contributions from pme in advective form using process based method
id_neut_rho_pbl_pme_pr = register_diag_field ('ocean_model', &
'neut_rho_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'process advective-form material time derivative from pme', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_pbl_pme_pr > 0) compute_watermass_diag=.true.
id_neut_rho_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'neut_rho_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'process advective-form material time derivative from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
id_wdian_rho_pbl_pme_pr = register_diag_field ('ocean_model', &
'wdian_rho_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'process advective-form dianeutral transport from pme', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_pbl_pme_pr > 0) compute_watermass_diag=.true.
id_wdian_rho_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'wdian_rho_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'process advective-form dianeutral transport from pme as binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
id_tform_rho_pbl_pme_pr = register_diag_field ('ocean_model', &
'tform_rho_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'process advective-form water mass transform from pme pre-binning to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_pbl_pme_pr > 0) compute_watermass_diag=.true.
id_tform_rho_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'tform_rho_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'process advective-form water mass transform from pme binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
! contributions from pme temp terms in advective form using process based method
id_neut_temp_pbl_pme_pr = register_diag_field ('ocean_model', &
'neut_temp_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'temp related process advective-form material time derivative from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_pbl_pme_pr > 0) compute_watermass_diag=.true.
id_neut_temp_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'neut_temp_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related process advective-form material time derivative from pme binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
idwdian_temp_pbl_pme_pr = register_diag_field ('ocean_model', &
'wdian_temp_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'temp related process advective-form dianeutral transport from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_temp_pbl_pme_pr > 0) compute_watermass_diag=.true.
idwdian_temp_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'wdian_temp_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related process advective-form dianeutral transport from pme binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_temp_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
idtform_temp_pbl_pme_pr = register_diag_field ('ocean_model', &
'tform_temp_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'temp related process advective-form water mass transform from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_temp_pbl_pme_pr > 0) compute_watermass_diag=.true.
idtform_temp_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'tform_temp_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related process advective-form water mass transform from pme binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_temp_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
! contributions from pme salinity terms in advective form using process based method
id_neut_salt_pbl_pme_pr = register_diag_field ('ocean_model', &
'neut_salt_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'salt related process advective-form material time derivative from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_pbl_pme_pr > 0) compute_watermass_diag=.true.
id_neut_salt_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'neut_salt_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related process advective-form material time derivative from pme binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
idwdian_salt_pbl_pme_pr = register_diag_field ('ocean_model', &
'wdian_salt_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'salt related process advective-form dianeutral transport from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_salt_pbl_pme_pr > 0) compute_watermass_diag=.true.
idwdian_salt_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'wdian_salt_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related process advective-form dianeutral transport from pme binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idwdian_salt_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
idtform_salt_pbl_pme_pr = register_diag_field ('ocean_model', &
'tform_salt_pbl_pme_pr', Grd%tracer_axes(1:3), Time%model_time, &
'salt related process advective-form water mass transform from pme',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_salt_pbl_pme_pr > 0) compute_watermass_diag=.true.
idtform_salt_pbl_pme_pr_on_nrho = register_diag_field ('ocean_model', &
'tform_salt_pbl_pme_pr_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related process advective-form water mass transform from pme binned to neutral density',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(idtform_salt_pbl_pme_pr_on_nrho > 0) compute_watermass_diag=.true.
if(vert_coordinate_class==DEPTH_BASED) then
id_mass_t_on_nrho = register_diag_field ('ocean_model', &
'mass_t_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'ocean t-cell volume (*rho0) as binned to neutral density layers',&
'kg', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_mass_t_on_nrho > 0) compute_watermass_diag=.true.
id_mass_t_tendency_on_nrho = register_diag_field ('ocean_model', &
'mass_t_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'tendency of ocean t-cell volume (*rho0) as binned to neutral density layers',&
'kg/s', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_mass_t_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_mass_nrho_layer = register_diag_field ('ocean_model', 'mass_nrho_layer', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'volume of a neutral rho layer (*rho0) computed via interpolation to get layer thickness', 'kg',&
missing_value=missing_value, range=(/0.0,1.e20/))
if(id_mass_nrho_layer > 0) compute_watermass_diag=.true.
id_mass_nrho_tendency_layer = register_diag_field ('ocean_model', &
'mass_nrho_tendency_layer', Dens%neutralrho_axes(1:3), Time%model_time, &
'tendency of volume of a neutral rho layer (*rho0) computed via interpolation to get layer thickness',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_mass_nrho_tendency_layer > 0) compute_watermass_diag=.true.
else
id_mass_t_on_nrho = register_diag_field ('ocean_model', &
'mass_t_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time,&
'ocean t-cell volume as binned to neutral density layers', &
'kg', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_mass_t_on_nrho > 0) compute_watermass_diag=.true.
id_mass_t_tendency_on_nrho = register_diag_field ('ocean_model', &
'mass_t_tendency_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'tendency of ocean t-cell volume as binned to neutral density layers', &
'kg/s', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_mass_t_tendency_on_nrho > 0) compute_watermass_diag=.true.
id_mass_nrho_layer = register_diag_field ('ocean_model', 'mass_nrho_layer', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'volume of a neutral rho layer computed via interpolation to get layer thickness', 'kg',&
missing_value=missing_value, range=(/0.0,1.e20/))
if(id_mass_nrho_layer > 0) compute_watermass_diag=.true.
id_mass_nrho_tendency_layer = register_diag_field ('ocean_model', &
'mass_nrho_tendency_layer', Dens%neutralrho_axes(1:3), Time%model_time, &
'tendency of volume of a neutral rho layer computed via interpolation to get layer thickness',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_mass_nrho_tendency_layer > 0) compute_watermass_diag=.true.
endif
if(compute_watermass_diag) then
write(stdoutunit,'(/a/)') &
'==>Note: running ocean_tracer_mod w/ compute_watermass_diag=.true.'
endif
end subroutine ocean_tracer_diagnostics_init
! NAME="ocean_tracer_diagnostics_init"
!#######################################################################
!
!
!
!
! For sending some tracer diagnostics
!
!
!
subroutine send_tracer_diagnostics(Time, T_prog, T_diag, Thickness, Dens, use_blobs)
type(ocean_time_type), intent(in) :: Time
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_diag_tracer_type), intent(in) :: T_diag(:)
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_density_type), intent(in) :: Dens
logical, intent(in) :: use_blobs
integer :: i,j,k,kbot,n
integer :: taum1,tau,taup1
real :: total_tracer
taum1 = Time%taum1
tau = Time%tau
taup1 = Time%taup1
do n=1,num_prog_tracers
if(id_prog(n) > 0) then
if( n==index_temp) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = T_prog(n)%field(i,j,k,tau) + Grd%tmask(i,j,k)*cmip_offset
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_prog(n), wrk1(:,:,:))
else
call diagnose_3d(Time, Grd, id_prog(n), T_prog(n)%field(:,:,:,tau))
endif
endif
if(id_prog_rhodzt(n) > 0) then
call diagnose_3d(Time, Grd, id_prog_rhodzt(n), Thickness%rho_dzt(:,:,:,tau)*T_prog(n)%field(:,:,:,tau))
endif
! vertical sum of tracer*rho_dzt
if(id_prog_int_rhodz(n) > 0) then
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = wrk1_2d(i,j) + Thickness%rho_dzt(i,j,k,tau)*T_prog(n)%field(i,j,k,tau)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_prog_int_rhodz(n), wrk1_2d(:,:))
endif
if (id_prog_on_depth(n) > 0) then
call remap_s_to_depth(Thickness, Time, T_prog(n)%field(:,:,:,tau), n)
endif
if (id_surf_tracer(n) > 0) then
if(n==index_temp ) then
wrk1_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = T_prog(n)%field(i,j,1,tau) + Grd%tmask(i,j,1)*cmip_offset
enddo
enddo
call diagnose_2d(Time, Grd, id_surf_tracer(n), wrk1_2d(:,:))
else
call diagnose_2d(Time, Grd, id_surf_tracer(n), T_prog(n)%field(:,:,1,tau))
endif
endif
if (id_surf_tracer_sq(n) > 0) then
if(n==index_temp) then
wrk1_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = T_prog(n)%field(i,j,1,tau) + Grd%tmask(i,j,1)*cmip_offset
enddo
enddo
call diagnose_2d(Time, Grd, id_surf_tracer_sq(n), wrk1_2d(:,:)**2)
else
call diagnose_2d(Time, Grd, id_surf_tracer_sq(n), T_prog(n)%field(:,:,1,tau)**2)
endif
endif
if (id_bott_tracer(n) > 0) then
wrk1_2d(:,:) = 0.0
do j = jsd,jed
do i = isd,ied
kbot = Grd%kmt(i,j)
if(kbot > 1) then
wrk1_2d(i,j) = T_prog(n)%field(i,j,kbot,tau)
endif
enddo
enddo
call diagnose_2d(Time, Grd, id_bott_tracer(n), wrk1_2d(:,:))
endif
! time tendency for tracer mass per horizontal area
if (id_tendency(n) > 0 .or. id_tendency_on_nrho(n) > 0) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtimer*T_prog(n)%conversion &
*(T_prog(n)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(n)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1))
enddo
enddo
enddo
if (id_tendency(n) > 0) then
call diagnose_3d(Time, Grd, id_tendency(n),wrk1(:,:,:))
endif
if (id_tendency_on_nrho(n) > 0) then
call diagnose_3d_rho(Time, Dens, id_tendency_on_nrho(n),wrk1)
endif
endif
! time tendency for tracer concentration
do k=1,nk
do j=jsd,jed
do i=isd,ied
T_prog(n)%tendency(i,j,k) = dtimer*(T_prog(n)%field(i,j,k,taup1)-T_prog(n)%field(i,j,k,taum1))
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_tendency_conc(n), T_prog(n)%tendency(:,:,:))
! send some diagnostics for the blob systems
if (use_blobs) then
if (id_progT(n) > 0) then
if(n==index_temp) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = T_prog(n)%fieldT(i,j,k) + Grd%tmask(i,j,k)*cmip_offset
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_progT(n), wrk1(:,:,:))
else
call diagnose_3d(Time, Grd, id_progT(n), T_prog(n)%fieldT(:,:,:))
endif
endif
if(id_prog_rhodztT(n) > 0) then
call diagnose_3d(Time, Grd, id_prog_rhodzt(n), Thickness%rho_dztT(:,:,:,tau)*T_prog(n)%fieldT(:,:,:))
endif
! vertical sum of tracer*rho_dzt
if(id_prog_int_rhodzT(n) > 0) then
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = wrk1_2d(i,j) + Thickness%rho_dztT(i,j,k,tau)*T_prog(n)%fieldT(i,j,k)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_prog_int_rhodzT(n), wrk1_2d(:,:))
endif
if(id_surf_tracerT(n) > 0) then
if(n==index_temp ) then
wrk1_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = T_prog(n)%fieldT(i,j,1) + Grd%tmask(i,j,1)*cmip_offset
enddo
enddo
call diagnose_2d(Time, Grd, id_surf_tracer_sq(n), wrk1_2d(:,:))
else
call diagnose_2d(Time, Grd, id_surf_tracer_sq(n), T_prog(n)%fieldT(:,:,1))
endif
endif
if (id_bott_tracerT(n) > 0) then
wrk1_2d(:,:) = 0.0
do j = jsd,jed
do i = isd,ied
kbot = Grd%kmt(i,j)
if(kbot > 1) then
wrk1_2d(i,j) = T_prog(n)%fieldT(i,j,kbot)*Grd%tmask(i,j,1)
endif
enddo
enddo
call diagnose_2d(Time, Grd, id_bott_tracerT(n),wrk1_2d(:,:))
endif
if (id_tendencyT(n) > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtimer*T_prog(n)%conversion &
*(T_prog(n)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(n)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1) &
+T_prog(n)%sum_blob(i,j,k,taup1)*Grd%datr(i,j) &
-T_prog(n)%sum_blob(i,j,k,taum1)*Grd%datr(i,j) )
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_tendencyT(n), wrk1(:,:,:))
endif
if (id_tendencyL(n) > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtimer*T_prog(n)%conversion &
*(T_prog(n)%sum_blob(i,j,k,taup1)*Grd%datr(i,j) &
-T_prog(n)%sum_blob(i,j,k,taum1)*Grd%datr(i,j) )
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_tendencyL(n),wrk1(:,:,:))
endif
do k=1,nk
do j=jsd,jed
do i=isd,ied
wrk1(i,j,k) = ( T_prog(n)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
+Grd%datr(i,j)*T_prog(n)%sum_blob(i,j,k,taup1) ) &
/(Thickness%rho_dztT(i,j,k,taup1)+epsln)
wrk2(i,j,k) = ( T_prog(n)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1) &
+Grd%datr(i,j)*T_prog(n)%sum_blob(i,j,k,taum1) ) &
/(Thickness%rho_dztT(i,j,k,taum1)+epsln)
T_prog(n)%tendency(i,j,k) = dtimer*(wrk1(i,j,k) - wrk2(i,j,k))
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_tendency_concT(n), T_prog(n)%tendency(:,:,:))
if (id_tendency_concL(n) > 0) then
do k=1,nk
wrk1(:,:,k) = Grd%datr(:,:)*T_prog(n)%sum_blob(:,:,k,taup1)/Thickness%rho_dztL(:,:,k,taup1)
wrk2(:,:,k) = Grd%datr(:,:)*T_prog(n)%sum_blob(:,:,k,taum1)/Thickness%rho_dztL(:,:,k,taum1)
enddo
call diagnose_3d(Time, Grd, id_tendency_concL(n), &
dtimer*(wrk1(:,:,:) - wrk2(:,:,:)))
endif
endif !use_blobs
call send_tracer_variance(Time, T_prog(n), Thickness, n)
enddo ! enddo for n-tracers
! send diagnostic tracers to diag manager
do n=1,num_diag_tracers
if(id_diag(n) > 0) then
if(n==index_temp) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = T_diag(n)%field(i,j,k) + Grd%tmask(i,j,k)*cmip_offset
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_diag(n), wrk1(:,:,:))
else
call diagnose_3d(Time, Grd, id_diag(n), T_diag(n)%field(:,:,:))
endif
endif
if (id_diag_surf_tracer(n) > 0) then
if(n==index_temp) then
wrk1_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = T_diag(n)%field(i,j,1) + Grd%tmask(i,j,1)*cmip_offset
enddo
enddo
call diagnose_2d(Time, Grd, id_diag_surf_tracer(n), wrk1_2d(:,:))
else
call diagnose_2d(Time, Grd, id_diag_surf_tracer(n), T_diag(n)%field(:,:,1))
endif
endif
if (id_diag_surf_tracer_sq(n) > 0) then
if(n==index_temp) then
wrk1_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = T_diag(n)%field(i,j,1) + Grd%tmask(i,j,1)*cmip_offset
enddo
enddo
call diagnose_2d(Time, Grd, id_diag_surf_tracer_sq(n), wrk1_2d(:,:)**2)
else
call diagnose_2d(Time, Grd, id_diag_surf_tracer_sq(n), T_diag(n)%field(:,:,1)**2)
endif
endif
if(id_diag_total(n) > 0) then
total_tracer=0.0
do k=1,nk
wrk1_2d(:,:) = T_diag(n)%conversion*Grd%tmask(:,:,k)*Grd%dat(:,:) &
*Thickness%rho_dzt(:,:,k,tau)*T_diag(n)%field(:,:,k)
if(have_obc) wrk1_2d(:,:) = wrk1_2d(:,:)*Grd%obc_tmask(:,:)
total_tracer = total_tracer + mpp_global_sum(Dom%domain2d,wrk1_2d(:,:), NON_BITWISE_EXACT_SUM)
enddo
used = send_data (id_diag_total(n), total_tracer*1e-18, Time%model_time)
endif
enddo ! enddo for n=1,num_diag_tracers
end subroutine send_tracer_diagnostics
! NAME="send_tracer_diagnostics"
!#######################################################################
!
!
!
! Diagnose effects from interior ocean vertical mixing of
! temp and salt on the watermass transformation diagnostics.
!
! This routine is called prior to implicit update of the tracer
! fields, so that taup1 value contains only explicit in-time
! tendencies. The implicit in-time tendencies are diagnosed
! in this routine by various calls to invtri using same methods
! as for the prognostic calculation.
!
!
!
subroutine watermass_diag(Time, T_prog, T_diag, Dens, Thickness, pme)
type(ocean_time_type), intent(in) :: Time
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
type(ocean_diag_tracer_type), intent(in) :: T_diag(:)
type(ocean_density_type), intent(in) :: Dens
type(ocean_thickness_type), intent(in) :: Thickness
real, dimension(isd:,jsd:), intent(in) :: pme
integer :: i,j,k,taum1,tau,taup1
real, dimension(isd:ied,jsd:jed) :: eta_tend
if(.not. compute_watermass_diag) return
taum1 = Time%taum1
tau = Time%tau
taup1 = Time%taup1
if (id_pot_rho_tendency > 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%dpotrhodT(i,j,k)* &
(T_prog(index_temp)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(index_temp)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1))
wrk2(i,j,k) = Dens%dpotrhodS(i,j,k)* &
(T_prog(index_salt)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(index_salt)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1))
wrk3(i,j,k) = Grd%tmask(i,j,k)*dtimer*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_pot_rho_tendency, wrk4(:,:,:))
endif
! full time tendency
if (id_neut_rho_tendency > 0 .or. id_neut_rho_tendency_on_nrho > 0 .or. &
id_wdian_rho_tendency > 0 .or. id_wdian_rho_tendency_on_nrho > 0 .or. &
id_tform_rho_tendency > 0 .or. id_tform_rho_tendency_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
wrk6(:,:,:) = 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)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(index_temp)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1))
wrk2(i,j,k) = Dens%drhodS(i,j,k)* &
(T_prog(index_salt)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(index_salt)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1))
wrk3(i,j,k) = Grd%tmask(i,j,k)*dtimer*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk5(i,j,k) = wrk4(i,j,k)*Dens%stratification_factor(i,j,k)
wrk6(i,j,k) = wrk3(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_rho_tendency, wrk4(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_tendency, wrk5(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_tendency, wrk6(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_tendency_on_nrho, wrk4)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_tendency_on_nrho, wrk5)
call diagnose_3d_rho(Time, Dens, id_tform_rho_tendency_on_nrho, wrk6)
endif
! time tendency arising from temperature effects
if (id_neut_temp_tendency > 0 .or. id_neut_temp_tendency_on_nrho > 0 .or. &
id_wdian_temp_tendency > 0 .or. id_wdian_temp_tendency_on_nrho > 0 .or. &
id_tform_temp_tendency > 0 .or. id_tform_temp_tendency_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
wrk6(:,:,:) = 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)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(index_temp)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1))
wrk2(i,j,k) = 0.0
wrk3(i,j,k) = Grd%tmask(i,j,k)*dtimer*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk5(i,j,k) = wrk4(i,j,k)*Dens%stratification_factor(i,j,k)
wrk6(i,j,k) = wrk3(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_temp_tendency,wrk4(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_temp_tendency, wrk5(:,:,:))
call diagnose_3d(Time, Grd, id_tform_temp_tendency, wrk6(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_tendency_on_nrho, wrk4)
call diagnose_3d_rho(Time, Dens, id_wdian_temp_tendency_on_nrho, wrk5)
call diagnose_3d_rho(Time, Dens, id_tform_temp_tendency_on_nrho, wrk6)
endif
! time tendency arising from salinity effects
if (id_neut_salt_tendency > 0 .or. id_neut_salt_tendency_on_nrho > 0 .or. &
id_wdian_salt_tendency > 0 .or. id_wdian_salt_tendency_on_nrho > 0 .or. &
id_tform_salt_tendency > 0 .or. id_tform_salt_tendency_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
wrk6(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = 0.0
wrk2(i,j,k) = Dens%drhodS(i,j,k)* &
(T_prog(index_salt)%field(i,j,k,taup1)*Thickness%rho_dzt(i,j,k,taup1) &
-T_prog(index_salt)%field(i,j,k,taum1)*Thickness%rho_dzt(i,j,k,taum1))
wrk3(i,j,k) = Grd%tmask(i,j,k)*dtimer*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk5(i,j,k) = wrk4(i,j,k)*Dens%stratification_factor(i,j,k)
wrk6(i,j,k) = wrk3(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_salt_tendency, wrk4(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_salt_tendency, wrk5(:,:,:))
call diagnose_3d(Time, Grd, id_tform_salt_tendency, wrk6(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_tendency_on_nrho, wrk4)
call diagnose_3d_rho(Time, Dens, id_wdian_salt_tendency_on_nrho, wrk5)
call diagnose_3d_rho(Time, Dens, id_tform_salt_tendency_on_nrho, wrk6)
endif
! effects from smooth eta or pbot
if (id_neut_rho_smooth > 0 .or. id_neut_rho_smooth_on_nrho > 0 .or. &
id_wdian_rho_smooth > 0 .or. id_wdian_rho_smooth_on_nrho > 0 .or. &
id_tform_rho_smooth > 0 .or. id_tform_rho_smooth_on_nrho > 0 .or. &
id_eta_tend_smooth > 0 .or. id_eta_tend_smooth_glob > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
wrk6(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k) &
*( Dens%drhodT(i,j,k)*T_prog(index_temp)%eta_smooth(i,j) &
+Dens%drhodS(i,j,k)*T_prog(index_salt)%eta_smooth(i,j))
enddo
enddo
do j=jsc,jec
do i=isc,iec
k=Grd%kmt(i,j)
if(k > 0) then
wrk2(i,j,k) = Grd%tmask(i,j,k) &
*( Dens%drhodT(i,j,k)*T_prog(index_temp)%pbot_smooth(i,j) &
+Dens%drhodS(i,j,k)*T_prog(index_salt)%pbot_smooth(i,j))
endif
enddo
enddo
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk3(i,j,k) = Grd%tmask(i,j,k)*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk5(i,j,k) = wrk4(i,j,k)*Dens%stratification_factor(i,j,k)
wrk6(i,j,k) = wrk3(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
wrk1(i,j,k) =-wrk3(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_smooth, wrk4(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_smooth, wrk5(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_smooth, wrk6(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_smooth_on_nrho, wrk4)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_smooth_on_nrho, wrk5)
call diagnose_3d_rho(Time, Dens, id_tform_rho_smooth_on_nrho, wrk6)
if(id_eta_tend_smooth > 0 .or. id_eta_tend_smooth_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) + wrk1(i,j,k)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_smooth, eta_tend(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_smooth_glob, eta_tend, cellarea_r)
endif
endif
! effects from smooth eta or pbot on potential density
if (id_pot_rho_smooth > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k) &
*( Dens%dpotrhodT(i,j,k)*T_prog(index_temp)%eta_smooth(i,j) &
+Dens%dpotrhodS(i,j,k)*T_prog(index_salt)%eta_smooth(i,j))
enddo
enddo
do j=jsc,jec
do i=isc,iec
k=Grd%kmt(i,j)
if(k > 0) then
wrk2(i,j,k) = Grd%tmask(i,j,k) &
*( Dens%dpotrhodT(i,j,k)*T_prog(index_temp)%pbot_smooth(i,j) &
+Dens%dpotrhodS(i,j,k)*T_prog(index_salt)%pbot_smooth(i,j))
endif
enddo
enddo
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk3(i,j,k) = Grd%tmask(i,j,k)*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_pot_rho_smooth, wrk4(:,:,:))
endif
! effects from smooth eta or pbot due to temperature contributions
if (id_neut_temp_smooth > 0 .or. id_neut_temp_smooth_on_nrho > 0 .or. &
id_wdian_temp_smooth > 0 .or. id_wdian_temp_smooth_on_nrho > 0 .or. &
id_tform_temp_smooth > 0 .or. id_tform_temp_smooth_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
wrk6(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k) &
*( Dens%drhodT(i,j,k)*T_prog(index_temp)%eta_smooth(i,j) &
+0.0)
enddo
enddo
do j=jsc,jec
do i=isc,iec
k=Grd%kmt(i,j)
if(k > 0) then
wrk2(i,j,k) = Grd%tmask(i,j,k) &
*( Dens%drhodT(i,j,k)*T_prog(index_temp)%pbot_smooth(i,j) &
+0.0)
endif
enddo
enddo
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk3(i,j,k) = Grd%tmask(i,j,k)*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk5(i,j,k) = wrk4(i,j,k)*Dens%stratification_factor(i,j,k)
wrk6(i,j,k) = wrk3(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_temp_smooth, wrk4(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_temp_smooth, wrk5(:,:,:))
call diagnose_3d(Time, Grd, id_tform_temp_smooth, wrk6(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_smooth_on_nrho, wrk4)
call diagnose_3d_rho(Time, Dens, id_wdian_temp_smooth_on_nrho, wrk5)
call diagnose_3d_rho(Time, Dens, id_tform_temp_smooth_on_nrho, wrk6)
endif
! effects from smooth eta or pbot due to salinty contributions
if (id_neut_salt_smooth > 0 .or. id_neut_salt_smooth_on_nrho > 0 .or. &
id_wdian_salt_smooth > 0 .or. id_wdian_salt_smooth_on_nrho > 0 .or. &
id_tform_salt_smooth > 0 .or. id_tform_salt_smooth_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
wrk6(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k) &
*( 0.0 &
+Dens%drhodS(i,j,k)*T_prog(index_salt)%eta_smooth(i,j))
enddo
enddo
do j=jsc,jec
do i=isc,iec
k=Grd%kmt(i,j)
if(k > 0) then
wrk2(i,j,k) = Grd%tmask(i,j,k) &
*( 0.0 &
+Dens%drhodS(i,j,k)*T_prog(index_salt)%pbot_smooth(i,j))
endif
enddo
enddo
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk3(i,j,k) = Grd%tmask(i,j,k)*(wrk1(i,j,k)+wrk2(i,j,k))
wrk4(i,j,k) = wrk3(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk5(i,j,k) = wrk4(i,j,k)*Dens%stratification_factor(i,j,k)
wrk6(i,j,k) = wrk3(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_salt_smooth,wrk4(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_salt_smooth, wrk5(:,:,:))
call diagnose_3d(Time, Grd, id_tform_salt_smooth, wrk6(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_smooth_on_nrho, wrk4)
call diagnose_3d_rho(Time, Dens, id_wdian_salt_smooth_on_nrho, wrk5)
call diagnose_3d_rho(Time, Dens, id_tform_salt_smooth_on_nrho, wrk6)
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to pme. formulated here according to flux form of material time derivative.
! note that river contribution is diagnosed in ocean_rivermix.F90.
if (id_neut_rho_pme > 0 .or. id_neut_rho_pme_on_nrho > 0 .or. &
id_wdian_rho_pme > 0 .or. id_wdian_rho_pme_on_nrho > 0 .or. &
id_tform_rho_pme > 0 .or. id_tform_rho_pme_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*( Dens%drhodT(i,j,k)*T_prog(index_temp)%tpme(i,j) &
+Dens%drhodS(i,j,k)*T_prog(index_salt)%tpme(i,j))
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
call diagnose_3d(Time, Grd, id_neut_rho_pme, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_pme, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_pme, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_pme_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_pme_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_pme_on_nrho, wrk4)
endif
! time tendency for potrho and dianeutral velocity component due
! to pme. formulated here according to flux form of material time derivative.
! note that river contribution is diagnosed in ocean_rivermix.F90.
if (id_pot_rho_pme > 0 ) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*( Dens%dpotrhodT(i,j,k)*T_prog(index_temp)%tpme(i,j) &
+Dens%dpotrhodS(i,j,k)*T_prog(index_salt)%tpme(i,j))
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
enddo
enddo
call diagnose_3d(Time, Grd, id_pot_rho_pme, wrk2(:,:,:))
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to pme induced buoyancy flux. written here according to advective form
! of material time derivative using the kinematic formulation.
! note that river contribution is diagnosed in ocean_rivermix.F90.
if (id_neut_rho_pbl_pme_kn > 0 .or. id_neut_rho_pbl_pme_kn_on_nrho > 0 .or. &
id_wdian_rho_pbl_pme_kn > 0 .or. id_wdian_rho_pbl_pme_kn_on_nrho > 0 .or. &
id_tform_rho_pbl_pme_kn > 0 .or. id_tform_rho_pbl_pme_kn_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*(Dens%drhodT(i,j,k)*(0.0-T_prog(index_temp)%field(i,j,k,tau)) &
+Dens%drhodS(i,j,k)*(0.0-T_prog(index_salt)%field(i,j,k,tau)))
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
call diagnose_3d(Time, Grd, id_neut_rho_pbl_pme_kn, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_pbl_pme_kn, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_pbl_pme_kn, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_pbl_pme_kn_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_pbl_pme_kn_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_pbl_pme_kn_on_nrho, wrk4)
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to pme induced heat flux. written here according to advective form
! of material time derivative using the kinematic formulation.
! note that river contribution is diagnosed in ocean_rivermix.F90.
if (id_neut_temp_pbl_pme_kn > 0 .or. id_neut_temp_pbl_pme_kn_on_nrho > 0 .or. &
idwdian_temp_pbl_pme_kn > 0 .or. idwdian_temp_pbl_pme_kn_on_nrho > 0 .or. &
idtform_temp_pbl_pme_kn > 0 .or. idtform_temp_pbl_pme_kn_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*(Dens%drhodT(i,j,k)*(0.0-T_prog(index_temp)%field(i,j,k,tau)) &
+0.0)
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
call diagnose_3d(Time, Grd, id_neut_temp_pbl_pme_kn, wrk2(:,:,:))
call diagnose_3d(Time, Grd, idwdian_temp_pbl_pme_kn, wrk3(:,:,:))
call diagnose_3d(Time, Grd, idtform_temp_pbl_pme_kn, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_pbl_pme_kn_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, idwdian_temp_pbl_pme_kn_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, idtform_temp_pbl_pme_kn_on_nrho, wrk4)
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to pme induced salt flux. written here according to advective form
! of material time derivative using the kinematic formulation.
! note that river contribution is diagnosed in ocean_rivermix.F90.
if (id_neut_salt_pbl_pme_kn > 0 .or. id_neut_salt_pbl_pme_kn_on_nrho > 0 .or. &
idwdian_salt_pbl_pme_kn > 0 .or. idwdian_salt_pbl_pme_kn_on_nrho > 0 .or. &
idtform_salt_pbl_pme_kn > 0 .or. idtform_salt_pbl_pme_kn_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*(0.0 &
+Dens%drhodS(i,j,k)*(0.0-T_prog(index_salt)%field(i,j,k,tau)))
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
call diagnose_3d(Time, Grd, id_neut_salt_pbl_pme_kn, wrk2(:,:,:))
call diagnose_3d(Time, Grd, idwdian_salt_pbl_pme_kn, wrk3(:,:,:))
call diagnose_3d(Time, Grd, idtform_salt_pbl_pme_kn, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_pbl_pme_kn_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, idwdian_salt_pbl_pme_kn_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, idtform_salt_pbl_pme_kn_on_nrho, wrk4)
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to pme induced buoyancy flux. written here according to advective form
! of material time derivative using the process formulation.
! note that river contribution is diagnosed in ocean_rivermix.F90.
if (id_neut_rho_pbl_pme_pr > 0 .or. id_neut_rho_pbl_pme_pr_on_nrho > 0 .or. &
id_wdian_rho_pbl_pme_pr > 0 .or. id_wdian_rho_pbl_pme_pr_on_nrho > 0 .or. &
id_tform_rho_pbl_pme_pr > 0 .or. id_tform_rho_pbl_pme_pr_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*(Dens%drhodT(i,j,k)*(T_prog(index_temp)%tpme(i,j)-T_prog(index_temp)%field(i,j,k,tau)) &
+Dens%drhodS(i,j,k)*(T_prog(index_salt)%tpme(i,j)-T_prog(index_salt)%field(i,j,k,tau)))
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
call diagnose_3d(Time, Grd, id_neut_rho_pbl_pme_pr, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_pbl_pme_pr, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_pbl_pme_pr, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_pbl_pme_pr_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_pbl_pme_pr_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_pbl_pme_pr_on_nrho, wrk4)
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to pme induced buoyancy flux. written here according to advective form
! of material time derivative using the process formulation.
! note that river contribution is diagnosed in ocean_rivermix.F90.
! temperature contributions are diagnosed here.
if (id_neut_temp_pbl_pme_pr > 0 .or. id_neut_temp_pbl_pme_pr_on_nrho > 0 .or. &
idwdian_temp_pbl_pme_pr > 0 .or. idwdian_temp_pbl_pme_pr_on_nrho > 0 .or. &
idtform_temp_pbl_pme_pr > 0 .or. idtform_temp_pbl_pme_pr_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*(Dens%drhodT(i,j,k)*(T_prog(index_temp)%tpme(i,j)-T_prog(index_temp)%field(i,j,k,tau)) &
+0.0)
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
call diagnose_3d(Time, Grd, id_neut_temp_pbl_pme_pr, wrk2(:,:,:))
call diagnose_3d(Time, Grd, idwdian_temp_pbl_pme_pr, wrk3(:,:,:))
call diagnose_3d(Time, Grd, idtform_temp_pbl_pme_pr, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_pbl_pme_pr_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, idwdian_temp_pbl_pme_pr_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, idtform_temp_pbl_pme_pr_on_nrho, wrk4)
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to pme induced buoyancy flux. written here according to advective form
! of material time derivative using the process formulation.
! note that river contribution is diagnosed in ocean_rivermix.F90.
! salinity contributions are diagnosed here.
if (id_neut_salt_pbl_pme_pr > 0 .or. id_neut_salt_pbl_pme_pr_on_nrho > 0 .or. &
idwdian_salt_pbl_pme_pr > 0 .or. idwdian_salt_pbl_pme_pr_on_nrho > 0 .or. &
idtform_salt_pbl_pme_pr > 0 .or. idtform_salt_pbl_pme_pr_on_nrho > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*pme(i,j) &
*(0.0 &
+Dens%drhodS(i,j,k)*(T_prog(index_salt)%tpme(i,j)-T_prog(index_salt)%field(i,j,k,tau)))
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
call diagnose_3d(Time, Grd, id_neut_salt_pbl_pme_pr, wrk2(:,:,:))
call diagnose_3d(Time, Grd, idwdian_salt_pbl_pme_pr, wrk3(:,:,:))
call diagnose_3d(Time, Grd, idtform_salt_pbl_pme_pr, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_pbl_pme_pr_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, idwdian_salt_pbl_pme_pr_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, idtform_salt_pbl_pme_pr_on_nrho, wrk4)
endif
! time tendency for locally ref potrho and dianeutral velocity component due
! to frazil formation, which is a heating of the liquid ocean.
if(index_frazil > 0) then
if (id_neut_rho_frazil > 0 .or. id_neut_rho_frazil_on_nrho > 0 .or. &
id_wdian_rho_frazil > 0 .or. id_wdian_rho_frazil_on_nrho > 0 .or. &
id_tform_rho_frazil > 0 .or. id_tform_rho_frazil_on_nrho > 0 .or. &
id_eta_tend_frazil > 0 .or. id_eta_tend_frazil_glob > 0 .or. &
id_frazil_on_nrho > 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)*dtimer*T_diag(index_frazil)%field(i,j,k)*cp_oceanr*Dens%drhodT(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_frazil, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_frazil, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_frazil, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_frazil_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_frazil_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_frazil_on_nrho, wrk4)
call diagnose_3d_rho(Time, Dens, id_frazil_on_nrho, T_diag(index_frazil)%field*dtimer)
if(id_eta_tend_frazil > 0 .or. id_eta_tend_frazil_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_frazil, eta_tend(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_frazil_glob, eta_tend, cellarea_r)
endif
endif
! frazil contribution to potential density tendency
if (id_pot_rho_frazil > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%tmask(i,j,k)*dtimer*T_diag(index_frazil)%field(i,j,k)*cp_oceanr*Dens%dpotrhodT(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_pot_rho_frazil, wrk2(:,:,:))
endif
endif ! endif for index_frazil > 0
! some further layer diagnostics for mass changes.
! these diagnostics do not lead to robust results as a function of
! layer binning. User should thus beware.
if(id_mass_t_on_nrho > 0 .or. id_mass_t_tendency_on_nrho > 0) then
! get mass of layers at time taup1 via binning mass_t.
nrho_work(:,:,:) = 0.0
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%dat(i,j)*Thickness%rho_dzt(i,j,k,taup1)
enddo
enddo
enddo
wrk3(:,:,:) = Grd%tmask(:,:,:)*neutral_density(T_prog(index_salt)%field(:,:,:,taup1), &
T_prog(index_temp)%field(:,:,:,taup1))
call rebin_onto_rho (Dens%neutralrho_bounds, wrk3, wrk1, nrho_work)
! we diagnose mass_t_on_nrho at taup1 in order to be
! consistent with mass_nrho_layer diagnosed below.
! The mass computed with mass_nrho_layer method can
! only get the mass(taup1) since that is the time
! where the level thicknesses are known.
if(id_mass_t_on_nrho > 0 ) then
used = send_data (id_mass_t_on_nrho, nrho_work(:,:,:), &
Time%model_time, &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=neutralrho_nk)
endif
! bin mass_t at two adjacent time steps to estimate the layer tendency
if(id_mass_t_tendency_on_nrho > 0) then
nrho_work2(:,:,:) = 0.0
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%dat(i,j)*Thickness%rho_dzt(i,j,k,tau)
enddo
enddo
enddo
! get mass of layers at time taup1 via binning mass_t.
! we already have the mass at time taup1 contained in nrho_work.
! note that if mass_t(tau)=mass_t(taup1), then only means for getting a
! layer tendency is if the layer interfaces evolve. Indeed, such is the
! dominant effect, since interior grid cells change their mass only a
! tiny amount, since interior ocean grid cells are nearly 3d non-divergent.
call rebin_onto_rho (Dens%neutralrho_bounds, Dens%neutralrho, wrk1, nrho_work2)
! compute time tendency for layer mass.
do k=1,neutralrho_nk
do j=jsc,jec
do i=isc,iec
nrho_work(i,j,k) = (nrho_work(i,j,k) - nrho_work2(i,j,k))*dtimer
enddo
enddo
enddo
used = send_data (id_mass_t_tendency_on_nrho, nrho_work(:,:,:),&
Time%model_time, &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=neutralrho_nk)
endif
endif ! endif corresponding to if(id_mass_t_on_nrho > 0 .or. id_mass_t_tendency_on_nrho > 0)
! estimate layer mass by estimating mass per area of the layer interfaces.
! note that we do not save nrho_mass in a restart file. So there will be
! a slight blip at the start of a new model segment. This blip is not
! so important, since we are only working with diagnostics here.
if(id_mass_nrho_layer > 0 .or. id_mass_nrho_tendency_layer > 0) then
nrho_work(:,:,:) = nrho_mass(:,:,:) ! save previous time step layer mass
nrho_mass(:,:,:) = 0.0 ! for the updated layer mass
wrk3(:,:,:) = Grd%tmask(:,:,:)*neutral_density(T_prog(index_salt)%field(:,:,:,taup1), &
T_prog(index_temp)%field(:,:,:,taup1))
! note that all thickness arrays are at taup1, so we
! compute mass(taup1) and mass_of_layer(taup1).
call diagnose_mass_of_layer (Grd%dat, &
Thickness%dzt, Thickness%dztlo, Thickness%dztup, &
Thickness%rho_dzt(:,:,:,taup1), &
wrk3, neutralrho_nk, Dens%neutralrho_bounds, nrho_mass)
if(id_mass_nrho_layer > 0) then
used = send_data (id_mass_nrho_layer, nrho_mass(:,:,:), &
Time%model_time, &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=neutralrho_nk)
endif
if(id_mass_nrho_tendency_layer > 0) then
! first time step need to compute mass_of_layer using the binning
! approach, since we do not have access to the tau values of all
! relevant Thickness arrays.
if(Time%itt <= 1) then
nrho_work(:,:,:) = 0.0
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Grd%dat(i,j)*Thickness%rho_dzt(i,j,k,tau)
enddo
enddo
enddo
call rebin_onto_rho (Dens%neutralrho_bounds, Dens%neutralrho, wrk1, nrho_work)
endif
do k=1,neutralrho_nk
do j=jsc,jec
do i=isc,iec
nrho_work(i,j,k) = (nrho_mass(i,j,k) - nrho_work(i,j,k))*dtimer
enddo
enddo
enddo
used = send_data (id_mass_nrho_tendency_layer, nrho_work(:,:,:),&
Time%model_time, &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=neutralrho_nk)
endif
endif
end subroutine watermass_diag
! NAME="watermass_diag"
end module ocean_tracer_mod