src/diagnostics/MOM_diagnostics.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> Calculates any requested diagnostic quantities
6		!! that are not calculated in the various subroutines.
7		!! Diagnostic quantities are requested by allocating them memory.
8		module MOM_diagnostics
9
10		use MOM_coms, only : reproducing_sum
11		use MOM_coupler_types, only : coupler_type_send_data
12		use MOM_density_integrals, only : int_density_dz
13		use MOM_diag_mediator, only : post_data, get_diag_time_end
14		use MOM_diag_mediator, only : post_product_u, post_product_sum_u
15		use MOM_diag_mediator, only : post_product_v, post_product_sum_v
16		use MOM_diag_mediator, only : register_diag_field, register_scalar_field
17		use MOM_diag_mediator, only : register_static_field, diag_register_area_ids
18		use MOM_diag_mediator, only : diag_ctrl, time_type, safe_alloc_ptr
19		use MOM_diag_mediator, only : diag_get_volume_cell_measure_dm_id
20		use MOM_diag_mediator, only : diag_grid_storage
21		use MOM_diag_mediator, only : diag_save_grids, diag_restore_grids, diag_copy_storage_to_diag
22		use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
23		use MOM_domains, only : To_North, To_East
24		use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain
25		use MOM_EOS, only : cons_temp_to_pot_temp, pot_temp_to_cons_temp
26		use MOM_EOS, only : prac_saln_to_abs_saln, abs_saln_to_prac_saln
27		use MOM_error_handler, only : MOM_error, FATAL, WARNING
28		use MOM_file_parser, only : get_param, log_version, param_file_type
29		use MOM_grid, only : ocean_grid_type
30		use MOM_interface_heights, only : find_eta, find_dz_for_eta, find_col_mass
31		use MOM_spatial_means, only : global_area_mean, global_layer_mean
32		use MOM_spatial_means, only : global_volume_mean, global_area_integral
33		use MOM_tracer_registry, only : tracer_registry_type, post_tracer_transport_diagnostics
34		use MOM_unit_scaling, only : unit_scale_type
35		use MOM_variables, only : thermo_var_ptrs, ocean_internal_state, p3d
36		use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, surface
37		use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units, get_flux_units
38		use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init
39		use Recon1d_EPPM_CWK, only : EPPM_CWK
40
41		implicit none ; private
42
43		#include <MOM_memory.h>
44
45		public calculate_diagnostic_fields, register_time_deriv, write_static_fields
46		public register_surface_diags, post_surface_dyn_diags, post_surface_thermo_diags
47		public register_transport_diags, post_transport_diagnostics
48		public MOM_diagnostics_init, MOM_diagnostics_end
49
50		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
51		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
52		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
53		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
54
55		!> The control structure for the MOM_diagnostics module
56		type, public :: diagnostics_CS ; private
57		logical :: initialized = .false. !< True if this control structure has been initialized.
58		real :: mono_N2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as
59		!! monotonic for the purposes of calculating the equivalent
60		!! barotropic wave speed [nondim].
61		real :: mono_N2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of
62		!! calculating the equivalent barotropic wave speed [H ~> m or kg m-2].
63		logical :: accurate_thick_cello !< If true, use the same careful integrals to find the diagnosed
64		!! non-Boussinesq layer thicknesses as are used to find the free
65		!! surface height, instead of using an approximate thickness
66		!! based on division by the mid-layer density.
67
68		type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
69		!! regulate the timing of diagnostic output.
70
71		! following arrays store diagnostics calculated here and unavailable outside.
72
73		! following fields have nz layers.
74		real, allocatable :: du_dt(:,:,:) !< net i-acceleration [L T-2 ~> m s-2]
75		real, allocatable :: dv_dt(:,:,:) !< net j-acceleration [L T-2 ~> m s-2]
76		real, allocatable :: dh_dt(:,:,:) !< thickness rate of change [H T-1 ~> m s-1 or kg m-2 s-1]
77
78		logical :: KE_term_on !< If true, at least one diagnostic term in the KE budget is in use.
79
80		!>@{ Diagnostic IDs
81		integer :: id_u = -1, id_v = -1, id_h = -1
82		integer :: id_usq = -1, id_vsq = -1, id_uv = -1
83		integer :: id_e = -1, id_e_D = -1
84		integer :: id_du_dt = -1, id_dv_dt = -1
85		! integer :: id_hf_du_dt = -1, id_hf_dv_dt = -1
86		integer :: id_h_du_dt = -1, id_h_dv_dt = -1
87		integer :: id_hf_du_dt_2d = -1, id_hf_dv_dt_2d = -1
88		integer :: id_col_ht = -1, id_dh_dt = -1
89		integer :: id_KE = -1, id_dKEdt = -1
90		integer :: id_PE_to_KE = -1, id_KE_BT = -1
91		integer :: id_KE_SAL = -1, id_KE_TIDES = -1
92		integer :: id_KE_BT_PF = -1, id_KE_BT_CF = -1
93		integer :: id_KE_BT_WD = -1
94		integer :: id_PE_to_KE_btbc = -1, id_KE_Coradv_btbc = -1
95		integer :: id_KE_Coradv = -1, id_KE_adv = -1
96		integer :: id_KE_visc = -1, id_KE_stress = -1
97		integer :: id_KE_visc_gl90 = -1
98		integer :: id_KE_horvisc = -1, id_KE_dia = -1
99		integer :: id_uh_Rlay = -1, id_vh_Rlay = -1
100		integer :: id_uhGM_Rlay = -1, id_vhGM_Rlay = -1
101		integer :: id_h_Rlay = -1, id_Rd1 = -1
102		integer :: id_Rml = -1, id_Rcv = -1
103		integer :: id_cg1 = -1, id_cfl_cg1 = -1
104		integer :: id_cfl_cg1_x = -1, id_cfl_cg1_y = -1
105		integer :: id_cg_ebt = -1, id_Rd_ebt = -1
106		integer :: id_p_ebt = -1
107		integer :: id_temp_int = -1, id_salt_int = -1
108		integer :: id_absscint = -1, id_pfscint = -1
109		integer :: id_scint = -1
110		integer :: id_chcint = -1, id_phcint = -1
111		integer :: id_mass_wt = -1, id_col_mass = -1
112		integer :: id_masscello = -1, id_masso = -1
113		integer :: id_volcello = -1
114		integer :: id_Tpot = -1, id_Sprac = -1
115		integer :: id_tob = -1, id_sob = -1
116		integer :: id_thetaoga = -1, id_soga = -1
117		integer :: id_bigthetaoga = -1, id_abssoga = -1
118		integer :: id_sosga = -1, id_tosga = -1
119		integer :: id_abssosga = -1, id_bigtosga = -1
120		integer :: id_temp_layer_ave = -1, id_salt_layer_ave = -1
121		integer :: id_bigtemp_layer_ave = -1, id_abssalt_layer_ave = -1
122		integer :: id_pbo = -1
123		integer :: id_thkcello = -1, id_rhoinsitu = -1
124		integer :: id_rhopot0 = -1, id_rhopot2 = -1
125		integer :: id_drho_dT = -1, id_drho_dS = -1
126		integer :: id_h_pre_sync = -1
127		integer :: id_tosq = -1, id_sosq = -1
128		integer :: id_t20d = -1, id_t17d = -1
129
130		!>@}
131		type(wave_speed_CS) :: wave_speed !< Wave speed control struct
132
133		type(p3d) :: var_ptr(MAX_FIELDS_) !< pointers to variables used in the calculation
134		!! of time derivatives
135		type(p3d) :: deriv(MAX_FIELDS_) !< Time derivatives of various fields
136		type(p3d) :: prev_val(MAX_FIELDS_) !< Previous values of variables used in the calculation
137		!! of time derivatives
138		!< previous values of variables used in calculation of time derivatives
139		integer :: nlay(MAX_FIELDS_) !< The number of layers in each diagnostics
140		integer :: num_time_deriv = 0 !< The number of time derivative diagnostics
141
142		type(group_pass_type) :: pass_KE_uv !< A handle used for group halo passes
143
144		end type diagnostics_CS
145
146
147		!> A structure with diagnostic IDs of the surface and integrated variables
148		type, public :: surface_diag_IDs ; private
149		!>@{ Diagnostic IDs for 2-d surface and bottom flux and state fields
150		!Diagnostic IDs for 2-d surface and bottom fields
151		integer :: id_zos = -1, id_zossq = -1
152		integer :: id_volo = -1, id_speed = -1
153		integer :: id_ssh = -1, id_ssh_ga = -1
154		integer :: id_sst = -1, id_sst_sq = -1, id_sstcon = -1
155		integer :: id_sss = -1, id_sss_sq = -1, id_sssabs = -1
156		integer :: id_ssu = -1, id_ssv = -1
157		integer :: id_ssu_east = -1, id_ssv_north = -1
158
159		! Diagnostic IDs for heat and salt flux fields
160		integer :: id_fraz = -1
161		integer :: id_salt_deficit = -1
162		integer :: id_Heat_PmE = -1
163		integer :: id_intern_heat = -1
164		!>@}
165		end type surface_diag_IDs
166
167
168		!> A structure with diagnostic IDs of mass transport related diagnostics
169		type, public :: transport_diag_IDs ; private
170		!>@{ Diagnostics for tracer horizontal transport
171		integer :: id_uhtr = -1, id_umo = -1, id_umo_2d = -1
172		integer :: id_vhtr = -1, id_vmo = -1, id_vmo_2d = -1
173		integer :: id_dynamics_h = -1, id_dynamics_h_tendency = -1
174		!>@}
175		end type transport_diag_IDs
176
177
178		contains
179		!> Diagnostics not more naturally calculated elsewhere are computed here.
180	12	subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
181		dt, diag_pre_sync, G, GV, US, CS)
182		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
183		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
184		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
185		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
186		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
187		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
188		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
189		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
190		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
191		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
192		intent(in) :: uh !< Transport through zonal faces = uhdy,
193		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
194		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
195		intent(in) :: vh !< Transport through meridional faces = vhdx,
196		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
197		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
198		!! thermodynamic variables.
199		type(accel_diag_ptrs), intent(in) :: ADp !< structure with pointers to
200		!! accelerations in momentum equation.
201		type(cont_diag_ptrs), intent(in) :: CDp !< structure with pointers to
202		!! terms in continuity equation.
203		real, dimension(:,:), pointer :: p_surf !< A pointer to the surface pressure [R L2 T-2 ~> Pa].
204		!! If p_surf is not associated, it is the same
205		!! as setting the surface pressure to 0.
206		real, intent(in) :: dt !< The time difference since the last
207		!! call to this subroutine [T ~> s].
208		type(diag_grid_storage), intent(in) :: diag_pre_sync !< Target grids from previous timestep
209		type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by a
210		!! previous call to diagnostics_init.
211
212		! Local variables
213	24	real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: uv ! u x v at h-points [L2 T-2 ~> m2 s-2]
214
215		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
216		integer i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
217
218	24	real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, either relative to a reference
219		! geopotential or the seafloor [Z ~> m].
220	24	real :: Rcv(SZI_(G),SZJ_(G),SZK_(GV)) ! Coordinate variable potential density [R ~> kg m-3].
221	24	real :: work_3d(SZI_(G),SZJ_(G),SZK_(GV)) ! A 3-d temporary work array in various units
222		! including [nondim] and [H ~> m or kg m-2].
223	24	real :: dz_lay(SZI_(G),SZJ_(G),SZK_(GV)) ! Height change across layers [Z ~> m]
224	24	real :: uh_tmp(SZIB_(G),SZJ_(G),SZK_(GV)) ! A temporary zonal transport [H L2 T-1 ~> m3 s-1 or kg s-1]
225	24	real :: vh_tmp(SZI_(G),SZJB_(G),SZK_(GV)) ! A temporary meridional transport [H L2 T-1 ~> m3 s-1 or kg s-1]
226	24	real :: mass_cell(SZI_(G),SZJ_(G)) ! The vertically integrated mass in a grid cell [R Z L2 ~> kg]
227	24	real :: rho_in_situ(SZI_(G)) ! In situ density [R ~> kg m-3]
228	24	real :: cg1(SZI_(G),SZJ_(G)) ! First baroclinic gravity wave speed [L T-1 ~> m s-1]
229	24	real :: Rd1(SZI_(G),SZJ_(G)) ! First baroclinic deformation radius [L ~> m]
230	24	real :: CFL_cg1(SZI_(G),SZJ_(G)) ! CFL for first baroclinic gravity wave speed, either based on the
231		! overall grid spacing or just one direction [nondim]
232
233		! tmp array for surface properties
234	24	real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS [R L2 T-2 ~> Pa]
235		real :: wt, wt_p ! The fractional weights of two successive values when interpolating from
236		! a list [nondim], scaled so that wt + wt_p = 1.
237		real :: f2_h ! Squared Coriolis parameter at to h-points [T-2 ~> s-2]
238		real :: mag_beta ! Magnitude of the gradient of f [T-1 L-1 ~> s-1 m-1]
239		real :: absurdly_small_freq2 ! Frequency squared used to avoid division by 0 [T-2 ~> s-2]
240
241		integer :: k_list
242
243	24	real, dimension(SZK_(GV)) :: temp_layer_ave ! The average temperature in a layer [C ~> degC]
244	24	real, dimension(SZK_(GV)) :: salt_layer_ave ! The average salinity in a layer [S ~> ppt]
245		real :: thetaoga ! The volume mean potential temperature [C ~> degC]
246		real :: soga ! The volume mean ocean salinity [S ~> ppt]
247		real :: masso ! The total mass of the ocean [R Z L2 ~> kg]
248		real :: tosga ! The area mean sea surface temperature [C ~> degC]
249		real :: sosga ! The area mean sea surface salinity [S ~> ppt]
250
251	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
252	12	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
253	12	nz = GV%ke ; nkmb = GV%nk_rho_varies
254
255		! This value is roughly (pi / (the age of the universe) )^2.
256	12	absurdly_small_freq2 = 1e-34US%T_to_s*2
257
258	12	if (.not. CS%initialized) call MOM_error(FATAL, &
259	0	"calculate_diagnostic_fields: Module must be initialized before used.")
260
261	12	call calculate_derivs(dt, G, CS)
262
263	12	if (dt > 0.0) then
264	12	call diag_save_grids(CS%diag)
265	12	call diag_copy_storage_to_diag(CS%diag, diag_pre_sync)
266
267	12	if (CS%id_h_pre_sync > 0) &
268	0	call post_data(CS%id_h_pre_sync, diag_pre_sync%h_state, CS%diag, alt_h=diag_pre_sync%h_state)
269
270	12	if (CS%id_du_dt>0) call post_data(CS%id_du_dt, CS%du_dt, CS%diag, alt_h=diag_pre_sync%h_state)
271
272	12	if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h=diag_pre_sync%h_state)
273
274	12	if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h=diag_pre_sync%h_state)
275
276		!! Diagnostics for terms multiplied by fractional thicknesses
277
278		! 3D diagnostics hf_du(dv)_dt are commented because there is no clarity on proper remapping grid option.
279		! The code is retained for debugging purposes in the future.
280		!if (CS%id_hf_du_dt > 0) then
281		! call post_product_u(CS%id_hf_du_dt, CS%du_dt, ADp%diag_hfrac_u, G, nz, CS%diag, alt_h=diag_pre_sync%h_state)
282		!if (CS%id_hf_dv_dt > 0) &
283		! call post_product_v(CS%id_hf_dv_dt, CS%dv_dt, ADp%diag_hfrac_v, G, nz, CS%diag, alt_h=diag_pre_sync%h_state)
284
285	12	if (CS%id_hf_du_dt_2d > 0) &
286	0	call post_product_sum_u(CS%id_hf_du_dt_2d, CS%du_dt, ADp%diag_hfrac_u, G, nz, CS%diag)
287	12	if (CS%id_hf_dv_dt_2d > 0) &
288	0	call post_product_sum_v(CS%id_hf_dv_dt_2d, CS%dv_dt, ADp%diag_hfrac_v, G, nz, CS%diag)
289
290	12	if (CS%id_h_du_dt > 0) &
291	0	call post_product_u(CS%id_h_du_dt, CS%du_dt, ADp%diag_hu, G, nz, CS%diag)
292	12	if (CS%id_h_dv_dt > 0) &
293	0	call post_product_v(CS%id_h_dv_dt, CS%dv_dt, ADp%diag_hv, G, nz, CS%diag)
294
295	12	call diag_restore_grids(CS%diag)
296
297	12	call calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS)
298		endif
299
300		! smg: is the following robust to ALE? It seems a bit opaque.
301		! If the model is NOT in isopycnal mode then nkmb=0. But we need all the
302		! following diagnostics to treat all layers as variable density, so we set
303		! nkmb = nz, on the expectation that loops nkmb+1,nz will not iterate.
304		! This behavior is ANSI F77 but some compiler options can force at least
305		! one iteration that would break the following one-line workaround!
306	12	if (nkmb==0 .and. nz > 1) nkmb = nz
307
308	12	if (CS%id_u > 0) call post_data(CS%id_u, u, CS%diag)
309
310	12	if (CS%id_v > 0) call post_data(CS%id_v, v, CS%diag)
311
312	12	if (CS%id_h > 0) call post_data(CS%id_h, h, CS%diag)
313
314	12	if (CS%id_usq > 0) call post_product_u(CS%id_usq, u, u, G, nz, CS%diag)
315
316	12	if (CS%id_vsq > 0) call post_product_v(CS%id_vsq, v, v, G, nz, CS%diag)
317
318	12	if (CS%id_uv > 0) then
319	0	do k=1,nz ; do j=js,je ; do i=is,ie
320		uv(i,j,k) = (0.5(u(I-1,j,k) + u(I,j,k))) &
321	0	(0.5*(v(i,J-1,k) + v(i,J,k)))
322		enddo ; enddo ; enddo
323	0	call post_data(CS%id_uv, uv, CS%diag)
324		endif
325
326		! Find the layer thicknesses in [Z ~> m] that can be used to determine interface heights
327	12	if ((CS%id_e > 0) .or. (CS%id_e_D > 0) .or. &
328		((CS%id_thkcello>0 .or. CS%id_volcello>0) .and. (CS%accurate_thick_cello))) &
329	0	call find_dz_for_eta(h, tv, G, GV, US, dz_lay)
330
331	12	if ((CS%id_e > 0) .or. (CS%id_e_D > 0)) then
332		! Find the interface heights, relative a reference height or to the bottom [Z ~> m]
333	0	do j=js,je ; do i=is,ie ; eta(i,j,nz+1) = -(G%bathyT(i,j) + G%Z_ref) ; enddo ; enddo
334	0	do k=nz,1,-1 ; do j=js,je ; do i=is,ie
335	0	eta(i,j,K) = eta(i,j,K+1) + dz_lay(i,j,K)
336		enddo ; enddo ; enddo
337	0	if (CS%id_e > 0) call post_data(CS%id_e, eta, CS%diag)
338
339	0	if (CS%id_e_D > 0) then
340		! Find the interface heights, relative to the bottom [Z ~> m]
341	0	do k=1,nz+1 ; do j=js,je ; do i=is,ie
342	0	eta(i,j,k) = eta(i,j,k) + (G%bathyT(i,j) + G%Z_ref)
343		enddo ; enddo ; enddo
344		! This is more accurate but changes answers in the e_D diagnostic:
345		! do j=js,je ; do i=is,ie ; eta(i,j,nz+1) = 0.0 ; enddo ; enddo
346		! do k=nz,1,-1 ; do j=js,je ; do i=is,ie
347		! eta(i,j,K) = eta(i,j,K+1) + dz_lay(i,j,K)
348		! enddo ; enddo ; enddo
349	0	call post_data(CS%id_e_D, eta, CS%diag)
350		endif
351		endif
352
353		! mass per area of grid cell (for Boussinesq, use Rho0)
354	12	if (CS%id_masscello > 0) then
355	0	call post_data(CS%id_masscello, h, CS%diag)
356		endif
357
358		! mass of liquid ocean (for Bouss, use Rho0) [R Z L2 ~> kg]
359	12	if (CS%id_masso > 0) then
360	0	mass_cell(:,:) = 0.0
361	0	do k=1,nz ; do j=js,je ; do i=is,ie
362	0	mass_cell(i,j) = mass_cell(i,j) + (GV%H_to_RZh(i,j,k)) G%areaT(i,j)
363		enddo ; enddo ; enddo
364	0	masso = reproducing_sum(mass_cell, unscale=US%RZL2_to_kg)
365	0	call post_data(CS%id_masso, masso, CS%diag)
366		endif
367
368		! diagnose thickness/volumes of grid cells [Z ~> m] and [m3]
369	12	if (CS%id_thkcello>0 .or. CS%id_volcello>0) then
370	0	if (GV%Boussinesq) then ! thkcello = h for Boussinesq
371	0	if (CS%id_thkcello > 0) then ; if (GV%H_to_Z == 1.0) then
372	0	call post_data(CS%id_thkcello, h, CS%diag)
373		else
374	0	do k=1,nz ; do j=js,je ; do i=is,ie
375	0	dz_lay(i,j,k) = GV%H_to_Z*h(i,j,k)
376		enddo ; enddo ; enddo
377	0	call post_data(CS%id_thkcello, dz_lay, CS%diag)
378		endif ; endif
379	0	if (CS%id_volcello > 0) then ! volcello = h*area for Boussinesq
380	0	do k=1,nz ; do j=js,je ; do i=is,ie
381	0	work_3d(i,j,k) = ( GV%H_to_Zh(i,j,k) ) US%Z_to_mUS%L_to_m2G%areaT(i,j)
382		enddo ; enddo ; enddo
383	0	call post_data(CS%id_volcello, work_3d, CS%diag)
384		endif
385		else ! thkcello is approximately dp/(rho*g) in non-Boussinesq mode.
386	0	if (.not.CS%accurate_thick_cello) then
387		! This is only an approximate calculation of dz_lay that does not use the careful integrals
388		! found in find_dz_for_eta that mirror what is done for the pressure gradient calculations.
389	0	EOSdom(:) = EOS_domain(G%HI)
390	0	do j=js,je
391	0	if (associated(p_surf)) then ! Pressure loading at top of surface layer [R L2 T-2 ~> Pa]
392	0	do i=is,ie
393	0	pressure_1d(i) = p_surf(i,j)
394		enddo
395		else
396	0	do i=is,ie
397	0	pressure_1d(i) = 0.0
398		enddo
399		endif
400	0	do k=1,nz ! Integrate vertically downward for pressure
401	0	do i=is,ie ! Pressure for EOS at the layer center [R L2 T-2 ~> Pa]
402	0	pressure_1d(i) = pressure_1d(i) + 0.5(GV%H_to_RZGV%g_Earth)*h(i,j,k)
403		enddo
404		! Store in-situ density [R ~> kg m-3] in work_3d
405		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, rho_in_situ, &
406	0	tv%eqn_of_state, EOSdom)
407	0	do i=is,ie ! Cell thickness = dz = dp/(g*rho) (meter); store in work_3d
408	0	dz_lay(i,j,k) = (GV%H_to_RZ*h(i,j,k)) / rho_in_situ(i)
409		enddo
410	0	do i=is,ie ! Pressure for EOS at the bottom interface [R L2 T-2 ~> Pa]
411	0	pressure_1d(i) = pressure_1d(i) + 0.5(GV%H_to_RZGV%g_Earth)*h(i,j,k)
412		enddo
413		enddo ! k
414		enddo ! j
415		endif ! Otherwise dz_lay is set in the call to find_dz_for_eta above.
416	0	if (CS%id_thkcello > 0) call post_data(CS%id_thkcello, dz_lay, CS%diag)
417	0	if (CS%id_volcello > 0) then
418	0	do k=1,nz ; do j=js,je ; do i=is,ie ! volcello = dp/(rhog)area for non-Boussinesq
419	0	work_3d(i,j,k) = US%Z_to_mUS%L_to_m2G%areaT(i,j) * dz_lay(i,j,k)
420		enddo ; enddo ; enddo
421	0	call post_data(CS%id_volcello, work_3d, CS%diag)
422		endif
423		endif
424		endif
425
426		! Calculate additional, potentially derived temperature diagnostics
427	12	if (tv%T_is_conT) then
428		! Internal T&S variables are conservative temperature & absolute salinity,
429		! so they need to converted to potential temperature and practical salinity
430		! for some diagnostics using TEOS-10 function calls.
431	0	if ((CS%id_Tpot > 0) .or. (CS%id_tob > 0) .or. (CS%id_tosq > 0)) then
432	0	EOSdom(:) = EOS_domain(G%HI)
433	0	do k=1,nz ; do j=js,je
434	0	call cons_temp_to_pot_temp(tv%T(:,j,k), tv%S(:,j,k), work_3d(:,j,k), tv%eqn_of_state, EOSdom)
435		enddo ; enddo
436	0	if (CS%id_Tpot > 0) call post_data(CS%id_Tpot, work_3d, CS%diag)
437	0	if (CS%id_tob > 0) call post_data(CS%id_tob, work_3d(:,:,nz), CS%diag)
438		! volume mean potential temperature
439	0	if (CS%id_thetaoga>0) then
440	0	thetaoga = global_volume_mean(work_3d, h, G, GV, tmp_scale=US%C_to_degC)
441	0	call post_data(CS%id_thetaoga, thetaoga, CS%diag)
442		endif
443		! volume mean conservative temperature
444	0	if (CS%id_bigthetaoga>0) then
445	0	thetaoga = global_volume_mean(tv%T, h, G, GV, tmp_scale=US%C_to_degC)
446	0	call post_data(CS%id_bigthetaoga, thetaoga, CS%diag)
447		endif
448		! area mean potential SST
449	0	if (CS%id_tosga > 0) then
450	0	tosga = global_area_mean(work_3d(:,:,1), G, tmp_scale=US%C_to_degC)
451	0	call post_data(CS%id_tosga, tosga, CS%diag)
452		endif
453		! area mean conservative SST
454	0	if (CS%id_bigtosga > 0) then
455	0	tosga = global_area_mean(tv%T(:,:,1), G, tmp_scale=US%C_to_degC)
456	0	call post_data(CS%id_bigtosga, tosga, CS%diag)
457		endif
458		! layer mean potential temperature
459	0	if (CS%id_temp_layer_ave>0) then
460	0	temp_layer_ave = global_layer_mean(work_3d, h, G, GV, tmp_scale=US%C_to_degC)
461	0	call post_data(CS%id_temp_layer_ave, temp_layer_ave, CS%diag)
462		endif
463		! layer mean conservative temperature
464	0	if (CS%id_bigtemp_layer_ave>0) then
465	0	temp_layer_ave = global_layer_mean(tv%T, h, G, GV, tmp_scale=US%C_to_degC)
466	0	call post_data(CS%id_bigtemp_layer_ave, temp_layer_ave, CS%diag)
467		endif
468	0	if (CS%id_tosq > 0) then
469	0	do k=1,nz ; do j=js,je ; do i=is,ie
470	0	work_3d(i,j,k) = work_3d(i,j,k)*work_3d(i,j,k)
471		enddo ; enddo ; enddo
472	0	call post_data(CS%id_tosq, work_3d, CS%diag)
473		endif
474		endif
475		else
476		! Internal T&S variables are potential temperature & practical salinity
477	12	if (CS%id_tob > 0) call post_data(CS%id_tob, tv%T(:,:,nz), CS%diag)
478	12	if (CS%id_tosq > 0) then
479	0	do k=1,nz ; do j=js,je ; do i=is,ie
480	0	work_3d(i,j,k) = tv%T(i,j,k)*tv%T(i,j,k)
481		enddo ; enddo ; enddo
482	0	call post_data(CS%id_tosq, work_3d, CS%diag)
483		endif
484		! volume mean potential temperature
485	12	if (CS%id_thetaoga>0) then
486	0	thetaoga = global_volume_mean(tv%T, h, G, GV, tmp_scale=US%C_to_degC)
487	0	call post_data(CS%id_thetaoga, thetaoga, CS%diag)
488		endif
489		! area mean SST
490	12	if (CS%id_tosga > 0) then
491	0	tosga = global_area_mean(tv%T(:,:,1), G, tmp_scale=US%C_to_degC)
492	0	call post_data(CS%id_tosga, tosga, CS%diag)
493		endif
494		! layer mean potential temperature
495	12	if (CS%id_temp_layer_ave>0) then
496	0	temp_layer_ave = global_layer_mean(tv%T, h, G, GV, tmp_scale=US%C_to_degC)
497	0	call post_data(CS%id_temp_layer_ave, temp_layer_ave, CS%diag)
498		endif
499		endif
500
501
502		! Calculate additional, potentially derived salinity diagnostics
503	12	if (tv%S_is_absS) then
504		! Internal T&S variables are conservative temperature & absolute salinity,
505		! so they need to converted to potential temperature and practical salinity
506		! for some diagnostics using TEOS-10 function calls.
507	0	if ((CS%id_Sprac > 0) .or. (CS%id_sob > 0) .or. (CS%id_sosq >0)) then
508	0	EOSdom(:) = EOS_domain(G%HI)
509	0	do k=1,nz ; do j=js,je
510	0	call abs_saln_to_prac_saln(tv%S(:,j,k), work_3d(:,j,k), tv%eqn_of_state, EOSdom)
511		enddo ; enddo
512	0	if (CS%id_Sprac > 0) call post_data(CS%id_Sprac, work_3d, CS%diag)
513	0	if (CS%id_sob > 0) call post_data(CS%id_sob, work_3d(:,:,nz), CS%diag)
514		! volume mean salinity
515	0	if (CS%id_soga>0) then
516	0	soga = global_volume_mean(work_3d, h, G, GV, tmp_scale=US%S_to_ppt)
517	0	call post_data(CS%id_soga, soga, CS%diag)
518		endif
519		! volume mean absolute salinity
520	0	if (CS%id_abssoga>0) then
521	0	soga = global_volume_mean(tv%S, h, G, GV, tmp_scale=US%S_to_ppt)
522	0	call post_data(CS%id_abssoga, soga, CS%diag)
523		endif
524		! area mean practical SSS
525	0	if (CS%id_sosga > 0) then
526	0	sosga = global_area_mean(work_3d(:,:,1), G, tmp_scale=US%S_to_ppt)
527	0	call post_data(CS%id_sosga, sosga, CS%diag)
528		endif
529		! area mean absolute SSS
530	0	if (CS%id_abssosga > 0) then
531	0	sosga = global_area_mean(tv%S(:,:,1), G, tmp_scale=US%S_to_ppt)
532	0	call post_data(CS%id_abssosga, sosga, CS%diag)
533		endif
534		! layer mean practical salinity
535	0	if (CS%id_salt_layer_ave>0) then
536	0	salt_layer_ave = global_layer_mean(work_3d, h, G, GV, tmp_scale=US%S_to_ppt)
537	0	call post_data(CS%id_salt_layer_ave, salt_layer_ave, CS%diag)
538		endif
539		! layer mean absolute salinity
540	0	if (CS%id_abssalt_layer_ave>0) then
541	0	salt_layer_ave = global_layer_mean(tv%S, h, G, GV, tmp_scale=US%S_to_ppt)
542	0	call post_data(CS%id_abssalt_layer_ave, salt_layer_ave, CS%diag)
543		endif
544	0	if (CS%id_sosq > 0) then
545	0	do k=1,nz ; do j=js,je ; do i=is,ie
546	0	work_3d(i,j,k) = work_3d(i,j,k)*work_3d(i,j,k)
547		enddo ; enddo ; enddo
548	0	call post_data(CS%id_sosq, work_3d, CS%diag)
549		endif
550		endif
551		else
552		! Internal T&S variables are potential temperature & practical salinity
553	12	if (CS%id_sob > 0) call post_data(CS%id_sob, tv%S(:,:,nz), CS%diag)
554	12	if (CS%id_sosq > 0) then
555	0	do k=1,nz ; do j=js,je ; do i=is,ie
556	0	work_3d(i,j,k) = tv%S(i,j,k)*tv%S(i,j,k)
557		enddo ; enddo ; enddo
558	0	call post_data(CS%id_sosq, work_3d, CS%diag)
559		endif
560		! volume mean salinity
561	12	if (CS%id_soga>0) then
562	0	soga = global_volume_mean(tv%S, h, G, GV, tmp_scale=US%S_to_ppt)
563	0	call post_data(CS%id_soga, soga, CS%diag)
564		endif
565		! area mean SSS
566	12	if (CS%id_sosga > 0) then
567	0	sosga = global_area_mean(tv%S(:,:,1), G, tmp_scale=US%S_to_ppt)
568	0	call post_data(CS%id_sosga, sosga, CS%diag)
569		endif
570		! layer mean salinity
571	12	if (CS%id_salt_layer_ave>0) then
572	0	salt_layer_ave = global_layer_mean(tv%S, h, G, GV, tmp_scale=US%S_to_ppt)
573	0	call post_data(CS%id_salt_layer_ave, salt_layer_ave, CS%diag)
574		endif
575		endif
576
577	12	call calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
578
579		if ((CS%id_Rml > 0) .or. (CS%id_Rcv > 0) .or. (CS%id_h_Rlay > 0) .or. &
580		(CS%id_uh_Rlay > 0) .or. (CS%id_vh_Rlay > 0) .or. &
581	12	(CS%id_uhGM_Rlay > 0) .or. (CS%id_vhGM_Rlay > 0)) then
582
583	0	if (associated(tv%eqn_of_state)) then
584	0	EOSdom(:) = EOS_domain(G%HI, halo=1)
585	0	pressure_1d(:) = tv%P_Ref
586		!$OMP parallel do default(shared)
587	0	do k=1,nz ; do j=js-1,je+1
588		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), tv%eqn_of_state, &
589	0	EOSdom)
590		enddo ; enddo
591		else ! Rcv should not be used much in this case, so fill in sensible values.
592	0	do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
593	0	Rcv(i,j,k) = GV%Rlay(k)
594		enddo ; enddo ; enddo
595		endif
596	0	if (CS%id_Rml > 0) call post_data(CS%id_Rml, Rcv, CS%diag)
597	0	if (CS%id_Rcv > 0) call post_data(CS%id_Rcv, Rcv, CS%diag)
598
599	0	if (CS%id_h_Rlay > 0) then
600		! Here work_3d is used for the layer thicknesses in potential density coordinates [H ~> m or kg m-2].
601	0	k_list = nz/2
602		!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
603	0	do j=js,je
604	0	do k=1,nkmb ; do i=is,ie
605	0	work_3d(i,j,k) = 0.0
606		enddo ; enddo
607	0	do k=nkmb+1,nz ; do i=is,ie
608	0	work_3d(i,j,k) = h(i,j,k)
609		enddo ; enddo
610	0	do k=1,nkmb ; do i=is,ie
611	0	call find_weights(GV%Rlay, Rcv(i,j,k), k_list, nz, wt, wt_p)
612	0	work_3d(i,j,k_list) = work_3d(i,j,k_list) + h(i,j,k)*wt
613	0	work_3d(i,j,k_list+1) = work_3d(i,j,k_list+1) + h(i,j,k)*wt_p
614		enddo ; enddo
615		enddo
616
617	0	call post_data(CS%id_h_Rlay, work_3d, CS%diag)
618		endif
619
620	0	if (CS%id_uh_Rlay > 0) then
621		! Calculate zonal transports in potential density coordinates [H L2 T-1 ~> m3 s-1 or kg s-1].
622	0	k_list = nz/2
623		!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
624	0	do j=js,je
625	0	do k=1,nkmb ; do I=Isq,Ieq
626	0	uh_tmp(I,j,k) = 0.0
627		enddo ; enddo
628	0	do k=nkmb+1,nz ; do I=Isq,Ieq
629	0	uh_tmp(I,j,k) = uh(I,j,k)
630		enddo ; enddo
631	0	k_list = nz/2
632	0	do k=1,nkmb ; do I=Isq,Ieq
633	0	call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
634	0	uh_tmp(I,j,k_list) = uh_tmp(I,j,k_list) + uh(I,j,k)*wt
635	0	uh_tmp(I,j,k_list+1) = uh_tmp(I,j,k_list+1) + uh(I,j,k)*wt_p
636		enddo ; enddo
637		enddo
638
639	0	call post_data(CS%id_uh_Rlay, uh_tmp, CS%diag)
640		endif
641
642	0	if (CS%id_vh_Rlay > 0) then
643		! Calculate meridional transports in potential density coordinates [H L2 T-1 ~> m3 s-1 or kg s-1].
644	0	k_list = nz/2
645		!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
646	0	do J=Jsq,Jeq
647	0	do k=1,nkmb ; do i=is,ie
648	0	vh_tmp(i,J,k) = 0.0
649		enddo ; enddo
650	0	do k=nkmb+1,nz ; do i=is,ie
651	0	vh_tmp(i,J,k) = vh(i,J,k)
652		enddo ; enddo
653	0	do k=1,nkmb ; do i=is,ie
654	0	call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
655	0	vh_tmp(i,J,k_list) = vh_tmp(i,J,k_list) + vh(i,J,k)*wt
656	0	vh_tmp(i,J,k_list+1) = vh_tmp(i,J,k_list+1) + vh(i,J,k)*wt_p
657		enddo ; enddo
658		enddo
659
660	0	call post_data(CS%id_vh_Rlay, vh_tmp, CS%diag)
661		endif
662
663	0	if ((CS%id_uhGM_Rlay > 0) .and. associated(CDp%uhGM)) then
664		! Calculate zonal Gent-McWilliams transports in potential density
665		! coordinates [H L2 T-1 ~> m3 s-1 or kg s-1].
666	0	k_list = nz/2
667		!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
668	0	do j=js,je
669	0	do k=1,nkmb ; do I=Isq,Ieq
670	0	uh_tmp(I,j,k) = 0.0
671		enddo ; enddo
672	0	do k=nkmb+1,nz ; do I=Isq,Ieq
673	0	uh_tmp(I,j,k) = CDp%uhGM(I,j,k)
674		enddo ; enddo
675	0	do k=1,nkmb ; do I=Isq,Ieq
676	0	call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
677	0	uh_tmp(I,j,k_list) = uh_tmp(I,j,k_list) + CDp%uhGM(I,j,k)*wt
678	0	uh_tmp(I,j,k_list+1) = uh_tmp(I,j,k_list+1) + CDp%uhGM(I,j,k)*wt_p
679		enddo ; enddo
680		enddo
681
682	0	if (CS%id_uhGM_Rlay > 0) call post_data(CS%id_uhGM_Rlay, uh_tmp, CS%diag)
683		endif
684
685	0	if ((CS%id_vhGM_Rlay > 0) .and. associated(CDp%vhGM)) then
686		! Calculate meridional Gent-McWilliams transports in potential density
687		! coordinates [H L2 T-1 ~> m3 s-1 or kg s-1].
688	0	k_list = nz/2
689		!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
690	0	do J=Jsq,Jeq
691	0	do k=1,nkmb ; do i=is,ie
692	0	vh_tmp(i,J,k) = 0.0
693		enddo ; enddo
694	0	do k=nkmb+1,nz ; do i=is,ie
695	0	vh_tmp(i,J,k) = CDp%vhGM(i,J,k)
696		enddo ; enddo
697	0	do k=1,nkmb ; do i=is,ie
698	0	call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
699	0	vh_tmp(i,J,k_list) = vh_tmp(i,J,k_list) + CDp%vhGM(i,J,k)*wt
700	0	vh_tmp(i,J,k_list+1) = vh_tmp(i,J,k_list+1) + CDp%vhGM(i,J,k)*wt_p
701		enddo ; enddo
702		enddo
703
704	0	if (CS%id_vhGM_Rlay > 0) call post_data(CS%id_vhGM_Rlay, vh_tmp, CS%diag)
705		endif
706		endif
707
708	12	if (associated(tv%eqn_of_state)) then
709	12	EOSdom(:) = EOS_domain(G%HI)
710	12	if (CS%id_rhopot0 > 0) then
711	0	pressure_1d(:) = 0.
712		!$OMP parallel do default(shared)
713	0	do k=1,nz ; do j=js,je
714		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), &
715	0	tv%eqn_of_state, EOSdom)
716		enddo ; enddo
717	0	if (CS%id_rhopot0 > 0) call post_data(CS%id_rhopot0, Rcv, CS%diag)
718		endif
719	12	if (CS%id_rhopot2 > 0) then
720	0	pressure_1d(:) = 2.0e7*US%Pa_to_RL2_T2 ! 2000 dbars
721		!$OMP parallel do default(shared)
722	0	do k=1,nz ; do j=js,je
723		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), &
724	0	tv%eqn_of_state, EOSdom)
725		enddo ; enddo
726	0	if (CS%id_rhopot2 > 0) call post_data(CS%id_rhopot2, Rcv, CS%diag)
727		endif
728	12	if (CS%id_rhoinsitu > 0) then
729		!$OMP parallel do default(shared) private(pressure_1d)
730	0	do j=js,je
731	0	pressure_1d(:) = 0. ! Start at p=0 Pa at surface
732	0	do k=1,nz
733	0	pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * (GV%H_to_RZ*GV%g_Earth) ! Pressure in middle of layer k
734		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), &
735	0	tv%eqn_of_state, EOSdom)
736	0	pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * (GV%H_to_RZ*GV%g_Earth) ! Pressure at bottom of layer k
737		enddo
738		enddo
739	0	if (CS%id_rhoinsitu > 0) call post_data(CS%id_rhoinsitu, Rcv, CS%diag)
740		endif
741
742	12	if (CS%id_drho_dT > 0 .or. CS%id_drho_dS > 0) then
743		!$OMP parallel do default(shared) private(pressure_1d)
744	0	do j=js,je
745	0	pressure_1d(:) = 0. ! Start at p=0 Pa at surface
746	0	do k=1,nz
747	0	pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * (GV%H_to_RZ*GV%g_Earth) ! Pressure in middle of layer k
748		! To avoid storing more arrays, put drho_dT into Rcv, and drho_dS into work3d
749		call calculate_density_derivs(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, &
750	0	Rcv(:,j,k), work_3d(:,j,k), tv%eqn_of_state, EOSdom)
751	0	pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * (GV%H_to_RZ*GV%g_Earth) ! Pressure at bottom of layer k
752		enddo
753		enddo
754	0	if (CS%id_drho_dT > 0) call post_data(CS%id_drho_dT, Rcv, CS%diag)
755	0	if (CS%id_drho_dS > 0) call post_data(CS%id_drho_dS, work_3d, CS%diag)
756		endif
757		endif
758
759		if ((CS%id_cg1>0) .or. (CS%id_Rd1>0) .or. (CS%id_cfl_cg1>0) .or. &
760	12	(CS%id_cfl_cg1_x>0) .or. (CS%id_cfl_cg1_y>0)) then
761	12	call wave_speed(h, tv, G, GV, US, cg1, CS%wave_speed)
762	12	if (CS%id_cg1>0) call post_data(CS%id_cg1, cg1, CS%diag)
763	12	if (CS%id_Rd1>0) then
764		!$OMP parallel do default(shared) private(f2_h,mag_beta)
765	87132	do j=js,je ; do i=is,ie
766		! Blend the equatorial deformation radius with the standard one.
767		f2_h = absurdly_small_freq2 + 0.25 * &
768		((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
769	86400	(G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1)))
770		mag_beta = sqrt(0.5 * ( &
771		((((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2) + &
772		(((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2)) + &
773		((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2) + &
774	86400	(((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2)) ))
775	87120	Rd1(i,j) = cg1(i,j) / sqrt(f2_h + cg1(i,j) * mag_beta)
776
777		enddo ; enddo
778	12	call post_data(CS%id_Rd1, Rd1, CS%diag)
779		endif
780	12	if (CS%id_cfl_cg1>0) then
781	0	do j=js,je ; do i=is,ie
782	0	CFL_cg1(i,j) = (dtcg1(i,j)) (G%IdxT(i,j) + G%IdyT(i,j))
783		enddo ; enddo
784	0	call post_data(CS%id_cfl_cg1, CFL_cg1, CS%diag)
785		endif
786	12	if (CS%id_cfl_cg1_x>0) then
787	0	do j=js,je ; do i=is,ie
788	0	CFL_cg1(i,j) = (dtcg1(i,j)) G%IdxT(i,j)
789		enddo ; enddo
790	0	call post_data(CS%id_cfl_cg1_x, CFL_cg1, CS%diag)
791		endif
792	12	if (CS%id_cfl_cg1_y>0) then
793	0	do j=js,je ; do i=is,ie
794	0	CFL_cg1(i,j) = (dtcg1(i,j)) G%IdyT(i,j)
795		enddo ; enddo
796	0	call post_data(CS%id_cfl_cg1_y, CFL_cg1, CS%diag)
797		endif
798		endif
799	12	if ((CS%id_cg_ebt>0) .or. (CS%id_Rd_ebt>0) .or. (CS%id_p_ebt>0)) then
800	0	if (CS%id_p_ebt>0) then
801		! Here work_3d is used for the equivalent barotropic modal structure [nondim].
802	0	work_3d(:,:,:) = 0.0
803		call wave_speed(h, tv, G, GV, US, cg1, CS%wave_speed, use_ebt_mode=.true., &
804		mono_N2_column_fraction=CS%mono_N2_column_fraction, &
805	0	mono_N2_depth=CS%mono_N2_depth, modal_structure=work_3d)
806	0	call post_data(CS%id_p_ebt, work_3d, CS%diag)
807		else
808		call wave_speed(h, tv, G, GV, US, cg1, CS%wave_speed, use_ebt_mode=.true., &
809		mono_N2_column_fraction=CS%mono_N2_column_fraction, &
810	0	mono_N2_depth=CS%mono_N2_depth)
811		endif
812	0	if (CS%id_cg_ebt>0) call post_data(CS%id_cg_ebt, cg1, CS%diag)
813	0	if (CS%id_Rd_ebt>0) then
814		!$OMP parallel do default(shared) private(f2_h,mag_beta)
815	0	do j=js,je ; do i=is,ie
816		! Blend the equatorial deformation radius with the standard one.
817		f2_h = absurdly_small_freq2 + 0.25 * &
818		((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
819	0	(G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1)))
820		mag_beta = sqrt(0.5 * ( &
821		((((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2) + &
822		(((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2)) + &
823		((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2) + &
824	0	(((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2)) ))
825	0	Rd1(i,j) = cg1(i,j) / sqrt(f2_h + cg1(i,j) * mag_beta)
826
827		enddo ; enddo
828	0	call post_data(CS%id_Rd_ebt, Rd1, CS%diag)
829		endif
830		endif
831
832	12	end subroutine calculate_diagnostic_fields
833
834		!> This subroutine finds the location of R_in in an increasing ordered
835		!! list, Rlist, returning as k the element such that
836		!! Rlist(k) <= R_in < Rlist(k+1), and where wt and wt_p are the linear
837		!! weights that should be assigned to elements k and k+1.
838	0	subroutine find_weights(Rlist, R_in, k, nz, wt, wt_p)
839		real, dimension(:), &
840		intent(in) :: Rlist !< The list of target densities [R ~> kg m-3]
841		real, intent(in) :: R_in !< The density being inserted into Rlist [R ~> kg m-3]
842		integer, intent(inout) :: k !< The value of k such that Rlist(k) <= R_in < Rlist(k+1)
843		!! The input value is a first guess
844		integer, intent(in) :: nz !< The number of layers in Rlist
845		real, intent(out) :: wt !< The weight of layer k for interpolation [nondim]
846		real, intent(out) :: wt_p !< The weight of layer k+1 for interpolation [nondim]
847
848		! This subroutine finds location of R_in in an increasing ordered
849		! list, Rlist, returning as k the element such that
850		! Rlist(k) <= R_in < Rlist(k+1), and where wt and wt_p are the linear
851		! weights that should be assigned to elements k and k+1.
852
853		integer :: k_upper, k_lower, k_new, inc
854
855		! First, bracket the desired point.
856	0	if ((k < 1) .or. (k > nz)) k = nz/2
857
858	0	k_upper = k ; k_lower = k ; inc = 1
859	0	if (R_in < Rlist(k)) then
860	0	do
861	0	k_lower = max(k_lower-inc, 1)
862	0	if ((k_lower == 1) .or. (R_in >= Rlist(k_lower))) exit
863	0	k_upper = k_lower
864	0	inc = inc*2
865		enddo
866		else
867	0	do
868	0	k_upper = min(k_upper+inc, nz)
869	0	if ((k_upper == nz) .or. (R_in < Rlist(k_upper))) exit
870	0	k_lower = k_upper
871	0	inc = inc*2
872		enddo
873		endif
874
875	0	if ((k_lower == 1) .and. (R_in <= Rlist(k_lower))) then
876	0	k = 1 ; wt = 1.0 ; wt_p = 0.0
877	0	elseif ((k_upper == nz) .and. (R_in >= Rlist(k_upper))) then
878	0	k = nz-1 ; wt = 0.0 ; wt_p = 1.0
879		else
880	0	do
881	0	if (k_upper <= k_lower+1) exit
882	0	k_new = (k_upper + k_lower) / 2
883	0	if (R_in < Rlist(k_new)) then
884	0	k_upper = k_new
885		else
886	0	k_lower = k_new
887		endif
888		enddo
889
890		! Uncomment this as a code check
891		! if ((R_in < Rlist(k_lower)) .or. (R_in >= Rlist(k_upper)) .or. (k_upper-k_lower /= 1)) &
892		! write (,) "Error: ",R_in," is not between R(",k_lower,") = ", &
893		! Rlist(k_lower)," and R(",k_upper,") = ",Rlist(k_upper),"."
894	0	k = k_lower
895	0	wt = (Rlist(k_upper) - R_in) / (Rlist(k_upper) - Rlist(k_lower))
896	0	wt_p = 1.0 - wt
897
898		endif
899
900	0	end subroutine find_weights
901
902		!> This subroutine calculates vertical integrals of several tracers, along
903		!! with the mass-weight of these tracers, the total column mass, and the
904		!! carefully calculated column height.
905	12	subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
906		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
907		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
908		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
909		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
910		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
911		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
912		!! thermodynamic variables.
913		real, dimension(:,:), pointer :: p_surf !< A pointer to the surface pressure [R L2 T-2 ~> Pa].
914		!! If p_surf is not associated, it is the same
915		!! as setting the surface pressure to 0.
916		type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by a
917		!! previous call to diagnostics_init.
918		! Local variables
919		real, dimension(SZI_(G),SZJ_(G)) :: &
920	24	z_top, & ! Height of the top of a layer or the ocean [Z ~> m].
921	24	z_bot, & ! Height of the bottom of a layer (for id_mass) or the
922		! (positive) depth of the ocean (for id_col_ht) [Z ~> m].
923	24	mass, & ! integrated mass of the water column [R Z ~> kg m-2]. For
924		! non-Boussinesq models this is rho*dz. For Boussinesq
925		! models, this is either the integral of in-situ density
926		! (rhodz for col_mass) or reference density (Rho_0dz for mass_wt).
927	24	btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'.
928		! This is the column mass multiplied by gravity plus the pressure
929		! at the ocean surface [R L2 T-2 ~> Pa].
930	24	tr_int,& ! vertical integral of a tracer times density,
931		! (Rho_0 in a Boussinesq model) [Conc R Z ~> Conc kg m-2].
932	24	d17,& ! Depth of 17 degC isotherm [Z ~> m]
933	24	d20 ! Depth of 20 degC isotherm [Z ~> m]
934	24	real :: tmp(SZI_(G),SZJ_(G),SZK_(GV)) ! Temporary array [defined at each usage]
935		real :: IG_Earth ! Inverse of gravitational acceleration [T2 Z L-2 ~> s2 m-1].
936		real :: Ttop, Tbot ! Temperature at top/bottom of cell [C ~> degC]
937	12	type(EPPM_CWK) :: PPM ! Class for reconstruction
938	24	real :: d_from_ssh(0:GV%ke) ! eta-z (Distance from surface) [Z ~> m]
939		real :: dz ! Layer thickness in Z [Z ~> m]
940
941		integer :: i, j, k, is, ie, js, je, nz
942		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
943	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
944
945	12	if (CS%id_mass_wt > 0) then
946	0	do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
947	0	do k=1,nz ; do j=js,je ; do i=is,ie
948	0	mass(i,j) = mass(i,j) + GV%H_to_RZ*h(i,j,k)
949		enddo ; enddo ; enddo
950	0	call post_data(CS%id_mass_wt, mass, CS%diag)
951		endif
952
953	12	if (CS%id_temp_int > 0) then
954	0	do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
955	0	do k=1,nz ; do j=js,je ; do i=is,ie
956	0	tr_int(i,j) = tr_int(i,j) + (GV%H_to_RZh(i,j,k))tv%T(i,j,k)
957		enddo ; enddo ; enddo
958	0	call post_data(CS%id_temp_int, tr_int, CS%diag)
959		endif
960
961	12	if (CS%id_salt_int > 0) then
962	0	do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
963	0	do k=1,nz ; do j=js,je ; do i=is,ie
964	0	tr_int(i,j) = tr_int(i,j) + (GV%H_to_RZh(i,j,k))tv%S(i,j,k)
965		enddo ; enddo ; enddo
966	0	call post_data(CS%id_salt_int, tr_int, CS%diag)
967		endif
968
969	12	if (CS%id_col_ht > 0) then
970		!$omp target update to(h)
971		!$omp target enter data map(alloc: z_top)
972	0	call find_eta(h, tv, G, GV, US, z_top)
973		!$omp target exit data map(from: z_top)
974	0	do j=js,je ; do i=is,ie
975	0	z_bot(i,j) = z_top(i,j) + G%bathyT(i,j)
976		enddo ; enddo
977	0	call post_data(CS%id_col_ht, z_bot, CS%diag)
978		endif
979
980	12	if (CS%id_col_mass > 0 .or. CS%id_pbo > 0) then
981	0	if (CS%id_pbo > 0) then
982	0	call find_col_mass(h, tv, G, GV, US, mass, btm_pres, p_surf)
983	0	call post_data(CS%id_pbo, btm_pres, CS%diag)
984		else
985	0	call find_col_mass(h, tv, G, GV, US, mass)
986		endif
987	0	if (CS%id_col_mass > 0) call post_data(CS%id_col_mass, mass, CS%diag)
988		endif
989	12	if (CS%id_t20d > 0 .or. CS%id_t17d > 0) then
990	0	call PPM%init(GV%ke, h_neglect=0.)
991	0	do j=js,je ; do i=is,ie
992		! Pre-calculate the interface depths relative to the surface
993	0	if (GV%Boussinesq) then
994	0	d_from_ssh(0) = 0.
995	0	do k=1,nz
996	0	d_from_ssh(k) = d_from_ssh(k-1) + h(i,j,k) * GV%H_to_Z
997		enddo
998		else
999		! Non-Boussinesq: use pre-computed layer-average specific volumes from tv%SpV_avg,
1000		! which are more accurate than cell-center specific volumes and correctly account
1001		! for surface pressure (including under ice-shelves).
1002	0	d_from_ssh(0) = 0.
1003	0	do k=1,nz
1004	0	d_from_ssh(k) = d_from_ssh(k-1) + ( h(i,j,k) * GV%H_to_RZ ) * tv%SpV_avg(i,j,k)
1005		enddo
1006		endif
1007	0	call PPM%reconstruct(h(i,j,:), tv%T(i,j,:))
1008	0	d17(i,j) = d_from_ssh(nz)
1009	0	d20(i,j) = d_from_ssh(nz)
1010	0	do k=nz,1,-1
1011	0	Ttop = PPM%f(k, 0.)
1012	0	Tbot = PPM%f(k, 1.)
1013	0	if ( Tbot>Ttop ) cycle ! The cell is inverted, skip to next
1014	0	if ( 20.<Tbot .and. Tbot<Ttop ) exit ! The whole remaining column is warmer than 20
1015	0	dz = d_from_ssh(k) - d_from_ssh(k-1) ! >=0
1016	0	if ( Tbot<=17. .and. 17.<=Ttop ) then
1017		! The 17 degC isotherm is within the cell which is non-negatively stratified
1018	0	d17(i,j) = d_from_ssh(k-1) + dz * PPM%x(k, 17.)
1019	0	elseif ( Ttop<17. ) then
1020		! The 17 degC isotherm is above the top of the cell
1021	0	d17(i,j) = d_from_ssh(k-1)
1022		endif
1023	0	if ( Tbot<=20. .and. 20.<=Ttop ) then
1024		! The 20 degC isotherm is within the cell which is non-negatively stratified
1025	0	d20(i,j) = d_from_ssh(k-1) + dz * PPM%x(k, 20.)
1026	0	elseif ( Ttop<20. ) then
1027		! The 20 degC isotherm is above the top of the cell
1028	0	d20(i,j) = d_from_ssh(k-1)
1029		endif
1030		enddo
1031		enddo ; enddo
1032	0	call PPM%destroy()
1033	0	if (CS%id_t17d > 0) call post_data(CS%id_t17d, d17, CS%diag)
1034	0	if (CS%id_t20d > 0) call post_data(CS%id_t20d, d20, CS%diag)
1035		endif
1036
1037		! Practical salinity expressed as salt mass content
1038	12	if (CS%id_scint > 0) then
1039	0	EOSdom(:) = EOS_domain(G%HI)
1040	0	if (tv%S_is_absS) then
1041	0	do k=1,nz ; do j=js,je
1042	0	call abs_saln_to_prac_saln(tv%S(:,j,k), tmp(:,j,k), tv%eqn_of_state, EOSdom) ! "tmp" [S ~> psu]
1043	0	do i=is,ie
1044	0	tmp(i,j,k) = ( GV%H_to_RZ * h(i,j,k) ) * tmp(i,j,k) ! "tmp" [R Z S ~> kg m-2]
1045		enddo
1046		enddo ; enddo
1047		else
1048	0	do k=1,nz ; do j=js,je ; do i=is,ie
1049	0	tmp(i,j,k) = ( GV%H_to_RZ * h(i,j,k) ) * tv%S(i,j,k) ! "tmp" [R Z S ~> kg m-2]
1050		enddo ; enddo ; enddo
1051		endif
1052	0	call post_data(CS%id_scint, tmp, CS%diag)
1053		endif
1054		! Absolute salinities expressed as salt mass content
1055	12	if (CS%id_absscint > 0 .or. CS%id_pfscint > 0) then
1056	0	EOSdom(:) = EOS_domain(G%HI)
1057	0	if (tv%S_is_absS) then
1058	0	do k=1,nz ; do j=js,je ; do i=is,ie
1059	0	tmp(i,j,k) = ( GV%H_to_RZ * h(i,j,k) ) * tv%S(i,j,k) ! "tmp" [R Z S ~> kg m-2]
1060		enddo ; enddo ; enddo
1061		else
1062	0	do k=1,nz ; do j=js,je
1063	0	call prac_saln_to_abs_saln(tv%S(:,j,k), tmp(:,j,k), tv%eqn_of_state, EOSdom) ! "tmp" [S ~> ppt]
1064	0	do i=is,ie
1065	0	tmp(i,j,k) = ( GV%H_to_RZ * h(i,j,k) ) * tmp(i,j,k) ! [R Z S ~> kg m-2]
1066		enddo
1067		enddo ; enddo
1068		endif
1069	0	if (CS%id_absscint > 0) call post_data(CS%id_absscint, tmp, CS%diag)
1070		! Based on the definitions in https://www.teos-10.org/pubs/gsw/pdf/TEOS-10_Manual.pdf
1071		! The preformed salinity, S*, is the conserved salinity used in models (page 8).
1072		! Although we appear to be labeling tv%S absolute salinity, we do not use the function
1073		! that calculates the "absolute salinity anomaly ratio" which accounts for the
1074		! geographic variations in the types of dissolved salts.
1075		! Hence, I think there is no difference between preformed and absolute salinity
1076		! for the current implementation of TEOS-10 and so we post the same data for
1077		! absscint and pfscint. -AJA
1078	0	if (CS%id_pfscint > 0) call post_data(CS%id_pfscint, tmp, CS%diag)
1079		endif
1080		! Potential temperature expressed as heat content
1081	12	if (CS%id_phcint > 0) then
1082	0	EOSdom(:) = EOS_domain(G%HI)
1083	0	if (tv%T_is_conT) then
1084	0	do k=1,nz ; do j=js,je
1085	0	call cons_temp_to_pot_temp(tv%T(:,j,k), tv%S(:,j,k), tmp(:,j,k), tv%eqn_of_state, EOSdom) ! "tmp" [C ~> degC]
1086	0	do i=is,ie
1087	0	tmp(i,j,k) = ( ( tv%C_p * GV%H_to_RZ ) * h(i,j,k) ) * tmp(i,j,k) ! "tmp" [ Q R Z ~> J m-2]
1088		enddo
1089		enddo ; enddo
1090		else
1091	0	do k=1,nz ; do j=js,je ; do i=is,ie
1092	0	tmp(i,j,k) = ( ( tv%C_p * GV%H_to_RZ ) * h(i,j,k) ) * tv%T(i,j,k) ! "tmp" [Q R Z ~> J m-2]
1093		enddo ; enddo ; enddo
1094		endif
1095	0	call post_data(CS%id_phcint, tmp, CS%diag)
1096		endif
1097		! Conservative temperature expressed as heat content
1098	12	if (CS%id_chcint > 0) then
1099	0	EOSdom(:) = EOS_domain(G%HI)
1100	0	if (tv%T_is_conT) then
1101	0	do k=1,nz ; do j=js,je ; do i=is,ie
1102	0	tmp(i,j,k) = ( ( tv%C_p * GV%H_to_RZ ) * h(i,j,k) ) * tv%T(i,j,k) ! "tmp" [Q R Z ~> J m-2]
1103		enddo ; enddo ; enddo
1104		else
1105	0	do k=1,nz ; do j=js,je
1106	0	call pot_temp_to_cons_temp(tv%T(:,j,k), tv%S(:,j,k), tmp(:,j,k), tv%eqn_of_state, EOSdom) ! "tmp" [C ~> degC]
1107	0	do i=is,ie
1108	0	tmp(i,j,k) = ( ( tv%C_p * GV%H_to_RZ ) * h(i,j,k) ) * tmp(i,j,k) ! "tmp" [ Q R Z ~> J m-2]
1109		enddo
1110		enddo ; enddo
1111		endif
1112	0	call post_data(CS%id_chcint, tmp, CS%diag)
1113		endif
1114
1115	12	end subroutine calculate_vertical_integrals
1116
1117		!> This subroutine calculates terms in the mechanical energy budget.
1118	12	subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS)
1119		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1120		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1121		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1122		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
1123		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1124		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
1125		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1126		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1127		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1128		intent(in) :: uh !< Transport through zonal faces=uhdy,
1129		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
1130		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1131		intent(in) :: vh !< Transport through merid faces=vhdx,
1132		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
1133		type(accel_diag_ptrs), intent(in) :: ADp !< Structure pointing to accelerations in momentum equation.
1134		type(cont_diag_ptrs), intent(in) :: CDp !< Structure pointing to terms in continuity equations.
1135		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1136		type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by a previous call to
1137		!! diagnostics_init.
1138
1139		! Local variables
1140	0	real :: KE(SZI_(G),SZJ_(G),SZK_(GV)) ! Kinetic energy per unit mass [L2 T-2 ~> m2 s-2]
1141	24	real :: KE_term(SZI_(G),SZJ_(G),SZK_(GV)) ! A term in the kinetic energy budget
1142		! [H L2 T-3 ~> m3 s-3 or W m-2]
1143	24	real :: KE_u(SZIB_(G),SZJ_(G)) ! The area integral of a KE term in a layer at u-points
1144		! [H L4 T-3 ~> m5 s-3 or W]
1145	24	real :: KE_v(SZI_(G),SZJB_(G)) ! The area integral of a KE term in a layer at v-points
1146		! [H L4 T-3 ~> m5 s-3 or W]
1147	12	real :: KE_h(SZI_(G),SZJ_(G)) ! A KE term contribution at tracer points
1148		! [H L2 T-3 ~> m3 s-3 or W m-2]
1149
1150		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1151	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1152	12	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1153
1154	12	if (.not.(CS%KE_term_on .or. (CS%id_KE > 0))) return
1155
1156	0	KE_u(:,:) = 0. ; KE_v(:,:) = 0.
1157
1158	0	do k=1,nz ; do j=js,je ; do i=is,ie
1159		KE(i,j,k) = (((u(I,j,k) * u(I,j,k)) + (u(I-1,j,k) * u(I-1,j,k))) &
1160	0	+ ((v(i,J,k) * v(i,J,k)) + (v(i,J-1,k) * v(i,J-1,k)))) * 0.25
1161		enddo ; enddo ; enddo
1162	0	if (CS%id_KE > 0) call post_data(CS%id_KE, KE, CS%diag)
1163
1164	0	if (CS%KE_term_on .and. .not.G%symmetric) then
1165	0	call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East)
1166		endif
1167
1168	0	if (CS%id_dKEdt > 0) then
1169		! Calculate the time derivative of the layer KE [H L2 T-3 ~> m3 s-3 or W m-2].
1170	0	do k=1,nz
1171	0	do j=js,je ; do I=Isq,Ieq
1172	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * CS%du_dt(I,j,k)
1173		enddo ; enddo
1174	0	do J=Jsq,Jeq ; do i=is,ie
1175	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * CS%dv_dt(i,J,k)
1176		enddo ; enddo
1177	0	do j=js,je ; do i=is,ie
1178	0	KE_h(i,j) = KE(i,j,k) * CS%dh_dt(i,j,k)
1179		enddo ; enddo
1180	0	if (.not.G%symmetric) &
1181	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1182	0	do j=js,je ; do i=is,ie
1183		KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) &
1184	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1185		enddo ; enddo
1186		enddo
1187	0	call post_data(CS%id_dKEdt, KE_term, CS%diag)
1188		endif
1189
1190	0	if (CS%id_PE_to_KE > 0) then
1191		! Calculate the potential energy to KE term [H L2 T-3 ~> m3 s-3 or W m-2].
1192	0	do k=1,nz
1193	0	do j=js,je ; do I=Isq,Ieq
1194	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu(I,j,k)
1195		enddo ; enddo
1196	0	do J=Jsq,Jeq ; do i=is,ie
1197	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv(i,J,k)
1198		enddo ; enddo
1199	0	if (.not.G%symmetric) &
1200	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1201	0	do j=js,je ; do i=is,ie
1202		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1203	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1204		enddo ; enddo
1205		enddo
1206	0	call post_data(CS%id_PE_to_KE, KE_term, CS%diag)
1207		endif
1208
1209	0	if (CS%id_KE_SAL > 0) then
1210		! Calculate the KE source from self-attraction and loading [H L2 T-3 ~> m3 s-3 or W m-2].
1211	0	do k=1,nz
1212	0	do j=js,je ; do I=Isq,Ieq
1213	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%sal_u(I,j,k)
1214		enddo ; enddo
1215	0	do J=Jsq,Jeq ; do i=is,ie
1216	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%sal_v(i,J,k)
1217		enddo ; enddo
1218	0	if (.not.G%symmetric) &
1219	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1220	0	do j=js,je ; do i=is,ie
1221		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1222	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1223		enddo ; enddo
1224		enddo
1225	0	call post_data(CS%id_KE_SAL, KE_term, CS%diag)
1226		endif
1227
1228	0	if (CS%id_KE_TIDES > 0) then
1229		! Calculate the KE source from astronomical tidal forcing [H L2 T-3 ~> m3 s-3 or W m-2].
1230	0	do k=1,nz
1231	0	do j=js,je ; do I=Isq,Ieq
1232	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%tides_u(I,j,k)
1233		enddo ; enddo
1234	0	do J=Jsq,Jeq ; do i=is,ie
1235	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%tides_v(i,J,k)
1236		enddo ; enddo
1237	0	if (.not.G%symmetric) &
1238	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1239	0	do j=js,je ; do i=is,ie
1240		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1241	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1242		enddo ; enddo
1243		enddo
1244	0	call post_data(CS%id_KE_TIDES, KE_term, CS%diag)
1245		endif
1246
1247	0	if (CS%id_KE_BT > 0) then
1248		! Calculate the barotropic contribution to KE term [H L2 T-3 ~> m3 s-3 or W m-2].
1249	0	do k=1,nz
1250	0	do j=js,je ; do I=Isq,Ieq
1251	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt(I,j,k)
1252		enddo ; enddo
1253	0	do J=Jsq,Jeq ; do i=is,ie
1254	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt(i,J,k)
1255		enddo ; enddo
1256	0	if (.not.G%symmetric) &
1257	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1258	0	do j=js,je ; do i=is,ie
1259		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1260	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1261		enddo ; enddo
1262		enddo
1263	0	call post_data(CS%id_KE_BT, KE_term, CS%diag)
1264		endif
1265
1266	0	if (CS%id_PE_to_KE_btbc > 0) then
1267		! Calculate the potential energy to KE term including barotropic solver contribution
1268		! [H L2 T-3 ~> m3 s-3 or W m-2].
1269	0	do k=1,nz
1270	0	do j=js,je ; do I=Isq,Ieq
1271	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * (ADp%PFu(I,j,k) + ADp%bt_pgf_u(I,j,k))
1272		enddo ; enddo
1273	0	do J=Jsq,Jeq ; do i=is,ie
1274	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * (ADp%PFv(i,J,k) + ADp%bt_pgf_v(i,J,k))
1275		enddo ; enddo
1276	0	if (.not.G%symmetric) &
1277	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1278	0	do j=js,je ; do i=is,ie
1279		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1280	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1281		enddo ; enddo
1282		enddo
1283	0	call post_data(CS%id_PE_to_KE_btbc, KE_term, CS%diag)
1284		endif
1285
1286	0	if (CS%id_KE_Coradv_btbc > 0) then
1287		! Calculate the KE source from Coriolis and advection terms including barotropic solver contribution
1288		! [H L2 T-3 ~> m3 s-3 or W m-2].
1289	0	do k=1,nz
1290	0	do j=js,je ; do I=Isq,Ieq
1291	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * (ADp%CAu(I,j,k) + ADp%bt_cor_u(I,j))
1292		enddo ; enddo
1293	0	do J=Jsq,Jeq ; do i=is,ie
1294	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * (ADp%CAv(i,J,k) + ADp%bt_cor_v(i,J))
1295		enddo ; enddo
1296	0	do j=js,je ; do i=is,ie
1297		KE_h(i,j) = -KE(i,j,k) * G%IareaT(i,j) &
1298	0	* ((uh(I,j,k) - uh(I-1,j,k)) + (vh(i,J,k) - vh(i,J-1,k)))
1299		enddo ; enddo
1300	0	if (.not.G%symmetric) &
1301	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1302	0	do j=js,je ; do i=is,ie
1303		KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) &
1304	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1305		enddo ; enddo
1306		enddo
1307	0	call post_data(CS%id_KE_Coradv_btbc, KE_term, CS%diag)
1308		endif
1309
1310	0	if (CS%id_KE_BT_PF > 0) then
1311		! Calculate the anomalous pressure gradient force contribution to KE term [H L2 T-3 ~> m3 s-3 or W m-2].
1312	0	do k=1,nz
1313	0	do j=js,je ; do I=Isq,Ieq
1314	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%bt_pgf_u(I,j,k)
1315		enddo ; enddo
1316	0	do J=Jsq,Jeq ; do i=is,ie
1317	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%bt_pgf_v(i,J,k)
1318		enddo ; enddo
1319	0	if (.not.G%symmetric) &
1320	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1321	0	do j=js,je ; do i=is,ie
1322		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1323	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1324		enddo ; enddo
1325		enddo
1326	0	call post_data(CS%id_KE_BT_PF, KE_term, CS%diag)
1327		endif
1328
1329	0	if (CS%id_KE_BT_CF > 0) then
1330		! Calculate the anomalous Coriolis force contribution to KE term [H L2 T-3 ~> m3 s-3 or W m-2].
1331	0	do k=1,nz
1332	0	do j=js,je ; do I=Isq,Ieq
1333	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%bt_cor_u(I,j)
1334		enddo ; enddo
1335	0	do J=Jsq,Jeq ; do i=is,ie
1336	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%bt_cor_v(i,J)
1337		enddo ; enddo
1338	0	if (.not.G%symmetric) &
1339	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1340	0	do j=js,je ; do i=is,ie
1341		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1342	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1343		enddo ; enddo
1344		enddo
1345	0	call post_data(CS%id_KE_BT_CF, KE_term, CS%diag)
1346		endif
1347
1348	0	if (CS%id_KE_BT_WD > 0) then
1349		! Calculate the barotropic linear wave drag contribution to KE term [H L2 T-3 ~> m3 s-3 or W m-2].
1350	0	do k=1,nz
1351	0	do j=js,je ; do I=Isq,Ieq
1352	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%bt_lwd_u(I,j)
1353		enddo ; enddo
1354	0	do J=Jsq,Jeq ; do i=is,ie
1355	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%bt_lwd_v(i,J)
1356		enddo ; enddo
1357	0	if (.not.G%symmetric) &
1358	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1359	0	do j=js,je ; do i=is,ie
1360		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1361	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1362		enddo ; enddo
1363		enddo
1364	0	call post_data(CS%id_KE_BT_WD, KE_term, CS%diag)
1365		endif
1366
1367	0	if (CS%id_KE_Coradv > 0) then
1368		! Calculate the KE source from the combined Coriolis and advection terms [H L2 T-3 ~> m3 s-3 or W m-2].
1369		! The Coriolis source should be zero, but is not due to truncation errors. There should be
1370		! near-cancellation of the global integral of this spurious Coriolis source.
1371	0	do k=1,nz
1372	0	do j=js,je ; do I=Isq,Ieq
1373	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%CAu(I,j,k)
1374		enddo ; enddo
1375	0	do J=Jsq,Jeq ; do i=is,ie
1376	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%CAv(i,J,k)
1377		enddo ; enddo
1378	0	do j=js,je ; do i=is,ie
1379		KE_h(i,j) = -KE(i,j,k) * G%IareaT(i,j) &
1380	0	* ((uh(I,j,k) - uh(I-1,j,k)) + (vh(i,J,k) - vh(i,J-1,k)))
1381		enddo ; enddo
1382	0	if (.not.G%symmetric) &
1383	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1384	0	do j=js,je ; do i=is,ie
1385		KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) &
1386	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1387		enddo ; enddo
1388		enddo
1389	0	call post_data(CS%id_KE_Coradv, KE_term, CS%diag)
1390		endif
1391
1392	0	if (CS%id_KE_adv > 0) then
1393		! Calculate the KE source from along-layer advection [H L2 T-3 ~> m3 s-3 or W m-2].
1394		! NOTE: All terms in KE_adv are multiplied by -1, which can easily produce
1395		! negative zeros and may signal a reproducibility issue over land.
1396		! We resolve this by re-initializing and only evaluating over water points.
1397	0	KE_u(:,:) = 0. ; KE_v(:,:) = 0.
1398	0	do k=1,nz
1399	0	do j=js,je ; do I=Isq,Ieq
1400	0	if (G%mask2dCu(i,j) /= 0.) &
1401	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%gradKEu(I,j,k)
1402		enddo ; enddo
1403	0	do J=Jsq,Jeq ; do i=is,ie
1404	0	if (G%mask2dCv(i,j) /= 0.) &
1405	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%gradKEv(i,J,k)
1406		enddo ; enddo
1407	0	do j=js,je ; do i=is,ie
1408		KE_h(i,j) = -KE(i,j,k) * G%IareaT(i,j) &
1409	0	* ((uh(I,j,k) - uh(I-1,j,k)) + (vh(i,J,k) - vh(i,J-1,k)))
1410		enddo ; enddo
1411	0	if (.not.G%symmetric) &
1412	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1413	0	do j=js,je ; do i=is,ie
1414		KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) &
1415	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1416		enddo ; enddo
1417		enddo
1418	0	call post_data(CS%id_KE_adv, KE_term, CS%diag)
1419		endif
1420
1421	0	if (CS%id_KE_visc > 0) then
1422		! Calculate the KE source from vertical viscosity [H L2 T-3 ~> m3 s-3 or W m-2].
1423	0	do k=1,nz
1424	0	do j=js,je ; do I=Isq,Ieq
1425	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc(I,j,k)
1426		enddo ; enddo
1427	0	do J=Jsq,Jeq ; do i=is,ie
1428	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc(i,J,k)
1429		enddo ; enddo
1430	0	if (.not.G%symmetric) &
1431	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1432	0	do j=js,je ; do i=is,ie
1433		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1434	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1435		enddo ; enddo
1436		enddo
1437	0	call post_data(CS%id_KE_visc, KE_term, CS%diag)
1438		endif
1439
1440	0	if (CS%id_KE_visc_gl90 > 0) then
1441		! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3 or W m-2].
1442	0	do k=1,nz
1443	0	do j=js,je ; do I=Isq,Ieq
1444	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc_gl90(I,j,k)
1445		enddo ; enddo
1446	0	do J=Jsq,Jeq ; do i=is,ie
1447	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc_gl90(i,J,k)
1448		enddo ; enddo
1449	0	if (.not.G%symmetric) &
1450	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1451	0	do j=js,je ; do i=is,ie
1452		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1453	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1454		enddo ; enddo
1455		enddo
1456	0	call post_data(CS%id_KE_visc_gl90, KE_term, CS%diag)
1457		endif
1458
1459	0	if (CS%id_KE_stress > 0) then
1460		! Calculate the KE source from surface stress (included in KE_visc) [H L2 T-3 ~> m3 s-3 or W m-2].
1461	0	do k=1,nz
1462	0	do j=js,je ; do I=Isq,Ieq
1463	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_str(I,j,k)
1464		enddo ; enddo
1465	0	do J=Jsq,Jeq ; do i=is,ie
1466	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_str(i,J,k)
1467		enddo ; enddo
1468	0	if (.not.G%symmetric) &
1469	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1470	0	do j=js,je ; do i=is,ie
1471		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) * &
1472	0	((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1473		enddo ; enddo
1474		enddo
1475	0	call post_data(CS%id_KE_stress, KE_term, CS%diag)
1476		endif
1477
1478	0	if (CS%id_KE_horvisc > 0) then
1479		! Calculate the KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3 or W m-2].
1480	0	do k=1,nz
1481	0	do j=js,je ; do I=Isq,Ieq
1482	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu(I,j,k)
1483		enddo ; enddo
1484	0	do J=Jsq,Jeq ; do i=is,ie
1485	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv(i,J,k)
1486		enddo ; enddo
1487	0	if (.not.G%symmetric) &
1488	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1489	0	do j=js,je ; do i=is,ie
1490		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1491	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1492		enddo ; enddo
1493		enddo
1494	0	call post_data(CS%id_KE_horvisc, KE_term, CS%diag)
1495		endif
1496
1497	0	if (CS%id_KE_dia > 0) then
1498		! Calculate the KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3 or W m-2].
1499	0	do k=1,nz
1500	0	do j=js,je ; do I=Isq,Ieq
1501	0	KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_dia(I,j,k)
1502		enddo ; enddo
1503	0	do J=Jsq,Jeq ; do i=is,ie
1504	0	KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_dia(i,J,k)
1505		enddo ; enddo
1506	0	do j=js,je ; do i=is,ie
1507	0	KE_h(i,j) = KE(i,j,k) * (CDp%diapyc_vel(i,j,k) - CDp%diapyc_vel(i,j,k+1))
1508		enddo ; enddo
1509	0	if (.not.G%symmetric) &
1510	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1511	0	do j=js,je ; do i=is,ie
1512		KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) &
1513	0	* ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1)))
1514		enddo ; enddo
1515		enddo
1516	0	call post_data(CS%id_KE_dia, KE_term, CS%diag)
1517		endif
1518
1519		end subroutine calculate_energy_diagnostics
1520
1521		!> This subroutine registers fields to calculate a diagnostic time derivative.
1522	0	subroutine register_time_deriv(lb, f_ptr, deriv_ptr, CS)
1523		integer, intent(in), dimension(3) :: lb !< Lower index bound of f_ptr
1524		real, dimension(lb(1):,lb(2):,:), target :: f_ptr
1525		!< Time derivative operand, in arbitrary units [A ~> a]
1526		real, dimension(lb(1):,lb(2):,:), target :: deriv_ptr
1527		!< Time derivative of f_ptr, in units derived from
1528		!! the arbitrary units of f_ptr [A T-1 ~> a s-1]
1529		type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by previous call to
1530		!! diagnostics_init.
1531
1532		! This subroutine registers fields to calculate a diagnostic time derivative.
1533		! NOTE: Lower bound is required for grid indexing in calculate_derivs().
1534		! We assume that the vertical axis is 1-indexed.
1535
1536		integer :: m !< New index of deriv_ptr in CS%deriv
1537		integer :: ub(3) !< Upper index bound of f_ptr, based on shape.
1538
1539	0	if (.not.CS%initialized) call MOM_error(FATAL, &
1540	0	"register_time_deriv: Module must be initialized before it is used.")
1541
1542	0	if (CS%num_time_deriv >= MAX_FIELDS_) then
1543		call MOM_error(WARNING,"MOM_diagnostics: Attempted to register more than " // &
1544	0	"MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.")
1545	0	return
1546		endif
1547
1548	0	m = CS%num_time_deriv+1 ; CS%num_time_deriv = m
1549
1550	0	ub(:) = lb(:) + shape(f_ptr) - 1
1551
1552	0	CS%nlay(m) = size(f_ptr, 3)
1553	0	CS%deriv(m)%p => deriv_ptr
1554	0	allocate(CS%prev_val(m)%p(lb(1):ub(1), lb(2):ub(2), CS%nlay(m)))
1555
1556	0	CS%var_ptr(m)%p => f_ptr
1557	0	CS%prev_val(m)%p(:,:,:) = f_ptr(:,:,:)
1558
1559		end subroutine register_time_deriv
1560
1561		!> This subroutine calculates all registered time derivatives.
1562	12	subroutine calculate_derivs(dt, G, CS)
1563		real, intent(in) :: dt !< The time interval over which differences occur [T ~> s].
1564		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1565		type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by previous call to
1566		!! diagnostics_init.
1567
1568		! This subroutine calculates all registered time derivatives.
1569		real :: Idt ! The inverse timestep [T-1 ~> s-1]
1570		integer :: i, j, k, m
1571
1572	12	if (dt > 0.0) then ; Idt = 1.0/dt
1573	0	else ; return ; endif
1574
1575		! Because the field is unknown, its grid index bounds are also unknown.
1576		! Additionally, two of the fields (dudt, dvdt) require calculation of spatial
1577		! derivatives when computing d(KE)/dt. This raises issues in non-symmetric
1578		! mode, where the symmetric boundaries (west, south) may not be updated.
1579
1580		! For this reason, we explicitly loop from isc-1:iec and jsc-1:jec, in order
1581		! to force boundary value updates, even though it may not be strictly valid
1582		! for all fields. Note this assumes a halo, and that it has been updated.
1583
1584	12	do m=1,CS%num_time_deriv
1585	0	do k=1,CS%nlay(m) ; do j=G%jsc-1,G%jec ; do i=G%isc-1,G%iec
1586	0	CS%deriv(m)%p(i,j,k) = (CS%var_ptr(m)%p(i,j,k) - CS%prev_val(m)%p(i,j,k)) * Idt
1587	0	CS%prev_val(m)%p(i,j,k) = CS%var_ptr(m)%p(i,j,k)
1588		enddo ; enddo ; enddo
1589		enddo
1590
1591		end subroutine calculate_derivs
1592
1593		!> This routine posts diagnostics of various dynamic ocean surface quantities,
1594		!! including velocities, speed and sea surface height, at the time the ocean
1595		!! state is reported back to the caller
1596	12	subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh)
1597		type(surface_diag_IDs), intent(in) :: IDs !< A structure with the diagnostic IDs.
1598		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1599		type(diag_ctrl), intent(in) :: diag !< regulates diagnostic output
1600		type(surface), intent(in) :: sfc_state !< structure describing the ocean surface state
1601		real, dimension(SZI_(G),SZJ_(G)), &
1602		intent(in) :: ssh !< Time mean surface height without corrections
1603		!! for ice displacement [Z ~> m]
1604
1605		! Local variables
1606	24	real, dimension(SZI_(G),SZJ_(G)) :: speed ! The surface speed [L T-1 ~> m s-1]
1607	24	real :: ssu_east(SZI_(G),SZJ_(G)) ! Surface velocity due east component [L T-1 ~> m s-1]
1608	24	real :: ssv_north(SZI_(G),SZJ_(G)) ! Surface velocity due north component [L T-1 ~> m s-1]
1609		integer :: i, j, is, ie, js, je
1610
1611	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1612
1613	12	if (IDs%id_ssh > 0) &
1614	0	call post_data(IDs%id_ssh, ssh, diag)
1615
1616	12	if (IDs%id_ssu > 0) &
1617	12	call post_data(IDs%id_ssu, sfc_state%u, diag)
1618
1619	12	if (IDs%id_ssv > 0) &
1620	12	call post_data(IDs%id_ssv, sfc_state%v, diag)
1621
1622	12	if (IDs%id_speed > 0) then
1623	0	do j=js,je ; do i=is,ie
1624		speed(i,j) = sqrt(0.5((sfc_state%u(I-1,j)2) + (sfc_state%u(I,j)*2)) + &
1625	0	0.5((sfc_state%v(i,J-1)2) + (sfc_state%v(i,J)*2)))
1626		enddo ; enddo
1627	0	call post_data(IDs%id_speed, speed, diag)
1628		endif
1629
1630	12	if (IDs%id_ssu_east > 0 .or. IDs%id_ssv_north > 0) then
1631	0	do j=js,je ; do i=is,ie
1632		ssu_east(i,j) = ((0.5(sfc_state%u(I-1,j) + sfc_state%u(I,j))) G%cos_rot(i,j)) + &
1633	0	((0.5(sfc_state%v(i,J-1) + sfc_state%v(i,J))) G%sin_rot(i,j))
1634		ssv_north(i,j) = ((0.5(sfc_state%v(i,J-1) + sfc_state%v(i,J))) G%cos_rot(i,j)) - &
1635	0	((0.5(sfc_state%u(I-1,j) + sfc_state%u(I,j))) G%sin_rot(i,j))
1636		enddo ; enddo
1637	0	if (IDs%id_ssu_east > 0 ) call post_data(IDs%id_ssu_east, ssu_east, diag)
1638	0	if (IDs%id_ssv_north > 0 ) call post_data(IDs%id_ssv_north, ssv_north, diag)
1639		endif
1640
1641	12	end subroutine post_surface_dyn_diags
1642
1643
1644		!> This routine posts diagnostics of various ocean surface and integrated
1645		!! quantities at the time the ocean state is reported back to the caller
1646	12	subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv, &
1647	12	ssh, ssh_ibc)
1648		type(surface_diag_IDs), intent(in) :: IDs !< A structure with the diagnostic IDs.
1649		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1650		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1651		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1652		type(diag_ctrl), intent(in) :: diag !< regulates diagnostic output
1653		real, intent(in) :: dt_int !< total time step associated with these diagnostics [T ~> s].
1654		type(surface), intent(in) :: sfc_state !< structure describing the ocean surface state
1655		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
1656		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh !< Time mean surface height without corrections
1657		!! for ice displacement [Z ~> m]
1658		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections
1659		!! for ice displacement and the inverse barometer [Z ~> m]
1660
1661	24	real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array [various]
1662		real, dimension(SZI_(G),SZJ_(G)) :: &
1663	24	zos ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [Z ~> m]
1664		real :: I_time_int ! The inverse of the time interval [T-1 ~> s-1].
1665		real :: zos_area_mean ! Global area mean sea surface height [Z ~> m]
1666		real :: volo ! Total volume of the ocean [Z L2 ~> m3]
1667		real :: ssh_ga ! Global ocean area weighted mean sea seaface height [Z ~> m]
1668		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1669		integer :: i, j, is, ie, js, je
1670
1671	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1672
1673		! area mean SSH
1674	12	if (IDs%id_ssh_ga > 0) then
1675	0	ssh_ga = global_area_mean(ssh, G, tmp_scale=US%Z_to_m)
1676	0	call post_data(IDs%id_ssh_ga, ssh_ga, diag)
1677		endif
1678
1679		! post the dynamic sea level, zos, and zossq.
1680		! zos is ave_ssh with sea ice inverse barometer removed, and with zero global area mean.
1681	12	if (IDs%id_zos > 0 .or. IDs%id_zossq > 0) then
1682	12	zos_area_mean = global_area_mean(ssh_ibc, G, tmp_scale=US%Z_to_m)
1683	87132	do j=js,je ; do i=is,ie
1684	87120	zos(i,j) = ssh_ibc(i,j) - G%mask2dT(i,j)*zos_area_mean
1685		enddo ; enddo
1686	12	if (IDs%id_zos > 0) call post_data(IDs%id_zos, zos, diag)
1687	12	if (IDs%id_zossq > 0) then
1688	0	do j=js,je ; do i=is,ie
1689	0	work_2d(i,j) = zos(i,j)*zos(i,j)
1690		enddo ; enddo
1691	0	call post_data(IDs%id_zossq, work_2d, diag)
1692		endif
1693		endif
1694
1695		! post total volume of the liquid ocean
1696	12	if (IDs%id_volo > 0) then
1697	0	do j=js,je ; do i=is,ie
1698	0	work_2d(i,j) = G%mask2dT(i,j) * (ssh(i,j) + G%bathyT(i,j))
1699		enddo ; enddo
1700	0	volo = global_area_integral(work_2d, G, tmp_scale=US%Z_to_m)
1701	0	call post_data(IDs%id_volo, volo, diag)
1702		endif
1703
1704		! Use Adcroft's rule of reciprocals; it does the right thing here.
1705	12	I_time_int = 0.0 ; if (dt_int > 0.0) I_time_int = 1.0 / dt_int
1706
1707		! post time-averaged rate of frazil formation
1708	12	if (associated(tv%frazil) .and. (IDs%id_fraz > 0)) then
1709	0	do j=js,je ; do i=is,ie
1710	0	work_2d(i,j) = tv%frazil(i,j) * I_time_int
1711		enddo ; enddo
1712	0	call post_data(IDs%id_fraz, work_2d, diag)
1713		endif
1714
1715		! post time-averaged salt deficit
1716	12	if (associated(tv%salt_deficit) .and. (IDs%id_salt_deficit > 0)) then
1717	0	do j=js,je ; do i=is,ie
1718	0	work_2d(i,j) = tv%salt_deficit(i,j) * I_time_int
1719		enddo ; enddo
1720	0	call post_data(IDs%id_salt_deficit, work_2d, diag)
1721		endif
1722
1723		! post temperature of P-E+R
1724	12	if (associated(tv%TempxPmE) .and. (IDs%id_Heat_PmE > 0)) then
1725	0	do j=js,je ; do i=is,ie
1726	0	work_2d(i,j) = tv%TempxPmE(i,j) * (tv%C_p * I_time_int)
1727		enddo ; enddo
1728	0	call post_data(IDs%id_Heat_PmE, work_2d, diag)
1729		endif
1730
1731		! post geothermal heating or internal heat source/sinks
1732	12	if (associated(tv%internal_heat) .and. (IDs%id_intern_heat > 0)) then
1733	0	do j=js,je ; do i=is,ie
1734	0	work_2d(i,j) = tv%internal_heat(i,j) * (tv%C_p * I_time_int)
1735		enddo ; enddo
1736	0	call post_data(IDs%id_intern_heat, work_2d, diag)
1737		endif
1738
1739	12	if (tv%T_is_conT) then
1740		! Internal T&S variables are conservative temperature & absolute salinity
1741	0	if (IDs%id_sstcon > 0) call post_data(IDs%id_sstcon, sfc_state%SST, diag)
1742		! Use TEOS-10 function calls convert T&S diagnostics from conservative temp
1743		! to potential temperature.
1744	0	EOSdom(:) = EOS_domain(G%HI)
1745	0	do j=js,je
1746	0	call cons_temp_to_pot_temp(sfc_state%SST(:,j), sfc_state%SSS(:,j), work_2d(:,j), tv%eqn_of_state, EOSdom)
1747		enddo
1748	0	if (IDs%id_sst > 0) call post_data(IDs%id_sst, work_2d, diag)
1749		else
1750		! Internal T&S variables are potential temperature & practical salinity
1751	12	if (IDs%id_sst > 0) call post_data(IDs%id_sst, sfc_state%SST, diag)
1752		endif
1753
1754	12	if (tv%S_is_absS) then
1755		! Internal T&S variables are conservative temperature & absolute salinity
1756	0	if (IDs%id_sssabs > 0) call post_data(IDs%id_sssabs, sfc_state%SSS, diag)
1757		! Use TEOS-10 function calls convert T&S diagnostics from absolute salinity
1758		! to practical salinity.
1759	0	EOSdom(:) = EOS_domain(G%HI)
1760	0	do j=js,je
1761	0	call abs_saln_to_prac_saln(sfc_state%SSS(:,j), work_2d(:,j), tv%eqn_of_state, EOSdom)
1762		enddo
1763	0	if (IDs%id_sss > 0) call post_data(IDs%id_sss, work_2d, diag)
1764		else
1765		! Internal T&S variables are potential temperature & practical salinity
1766	12	if (IDs%id_sss > 0) call post_data(IDs%id_sss, sfc_state%SSS, diag)
1767		endif
1768
1769	12	if (IDs%id_sst_sq > 0) then
1770	0	do j=js,je ; do i=is,ie
1771	0	work_2d(i,j) = sfc_state%SST(i,j)*sfc_state%SST(i,j)
1772		enddo ; enddo
1773	0	call post_data(IDs%id_sst_sq, work_2d, diag)
1774		endif
1775	12	if (IDs%id_sss_sq > 0) then
1776	0	do j=js,je ; do i=is,ie
1777	0	work_2d(i,j) = sfc_state%SSS(i,j)*sfc_state%SSS(i,j)
1778		enddo ; enddo
1779	0	call post_data(IDs%id_sss_sq, work_2d, diag)
1780		endif
1781
1782	12	call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag))
1783
1784	12	end subroutine post_surface_thermo_diags
1785
1786
1787		!> This routine posts diagnostics of the transports, including the subgridscale
1788		!! contributions.
1789	12	subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dyn, &
1790		diag, dt_trans, Reg)
1791		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1792		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1793		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1794		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: uhtr !< Accumulated zonal thickness fluxes
1795		!! used to advect tracers [H L2 ~> m3 or kg]
1796		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vhtr !< Accumulated meridional thickness fluxes
1797		!! used to advect tracers [H L2 ~> m3 or kg]
1798		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1799		intent(in) :: h !< The updated layer thicknesses [H ~> m or kg m-2]
1800		type(transport_diag_IDs), intent(in) :: IDs !< A structure with the diagnostic IDs.
1801		type(diag_grid_storage), intent(inout) :: diag_pre_dyn !< Stored grids from before dynamics
1802		type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
1803		real, intent(in) :: dt_trans !< total time step associated with the transports [T ~> s].
1804		type(tracer_registry_type), pointer :: Reg !< Pointer to the tracer registry
1805
1806		! Local variables
1807	24	real, dimension(SZIB_(G),SZJ_(G)) :: umo2d ! Diagnostics of integrated mass transport [R Z L2 T-1 ~> kg s-1]
1808	24	real, dimension(SZI_(G),SZJB_(G)) :: vmo2d ! Diagnostics of integrated mass transport [R Z L2 T-1 ~> kg s-1]
1809	24	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: umo ! Diagnostics of layer mass transport [R Z L2 T-1 ~> kg s-1]
1810	24	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vmo ! Diagnostics of layer mass transport [R Z L2 T-1 ~> kg s-1]
1811	24	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_tend ! Change in layer thickness due to dynamics
1812		! [H T-1 ~> m s-1 or kg m-2 s-1].
1813		real :: Idt ! The inverse of the time interval [T-1 ~> s-1]
1814		real :: H_to_RZ_dt ! A conversion factor from accumulated transports to fluxes
1815		! [R Z H-1 T-1 ~> kg m-3 s-1 or s-1].
1816		integer :: i, j, k, is, ie, js, je, nz
1817	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1818
1819	12	Idt = 1. / dt_trans
1820	12	H_to_RZ_dt = GV%H_to_RZ * Idt
1821
1822	12	call diag_save_grids(diag)
1823	12	call diag_copy_storage_to_diag(diag, diag_pre_dyn)
1824
1825		!$omp target update from(uhtr) if (any([IDs%id_umo_2d, IDs%id_umo, IDs%id_uhtr] > 0))
1826		!$omp target update from(vhtr) if (any([IDs%id_vmo_2d, IDs%id_vmo, IDs%id_vhtr] > 0))
1827		!$omp target update from(h) if (IDs%id_dynamics_h_tendency > 0)
1828
1829	12	if (IDs%id_umo_2d > 0) then
1830	0	umo2d(:,:) = 0.0
1831	0	do k=1,nz ; do j=js,je ; do I=is-1,ie
1832	0	umo2d(I,j) = umo2d(I,j) + uhtr(I,j,k) * H_to_RZ_dt
1833		enddo ; enddo ; enddo
1834	0	call post_data(IDs%id_umo_2d, umo2d, diag)
1835		endif
1836	12	if (IDs%id_umo > 0) then
1837		! Convert to kg/s.
1838	0	do k=1,nz ; do j=js,je ; do I=is-1,ie
1839	0	umo(I,j,k) = uhtr(I,j,k) * H_to_RZ_dt
1840		enddo ; enddo ; enddo
1841	0	call post_data(IDs%id_umo, umo, diag, alt_h=diag_pre_dyn%h_state)
1842		endif
1843	12	if (IDs%id_vmo_2d > 0) then
1844	0	vmo2d(:,:) = 0.0
1845	0	do k=1,nz ; do J=js-1,je ; do i=is,ie
1846	0	vmo2d(i,J) = vmo2d(i,J) + vhtr(i,J,k) * H_to_RZ_dt
1847		enddo ; enddo ; enddo
1848	0	call post_data(IDs%id_vmo_2d, vmo2d, diag)
1849		endif
1850	12	if (IDs%id_vmo > 0) then
1851		! Convert to kg/s.
1852	0	do k=1,nz ; do J=js-1,je ; do i=is,ie
1853	0	vmo(i,J,k) = vhtr(i,J,k) * H_to_RZ_dt
1854		enddo ; enddo ; enddo
1855	0	call post_data(IDs%id_vmo, vmo, diag, alt_h=diag_pre_dyn%h_state)
1856		endif
1857
1858	12	if (IDs%id_uhtr > 0) call post_data(IDs%id_uhtr, uhtr, diag, alt_h=diag_pre_dyn%h_state)
1859	12	if (IDs%id_vhtr > 0) call post_data(IDs%id_vhtr, vhtr, diag, alt_h=diag_pre_dyn%h_state)
1860	12	if (IDs%id_dynamics_h > 0) call post_data(IDs%id_dynamics_h, diag_pre_dyn%h_state, diag, &
1861	0	alt_h=diag_pre_dyn%h_state)
1862		! Post the change in thicknesses
1863	12	if (IDs%id_dynamics_h_tendency > 0) then
1864	0	h_tend(:,:,:) = 0.
1865	0	do k=1,nz ; do j=js,je ; do i=is,ie
1866	0	h_tend(i,j,k) = (h(i,j,k) - diag_pre_dyn%h_state(i,j,k))*Idt
1867		enddo ; enddo ; enddo
1868	0	call post_data(IDs%id_dynamics_h_tendency, h_tend, diag, alt_h=diag_pre_dyn%h_state)
1869		endif
1870
1871	12	call post_tracer_transport_diagnostics(G, GV, Reg, diag_pre_dyn%h_state, diag)
1872
1873	12	call diag_restore_grids(diag)
1874
1875	12	end subroutine post_transport_diagnostics
1876
1877		!> This subroutine registers various diagnostics and allocates space for fields
1878		!! that other diagnostics depend upon.
1879	1	subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag, CS, tv)
1880		type(ocean_internal_state), intent(in) :: MIS !< For "MOM Internal State" a set of pointers to
1881		!! the fields and accelerations that make up the
1882		!! ocean's internal physical state.
1883		type(accel_diag_ptrs), intent(inout) :: ADp !< Structure with pointers to momentum equation
1884		!! terms.
1885		type(cont_diag_ptrs), intent(inout) :: CDp !< Structure with pointers to continuity
1886		!! equation terms.
1887		type(time_type), intent(in) :: Time !< Current model time.
1888		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1889		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1890		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1891		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1892		!! parameters.
1893		type(diag_ctrl), target, intent(inout) :: diag !< Structure to regulate diagnostic output.
1894		type(diagnostics_CS), intent(inout) :: CS !< Diagnostic control struct
1895		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
1896		!! thermodynamic variables.
1897
1898		! Local variables
1899		real :: wave_speed_min ! A floor in the first mode speed below which 0 is returned [L T-1 ~> m s-1]
1900		real :: wave_speed_tol ! The fractional tolerance for finding the wave speeds [nondim]
1901		real :: convert_H ! A conversion factor from internal thickness units to the appropriate
1902		! MKS units (m or kg m-2) for thicknesses depending on whether the
1903		! Boussinesq approximation is being made [m H-1 ~> 1] or [kg m-2 H-1 ~> 1]
1904		logical :: better_speed_est ! If true, use a more robust estimate of the first
1905		! mode wave speed as the starting point for iterations.
1906		logical :: split ! True if using the barotropic-baroclinic split algorithm
1907		logical :: calc_tides ! True if using tidal forcing
1908		logical :: calc_sal ! True if using self-attraction and loading
1909		logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure
1910		! This include declares and sets the variable "version".
1911		# include "version_variable.h"
1912		character(len=40) :: mdl = "MOM_diagnostics" ! This module's name.
1913		character(len=48) :: thickness_units, flux_units
1914		logical :: use_temperature, adiabatic
1915		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1916		integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
1917		! for remapping. Values below 20190101 recover the remapping
1918		! answers from 2018, while higher values use more robust
1919		! forms of the same remapping expressions.
1920
1921	1	CS%initialized = .true.
1922
1923	1	CS%diag => diag
1924	1	use_temperature = associated(tv%T)
1925	1	call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., do_not_log=.true.)
1926
1927		! Read all relevant parameters and write them to the model log.
1928	1	call log_version(param_file, mdl, version, "")
1929		call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_COLUMN_FRACTION", CS%mono_N2_column_fraction, &
1930		"The lower fraction of water column over which N2 is limited as monotonic "// &
1931		"for the purposes of calculating the equivalent barotropic wave speed.", &
1932	1	units='nondim', default=0.)
1933		call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_DEPTH", CS%mono_N2_depth, &
1934		"The depth below which N2 is limited as monotonic for the "// &
1935		"purposes of calculating the equivalent barotropic wave speed.", &
1936	1	units='m', scale=GV%m_to_H, default=-1.)
1937		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
1938		"The fractional tolerance for finding the wave speeds.", &
1939	1	units="nondim", default=0.001)
1940		!### Set defaults so that wave_speed_min*wave_speed_tol >= 1e-9 m s-1
1941		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_MIN", wave_speed_min, &
1942		"A floor in the first mode speed below which 0 used instead.", &
1943	1	units="m s-1", default=0.0, scale=US%m_s_to_L_T)
1944		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
1945		"If true, use a more robust estimate of the first mode wave speed as the "//&
1946	1	"starting point for iterations.", default=.true.)
1947		call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
1948	1	do_not_log=.true., default=.true.)
1949
1950		call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
1951		"If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//&
1952		"See REMAPPING_USE_OM4_SUBCELLS for details. "//&
1953	1	"We recommend setting this option to false.", default=om4_remap_via_sub_cells)
1954		call get_param(param_file, mdl, "ACCURATE_NONBOUS_THICK_CELLO", CS%accurate_thick_cello, &
1955		"If true, use the same careful integrals to find the diagnosed non-Boussinesq "//&
1956		"layer thicknesses as are used to find the free surface height, instead of "//&
1957		"using an approximate thickness based on division by the mid-layer density.", &
1958	1	default=.false., do_not_log=GV%Boussinesq)
1959	1	if (GV%Boussinesq) CS%accurate_thick_cello = .false.
1960		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
1961		"This sets the default value for the various _ANSWER_DATE parameters.", &
1962	1	default=99991231)
1963		call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
1964		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
1965		"Values below 20190101 result in the use of older, less accurate expressions "//&
1966		"that were in use at the end of 2018. Higher values result in the use of more "//&
1967		"robust and accurate forms of mathematically equivalent expressions.", &
1968	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
1969	1	if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
1970
1971	1	call get_param(param_file, mdl, "SPLIT", split, default=.true., do_not_log=.true.)
1972	1	call get_param(param_file, mdl, "TIDES", calc_tides, default=.false., do_not_log=.true.)
1973	1	call get_param(param_file, mdl, "CALCULATE_SAL", calc_sal, default=calc_tides, do_not_log=.true.)
1974
1975	1	thickness_units = get_thickness_units(GV)
1976	1	flux_units = get_flux_units(GV)
1977	1	convert_H = GV%H_to_MKS
1978
1979		CS%id_masscello = register_diag_field('ocean_model', 'masscello', diag%axesTL, &
1980		Time, 'Mass per unit area of liquid ocean grid cell', 'kg m-2', conversion=GV%H_to_kg_m2, &
1981	1	standard_name='sea_water_mass_per_unit_area', v_extensive=.true.)
1982
1983		CS%id_masso = register_scalar_field('ocean_model', 'masso', Time, diag, &
1984		'Mass of liquid ocean', units='kg', conversion=US%RZL2_to_kg, &
1985	1	standard_name='sea_water_mass')
1986
1987		CS%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, Time, &
1988		long_name='Cell Thickness', standard_name='cell_thickness', &
1989	1	units='m', conversion=US%Z_to_m, v_extensive=.true.)
1990		CS%id_h_pre_sync = register_diag_field('ocean_model', 'h_pre_sync', diag%axesTL, Time, &
1991		long_name='Cell thickness from the previous timestep', &
1992	1	units=thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
1993
1994		! Note that CS%id_volcello would normally be registered here but because it is a "cell measure" and
1995		! must be registered first. We earlier stored the handle of volcello but need it here for posting
1996		! by this module.
1997	1	CS%id_volcello = diag_get_volume_cell_measure_dm_id(diag)
1998
1999	1	if (use_temperature) then
2000	1	if (tv%T_is_conT) then
2001		CS%id_Tpot = register_diag_field('ocean_model', 'temp', diag%axesTL, &
2002	0	Time, 'Potential Temperature', 'degC', conversion=US%C_to_degC, cmor_field_name="thetao")
2003		endif
2004	1	if (tv%S_is_absS) then
2005		CS%id_Sprac = register_diag_field('ocean_model', 'salt', diag%axesTL, &
2006	0	Time, 'Salinity', 'psu', conversion=US%S_to_ppt, cmor_field_name='so')
2007		endif
2008
2009		CS%id_tob = register_diag_field('ocean_model','tob', diag%axesT1, Time, &
2010		long_name='Sea Water Potential Temperature at Sea Floor', &
2011		standard_name='sea_water_potential_temperature_at_sea_floor', &
2012	1	units='degC', conversion=US%C_to_degC)
2013		CS%id_sob = register_diag_field('ocean_model','sob',diag%axesT1, Time, &
2014		long_name='Sea Water Salinity at Sea Floor', &
2015		standard_name='sea_water_salinity_at_sea_floor', &
2016	1	units='psu', conversion=US%S_to_ppt)
2017
2018		CS%id_tosq = register_diag_field('ocean_model', 'tosq', diag%axesTL, &
2019		Time, 'Square of Potential Temperature', 'degC2', conversion=US%C_to_degC**2, &
2020	1	standard_name='Potential Temperature Squared')
2021		CS%id_sosq = register_diag_field('ocean_model', 'sosq', diag%axesTL, &
2022		Time, 'Square of Salinity', 'psu2', conversion=US%S_to_ppt**2, &
2023	1	standard_name='Salinity Squared')
2024
2025		CS%id_temp_layer_ave = register_diag_field('ocean_model', 'temp_layer_ave', &
2026	1	diag%axesZL, Time, 'Layer Average Ocean Temperature', units='degC', conversion=US%C_to_degC)
2027		CS%id_bigtemp_layer_ave = register_diag_field('ocean_model', 'contemp_layer_ave', &
2028	1	diag%axesZL, Time, 'Layer Average Ocean Conservative Temperature', units='Celsius', conversion=US%C_to_degC)
2029		CS%id_salt_layer_ave = register_diag_field('ocean_model', 'salt_layer_ave', &
2030	1	diag%axesZL, Time, 'Layer Average Ocean Salinity', units='psu', conversion=US%S_to_ppt)
2031		CS%id_abssalt_layer_ave = register_diag_field('ocean_model', 'abssalt_layer_ave', &
2032	1	diag%axesZL, Time, 'Layer Average Ocean Absolute Salinity', units='g kg-1', conversion=US%S_to_ppt)
2033
2034		CS%id_thetaoga = register_scalar_field('ocean_model', 'thetaoga', &
2035		Time, diag, 'Global Mean Ocean Potential Temperature', units='degC', conversion=US%C_to_degC, &
2036	1	standard_name='sea_water_potential_temperature')
2037		CS%id_bigthetaoga = register_scalar_field('ocean_model', 'bigthetaoga', &
2038		Time, diag, 'Global Mean Ocean Conservative Temperature', units='Celsius', conversion=US%C_to_degC, &
2039	1	standard_name='sea_water_conservative_temperature')
2040		CS%id_soga = register_scalar_field('ocean_model', 'soga', &
2041		Time, diag, 'Global Mean Ocean Salinity', units='psu', conversion=US%S_to_ppt, &
2042	1	standard_name='sea_water_salinity')
2043		CS%id_abssoga = register_scalar_field('ocean_model', 'abssoga', &
2044		Time, diag, 'Global Mean Ocean Absolute Salinity', units='g kg-1', conversion=US%S_to_ppt, &
2045	1	standard_name='sea_water_absolute_salinity')
2046
2047		! The CMIP convention is potential temperature, but not indicated in the CMIP long name.
2048		CS%id_tosga = register_scalar_field('ocean_model', 'sst_global', Time, diag, &
2049		long_name='Global Area Average Sea Surface Temperature', &
2050		units='degC', conversion=US%C_to_degC, standard_name='sea_surface_temperature', &
2051		cmor_field_name='tosga', cmor_standard_name='sea_surface_temperature', &
2052	1	cmor_long_name='Sea Surface Temperature')
2053		CS%id_bigtosga = register_scalar_field('ocean_model', 'sscont_global', Time, diag, &
2054		long_name='Global Area Average Sea Surface Conservative Temperature', &
2055	1	units='Celsius', conversion=US%C_to_degC, standard_name='sea_surface_temperature')
2056		! The CMIP convention is practical salinity, but not indicated in the CMIP long name.
2057		CS%id_sosga = register_scalar_field('ocean_model', 'sss_global', Time, diag, &
2058		long_name='Global Area Average Sea Surface Salinity', &
2059		units='psu', conversion=US%S_to_ppt, standard_name='sea_surface_salinity', &
2060		cmor_field_name='sosga', cmor_standard_name='sea_surface_salinity', &
2061	1	cmor_long_name='Sea Surface Salinity')
2062		CS%id_abssosga = register_scalar_field('ocean_model', 'ssabss_global', Time, diag, &
2063		long_name='Global Area Average Sea Surface Absolute Salinity', &
2064	1	units='psu', conversion=US%S_to_ppt, standard_name='sea_surface_absolute_salinity')
2065
2066		! 2d column integrated
2067		CS%id_temp_int = register_diag_field('ocean_model', 'temp_int', diag%axesT1, Time, &
2068		'Density weighted column integrated potential temperature', &
2069		'degC kg m-2', conversion=US%C_to_degC*US%RZ_to_kg_m2, &
2070		cmor_field_name='opottempmint', &
2071		cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_potential_temperature', &
2072	1	cmor_standard_name='Depth integrated density times potential temperature')
2073		CS%id_salt_int = register_diag_field('ocean_model', 'salt_int', diag%axesT1, Time, &
2074		'Density weighted column integrated salinity', &
2075		'psu kg m-2', conversion=US%S_to_ppt*US%RZ_to_kg_m2, v_extensive=.true., &
2076		cmor_field_name='somint', &
2077		cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_salinity', &
2078	1	cmor_standard_name='Depth integrated density times salinity')
2079
2080		! 3d vertically integrated
2081		CS%id_absscint = register_diag_field('ocean_model', 'absscint', diag%axesTL, Time, &
2082		'Integral wrt depth of seawater absolute salinity expressed as salt mass content', &
2083		units='kg m-2', conversion=US%S_to_ppt*US%RZ_to_kg_m2, v_extensive=.true., &
2084	1	standard_name='integral_wrt_depth_of_sea_water_absolute_salinity_expressed_as_salt_mass_content')
2085		CS%id_pfscint = register_diag_field('ocean_model', 'pfscint', diag%axesTL, Time, &
2086		' Integral wrt depth of seawater preformed salinity expressed as salt mass content', &
2087		units='kg m-2', conversion=US%S_to_ppt*US%RZ_to_kg_m2, v_extensive=.true., &
2088	1	standard_name='integral_wrt_depth_of_sea_water_preformed_salinity_expressed_as_salt_mass_content')
2089		CS%id_scint = register_diag_field('ocean_model', 'scint', diag%axesTL, Time, &
2090		'Integral wrt depth of seawater practical salinity expressed as salt mass content', &
2091		units='kg m-2', conversion=US%S_to_ppt*US%RZ_to_kg_m2, v_extensive=.true., &
2092	1	standard_name='integral_wrt_depth_of_sea_water_practical_salinity_expressed_as_salt_mass_content')
2093		CS%id_chcint = register_diag_field('ocean_model', 'chcint', diag%axesTL, Time, &
2094		'Depth Integrated Seawater Conservative Temperature Expressed As Heat Content', &
2095		units='J m-2', conversion=US%Q_to_J_kg*US%RZ_to_kg_m2, v_extensive=.true., &
2096	1	standard_name='integral_wrt_depth_of_sea_water_conservative_temperature_expressed_as_heat_content')
2097		CS%id_phcint = register_diag_field('ocean_model', 'phcint', diag%axesTL, Time, &
2098		'Integrated Ocean Heat Content from Potential Temperature', &
2099		units='J m-2', conversion=US%Q_to_J_kg*US%RZ_to_kg_m2, v_extensive=.true., &
2100	1	standard_name='integral_wrt_depth_of_sea_water_potential_temperature_expressed_as_heat_content')
2101
2102		CS%id_t20d = register_diag_field('ocean_model', 't20d', diag%axesT1, Time, &
2103		'Depth of 20 degree Celsius Isotherm', &
2104		units='m', conversion=US%Z_to_m, &
2105	1	standard_name='depth_of_isosurface_of_sea_water_potential_temperature')
2106		CS%id_t17d = register_diag_field('ocean_model', 't17d', diag%axesT1, Time, &
2107		'Depth of 17 degree Celsius Isotherm', &
2108		units='m', conversion=US%Z_to_m, &
2109	1	standard_name='depth_of_isosurface_of_sea_water_potential_temperature')
2110		endif
2111
2112		CS%id_u = register_diag_field('ocean_model', 'u', diag%axesCuL, Time, &
2113		'Zonal velocity', 'm s-1', conversion=US%L_T_to_m_s, cmor_field_name='uo', &
2114	1	cmor_standard_name='sea_water_x_velocity', cmor_long_name='Sea Water X Velocity')
2115		CS%id_v = register_diag_field('ocean_model', 'v', diag%axesCvL, Time, &
2116		'Meridional velocity', 'm s-1', conversion=US%L_T_to_m_s, cmor_field_name='vo', &
2117	1	cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity')
2118		CS%id_usq = register_diag_field('ocean_model', 'usq', diag%axesCuL, Time, &
2119	1	'Zonal velocity squared', 'm2 s-2', conversion=US%L_T_to_m_s**2)
2120		CS%id_vsq = register_diag_field('ocean_model', 'vsq', diag%axesCvL, Time, &
2121	1	'Meridional velocity squared', 'm2 s-2', conversion=US%L_T_to_m_s**2)
2122		CS%id_uv = register_diag_field('ocean_model', 'uv', diag%axesTL, Time, &
2123		'Product between zonal and meridional velocities at h-points', &
2124	1	'm2 s-2', conversion=US%L_T_to_m_s**2)
2125		CS%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, Time, &
2126	1	'Layer Thickness', thickness_units, v_extensive=.true., conversion=convert_H)
2127
2128		CS%id_e = register_diag_field('ocean_model', 'e', diag%axesTi, Time, &
2129	1	'Interface Height Relative to Mean Sea Level', 'm', conversion=US%Z_to_m)
2130		CS%id_e_D = register_diag_field('ocean_model', 'e_D', diag%axesTi, Time, &
2131	1	'Interface Height above the Seafloor', 'm', conversion=US%Z_to_m)
2132
2133		CS%id_Rml = register_diag_field('ocean_model', 'Rml', diag%axesTL, Time, &
2134	1	'Mixed Layer Coordinate Potential Density', 'kg m-3', conversion=US%R_to_kg_m3)
2135
2136		CS%id_Rcv = register_diag_field('ocean_model', 'Rho_cv', diag%axesTL, Time, &
2137	1	'Coordinate Potential Density', 'kg m-3', conversion=US%R_to_kg_m3)
2138
2139		CS%id_rhopot0 = register_diag_field('ocean_model', 'rhopot0', diag%axesTL, Time, &
2140	1	'Potential density referenced to surface', 'kg m-3', conversion=US%R_to_kg_m3)
2141		CS%id_rhopot2 = register_diag_field('ocean_model', 'rhopot2', diag%axesTL, Time, &
2142	1	'Potential density referenced to 2000 dbar', 'kg m-3', conversion=US%R_to_kg_m3)
2143		CS%id_rhoinsitu = register_diag_field('ocean_model', 'rhoinsitu', diag%axesTL, Time, &
2144	1	'In situ density', 'kg m-3', conversion=US%R_to_kg_m3)
2145		CS%id_drho_dT = register_diag_field('ocean_model', 'drho_dT', diag%axesTL, Time, &
2146		'Partial derivative of rhoinsitu with respect to temperature (alpha)', &
2147	1	'kg m-3 degC-1', conversion=US%R_to_kg_m3*US%degC_to_C)
2148		CS%id_drho_dS = register_diag_field('ocean_model', 'drho_dS', diag%axesTL, Time, &
2149		'Partial derivative of rhoinsitu with respect to salinity (beta)', &
2150	1	'kg^2 g-1 m-3', conversion=US%R_to_kg_m3*US%ppt_to_S)
2151
2152		CS%id_du_dt = register_diag_field('ocean_model', 'dudt', diag%axesCuL, Time, &
2153	1	'Zonal Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
2154
2155		CS%id_dv_dt = register_diag_field('ocean_model', 'dvdt', diag%axesCvL, Time, &
2156	1	'Meridional Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
2157
2158		CS%id_dh_dt = register_diag_field('ocean_model', 'dhdt', diag%axesTL, Time, &
2159	1	'Thickness tendency', trim(thickness_units)//" s-1", conversion=convert_H*US%s_to_T, v_extensive=.true.)
2160
2161		!CS%id_hf_du_dt = register_diag_field('ocean_model', 'hf_dudt', diag%axesCuL, Time, &
2162		! 'Fractional Thickness-weighted Zonal Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2, &
2163		! v_extensive=.true.)
2164
2165		!CS%id_hf_dv_dt = register_diag_field('ocean_model', 'hf_dvdt', diag%axesCvL, Time, &
2166		! 'Fractional Thickness-weighted Meridional Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2, &
2167		! v_extensive=.true.)
2168
2169		CS%id_hf_du_dt_2d = register_diag_field('ocean_model', 'hf_dudt_2d', diag%axesCu1, Time, &
2170	1	'Depth-sum Fractional Thickness-weighted Zonal Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
2171
2172		CS%id_hf_dv_dt_2d = register_diag_field('ocean_model', 'hf_dvdt_2d', diag%axesCv1, Time, &
2173	1	'Depth-sum Fractional Thickness-weighted Meridional Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
2174
2175		CS%id_h_du_dt = register_diag_field('ocean_model', 'h_du_dt', diag%axesCuL, Time, &
2176	1	'Thickness Multiplied Zonal Acceleration', 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2177
2178		CS%id_h_dv_dt = register_diag_field('ocean_model', 'h_dv_dt', diag%axesCvL, Time, &
2179	1	'Thickness Multiplied Meridional Acceleration', 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2180
2181		! layer thickness variables
2182		!if (GV%nk_rho_varies > 0) then
2183		CS%id_h_Rlay = register_diag_field('ocean_model', 'h_rho', diag%axesTL, Time, &
2184		'Layer thicknesses in pure potential density coordinates', &
2185	1	thickness_units, conversion=convert_H)
2186
2187		CS%id_uh_Rlay = register_diag_field('ocean_model', 'uh_rho', diag%axesCuL, Time, &
2188		'Zonal volume transport in pure potential density coordinates', &
2189	1	flux_units, conversion=US%L_to_m*2US%s_to_T*convert_H)
2190
2191		CS%id_vh_Rlay = register_diag_field('ocean_model', 'vh_rho', diag%axesCvL, Time, &
2192		'Meridional volume transport in pure potential density coordinates', &
2193	1	flux_units, conversion=US%L_to_m*2US%s_to_T*convert_H)
2194
2195		CS%id_uhGM_Rlay = register_diag_field('ocean_model', 'uhGM_rho', diag%axesCuL, Time, &
2196		'Zonal volume transport due to interface height diffusion in pure potential '//&
2197	1	'density coordinates', flux_units, conversion=US%L_to_m*2US%s_to_T*convert_H)
2198
2199		CS%id_vhGM_Rlay = register_diag_field('ocean_model', 'vhGM_rho', diag%axesCvL, Time, &
2200		'Meridional volume transport due to interface height diffusion in pure potential '//&
2201	1	'density coordinates', flux_units, conversion=US%L_to_m*2US%s_to_T*convert_H)
2202		!endif
2203
2204
2205		! terms in the kinetic energy budget
2206		CS%id_KE = register_diag_field('ocean_model', 'KE', diag%axesTL, Time, &
2207		'Layer kinetic energy per unit mass', &
2208	1	'm2 s-2', conversion=US%L_T_to_m_s**2)
2209		CS%id_dKEdt = register_diag_field('ocean_model', 'dKE_dt', diag%axesTL, Time, &
2210		'Kinetic Energy Tendency of Layer', &
2211	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2212		CS%id_PE_to_KE = register_diag_field('ocean_model', 'PE_to_KE', diag%axesTL, Time, &
2213		'Potential to Kinetic Energy Conversion of Layer', &
2214	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2215	1	if (calc_sal) &
2216		CS%id_KE_SAL = register_diag_field('ocean_model', 'KE_SAL', diag%axesTL, Time, &
2217		'Kinetic Energy Source from Self-Attraction and Loading', &
2218	0	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2219	1	if (calc_tides) &
2220		CS%id_KE_TIDES = register_diag_field('ocean_model', 'KE_tides', diag%axesTL, Time, &
2221		'Kinetic Energy Source from Astronomical Tidal Forcing', &
2222	0	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2223	1	if (split) then
2224		CS%id_KE_BT = register_diag_field('ocean_model', 'KE_BT', diag%axesTL, Time, &
2225		'Barotropic contribution to Kinetic Energy', &
2226	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2227		CS%id_PE_to_KE_btbc = register_diag_field('ocean_model', 'PE_to_KE_btbc', diag%axesTL, Time, &
2228		'Potential to Kinetic Energy Conversion of Layer (including barotropic solver contribution)', &
2229	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2230		CS%id_KE_Coradv_btbc = register_diag_field('ocean_model', 'KE_Coradv_btbc', diag%axesTL, Time, &
2231		'Kinetic Energy Source from Coriolis and Advection (including barotropic solver contribution)', &
2232	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2233		CS%id_KE_BT_PF = register_diag_field('ocean_model', 'KE_BTPF', diag%axesTL, Time, &
2234		'Kinetic Energy Source from Barotropic Pressure Gradient Force.', &
2235	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2236		CS%id_KE_BT_CF = register_diag_field('ocean_model', 'KE_BTCF', diag%axesTL, Time, &
2237		'Kinetic Energy Source from Barotropic Coriolis Force.', &
2238	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2239		CS%id_KE_BT_WD = register_diag_field('ocean_model', 'KE_BTWD', diag%axesTL, Time, &
2240		'Kinetic Energy Source from Barotropic Linear Wave Drag.', &
2241	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2242		endif
2243		CS%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, Time, &
2244		'Kinetic Energy Source from Coriolis and Advection', &
2245	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2246		CS%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, Time, &
2247		'Kinetic Energy Source from Advection', &
2248	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2249		CS%id_KE_visc = register_diag_field('ocean_model', 'KE_visc', diag%axesTL, Time, &
2250		'Kinetic Energy Source from Vertical Viscosity and Stresses', &
2251	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2252		CS%id_KE_visc_gl90 = register_diag_field('ocean_model', 'KE_visc_gl90', diag%axesTL, Time, &
2253		'Kinetic Energy Source from GL90 Vertical Viscosity', &
2254	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2255		CS%id_KE_stress = register_diag_field('ocean_model', 'KE_stress', diag%axesTL, Time, &
2256		'Kinetic Energy Source from Surface Stresses or Body Wind Stress', &
2257	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2258		CS%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, Time, &
2259		'Kinetic Energy Source from Horizontal Viscosity', &
2260	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2261	1	if (.not. adiabatic) then
2262		CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, &
2263		'Kinetic Energy Source from Diapycnal Diffusion', &
2264	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
2265		endif
2266
2267		! gravity wave CFLs
2268		CS%id_cg1 = register_diag_field('ocean_model', 'cg1', diag%axesT1, Time, &
2269	1	'First baroclinic gravity wave speed', 'm s-1', conversion=US%L_T_to_m_s)
2270		CS%id_Rd1 = register_diag_field('ocean_model', 'Rd1', diag%axesT1, Time, &
2271	1	'First baroclinic deformation radius', 'm', conversion=US%L_to_m)
2272		CS%id_cfl_cg1 = register_diag_field('ocean_model', 'CFL_cg1', diag%axesT1, Time, &
2273	1	'CFL of first baroclinic gravity wave = dtcg1(1/dx+1/dy)', 'nondim')
2274		CS%id_cfl_cg1_x = register_diag_field('ocean_model', 'CFL_cg1_x', diag%axesT1, Time, &
2275	1	'i-component of CFL of first baroclinic gravity wave = dtcg1/dx', 'nondim')
2276		CS%id_cfl_cg1_y = register_diag_field('ocean_model', 'CFL_cg1_y', diag%axesT1, Time, &
2277	1	'j-component of CFL of first baroclinic gravity wave = dtcg1/dy', 'nondim')
2278		CS%id_cg_ebt = register_diag_field('ocean_model', 'cg_ebt', diag%axesT1, Time, &
2279	1	'Equivalent barotropic gravity wave speed', 'm s-1', conversion=US%L_T_to_m_s)
2280		CS%id_Rd_ebt = register_diag_field('ocean_model', 'Rd_ebt', diag%axesT1, Time, &
2281	1	'Equivalent barotropic deformation radius', 'm', conversion=US%L_to_m)
2282		CS%id_p_ebt = register_diag_field('ocean_model', 'p_ebt', diag%axesTL, Time, &
2283	1	'Equivalent barotropic modal strcuture', 'nondim')
2284
2285		if ((CS%id_cg1>0) .or. (CS%id_Rd1>0) .or. (CS%id_cfl_cg1>0) .or. &
2286		(CS%id_cfl_cg1_x>0) .or. (CS%id_cfl_cg1_y>0) .or. &
2287	1	(CS%id_cg_ebt>0) .or. (CS%id_Rd_ebt>0) .or. (CS%id_p_ebt>0)) then
2288		call wave_speed_init(CS%wave_speed, GV, remap_answer_date=remap_answer_date, &
2289		better_speed_est=better_speed_est, min_speed=wave_speed_min, &
2290	1	wave_speed_tol=wave_speed_tol, om4_remap_via_sub_cells=om4_remap_via_sub_cells)
2291		endif
2292
2293		CS%id_mass_wt = register_diag_field('ocean_model', 'mass_wt', diag%axesT1, Time, &
2294	1	'The column mass for calculating mass-weighted average properties', 'kg m-2', conversion=US%RZ_to_kg_m2)
2295
2296		CS%id_col_mass = register_diag_field('ocean_model', 'col_mass', diag%axesT1, Time, &
2297	1	'The column integrated in situ density', 'kg m-2', conversion=US%RZ_to_kg_m2)
2298
2299		CS%id_col_ht = register_diag_field('ocean_model', 'col_height', diag%axesT1, Time, &
2300	1	'The height of the water column', 'm', conversion=US%Z_to_m)
2301		CS%id_pbo = register_diag_field('ocean_model', 'pbo', diag%axesT1, Time, &
2302		long_name='Sea Water Pressure at Sea Floor', standard_name='sea_water_pressure_at_sea_floor', &
2303	1	units='Pa', conversion=US%RL2_T2_to_Pa)
2304
2305		! Register time derivatives and allocate memory for diagnostics that need
2306		! access from across several modules.
2307	1	call set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
2308
2309	1	end subroutine MOM_diagnostics_init
2310
2311
2312		!> Register diagnostics of the surface state and integrated quantities
2313	1	subroutine register_surface_diags(Time, G, US, IDs, diag, tv)
2314		type(time_type), intent(in) :: Time !< current model time
2315		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2316		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2317		type(surface_diag_IDs), intent(inout) :: IDs !< A structure with the diagnostic IDs.
2318		type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
2319		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
2320
2321		! Vertically integrated, budget, and surface state diagnostics
2322		IDs%id_volo = register_scalar_field('ocean_model', 'volo', Time, diag, &
2323		long_name='Total volume of liquid ocean', units='m3', conversion=US%Z_to_mUS%L_to_m*2, &
2324	1	standard_name='sea_water_volume')
2325		IDs%id_zos = register_diag_field('ocean_model', 'zos', diag%axesT1, Time, &
2326		standard_name = 'sea_surface_height_above_geoid', &
2327	1	long_name= 'Sea surface height above geoid', units='m', conversion=US%Z_to_m)
2328		IDs%id_zossq = register_diag_field('ocean_model', 'zossq', diag%axesT1, Time, &
2329		standard_name='square_of_sea_surface_height_above_geoid', &
2330	1	long_name='Square of sea surface height above geoid', units='m2', conversion=US%Z_to_m**2)
2331		IDs%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, Time, &
2332	1	'Sea Surface Height', 'm', conversion=US%Z_to_m)
2333		IDs%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', Time, diag, &
2334		long_name='Area averaged sea surface height', units='m', conversion=US%Z_to_m, &
2335	1	standard_name='area_averaged_sea_surface_height')
2336		IDs%id_ssu = register_diag_field('ocean_model', 'SSU', diag%axesCu1, Time, &
2337	1	'Sea Surface Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s)
2338		IDs%id_ssv = register_diag_field('ocean_model', 'SSV', diag%axesCv1, Time, &
2339	1	'Sea Surface Meridional Velocity', 'm s-1', conversion=US%L_T_to_m_s)
2340		IDs%id_speed = register_diag_field('ocean_model', 'speed', diag%axesT1, Time, &
2341	1	'Sea Surface Speed', 'm s-1', conversion=US%L_T_to_m_s)
2342		IDs%id_ssu_east = register_diag_field('ocean_model', 'ssu_east', diag%axesT1, Time, &
2343	1	'Eastward velocity', 'm s-1', conversion=US%L_T_to_m_s)
2344		IDs%id_ssv_north = register_diag_field('ocean_model', 'ssv_north', diag%axesT1, Time, &
2345	1	'Northward velocity', 'm s-1', conversion=US%L_T_to_m_s)
2346
2347	1	if (associated(tv%T)) then
2348		IDs%id_sst = register_diag_field('ocean_model', 'SST', diag%axesT1, Time, &
2349		'Sea Surface Temperature', 'degC', conversion=US%C_to_degC, &
2350		cmor_field_name='tos', cmor_long_name='Sea Surface Temperature', &
2351	1	cmor_standard_name='sea_surface_temperature')
2352		IDs%id_sst_sq = register_diag_field('ocean_model', 'SST_sq', diag%axesT1, Time, &
2353		'Sea Surface Temperature Squared', 'degC2', conversion=US%C_to_degC**2, &
2354		cmor_field_name='tossq', cmor_long_name='Square of Sea Surface Temperature ', &
2355	1	cmor_standard_name='square_of_sea_surface_temperature')
2356		IDs%id_sss = register_diag_field('ocean_model', 'SSS', diag%axesT1, Time, &
2357		'Sea Surface Salinity', 'psu', conversion=US%S_to_ppt, &
2358		cmor_field_name='sos', cmor_long_name='Sea Surface Salinity', &
2359	1	cmor_standard_name='sea_surface_salinity')
2360		IDs%id_sss_sq = register_diag_field('ocean_model', 'SSS_sq', diag%axesT1, Time, &
2361		'Sea Surface Salinity Squared', 'psu2', conversion=US%S_to_ppt**2, &
2362		cmor_field_name='sossq', cmor_long_name='Square of Sea Surface Salinity ', &
2363	1	cmor_standard_name='square_of_sea_surface_salinity')
2364	1	if (tv%T_is_conT) then
2365		IDs%id_sstcon = register_diag_field('ocean_model', 'conSST', diag%axesT1, Time, &
2366	0	'Sea Surface Conservative Temperature', 'Celsius', conversion=US%C_to_degC)
2367		endif
2368	1	if (tv%S_is_absS) then
2369		IDs%id_sssabs = register_diag_field('ocean_model', 'absSSS', diag%axesT1, Time, &
2370	0	'Sea Surface Absolute Salinity', 'g kg-1', conversion=US%S_to_ppt)
2371		endif
2372	1	if (associated(tv%frazil)) then
2373		IDs%id_fraz = register_diag_field('ocean_model', 'frazil', diag%axesT1, Time, &
2374		'Heat from frazil formation', 'W m-2', conversion=US%QRZ_T_to_W_m2, &
2375		cmor_field_name='hfsifrazil', &
2376		cmor_standard_name='heat_flux_into_sea_water_due_to_frazil_ice_formation', &
2377	1	cmor_long_name='Heat Flux into Sea Water due to Frazil Ice Formation')
2378		endif
2379		endif
2380
2381		IDs%id_salt_deficit = register_diag_field('ocean_model', 'salt_deficit', diag%axesT1, Time, &
2382		'Salt source in ocean required to supply excessive ice salt fluxes', &
2383	1	'ppt kg m-2 s-1', conversion=US%S_to_ppt*US%RZ_T_to_kg_m2s)
2384		IDs%id_Heat_PmE = register_diag_field('ocean_model', 'Heat_PmE', diag%axesT1, Time, &
2385		'Heat flux into ocean from mass flux into ocean', &
2386	1	'W m-2', conversion=US%QRZ_T_to_W_m2)
2387		IDs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, Time, &
2388		'Heat flux into ocean from geothermal or other internal sources', &
2389	1	'W m-2', conversion=US%QRZ_T_to_W_m2)
2390
2391	1	end subroutine register_surface_diags
2392
2393		!> Register certain diagnostics related to transports
2394	1	subroutine register_transport_diags(Time, G, GV, US, IDs, diag)
2395		type(time_type), intent(in) :: Time !< current model time
2396		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2397		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2398		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2399		type(transport_diag_IDs), intent(inout) :: IDs !< A structure with the diagnostic IDs.
2400		type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
2401
2402		character(len=48) :: thickness_units, accum_flux_units
2403
2404	1	thickness_units = get_thickness_units(GV)
2405	1	if (GV%Boussinesq) then
2406	1	accum_flux_units = "m3"
2407		else
2408	0	accum_flux_units = "kg"
2409		endif
2410
2411		! Diagnostics related to tracer and mass transport
2412		IDs%id_uhtr = register_diag_field('ocean_model', 'uhtr', diag%axesCuL, Time, &
2413		'Accumulated zonal thickness fluxes to advect tracers', &
2414	1	accum_flux_units, y_cell_method='sum', v_extensive=.true., conversion=GV%H_to_MKSUS%L_to_m*2)
2415		IDs%id_vhtr = register_diag_field('ocean_model', 'vhtr', diag%axesCvL, Time, &
2416		'Accumulated meridional thickness fluxes to advect tracers', &
2417	1	accum_flux_units, x_cell_method='sum', v_extensive=.true., conversion=GV%H_to_MKSUS%L_to_m*2)
2418		IDs%id_umo = register_diag_field('ocean_model', 'umo', &
2419		diag%axesCuL, Time, 'Ocean Mass X Transport', &
2420		'kg s-1', conversion=US%RZ_T_to_kg_m2sUS%L_to_m*2, &
2421	1	standard_name='ocean_mass_x_transport', y_cell_method='sum', v_extensive=.true.)
2422		IDs%id_vmo = register_diag_field('ocean_model', 'vmo', &
2423		diag%axesCvL, Time, 'Ocean Mass Y Transport', &
2424		'kg s-1', conversion=US%RZ_T_to_kg_m2sUS%L_to_m*2, &
2425	1	standard_name='ocean_mass_y_transport', x_cell_method='sum', v_extensive=.true.)
2426		IDs%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', &
2427		diag%axesCu1, Time, 'Ocean Mass X Transport Vertical Sum', &
2428		'kg s-1', conversion=US%RZ_T_to_kg_m2sUS%L_to_m*2, &
2429	1	standard_name='ocean_mass_x_transport_vertical_sum', y_cell_method='sum')
2430		IDs%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d', &
2431		diag%axesCv1, Time, 'Ocean Mass Y Transport Vertical Sum', &
2432		'kg s-1', conversion=US%RZ_T_to_kg_m2sUS%L_to_m*2, &
2433	1	standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum')
2434		IDs%id_dynamics_h = register_diag_field('ocean_model','dynamics_h', &
2435		diag%axesTl, Time, 'Layer thicknesses prior to horizontal dynamics', &
2436	1	thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
2437		IDs%id_dynamics_h_tendency = register_diag_field('ocean_model','dynamics_h_tendency', &
2438		diag%axesTl, Time, 'Change in layer thicknesses due to horizontal dynamics', &
2439	1	trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)
2440
2441	1	end subroutine register_transport_diags
2442
2443		!> Offers the static fields in the ocean grid type for output via the diag_manager.
2444	1	subroutine write_static_fields(G, GV, US, tv, diag)
2445		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2446		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2447		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2448		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
2449		type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
2450
2451		! Local variables
2452	2	real :: work_2d(SZI_(G),SZJ_(G)) ! A 2-d temporary work array [Z ~> m]
2453		integer :: id, i, j
2454		logical :: use_temperature
2455
2456		id = register_static_field('ocean_model', 'geolat', diag%axesT1, &
2457	1	'Latitude of tracer (T) points', 'degrees_north')
2458	1	if (id > 0) call post_data(id, G%geoLatT, diag, .true.)
2459
2460		id = register_static_field('ocean_model', 'geolon', diag%axesT1, &
2461	1	'Longitude of tracer (T) points', 'degrees_east')
2462	1	if (id > 0) call post_data(id, G%geoLonT, diag, .true.)
2463
2464		id = register_static_field('ocean_model', 'geolat_c', diag%axesB1, &
2465	1	'Latitude of corner (Bu) points', 'degrees_north', interp_method='none')
2466	1	if (id > 0) call post_data(id, G%geoLatBu, diag, .true.)
2467
2468		id = register_static_field('ocean_model', 'geolon_c', diag%axesB1, &
2469	1	'Longitude of corner (Bu) points', 'degrees_east', interp_method='none')
2470	1	if (id > 0) call post_data(id, G%geoLonBu, diag, .true.)
2471
2472		id = register_static_field('ocean_model', 'geolat_v', diag%axesCv1, &
2473	1	'Latitude of meridional velocity (Cv) points', 'degrees_north', interp_method='none')
2474	1	if (id > 0) call post_data(id, G%geoLatCv, diag, .true.)
2475
2476		id = register_static_field('ocean_model', 'geolon_v', diag%axesCv1, &
2477	1	'Longitude of meridional velocity (Cv) points', 'degrees_east', interp_method='none')
2478	1	if (id > 0) call post_data(id, G%geoLonCv, diag, .true.)
2479
2480		id = register_static_field('ocean_model', 'geolat_u', diag%axesCu1, &
2481	1	'Latitude of zonal velocity (Cu) points', 'degrees_north', interp_method='none')
2482	1	if (id > 0) call post_data(id, G%geoLatCu, diag, .true.)
2483
2484		id = register_static_field('ocean_model', 'geolon_u', diag%axesCu1, &
2485	1	'Longitude of zonal velocity (Cu) points', 'degrees_east', interp_method='none')
2486	1	if (id > 0) call post_data(id, G%geoLonCu, diag, .true.)
2487
2488		id = register_static_field('ocean_model', 'area_t', diag%axesT1, &
2489		'Surface area of tracer (T) cells', 'm2', conversion=US%L_to_m**2, &
2490		cmor_field_name='areacello', cmor_standard_name='cell_area', &
2491		cmor_long_name='Ocean Grid-Cell Area', &
2492	1	x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
2493	1	if (id > 0) then
2494	1	call post_data(id, G%areaT, diag, .true.)
2495	1	call diag_register_area_ids(diag, id_area_t=id)
2496		endif
2497
2498		id = register_static_field('ocean_model', 'area_u', diag%axesCu1, &
2499		'Surface area of x-direction flow (U) cells', 'm2', conversion=US%L_to_m**2, &
2500		cmor_field_name='areacello_cu', cmor_standard_name='cell_area', &
2501		cmor_long_name='Ocean Grid-Cell Area', &
2502	1	x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
2503	1	if (id > 0) call post_data(id, G%areaCu, diag, .true.)
2504
2505		id = register_static_field('ocean_model', 'area_v', diag%axesCv1, &
2506		'Surface area of y-direction flow (V) cells', 'm2', conversion=US%L_to_m**2, &
2507		cmor_field_name='areacello_cv', cmor_standard_name='cell_area', &
2508		cmor_long_name='Ocean Grid-Cell Area', &
2509	1	x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
2510	1	if (id > 0) call post_data(id, G%areaCv, diag, .true.)
2511
2512		id = register_static_field('ocean_model', 'area_q', diag%axesB1, &
2513		'Surface area of B-grid flow (Q) cells', 'm2', conversion=US%L_to_m**2, &
2514		cmor_field_name='areacello_bu', cmor_standard_name='cell_area', &
2515		cmor_long_name='Ocean Grid-Cell Area', &
2516	1	x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
2517	1	if (id > 0) call post_data(id, G%areaBu, diag, .true.)
2518
2519		id = register_static_field('ocean_model', 'depth_ocean', diag%axesT1, &
2520		'Depth of the ocean at tracer points', 'm', conversion=US%Z_to_m, &
2521		standard_name='sea_floor_depth_below_geoid', &
2522		cmor_field_name='deptho', cmor_long_name='Sea Floor Depth', &
2523		cmor_standard_name='sea_floor_depth_below_geoid', area=diag%axesT1%id_area, &
2524	1	x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
2525	1	if (id > 0) then
2526	7261	do j=G%jsc,G%jec ; do i=G%isc,G%iec ; work_2d(i,j) = G%bathyT(i,j)+G%Z_ref ; enddo ; enddo
2527		! A mask argument is required here because masks are not applied to static fields by default.
2528	1	call post_data(id, work_2d, diag, .true., mask=G%mask2dT)
2529		endif
2530
2531		id = register_static_field('ocean_model', 'wet', diag%axesT1, &
2532	1	'0 if land, 1 if ocean at tracer points', 'none', area=diag%axesT1%id_area)
2533	1	if (id > 0) call post_data(id, G%mask2dT, diag, .true.)
2534
2535		id = register_static_field('ocean_model', 'wet_c', diag%axesB1, &
2536	1	'0 if land, 1 if ocean at corner (Bu) points', 'none', interp_method='none')
2537	1	if (id > 0) call post_data(id, G%mask2dBu, diag, .true.)
2538
2539		id = register_static_field('ocean_model', 'wet_u', diag%axesCu1, &
2540	1	'0 if land, 1 if ocean at zonal velocity (Cu) points', 'none', interp_method='none')
2541	1	if (id > 0) call post_data(id, G%mask2dCu, diag, .true.)
2542
2543		id = register_static_field('ocean_model', 'wet_v', diag%axesCv1, &
2544	1	'0 if land, 1 if ocean at meridional velocity (Cv) points', 'none', interp_method='none')
2545	1	if (id > 0) call post_data(id, G%mask2dCv, diag, .true.)
2546
2547		id = register_static_field('ocean_model', 'Coriolis', diag%axesB1, &
2548	1	'Coriolis parameter at corner (Bu) points', 's-1', interp_method='none', conversion=US%s_to_T)
2549	1	if (id > 0) call post_data(id, G%CoriolisBu, diag, .true.)
2550
2551		id = register_static_field('ocean_model', 'dxt', diag%axesT1, &
2552	1	'Delta(x) at thickness/tracer points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2553	1	if (id > 0) call post_data(id, G%dxT, diag, .true.)
2554
2555		id = register_static_field('ocean_model', 'dyt', diag%axesT1, &
2556	1	'Delta(y) at thickness/tracer points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2557	1	if (id > 0) call post_data(id, G%dyT, diag, .true.)
2558
2559		id = register_static_field('ocean_model', 'dxCu', diag%axesCu1, &
2560	1	'Delta(x) at u points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2561	1	if (id > 0) call post_data(id, G%dxCu, diag, .true.)
2562
2563		id = register_static_field('ocean_model', 'dyCu', diag%axesCu1, &
2564	1	'Delta(y) at u points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2565	1	if (id > 0) call post_data(id, G%dyCu, diag, .true.)
2566
2567		id = register_static_field('ocean_model', 'dxCv', diag%axesCv1, &
2568	1	'Delta(x) at v points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2569	1	if (id > 0) call post_data(id, G%dxCv, diag, .true.)
2570
2571		id = register_static_field('ocean_model', 'dyCv', diag%axesCv1, &
2572	1	'Delta(y) at v points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2573	1	if (id > 0) call post_data(id, G%dyCv, diag, .true.)
2574
2575		id = register_static_field('ocean_model', 'dyCuo', diag%axesCu1, &
2576	1	'Open meridional grid spacing at u points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2577	1	if (id > 0) call post_data(id, G%dy_Cu, diag, .true.)
2578
2579		id = register_static_field('ocean_model', 'dxCvo', diag%axesCv1, &
2580	1	'Open zonal grid spacing at v points (meter)', 'm', interp_method='none', conversion=US%L_to_m)
2581	1	if (id > 0) call post_data(id, G%dx_Cv, diag, .true.)
2582
2583		id = register_static_field('ocean_model', 'sin_rot', diag%axesT1, &
2584	1	'sine of the clockwise angle of the ocean grid north to true north', 'none')
2585	1	if (id > 0) call post_data(id, G%sin_rot, diag, .true.)
2586
2587		id = register_static_field('ocean_model', 'cos_rot', diag%axesT1, &
2588	1	'cosine of the clockwise angle of the ocean grid north to true north', 'none')
2589	1	if (id > 0) call post_data(id, G%cos_rot, diag, .true.)
2590
2591
2592		! This static diagnostic is from CF 1.8, and is the fraction of a cell
2593		! covered by ocean, given as a percentage (poorly named).
2594		id = register_static_field('ocean_model', 'area_t_percent', diag%axesT1, &
2595		'Percentage of cell area covered by ocean', '%', conversion=100.0, &
2596		cmor_field_name='sftof', cmor_standard_name='SeaAreaFraction', &
2597		cmor_long_name='Sea Area Fraction', &
2598	1	x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
2599	1	if (id > 0) call post_data(id, G%mask2dT, diag, .true.)
2600
2601
2602		id = register_static_field('ocean_model','Rho_0', diag%axesNull, &
2603		'mean ocean density used with the Boussinesq approximation', &
2604		'kg m-3', conversion=US%R_to_kg_m3, cmor_field_name='rhozero', &
2605		cmor_standard_name='reference_sea_water_density_for_boussinesq_approximation', &
2606	1	cmor_long_name='reference sea water density for boussinesq approximation')
2607	1	if (id > 0) call post_data(id, GV%Rho0, diag, .true.)
2608
2609	1	use_temperature = associated(tv%T)
2610	1	if (use_temperature) then
2611		id = register_static_field('ocean_model','C_p', diag%axesNull, &
2612		'heat capacity of sea water', 'J kg-1 K-1', conversion=US%Q_to_J_kg*US%degC_to_C, &
2613		cmor_field_name='cpocean', &
2614		cmor_standard_name='specific_heat_capacity_of_sea_water', &
2615	1	cmor_long_name='specific_heat_capacity_of_sea_water')
2616	1	if (id > 0) call post_data(id, tv%C_p, diag, .true.)
2617		endif
2618
2619	1	end subroutine write_static_fields
2620
2621
2622		!> This subroutine sets up diagnostics upon which other diagnostics depend.
2623	1	subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
2624		type(ocean_internal_state), intent(in) :: MIS !< For "MOM Internal State" a set of pointers to
2625		!! the fields and accelerations making up ocean
2626		!! internal physical state.
2627		type(accel_diag_ptrs), intent(inout) :: ADp !< Structure pointing to accelerations in
2628		!! momentum equation.
2629		type(cont_diag_ptrs), intent(inout) :: CDp !< Structure pointing to terms in continuity
2630		!! equation.
2631		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2632		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2633		type(diagnostics_CS), intent(inout) :: CS !< Pointer to the control structure for this
2634		!! module.
2635
2636		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
2637	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
2638	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
2639
2640		! Allocate and register time derivatives.
2641		if ( ( (CS%id_du_dt>0) .or. (CS%id_dKEdt > 0) .or. &
2642		! (CS%id_hf_du_dt > 0) .or. &
2643	1	(CS%id_h_du_dt > 0) .or. (CS%id_hf_du_dt_2d > 0) ) .and. &
2644		(.not. allocated(CS%du_dt)) ) then
2645	0	allocate(CS%du_dt(IsdB:IedB,jsd:jed,nz), source=0.)
2646	0	call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
2647		endif
2648		if ( ( (CS%id_dv_dt>0) .or. (CS%id_dKEdt > 0) .or. &
2649		! (CS%id_hf_dv_dt > 0) .or. &
2650	1	(CS%id_h_dv_dt > 0) .or. (CS%id_hf_dv_dt_2d > 0) ) .and. &
2651		(.not. allocated(CS%dv_dt)) ) then
2652	0	allocate(CS%dv_dt(isd:ied,JsdB:JedB,nz), source=0.)
2653	0	call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
2654		endif
2655	1	if ( ( (CS%id_dh_dt>0) .or. (CS%id_dKEdt > 0) ) .and. &
2656		(.not. allocated(CS%dh_dt)) ) then
2657	0	allocate(CS%dh_dt(isd:ied,jsd:jed,nz), source=0.)
2658	0	call register_time_deriv(lbound(MIS%h), MIS%h, CS%dh_dt, CS)
2659		endif
2660
2661		! Allocate memory for other dependent diagnostics.
2662	1	if (CS%id_KE_adv > 0) then
2663	0	call safe_alloc_ptr(ADp%gradKEu,IsdB,IedB,jsd,jed,nz)
2664	0	call safe_alloc_ptr(ADp%gradKEv,isd,ied,JsdB,JedB,nz)
2665		endif
2666	1	if (CS%id_KE_visc > 0) then
2667	0	call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
2668	0	call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
2669		endif
2670	1	if (CS%id_KE_visc_gl90 > 0) then
2671	0	call safe_alloc_ptr(ADp%du_dt_visc_gl90,IsdB,IedB,jsd,jed,nz)
2672	0	call safe_alloc_ptr(ADp%dv_dt_visc_gl90,isd,ied,JsdB,JedB,nz)
2673		endif
2674	1	if (CS%id_KE_stress > 0) then
2675	0	call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz)
2676	0	call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz)
2677		endif
2678
2679	1	if (CS%id_KE_dia > 0) then
2680	0	call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz)
2681	0	call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz)
2682	0	call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1)
2683		endif
2684
2685	1	if ((CS%id_PE_to_KE_btbc > 0) .or. (CS%id_KE_BT_PF > 0)) then
2686	0	call safe_alloc_ptr(ADp%bt_pgf_u, IsdB, IedB, jsd, jed, nz)
2687	0	call safe_alloc_ptr(ADp%bt_pgf_v, isd, ied, JsdB, JedB, nz)
2688		endif
2689
2690	1	if ((CS%id_KE_Coradv_btbc > 0) .or. (CS%id_KE_BT_CF > 0)) then
2691	0	call safe_alloc_ptr(ADp%bt_cor_u, IsdB, IedB, jsd, jed)
2692	0	call safe_alloc_ptr(ADp%bt_cor_v, isd, ied, JsdB, JedB)
2693		endif
2694
2695	1	if (CS%id_KE_BT_WD > 0) then
2696	0	call safe_alloc_ptr(ADp%bt_lwd_u, IsdB, IedB, jsd, jed)
2697	0	call safe_alloc_ptr(ADp%bt_lwd_v, isd, ied, JsdB, JedB)
2698		endif
2699
2700	1	if (CS%id_KE_SAL > 0) then
2701	0	call safe_alloc_ptr(ADp%sal_u, IsdB, IedB, jsd, jed, nz)
2702	0	call safe_alloc_ptr(ADp%sal_v, isd, ied, JsdB, JedB, nz)
2703		endif
2704
2705	1	if (CS%id_KE_TIDES > 0) then
2706	0	call safe_alloc_ptr(ADp%tides_u, IsdB, IedB, jsd, jed, nz)
2707	0	call safe_alloc_ptr(ADp%tides_v, isd, ied, JsdB, JedB, nz)
2708		endif
2709
2710		CS%KE_term_on = ((CS%id_dKEdt > 0) .or. (CS%id_PE_to_KE > 0) .or. (CS%id_KE_BT > 0) .or. &
2711		(CS%id_KE_Coradv > 0) .or. (CS%id_KE_adv > 0) .or. (CS%id_KE_visc > 0) .or. &
2712		(CS%id_KE_visc_gl90 > 0) .or. (CS%id_KE_stress > 0) .or. (CS%id_KE_horvisc > 0) .or. &
2713		(CS%id_KE_dia > 0) .or. (CS%id_PE_to_KE_btbc > 0) .or. (CS%id_KE_BT_PF > 0) .or. &
2714		(CS%id_KE_Coradv_btbc > 0) .or. (CS%id_KE_BT_CF > 0) .or. (CS%id_KE_BT_WD > 0) .or. &
2715	1	(CS%id_KE_SAL > 0) .or. (CS%id_KE_TIDES > 0))
2716
2717	1	if (CS%id_h_du_dt > 0) call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
2718	1	if (CS%id_h_dv_dt > 0) call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
2719
2720	1	if (CS%id_hf_du_dt_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
2721	1	if (CS%id_hf_dv_dt_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
2722		! if (CS%id_hf_du_dt > 0) call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
2723		! if (CS%id_hf_dv_dt > 0) call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
2724
2725	1	if (CS%id_uhGM_Rlay > 0) call safe_alloc_ptr(CDp%uhGM,IsdB,IedB,jsd,jed,nz)
2726	1	if (CS%id_vhGM_Rlay > 0) call safe_alloc_ptr(CDp%vhGM,isd,ied,JsdB,JedB,nz)
2727
2728	1	end subroutine set_dependent_diagnostics
2729
2730		!> Deallocate memory associated with the diagnostics module
2731	1	subroutine MOM_diagnostics_end(CS, ADp, CDp)
2732		type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by a
2733		!! previous call to diagnostics_init.
2734		type(accel_diag_ptrs), intent(inout) :: ADp !< structure with pointers to
2735		!! accelerations in momentum equation.
2736		type(cont_diag_ptrs), intent(inout) :: CDp !< Structure pointing to terms in continuity
2737		!! equation.
2738		integer :: m
2739
2740	1	if (allocated(CS%dh_dt)) deallocate(CS%dh_dt)
2741	1	if (allocated(CS%dv_dt)) deallocate(CS%dv_dt)
2742	1	if (allocated(CS%du_dt)) deallocate(CS%du_dt)
2743
2744	1	if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu)
2745	1	if (associated(ADp%gradKEv)) deallocate(ADp%gradKEv)
2746	1	if (associated(ADp%du_dt_visc)) deallocate(ADp%du_dt_visc)
2747	1	if (associated(ADp%dv_dt_visc)) deallocate(ADp%dv_dt_visc)
2748	1	if (associated(ADp%du_dt_str)) deallocate(ADp%du_dt_str)
2749	1	if (associated(ADp%dv_dt_str)) deallocate(ADp%dv_dt_str)
2750	1	if (associated(ADp%du_dt_dia)) deallocate(ADp%du_dt_dia)
2751	1	if (associated(ADp%dv_dt_dia)) deallocate(ADp%dv_dt_dia)
2752	1	if (associated(ADp%du_other)) deallocate(ADp%du_other)
2753	1	if (associated(ADp%dv_other)) deallocate(ADp%dv_other)
2754
2755	1	if (associated(ADp%bt_pgf_u)) deallocate(ADp%bt_pgf_u)
2756	1	if (associated(ADp%bt_pgf_v)) deallocate(ADp%bt_pgf_v)
2757	1	if (associated(ADp%bt_cor_u)) deallocate(ADp%bt_cor_u)
2758	1	if (associated(ADp%bt_cor_v)) deallocate(ADp%bt_cor_v)
2759	1	if (associated(ADp%bt_lwd_u)) deallocate(ADp%bt_lwd_u)
2760	1	if (associated(ADp%bt_lwd_v)) deallocate(ADp%bt_lwd_v)
2761
2762		! NOTE: sal_[uv] and tide_[uv] may be allocated either here (KE budget diagnostics) or
2763		! PressureForce module (momentum acceleration diagnostics)
2764	1	if (associated(ADp%sal_u)) deallocate(ADp%sal_u)
2765	1	if (associated(ADp%sal_v)) deallocate(ADp%sal_v)
2766	1	if (associated(ADp%tides_u)) deallocate(ADp%tides_u)
2767	1	if (associated(ADp%tides_v)) deallocate(ADp%tides_v)
2768
2769	1	if (associated(ADp%diag_hfrac_u)) deallocate(ADp%diag_hfrac_u)
2770	1	if (associated(ADp%diag_hfrac_v)) deallocate(ADp%diag_hfrac_v)
2771
2772		! NOTE: [uv]hGM may be allocated either here or the thickness diffuse module
2773	1	if (associated(CDp%uhGM)) deallocate(CDp%uhGM)
2774	1	if (associated(CDp%vhGM)) deallocate(CDp%vhGM)
2775	1	if (associated(CDp%diapyc_vel)) deallocate(CDp%diapyc_vel)
2776
2777	1	do m=1,CS%num_time_deriv ; deallocate(CS%prev_val(m)%p) ; enddo
2778	1	end subroutine MOM_diagnostics_end
2779
2780	0	end module MOM_diagnostics