MOM_diagnostics.F90

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