!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module ocean_tracer_mod
!chd Dirty fix to spinup baltic
!chd temperature and salinity are constrained to
!chd extreme values; tracer conservation is violated
!
! Stephen Griffies
!
!
! Matt Harrison
!
!
! Richard D. Slater (initialization)
!
!
! Ron Pacanowski
!
!
!
! 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=.true.
!
!
!
! 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 mom4 preprocessing
! schemes, onto s-coordinates. This method is of use for initializing
! terrain following coordinate simulations with mom4.
!
!
!
! 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.
!
!
!
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, write_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: 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
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_convect_mod, only: convection
use ocean_domains_mod, only: get_local_indices
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
#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
use ocean_parameters_mod, only: missing_value, sec_in_yr_r, rho0
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_passive_mod, only: passive_tracer_init, update_tracer_passive
use ocean_polar_filter_mod, only: polar_filter_tracers
use ocean_shortwave_mod, only: ocean_irradiance_init
use ocean_tempsalt_mod, only: contemp_from_pottemp, pottemp_from_contemp
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_util_mod, only: tracer_prog_chksum, tracer_diag_chksum, tracer_min_max
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
use ocean_types_mod, only: ocean_public_type, ocean_density_type, ocean_options_type
use ocean_util_mod, only: write_timestamp
use ocean_vert_mix_mod, only: vert_diffuse, vert_diffuse_implicit
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, 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 1.1.2.2.2.43.20.6.18.3.2.3.10.6.20.1.4.2 2009/12/04 14:39:12 smg Exp $'
character(len=256) :: tagname='Tag $Name: mom4p1_pubrel_dec2009_nnz $'
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
integer :: index_temp =-1
integer :: index_salt =-1
integer :: index_temp_sq =-1
integer :: index_salt_sq =-1
integer :: index_diag_temp =-1
! for obc
logical :: have_obc=.false.
! for setting the temperature variable
integer :: prog_temp_variable=0
! for CMIP units (deg K rather than degC)
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_polar_filter
! 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(:)
! for diagnostics
logical :: used
integer, allocatable, dimension(:) :: id_eta_smooth
integer, allocatable, dimension(:) :: id_pbot_smooth
integer, allocatable, dimension(:) :: id_prog
integer, allocatable, dimension(:) :: id_prog_explicit
integer, allocatable, dimension(:) :: id_prog_rhodzt
integer, allocatable, dimension(:) :: id_prog_int_rhodz
integer, allocatable, dimension(:) :: id_prog_on_depth
integer, allocatable, dimension(:) :: id_tendency_conc
integer, allocatable, dimension(:) :: id_tendency
integer, allocatable, dimension(:) :: id_tendency_expl
integer, allocatable, dimension(:) :: id_surf_tracer
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_diag
integer, allocatable, dimension(:) :: id_diag_total
integer, allocatable, dimension(:) :: id_tmask_limit
integer :: id_neutral_rho_tendency
integer :: id_wdian_rho_tendency
integer :: id_neutral_rho_smooth
integer :: id_wdian_rho_smooth
integer :: id_neutral_rho_pme
integer :: id_wdian_rho_pme
! for ascii output
integer :: unit=6
public update_ocean_tracer
public ocean_prog_tracer_init
public ocean_diag_tracer_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
!---------------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 = .true.
logical :: frazil_heating_before_vphysics = .false.
logical :: frazil_heating_after_vphysics = .false.
logical :: compute_tmask_limit_on = .true.
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
contains
!#######################################################################
!
!
!
! Initialization code for prognostic tracers, returning a pointer to
! the T_prog array.
!
!
function ocean_prog_tracer_init (Grid, Thickness, Ocean, Ocean_options, Domain, Time, Time_steps, &
num_prog, vert_coordinate_type, obc, cmip_units, debug) &
result (T_prog)
type(ocean_grid_type), intent(in), target :: Grid
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_public_type), intent(inout) :: Ocean
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
logical, intent(in), optional :: debug
! return value
type(ocean_prog_tracer_type), dimension(:), pointer :: T_prog
integer :: i, j, k, n, kb, length, 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, longname, units, type
character(len=128) :: filename
logical :: initialize_as_a_passive_tracer=.false.
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
! provide for namelist over-ride
ioun = open_namelist_file()
read (ioun, ocean_tracer_nml,iostat=io_status)
write (stdoutunit,'(/)')
write (stdoutunit, ocean_tracer_nml)
write (stdlogunit, ocean_tracer_nml)
ierr = check_nml_error(io_status,'ocean_tracer_nml')
call close_file (ioun)
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) then
cmip_offset = kelvin
temp_units='degrees 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_polar_filter = mpp_clock_id('(Ocean tracer: polar filter) ' ,grain=CLOCK_MODULE)
! 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_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_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) )
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_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_expl(:) = -1
id_surf_tracer(:) = -1
id_surf_tracer_sq(:)= -1
id_bott_tracer(:) = -1
! 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 MOM4_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 MOM4_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))
#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.
enddo
if(inflow_nboundary) then
call inflow_nboundary_init
endif
! 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 ('dst_linear')
T_prog(n)%horz_advect_scheme = ADVECT_DST_LINEAR
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 ('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 default
call mpp_error(FATAL, trim(error_header) // ' Invalid horz-advect-scheme 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 ('dst_linear')
T_prog(n)%vert_advect_scheme = ADVECT_DST_LINEAR
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 default
call mpp_error(FATAL, trim(error_header) // ' Invalid vert-advect scheme 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)%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_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_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_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/))
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/))
range_array(1) = T_prog(n)%min_range
range_array(2) = T_prog(n)%max_range
if(T_prog(n)%longname=='Potential temperature') then
id_prog(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_tracer(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')
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_prog(n)%name), &
'squared '//trim(temp_units), &
missing_value=missing_value, range=(/0.0,1e10/), &
standard_name='square_of_sea_surface_temperature')
elseif(T_prog(n)%name=='salt') then
id_prog(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_tracer(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')
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_prog(n)%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(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_tracer(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)
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_prog(n)%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(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_tracer(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)
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_prog(n)%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(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_tracer(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)
id_surf_tracer_sq(n) = register_diag_field ('ocean_model', &
'squared_surface_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:2), &
Time%model_time, 'squared '//trim(T_prog(n)%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(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)
! 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(T_prog(n)%name)//'_tendency_conc', Grd%tracer_axes(1:3), &
Time%model_time, 'time tendency of tracer concentration for '//trim(T_prog(n)%longname),&
trim(T_prog(n)%units)//' per second', missing_value=missing_value, range=range_array)
id_tendency(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, range=range_array)
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)
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, range=range_array)
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, range=range_array)
else
id_tendency_conc(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(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_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)
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)
endif
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)%longname == 'Conservative temperature') prog_temp_variable=CONSERVATIVE_TEMP
if (T_prog(n)%longname == 'Potential temperature') prog_temp_variable=POTENTIAL_TEMP
enddo !} n-loop
id_neutral_rho_tendency = register_diag_field ('ocean_model','neutral_rho_tendency', &
Grd%tracer_axes(1:3), Time%model_time,'time tendency of neutral density', &
'rho*rho_dz/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
id_wdian_rho_tendency = register_diag_field ('ocean_model','wdian_rho_tendency', &
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral velocity component due to time tendency of neutral density',&
'm/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
id_neutral_rho_smooth = register_diag_field ('ocean_model','neutral_rho_smooth', &
Grd%tracer_axes(1:3), Time%model_time,'time tendency of neutral density from eta/pbot smoother',&
'rho*rho_dz/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
id_wdian_rho_smooth = register_diag_field ('ocean_model','wdian_rho_smooth', &
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral velocity component due to eta/pbot smoother', &
'm/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
id_neutral_rho_pme = register_diag_field ('ocean_model','neutral_rho_pme', &
Grd%tracer_axes(1:3), Time%model_time,'time tendency of neutral density from pme', &
'rho*rho_dz/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
id_wdian_rho_pme = register_diag_field ('ocean_model','wdian_rho_pme', &
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral velocity component due to pme', &
'm/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
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. mom4 needs both.')
endif
#ifdef USE_OCEAN_BGC
!Get the %filed for "generic" tracers as it might have already been set.
!nnz: find a way to use their already allocated field pointer directly.
do n=1,num_prog_tracers
if(T_prog(n)%type .eq. 'generic') then
call ocean_generic_get_field(T_prog(n)%name,T_prog(n)%field)
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 mom4 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
fact = Thickness%dzwt(i,j,kb-1)/Grd%dzw(kb-1)
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, cmip_units, debug) &
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) :: cmip_units
logical, intent(in), optional :: debug
!
! Return type
!
type(ocean_diag_tracer_type), dimension(:), pointer :: T_diag
integer :: kb, i, j, k, n, length, l
character(len=32) :: name, longname, units, type
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
integer :: index_diag_tracers
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
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 MOM4_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(T_diag(n)%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(T_diag(n)%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(T_diag(n)%units), &
missing_value=missing_value, range=range_array, &
standard_name='square_of_sea_surface_temperature')
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=range_array)
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)
!nnz: find a way to use their already allocated field pointer directly.
!#ifndef MOM4_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 !{
fact = Thickness%dzwt(i,j,kb-1)/Grd%dzw(kb-1)
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. mom4p1 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, &
pressure_at_depth, T_prog, T_diag)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
type(ocean_adv_vel_type), intent(in) :: Adv_vel
type(ocean_thickness_type), intent(in) :: Thickness
real, dimension(isd:,jsd:), intent(in) :: pme
real, dimension(isd:,jsd:,:,:), intent(in) :: diff_cbt
real, dimension(isd:,jsd:,:), intent(in) :: pressure_at_depth
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_diag_tracer_type), intent(inout) :: T_diag(:)
type(time_type) :: next_time
type(time_type) :: time_step
real, dimension(isd:ied,jsd:jed) :: tmp
real :: temporary
real :: total_tracer
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, &
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, Dens, 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
! 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
!chd dirty fix constaining salinity and temperature
!chd prevent SST from dropping below -2degC by resetting SST to -2degC
!chd this is crude!
!chd do k=1,nk !nk
!chd do j=jsc,jec
!chd do i=isc,iec
!chd if(T_prog(index_temp)%field(i,j,k,taup1) < -1.0) then
!chd T_prog(index_temp)%field(i,j,k,taup1) = -1.0
!chd endif
!chd enddo
!chd enddo
!chd enddo
!chd prevent SST from rising above 38 degC by resetting SST to 38degC
!chd this is crude!
!chd do k=1,nk
!chd do j=jsc,jec
!chd do i=isc,iec
!chd if(T_prog(index_temp)%field(i,j,k,taup1) > 39.0) then
!chd T_prog(index_temp)%field(i,j,k,taup1) = 39.0
!chd endif
!chd enddo
!chd enddo
!chd enddo
!chd prevent salinity from dropping below 1 PSU by resetting sal to 1
!chd this is crude!
!chd do k=1,nk
!chd do j=jsc,jec
!chd do i=isc,iec
!chd if(T_prog(index_salt)%field(i,j,k,taup1) < 1e-1) then
!chd T_prog(index_salt)%field(i,j,k,taup1) = 1e-1
!chd endif
!chd enddo
!chd enddo
!chd enddo
!chd prevent salinity from rising above 39 psu by resetting salinity to 40degC
!chd this is crude!
!chd do k=1,nk
!chd do j=jsc,jec
!chd do i=isc,iec
!chd if(T_prog(index_salt)%field(i,j,k,taup1) > 39.) then
!chd T_prog(index_salt)%field(i,j,k,taup1) = 39.
!chd endif
!chd enddo
!chd enddo
!chd enddo
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, T_prog(n)%name, 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
used = send_data(id_tendency_expl(n), T_prog(n)%th_tendency(:,:,:)*T_prog(n)%conversion, &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if (id_eta_smooth(n)> 0) then
used = send_data(id_eta_smooth(n), T_prog(n)%eta_smooth(:,:)*T_prog(n)%conversion, &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_pbot_smooth(n)> 0) then
used = send_data(id_pbot_smooth(n), T_prog(n)%pbot_smooth(:,:)*T_prog(n)%conversion, &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if(id_prog_explicit(n) > 0) then
used = send_data (id_prog_explicit(n), T_prog(n)%field(:,:,:,taup1),&
Time%model_time,rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
enddo ! 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 compute_frazil_heating(Time, Thickness, pressure_at_depth, &
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)
! 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 compute_frazil_heating(Time, Thickness, pressure_at_depth, &
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))
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
! apply polar filter to taup1 tracer concentrations
call mpp_clock_begin(id_clock_polar_filter)
call polar_filter_tracers(Time, Thickness, T_prog(1:num_prog_tracers), index_temp)
call mpp_clock_end(id_clock_polar_filter)
! 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, taup1)
endif ! endif for zero_tendency
! update the diagnostic temperature variable to taup1
if(prog_temp_variable==CONSERVATIVE_TEMP) then
T_diag(index_diag_temp)%field(:,:,:) = pottemp_from_contemp(T_prog(index_salt)%field(:,:,:,taup1), &
T_prog(index_temp)%field(:,:,:,taup1))
else
T_diag(index_diag_temp)%field(:,:,:) = contemp_from_pottemp(T_prog(index_salt)%field(:,:,:,taup1), &
T_prog(index_temp)%field(:,:,:,taup1))
endif
! send diagnostics to diag_manager
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
used = send_data (id_prog(n), wrk1(:,:,:), &
Time%model_time,rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
else
used = send_data (id_prog(n), T_prog(n)%field(:,:,:,tau),&
Time%model_time,rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
endif
if(id_prog_rhodzt(n) > 0) then
used = send_data (id_prog_rhodzt(n), Thickness%rho_dzt(:,:,:,tau)*T_prog(n)%field(:,:,:,tau), &
Time%model_time,rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
! 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
used = send_data (id_prog_int_rhodz(n), wrk1_2d(:,:), &
Time%model_time,rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
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
used = send_data (id_surf_tracer(n), wrk1_2d(:,:),&
Time%model_time,rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
else
used = send_data (id_surf_tracer(n), T_prog(n)%field(:,:,1,tau),&
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
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
used = send_data (id_surf_tracer_sq(n), wrk1_2d(:,:)**2,&
Time%model_time,rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
else
used = send_data (id_surf_tracer_sq(n), T_prog(n)%field(:,:,1,tau)**2,&
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
endif
if (id_bott_tracer(n) > 0) then
tmp(:,:) = 0.0
do j = jsd,jed
do i = isd,ied
kbot = Grd%kmt(i,j)
if(kbot > 1) then
tmp(i,j) = T_prog(n)%field(i,j,kbot,tau)
endif
enddo
enddo
used = send_data (id_bott_tracer(n), tmp(:,:), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
! total time tendency for tracer mass per horizontal area
if (id_tendency(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))
enddo
enddo
enddo
used = send_data(id_tendency(n),wrk1(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
! total time tendency for the 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
if (id_tendency_conc(n) > 0) then
used = send_data(id_tendency_conc(n),T_prog(n)%tendency(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
enddo ! enddo for n-tracers
! time tendency for neutral density and dianeutral velocity component
if (id_neutral_rho_tendency > 0 .or. id_wdian_rho_tendency > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtimer*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) = dtimer*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))
enddo
enddo
enddo
if(id_neutral_rho_tendency > 0) then
used = send_data(id_neutral_rho_tendency,wrk1(:,:,:)+wrk2(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_wdian_rho_tendency > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
temporary = Dens%drhodz_zt(i,j,k)*Thickness%rho_dzt(i,j,k,tau)
wrk3(i,j,k) = Grd%tmask(i,j,k)*(wrk1(i,j,k)+wrk2(i,j,k))/(epsln+temporary)
enddo
enddo
enddo
used = send_data(id_wdian_rho_tendency, wrk3(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
endif
! time tendency for neutral density and dianeutral velocity component due to smoothing
if (id_neutral_rho_smooth > 0 .or. id_wdian_rho_smooth > 0) then
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 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
if(id_neutral_rho_smooth > 0) then
used = send_data(id_neutral_rho_smooth,wrk1(:,:,:)+wrk2(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_wdian_rho_smooth > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
temporary = Dens%drhodz_zt(i,j,k)*Thickness%rho_dzt(i,j,k,tau)
wrk3(i,j,k) = Grd%tmask(i,j,k)*(wrk1(i,j,k)+wrk2(i,j,k))/(epsln+temporary)
enddo
enddo
enddo
used = send_data(id_wdian_rho_smooth, wrk3(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
endif
! time tendency for neutral density and dianeutral velocity component due.
! to pme. note that river contribution is diagnosed in ocean_rivermix.F90.
if (id_neutral_rho_pme > 0 .or. id_wdian_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) &
*( Dens%drhodT(i,j,k)*pme(i,j)*T_prog(index_temp)%tpme(i,j) &
+Dens%drhodS(i,j,k)*pme(i,j)*T_prog(index_salt)%tpme(i,j))
enddo
enddo
if(id_neutral_rho_pme > 0) then
used = send_data(id_neutral_rho_pme,wrk1(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_wdian_rho_pme > 0) then
k=1
do j=jsc,jec
do i=isc,iec
temporary = Dens%drhodz_zt(i,j,k)*Thickness%rho_dzt(i,j,k,tau)
wrk2(i,j,k) = Grd%tmask(i,j,k)*wrk1(i,j,k)/(epsln+temporary)
enddo
enddo
used = send_data(id_wdian_rho_pme, wrk2(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
endif
! send diagnostic tracers to diag manager
do n=1,num_diag_tracers
if (id_diag(n) > 0) used = send_data (id_diag(n), T_diag(n)%field(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_diag_surf_tracer(n) > 0) used = send_data (id_diag_surf_tracer(n), &
T_diag(n)%field(:,:,1), &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
if (id_diag_surf_tracer_sq(n) > 0) used = send_data (id_diag_surf_tracer_sq(n), &
T_diag(n)%field(:,:,1)**2, &
Time%model_time, rmask=Grd%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
if(id_diag_total(n) > 0) then
total_tracer=0.0
do k=1,nk
tmp(:,:) = T_diag(n)%conversion*Grd%tmask(:,:,k)*Grd%dat(:,:) &
*Thickness%rho_dzt(:,:,k,tau)*T_diag(n)%field(:,:,k)
if(have_obc) tmp(:,:) = tmp(:,:)*Grd%obc_tmask(:,:)
total_tracer = total_tracer + mpp_global_sum(Dom%domain2d,tmp(:,:), NON_BITWISE_EXACT_SUM)
enddo
used = send_data (id_diag_total(n), total_tracer*1e-18, Time%model_time)
endif
enddo
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, &
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) )
end if
end do
do n=1, num_tracer_restart
call save_restart(Tracer_restart(n), time_stamp)
end do
end subroutine ocean_tracer_restart
! NAME="ocean_tracer_restart"
!#######################################################################
!
!
!
! Write ocean tracer restarts
!
!
subroutine ocean_tracer_end(Time, T_prog, T_diag)
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(:)
integer :: n, tau, taup1, unit
integer(i8_kind) :: chksum
character(len=128) :: filename
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
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 in mom4. 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
! debugging and diagnostics
do n=1,num_prog_tracers
if(debug_this_module) then
write(stdoutunit,*) ' '
call write_timestamp(Time%model_time)
write(stdoutunit,*) 'tmask_limit chksum for ',trim(T_prog(n)%name),' = ',&
mpp_chksum(T_prog(n)%tmask_limit(isc:iec,jsc:jec,:))
endif
if (id_tmask_limit(n) > 0) then
used = send_data(id_tmask_limit(n), T_prog(n)%tmask_limit(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
enddo
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@noaa.gov
!
!
!
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@noaa.gov
!
!
!
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"
end module ocean_tracer_mod