MOM_tracer_registry.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!> This module contains subroutines that handle registration of tracers
6!! and related subroutines. The primary subroutine, register_tracer, is
7!! called to indicate the tracers advected and diffused.
8!! It also makes public the types defined in MOM_tracer_types.
10
11! use MOM_diag_mediator, only : diag_ctrl
12use mom_coms, only : reproducing_sum
13use mom_debugging, only : hchksum
14use mom_diag_mediator, only : diag_ctrl, register_diag_field, post_data, safe_alloc_ptr
17use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
18use mom_file_parser, only : get_param, log_version, param_file_type
20use mom_grid, only : ocean_grid_type
21use mom_interface_heights, only : thickness_to_dz
23use mom_restart, only : register_restart_field, mom_restart_cs
25use mom_time_manager, only : time_type
30
31implicit none ; private
32
33#include <MOM_memory.h>
34
35public register_tracer
36public mom_tracer_chksum, mom_tracer_chkinv
44
45!> Write out checksums for registered tracers
46interface mom_tracer_chksum
48end interface mom_tracer_chksum
49
50!> Calculate and print the global inventories of registered tracers
51interface mom_tracer_chkinv
53end interface mom_tracer_chkinv
54
55contains
56
57!> This subroutine registers a tracer to be advected and laterally diffused.
58subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, units, &
59 cmor_name, cmor_units, cmor_longname, net_surfflux_name, &
60 NLT_budget_name, net_surfflux_longname, tr_desc, OBC_inflow, &
61 OBC_in_u, OBC_in_v, ad_x, ad_y, df_x, df_y, ad_2d_x, ad_2d_y, &
62 df_2d_x, df_2d_y, advection_xy, registry_diags, &
63 conc_scale, flux_nameroot, flux_longname, flux_units, flux_scale, &
64 convergence_units, convergence_scale, cmor_tendprefix, diag_form, &
65 restart_CS, mandatory, underflow_conc, Tr_out, advect_scheme)
66 type(hor_index_type), intent(in) :: hi !< horizontal index type
67 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
68 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
69 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
70 target :: tr_ptr !< target or pointer to the tracer array [CU ~> conc]
71 type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values
72 character(len=*), optional, intent(in) :: name !< Short tracer name
73 character(len=*), optional, intent(in) :: longname !< The long tracer name
74 character(len=*), optional, intent(in) :: units !< The units of this tracer
75 character(len=*), optional, intent(in) :: cmor_name !< CMOR name
76 character(len=*), optional, intent(in) :: cmor_units !< CMOR physical dimensions of variable
77 character(len=*), optional, intent(in) :: cmor_longname !< CMOR long name
78 character(len=*), optional, intent(in) :: net_surfflux_name !< Name for net_surfflux diag
79 character(len=*), optional, intent(in) :: nlt_budget_name !< Name for NLT_budget diag
80 character(len=*), optional, intent(in) :: net_surfflux_longname !< Long name for net_surfflux diag
81 type(vardesc), optional, intent(in) :: tr_desc !< A structure with metadata about the tracer
82
83 real, optional, intent(in) :: obc_inflow !< the tracer for all inflows via OBC for which OBC_in_u
84 !! or OBC_in_v are not specified [CU ~> conc]
85 real, dimension(:,:,:), optional, pointer :: obc_in_u !< tracer at inflows through u-faces of
86 !! tracer cells [CU ~> conc]
87 real, dimension(:,:,:), optional, pointer :: obc_in_v !< tracer at inflows through v-faces of
88 !! tracer cells [CU ~> conc]
89
90 ! The following are probably not necessary if registry_diags is present and true.
91 real, dimension(:,:,:), optional, pointer :: ad_x !< diagnostic x-advective flux
92 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
93 real, dimension(:,:,:), optional, pointer :: ad_y !< diagnostic y-advective flux
94 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
95 real, dimension(:,:,:), optional, pointer :: df_x !< diagnostic x-diffusive flux
96 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
97 real, dimension(:,:,:), optional, pointer :: df_y !< diagnostic y-diffusive flux
98 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
99 real, dimension(:,:), optional, pointer :: ad_2d_x !< vert sum of diagnostic x-advect flux
100 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
101 real, dimension(:,:), optional, pointer :: ad_2d_y !< vert sum of diagnostic y-advect flux
102 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
103 real, dimension(:,:), optional, pointer :: df_2d_x !< vert sum of diagnostic x-diffuse flux
104 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
105 real, dimension(:,:), optional, pointer :: df_2d_y !< vert sum of diagnostic y-diffuse flux
106 !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
107
108 real, dimension(:,:,:), optional, pointer :: advection_xy !< convergence of lateral advective tracer fluxes
109 !! [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1]
110 logical, optional, intent(in) :: registry_diags !< If present and true, use the registry for
111 !! the diagnostics of this tracer.
112 real, optional, intent(in) :: conc_scale !< A scaling factor used to convert the concentration
113 !! of this tracer to its desired units [conc CU-1 ~> 1]
114 character(len=*), optional, intent(in) :: flux_nameroot !< Short tracer name snippet used construct the
115 !! names of flux diagnostics.
116 character(len=*), optional, intent(in) :: flux_longname !< A word or phrase used construct the long
117 !! names of flux diagnostics.
118 character(len=*), optional, intent(in) :: flux_units !< The units for the fluxes of this tracer.
119 real, optional, intent(in) :: flux_scale !< A scaling factor used to convert the fluxes
120 !! of this tracer to its desired units
121 !! [conc m CU-1 H-1 ~> 1] or [conc kg m-2 CU-1 H-1 ~> 1]
122 character(len=*), optional, intent(in) :: convergence_units !< The units for the flux convergence of
123 !! this tracer.
124 real, optional, intent(in) :: convergence_scale !< A scaling factor used to convert the flux
125 !! convergence of this tracer to its desired units.
126 !! [conc m CU-1 H-1 ~> 1] or [conc kg m-2 CU-1 H-1 ~> 1]
127 character(len=*), optional, intent(in) :: cmor_tendprefix !< The CMOR name for the layer-integrated
128 !! tendencies of this tracer.
129 integer, optional, intent(in) :: diag_form !< An integer (1 or 2, 1 by default) indicating the
130 !! character string template to use in
131 !! labeling diagnostics
132 type(mom_restart_cs), optional, intent(inout) :: restart_cs !< MOM restart control struct
133 logical, optional, intent(in) :: mandatory !< If true, this tracer must be read
134 !! from a restart file.
135 real, optional, intent(in) :: underflow_conc !< A tiny concentration, below which the tracer
136 !! concentration underflows to 0 [CU ~> conc].
137 type(tracer_type), optional, pointer :: tr_out !< If present, returns pointer into registry
138
139 integer, optional, intent(in) :: advect_scheme !< Advection scheme for this tracer, the default is -1
140 !! indicating to use the scheme from MOM_tracer_advect
141
142 logical :: mand
143 type(tracer_type), pointer :: tr=>null()
144 character(len=256) :: mesg ! Message for error messages.
145
146 if (.not. associated(reg)) call tracer_registry_init(param_file, reg)
147
148 if (reg%ntr>=max_fields_) then
149 write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I0," to allow for &
150 &all the tracers being registered via register_tracer.")') reg%ntr+1
151 call mom_error(fatal,"MOM register_tracer: "//mesg)
152 endif
153 reg%ntr = reg%ntr + 1
154
155 tr => reg%Tr(reg%ntr)
156 if (present(tr_out)) tr_out => reg%Tr(reg%ntr)
157
158 if (present(name)) then
159 tr%name = name
160 tr%longname = name ; if (present(longname)) tr%longname = longname
161 tr%units = "Conc" ; if (present(units)) tr%units = units
162
163 tr%cmor_name = ""
164 if (present(cmor_name)) tr%cmor_name = cmor_name
165
166 tr%cmor_units = tr%units
167 if (present(cmor_units)) tr%cmor_units = cmor_units
168
169 tr%cmor_longname = ""
170 if (present(cmor_longname)) tr%cmor_longname = cmor_longname
171
172 if (present(tr_desc)) call mom_error(warning, "MOM register_tracer: "//&
173 "It is a bad idea to use both name and tr_desc when registring "//trim(name))
174 elseif (present(tr_desc)) then
175 call query_vardesc(tr_desc, name=tr%name, units=tr%units, &
176 longname=tr%longname, cmor_field_name=tr%cmor_name, &
177 cmor_longname=tr%cmor_longname, caller="register_tracer")
178 tr%cmor_units = tr%units
179 else
180 call mom_error(fatal,"MOM register_tracer: Either name or "//&
181 "tr_desc must be present when registering a tracer.")
182 endif
183
184 if (reg%locked) call mom_error(fatal, &
185 "MOM register_tracer was called for variable "//trim(tr%name)//&
186 " with a locked tracer registry.")
187
188 tr%conc_scale = 1.0
189 if (present(conc_scale)) tr%conc_scale = conc_scale
190
191 tr%conc_underflow = 0.0
192 if (present(underflow_conc)) tr%conc_underflow = underflow_conc
193
194 tr%flux_nameroot = tr%name
195 if (present(flux_nameroot)) then
196 if (len_trim(flux_nameroot) > 0) tr%flux_nameroot = flux_nameroot
197 endif
198
199 tr%flux_longname = tr%longname
200 if (present(flux_longname)) then
201 if (len_trim(flux_longname) > 0) tr%flux_longname = flux_longname
202 endif
203
204 tr%net_surfflux_name = "KPP_net"//trim(tr%name)
205 if (present(net_surfflux_name)) then
206 tr%net_surfflux_name = net_surfflux_name
207 endif
208
209 tr%NLT_budget_name = 'KPP_NLT_'//trim(tr%flux_nameroot)//'_budget'
210 if (present(nlt_budget_name)) then
211 tr%NLT_budget_name = nlt_budget_name
212 endif
213
214 tr%net_surfflux_longname = 'Effective net surface '//trim(lowercase(tr%flux_longname))//&
215 ' flux, as used by [CVMix] KPP'
216 if (present(net_surfflux_longname)) then
217 tr%net_surfflux_longname = net_surfflux_longname
218 endif
219
220 tr%flux_units = ""
221 if (present(flux_units)) tr%flux_units = flux_units
222
223 tr%flux_scale = gv%H_to_MKS*tr%conc_scale
224 if (present(flux_scale)) tr%flux_scale = flux_scale
225
226 tr%conv_units = ""
227 if (present(convergence_units)) tr%conv_units = convergence_units
228
229 tr%cmor_tendprefix = ""
230 if (present(cmor_tendprefix)) tr%cmor_tendprefix = cmor_tendprefix
231
232 tr%conv_scale = gv%H_to_MKS*tr%conc_scale
233 if (present(convergence_scale)) then
234 tr%conv_scale = convergence_scale
235 elseif (present(flux_scale)) then
236 tr%conv_scale = flux_scale
237 endif
238
239 tr%diag_form = 1
240 if (present(diag_form)) tr%diag_form = diag_form
241
242 tr%advect_scheme = -1
243 if (present(advect_scheme)) tr%advect_scheme = advect_scheme
244
245 tr%t => tr_ptr
246
247 if (present(registry_diags)) tr%registry_diags = registry_diags
248
249 if (present(ad_x)) then ; if (associated(ad_x)) tr%ad_x => ad_x ; endif
250 if (present(ad_y)) then ; if (associated(ad_y)) tr%ad_y => ad_y ; endif
251 if (present(df_x)) then ; if (associated(df_x)) tr%df_x => df_x ; endif
252 if (present(df_y)) then ; if (associated(df_y)) tr%df_y => df_y ; endif
253! if (present(OBC_in_u)) then ; if (associated(OBC_in_u)) Tr%OBC_in_u => OBC_in_u ; endif
254! if (present(OBC_in_v)) then ; if (associated(OBC_in_v)) Tr%OBC_in_v => OBC_in_v ; endif
255 if (present(ad_2d_x)) then ; if (associated(ad_2d_x)) tr%ad2d_x => ad_2d_x ; endif
256 if (present(ad_2d_y)) then ; if (associated(ad_2d_y)) tr%ad2d_y => ad_2d_y ; endif
257 if (present(df_2d_x)) then ; if (associated(df_2d_x)) tr%df2d_x => df_2d_x ; endif
258 if (present(df_2d_y)) then ; if (associated(df_2d_y)) tr%df2d_y => df_2d_y ; endif
259
260 if (present(advection_xy)) then
261 if (associated(advection_xy)) tr%advection_xy => advection_xy
262 endif
263
264 if (present(restart_cs)) then
265 ! Register this tracer to be read from and written to restart files.
266 mand = .true. ; if (present(mandatory)) mand = mandatory
267
268 call register_restart_field(tr_ptr, tr%name, mand, restart_cs, &
269 longname=tr%longname, units=tr%units, conversion=conc_scale)
270 endif
271end subroutine register_tracer
272
273
274!> This subroutine locks the tracer registry to prevent the addition of more
275!! tracers. After locked=.true., can then register common diagnostics.
276subroutine lock_tracer_registry(Reg)
277 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
278
279 if (.not. associated(reg)) call mom_error(warning, &
280 "lock_tracer_registry called with an unassociated registry.")
281
282 reg%locked = .true.
283
284end subroutine lock_tracer_registry
285
286!> register_tracer_diagnostics does a set of register_diag_field calls for any previously
287!! registered in a tracer registry with a value of registry_diags set to .true.
288subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, use_KPP)
289 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
290 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
291 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
292 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
293 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
294 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
295 type(time_type), intent(in) :: time !< current model time
296 type(diag_ctrl), intent(in) :: diag !< structure to regulate diagnostic output
297 logical, intent(in) :: use_ale !< If true active diagnostics that only
298 !! apply to ALE configurations
299 logical, intent(in) :: use_kpp !< If true active diagnostics that only
300 !! apply to CVMix KPP mixings
301
302 character(len=24) :: name ! A variable's name in a NetCDF file.
303 character(len=24) :: shortnm ! A shortened version of a variable's name for
304 ! creating additional diagnostics.
305 character(len=72) :: longname ! The long name of that tracer variable.
306 character(len=72) :: flux_longname ! The tracer name in the long names of fluxes.
307 character(len=48) :: units ! The dimensions of the tracer.
308 character(len=48) :: flux_units ! The units for fluxes, either
309 ! [units] m3 s-1 or [units] kg s-1.
310 character(len=48) :: conv_units ! The units for flux convergences, either
311 ! [units] m s-1 or [units] kg m-2 s-1.
312 character(len=48) :: unit2 ! The dimensions of the tracer squared
313 character(len=72) :: cmorname ! The CMOR name of this tracer.
314 character(len=120) :: cmor_longname ! The CMOR long name of that variable.
315 character(len=120) :: var_lname ! A temporary longname for a diagnostic.
316 character(len=120) :: cmor_var_lname ! The temporary CMOR long name for a diagnostic
317 real :: conversion ! Temporary term while we address a bug [conc m CU-1 H-1 ~> 1] or [conc kg m-2 CU-1 H-1 ~> 1]
318 type(tracer_type), pointer :: tr=>null()
319 integer :: i, j, k, is, ie, js, je, nz, m, m2, ntr_in
320 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
321 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
322 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
323 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
324
325 if (.not. associated(reg)) call mom_error(fatal, "register_tracer_diagnostics: "//&
326 "register_tracer must be called before register_tracer_diagnostics")
327
328 ntr_in = reg%ntr
329
330 do m=1,ntr_in ; if (reg%Tr(m)%registry_diags) then
331 tr => reg%Tr(m)
332! call query_vardesc(Tr%vd, name, units=units, longname=longname, &
333! cmor_field_name=cmorname, cmor_longname=cmor_longname, &
334! caller="register_tracer_diagnostics")
335 name = tr%name ; units=adjustl(tr%units) ; longname = tr%longname
336 cmorname = tr%cmor_name ; cmor_longname = tr%cmor_longname
337 shortnm = tr%flux_nameroot
338 flux_longname = tr%flux_longname
339 if (len_trim(cmor_longname) == 0) cmor_longname = longname
340
341 if (len_trim(tr%flux_units) > 0) then ; flux_units = tr%flux_units
342 elseif (gv%Boussinesq) then ; flux_units = trim(units)//" m3 s-1"
343 else ; flux_units = trim(units)//" kg s-1" ; endif
344
345 if (len_trim(tr%conv_units) > 0) then ; conv_units = tr%conv_units
346 elseif (gv%Boussinesq) then ; conv_units = trim(units)//" m s-1"
347 else ; conv_units = trim(units)//" kg m-2 s-1" ; endif
348
349 if (len_trim(cmorname) == 0) then
350 tr%id_tr = register_diag_field("ocean_model", trim(name), diag%axesTL, &
351 time, trim(longname), trim(units), conversion=tr%conc_scale)
352 else
353 tr%id_tr = register_diag_field("ocean_model", trim(name), diag%axesTL, &
354 time, trim(longname), trim(units), conversion=tr%conc_scale, &
355 cmor_field_name=cmorname, cmor_long_name=cmor_longname, &
356 cmor_units=tr%cmor_units, cmor_standard_name=cmor_long_std(cmor_longname))
357 endif
358 tr%id_tr_post_horzn = register_diag_field("ocean_model", &
359 trim(name)//"_post_horzn", diag%axesTL, time, &
360 trim(longname)//" after horizontal transport (advection/diffusion) has occurred", &
361 trim(units), conversion=tr%conc_scale)
362 if (tr%diag_form == 1) then
363 tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", &
364 diag%axesCuL, time, trim(flux_longname)//" advective zonal flux" , &
365 trim(flux_units), v_extensive=.true., y_cell_method='sum', &
366 conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
367 tr%id_ady = register_diag_field("ocean_model", trim(shortnm)//"_ady", &
368 diag%axesCvL, time, trim(flux_longname)//" advective meridional flux" , &
369 trim(flux_units), v_extensive=.true., x_cell_method='sum', &
370 conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
371 tr%id_adx_resolved = register_diag_field("ocean_model", trim(shortnm)//"_adx_resolved", &
372 diag%axesCuL, time, trim(flux_longname)//" resolved advective zonal flux" , &
373 trim(flux_units), v_extensive=.true., y_cell_method='sum', &
374 conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
375 tr%id_ady_resolved = register_diag_field("ocean_model", trim(shortnm)//"_ady_resolved", &
376 diag%axesCvL, time, trim(flux_longname)//" resolved advective meridional flux" , &
377 trim(flux_units), v_extensive=.true., x_cell_method='sum', &
378 conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
379 tr%id_adx_param = register_diag_field("ocean_model", trim(shortnm)//"_adx_param", &
380 diag%axesCuL, time, trim(flux_longname)//" parameterized advective zonal flux" , &
381 trim(flux_units), v_extensive=.true., y_cell_method='sum', &
382 conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
383 tr%id_ady_param = register_diag_field("ocean_model", trim(shortnm)//"_ady_param", &
384 diag%axesCvL, time, trim(flux_longname)//" resolved parameterized meridional flux" , &
385 trim(flux_units), v_extensive=.true., x_cell_method='sum', &
386 conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
387 tr%id_dfx = register_diag_field("ocean_model", trim(shortnm)//"_dfx", &
388 diag%axesCuL, time, trim(flux_longname)//" diffusive zonal flux" , &
389 trim(flux_units), v_extensive=.true., y_cell_method='sum', &
390 conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
391 tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_dfy", &
392 diag%axesCvL, time, trim(flux_longname)//" diffusive meridional flux" , &
393 trim(flux_units), v_extensive=.true., x_cell_method='sum', &
394 conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
395 tr%id_hbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx", &
396 diag%axesCuL, time, trim(flux_longname)//" diffusive zonal flux " //&
397 "from the horizontal boundary diffusion scheme", trim(flux_units), v_extensive=.true., &
398 y_cell_method='sum', conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
399 tr%id_hbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy", &
400 diag%axesCvL, time, trim(flux_longname)//" diffusive meridional " //&
401 "flux from the horizontal boundary diffusion scheme", &
402 trim(flux_units), v_extensive=.true., x_cell_method='sum', &
403 conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
404 else
405 tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", &
406 diag%axesCuL, time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), &
407 flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, y_cell_method='sum')
408 tr%id_ady = register_diag_field("ocean_model", trim(shortnm)//"_ady", &
409 diag%axesCvL, time, "Advective (by residual mean) Meridional Flux of "//trim(flux_longname), &
410 flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, x_cell_method='sum')
411 tr%id_adx_resolved = register_diag_field("ocean_model", trim(shortnm)//"_adx_resolved", &
412 diag%axesCuL, time, "Advective (by resolved flow) Zonal Flux of "//trim(flux_longname), &
413 flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, y_cell_method='sum')
414 tr%id_ady_resolved = register_diag_field("ocean_model", trim(shortnm)//"_ady_resolved", &
415 diag%axesCvL, time, "Advective (by resolved flow) Meridional Flux of "//trim(flux_longname), &
416 flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, x_cell_method='sum')
417 tr%id_adx_param = register_diag_field("ocean_model", trim(shortnm)//"_adx_param", &
418 diag%axesCuL, time, "Advective (by parameterized flow) Zonal Flux of "//trim(flux_longname), &
419 flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, y_cell_method='sum')
420 tr%id_ady_param = register_diag_field("ocean_model", trim(shortnm)//"_ady_param", &
421 diag%axesCvL, time, "Advective (by parameterized flow) Meridional Flux of "//trim(flux_longname), &
422 flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, x_cell_method='sum')
423 tr%id_dfx = register_diag_field("ocean_model", trim(shortnm)//"_diffx", &
424 diag%axesCuL, time, "Diffusive Zonal Flux of "//trim(flux_longname), &
425 flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
426 y_cell_method='sum')
427 tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_diffy", &
428 diag%axesCvL, time, "Diffusive Meridional Flux of "//trim(flux_longname), &
429 flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
430 x_cell_method='sum')
431 tr%id_hbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx", &
432 diag%axesCuL, time, &
433 "Horizontal Boundary Diffusive Zonal Flux of "//trim(flux_longname), &
434 flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
435 y_cell_method='sum')
436 tr%id_hbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy", &
437 diag%axesCvL, time, &
438 "Horizontal Boundary Diffusive Meridional Flux of "//trim(flux_longname), &
439 flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
440 x_cell_method='sum')
441 endif
442 tr%id_zint = register_diag_field("ocean_model", trim(shortnm)//"_zint", &
443 diag%axesT1, time, &
444 "Thickness-weighted integral of " // trim(longname), &
445 trim(units) // " m", conversion=tr%conc_scale*us%Z_to_m)
446 tr%id_zint_100m = register_diag_field("ocean_model", trim(shortnm)//"_zint_100m", &
447 diag%axesT1, time, &
448 "Thickness-weighted integral of "// trim(longname) // " over top 100m", &
449 trim(units) // " m", conversion=tr%conc_scale*us%Z_to_m)
450 tr%id_surf = register_diag_field("ocean_model", trim(shortnm)//"_SURF", &
451 diag%axesT1, time, "Surface values of "// trim(longname), trim(units), conversion=tr%conc_scale)
452 if (tr%id_adx > 0) call safe_alloc_ptr(tr%ad_x,isdb,iedb,jsd,jed,nz)
453 if (tr%id_ady > 0) call safe_alloc_ptr(tr%ad_y,isd,ied,jsdb,jedb,nz)
454 if (tr%id_adx_resolved > 0) call safe_alloc_ptr(tr%ad_x_resolved,isdb,iedb,jsd,jed,nz)
455 if (tr%id_ady_resolved > 0) call safe_alloc_ptr(tr%ad_y_resolved,isd,ied,jsdb,jedb,nz)
456 if (tr%id_adx_param > 0) call safe_alloc_ptr(tr%ad_x_param,isdb,iedb,jsd,jed,nz)
457 if (tr%id_ady_param > 0) call safe_alloc_ptr(tr%ad_y_param,isd,ied,jsdb,jedb,nz)
458 if (tr%id_dfx > 0) call safe_alloc_ptr(tr%df_x,isdb,iedb,jsd,jed,nz)
459 if (tr%id_dfy > 0) call safe_alloc_ptr(tr%df_y,isd,ied,jsdb,jedb,nz)
460 if (tr%id_hbd_dfx > 0) call safe_alloc_ptr(tr%hbd_dfx,isdb,iedb,jsd,jed,nz)
461 if (tr%id_hbd_dfy > 0) call safe_alloc_ptr(tr%hbd_dfy,isd,ied,jsdb,jedb,nz)
462
463 tr%id_adx_2d = register_diag_field("ocean_model", trim(shortnm)//"_adx_2d", &
464 diag%axesCu1, time, &
465 "Vertically Integrated Advective Zonal Flux of "//trim(flux_longname), &
466 flux_units, conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, y_cell_method='sum')
467 tr%id_ady_2d = register_diag_field("ocean_model", trim(shortnm)//"_ady_2d", &
468 diag%axesCv1, time, &
469 "Vertically Integrated Advective Meridional Flux of "//trim(flux_longname), &
470 flux_units, conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, x_cell_method='sum')
471 tr%id_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_diffx_2d", &
472 diag%axesCu1, time, &
473 "Vertically Integrated Diffusive Zonal Flux of "//trim(flux_longname), &
474 flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
475 y_cell_method='sum')
476 tr%id_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_diffy_2d", &
477 diag%axesCv1, time, &
478 "Vertically Integrated Diffusive Meridional Flux of "//trim(flux_longname), &
479 flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
480 x_cell_method='sum')
481 tr%id_hbd_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx_2d", &
482 diag%axesCu1, time, "Vertically-integrated zonal diffusive flux from the horizontal boundary diffusion "//&
483 "scheme for "//trim(flux_longname), flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
484 y_cell_method='sum')
485 tr%id_hbd_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy_2d", &
486 diag%axesCv1, time, "Vertically-integrated meridional diffusive flux from the horizontal boundary diffusion "//&
487 "scheme for "//trim(flux_longname), flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
488 x_cell_method='sum')
489
490 if (tr%id_adx_2d > 0) call safe_alloc_ptr(tr%ad2d_x,isdb,iedb,jsd,jed)
491 if (tr%id_ady_2d > 0) call safe_alloc_ptr(tr%ad2d_y,isd,ied,jsdb,jedb)
492 if (tr%id_dfx_2d > 0) call safe_alloc_ptr(tr%df2d_x,isdb,iedb,jsd,jed)
493 if (tr%id_dfy_2d > 0) call safe_alloc_ptr(tr%df2d_y,isd,ied,jsdb,jedb)
494 if (tr%id_hbd_dfx_2d > 0) call safe_alloc_ptr(tr%hbd_dfx_2d,isdb,iedb,jsd,jed)
495 if (tr%id_hbd_dfy_2d > 0) call safe_alloc_ptr(tr%hbd_dfy_2d,isd,ied,jsdb,jedb)
496
497 tr%id_adv_xy = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy", &
498 diag%axesTL, time, &
499 'Horizontal convergence of residual mean advective fluxes of '//&
500 trim(lowercase(flux_longname)), &
501 conv_units, v_extensive=.true., conversion=tr%conv_scale*us%s_to_T)
502 tr%id_adv_xy_2d = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy_2d", &
503 diag%axesT1, time, &
504 'Vertical sum of horizontal convergence of residual mean advective fluxes of '//&
505 trim(lowercase(flux_longname)), conv_units, conversion=tr%conv_scale*us%s_to_T)
506 if ((tr%id_adv_xy > 0) .or. (tr%id_adv_xy_2d > 0)) &
507 call safe_alloc_ptr(tr%advection_xy,isd,ied,jsd,jed,nz)
508
509 tr%id_tendency = register_diag_field('ocean_model', trim(shortnm)//'_tendency', &
510 diag%axesTL, time, &
511 'Net time tendency for '//trim(lowercase(longname)), &
512 trim(units)//' s-1', conversion=tr%conc_scale*us%s_to_T)
513
514 if (tr%id_tendency > 0) then
515 call safe_alloc_ptr(tr%t_prev,isd,ied,jsd,jed,nz)
516 do k=1,nz ; do j=js,je ; do i=is,ie
517 tr%t_prev(i,j,k) = tr%t(i,j,k)
518 enddo ; enddo ; enddo
519 endif
520
521 ! Neutral/Horizontal diffusion convergence tendencies
522 if (tr%diag_form == 1) then
523 tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', &
524 diag%axesTL, time, "Neutral diffusion tracer content tendency for "//trim(shortnm), &
525 conv_units, conversion=tr%conv_scale*us%s_to_T, v_extensive=.true.)
526
527 tr%id_dfxy_cont_2d = register_diag_field("ocean_model", &
528 trim(shortnm)//'_dfxy_cont_tendency_2d', &
529 diag%axesT1, time, "Depth integrated neutral diffusion tracer content "//&
530 "tendency for "//trim(shortnm), conv_units, conversion=tr%conv_scale*us%s_to_T)
531
532 tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', &
533 diag%axesTL, time, "Horizontal boundary diffusion tracer content tendency for "//&
534 trim(shortnm), &
535 conv_units, conversion=tr%conv_scale*us%s_to_T, v_extensive=.true.)
536
537 tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", &
538 trim(shortnm)//'_hbdxy_cont_tendency_2d', &
539 diag%axesT1, time, "Depth integrated horizontal boundary diffusion tracer content "//&
540 "tendency for "//trim(shortnm), conv_units, conversion=tr%conv_scale*us%s_to_T)
541 else
542 cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//&
543 trim(lowercase(flux_longname))//&
544 ' content due to parameterized mesoscale neutral diffusion'
545 tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', &
546 diag%axesTL, time, "Neutral diffusion tracer content tendency for "//trim(shortnm), &
547 conv_units, conversion=tr%conv_scale*us%s_to_T, v_extensive=.true., &
548 cmor_field_name=trim(tr%cmor_tendprefix)//'pmdiff', &
549 cmor_long_name=trim(cmor_var_lname), &
550 cmor_standard_name=trim(cmor_long_std(cmor_var_lname)))
551
552 cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//&
553 trim(lowercase(flux_longname))//&
554 ' content due to parameterized mesoscale neutral diffusion'
555 tr%id_dfxy_cont_2d = register_diag_field("ocean_model", &
556 trim(shortnm)//'_dfxy_cont_tendency_2d', &
557 diag%axesT1, time, "Depth integrated neutral diffusion tracer "//&
558 "content tendency for "//trim(shortnm), conv_units, conversion=tr%conv_scale*us%s_to_T, &
559 cmor_field_name=trim(tr%cmor_tendprefix)//'pmdiff_2d', &
560 cmor_long_name=trim(cmor_var_lname), &
561 cmor_standard_name=trim(cmor_long_std(cmor_var_lname)))
562
563 tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', &
564 diag%axesTL, time, &
565 "Horizontal boundary diffusion tracer content tendency for "//trim(shortnm), &
566 conv_units, conversion=tr%conv_scale*us%s_to_T, v_extensive=.true.)
567
568 tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", &
569 trim(shortnm)//'_hbdxy_cont_tendency_2d', &
570 diag%axesT1, time, "Depth integrated horizontal boundary diffusion of tracer "//&
571 "content tendency for "//trim(shortnm), conv_units, conversion=tr%conv_scale*us%s_to_T)
572 endif
573 tr%id_dfxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_conc_tendency', &
574 diag%axesTL, time, "Neutral diffusion tracer concentration tendency for "//trim(shortnm), &
575 trim(units)//' s-1', conversion=tr%conc_scale*us%s_to_T)
576
577 tr%id_hbdxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_conc_tendency', &
578 diag%axesTL, time, &
579 "Horizontal diffusion tracer concentration tendency for "//trim(shortnm), &
580 trim(units)//' s-1', conversion=tr%conc_scale*us%s_to_T)
581
582 var_lname = "Net time tendency for "//lowercase(flux_longname)
583 if (len_trim(tr%cmor_tendprefix) == 0) then
584 tr%id_trxh_tendency = register_diag_field('ocean_model', trim(shortnm)//'h_tendency', &
585 diag%axesTL, time, var_lname, conv_units, conversion=tr%conv_scale*us%s_to_T, &
586 v_extensive=.true.)
587 tr%id_trxh_tendency_2d = register_diag_field('ocean_model', trim(shortnm)//'h_tendency_2d', &
588 diag%axesT1, time, "Vertical sum of "//trim(lowercase(var_lname)), &
589 conv_units, conversion=tr%conv_scale*us%s_to_T)
590 else
591 cmor_var_lname = "Tendency of "//trim(cmor_longname)//" Expressed as "//&
592 trim(flux_longname)//" Content"
593 tr%id_trxh_tendency = register_diag_field('ocean_model', trim(shortnm)//'h_tendency', &
594 diag%axesTL, time, var_lname, conv_units, conversion=tr%conv_scale*us%s_to_T, &
595 cmor_field_name=trim(tr%cmor_tendprefix)//"tend", &
596 cmor_standard_name=cmor_long_std(cmor_var_lname), cmor_long_name=cmor_var_lname, &
597 v_extensive=.true.)
598 cmor_var_lname = trim(cmor_var_lname)//" Vertical Sum"
599 tr%id_trxh_tendency_2d = register_diag_field('ocean_model', trim(shortnm)//'h_tendency_2d', &
600 diag%axesT1, time, "Vertical sum of "//trim(lowercase(var_lname)), &
601 conv_units, conversion=tr%conv_scale*us%s_to_T, &
602 cmor_field_name=trim(tr%cmor_tendprefix)//"tend_2d", &
603 cmor_standard_name=cmor_long_std(cmor_var_lname), cmor_long_name=cmor_var_lname)
604 endif
605 if ((tr%id_trxh_tendency > 0) .or. (tr%id_trxh_tendency_2d > 0)) then
606 call safe_alloc_ptr(tr%Trxh_prev,isd,ied,jsd,jed,nz)
607 do k=1,nz ; do j=js,je ; do i=is,ie
608 tr%Trxh_prev(i,j,k) = tr%t(i,j,k) * h(i,j,k)
609 enddo ; enddo ; enddo
610 endif
611
612 ! Vertical regridding/remapping tendencies
613 if (use_ale .and. tr%remap_tr) then
614 var_lname = "Vertical remapping tracer concentration tendency for "//trim(reg%Tr(m)%name)
615 tr%id_remap_conc= register_diag_field('ocean_model', &
616 trim(tr%flux_nameroot)//'_tendency_vert_remap', diag%axesTL, time, var_lname, &
617 trim(units)//' s-1', conversion=tr%conc_scale*us%s_to_T)
618
619 var_lname = "Vertical remapping tracer content tendency for "//trim(reg%Tr(m)%flux_longname)
620 tr%id_remap_cont = register_diag_field('ocean_model', &
621 trim(tr%flux_nameroot)//'h_tendency_vert_remap', &
622 diag%axesTL, time, var_lname, conv_units, v_extensive=.true., conversion=tr%conv_scale*us%s_to_T)
623
624 var_lname = "Vertical sum of vertical remapping tracer content tendency for "//&
625 trim(reg%Tr(m)%flux_longname)
626 tr%id_remap_cont_2d = register_diag_field('ocean_model', &
627 trim(tr%flux_nameroot)//'h_tendency_vert_remap_2d', &
628 diag%axesT1, time, var_lname, conv_units, conversion=tr%conv_scale*us%s_to_T)
629
630 endif
631
632 if (use_ale .and. (reg%ntr<max_fields_) .and. tr%remap_tr) then
633 unit2 = trim(units)//"2"
634 if (index(units(1:len_trim(units))," ") > 0) unit2 = "("//trim(units)//")2"
635 tr%id_tr_vardec = register_diag_field('ocean_model', trim(shortnm)//"_vardec", diag%axesTL, &
636 time, "ALE variance decay for "//lowercase(longname), &
637 trim(unit2)//" s-1", conversion=tr%conc_scale**2*us%s_to_T)
638 if (tr%id_tr_vardec > 0) then
639 ! Set up a new tracer for this tracer squared
640 m2 = reg%ntr+1
641 tr%ind_tr_squared = m2
642 call safe_alloc_ptr(reg%Tr(m2)%t,isd,ied,jsd,jed,nz) ; reg%Tr(m2)%t(:,:,:) = 0.0
643 reg%Tr(m2)%name = trim(shortnm)//"2"
644 reg%Tr(m2)%longname = "Squared "//trim(longname)
645 reg%Tr(m2)%units = unit2
646 reg%Tr(m2)%registry_diags = .false.
647 reg%Tr(m2)%ind_tr_squared = -1
648 ! Augment the total number of tracers, including the squared tracers.
649 reg%ntr = reg%ntr + 1
650 endif
651 endif
652
653 ! KPP nonlocal term diagnostics
654 if (use_kpp) then
655 tr%id_net_surfflux = register_diag_field('ocean_model', tr%net_surfflux_name, diag%axesT1, time, &
656 tr%net_surfflux_longname, trim(units)//' m s-1', conversion=tr%conc_scale*gv%H_to_m*us%s_to_T)
657 tr%id_NLT_tendency = register_diag_field('ocean_model', "KPP_NLT_d"//trim(shortnm)//"dt", &
658 diag%axesTL, time, &
659 trim(longname)//' tendency due to non-local transport of '//trim(lowercase(flux_longname))//&
660 ', as calculated by [CVMix] KPP', trim(units)//' s-1', conversion=tr%conc_scale*us%s_to_T)
661 if (tr%conv_scale == 0.001*gv%H_to_kg_m2) then
662 conversion = gv%H_to_kg_m2
663 else
664 conversion = tr%conv_scale
665 endif
666 ! We actually want conversion=Tr%conv_scale for all tracers, but introducing the local variable
667 ! 'conversion' and setting it to GV%H_to_kg_m2 instead of 0.001*GV%H_to_kg_m2 for salt tracers
668 ! keeps changes introduced by this refactoring limited to round-off level; as it turns out,
669 ! there is a bug in the code and the NLT budget term for salinity is off by a factor of 10^3
670 ! so introducing the 0.001 here will fix that bug.
671 tr%id_NLT_budget = register_diag_field('ocean_model', tr%NLT_budget_name, &
672 diag%axesTL, time, &
673 trim(flux_longname)//&
674 ' content change due to non-local transport, as calculated by [CVMix] KPP', &
675 conv_units, conversion=conversion*us%s_to_T, v_extensive=.true.)
676 endif
677
678 endif ; enddo
679
680end subroutine register_tracer_diagnostics
681
682subroutine preale_tracer_diagnostics(Reg, G, GV)
683 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
684 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
685 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
686
687 integer :: i, j, k, is, ie, js, je, nz, m, m2
688 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
689
690 do m=1,reg%ntr ; if (reg%Tr(m)%ind_tr_squared > 0) then
691 m2 = reg%Tr(m)%ind_tr_squared
692 ! Update squared quantities
693 do k=1,nz ; do j=js,je ; do i=is,ie
694 reg%Tr(m2)%T(i,j,k) = reg%Tr(m)%T(i,j,k)**2
695 enddo ; enddo ; enddo
696 endif ; enddo
697
698end subroutine preale_tracer_diagnostics
699
700subroutine postale_tracer_diagnostics(Reg, G, GV, diag, dt)
701 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
702 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
703 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
704 type(diag_ctrl), intent(in) :: diag !< regulates diagnostic output
705 real, intent(in) :: dt !< total time interval for these diagnostics [T ~> s]
706
707 real :: work(szi_(g),szj_(g),szk_(gv)) ! Variance decay [CU2 T-1 ~> conc2 s-1]
708 real :: idt ! The inverse of the time step [T-1 ~> s-1]
709 integer :: i, j, k, is, ie, js, je, nz, m, m2
710 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
711
712 ! The "if" is to avoid NaNs if the diagnostic is called for a zero length interval
713 idt = 0.0 ; if (dt /= 0.0) idt = 1.0 / dt
714
715 do m=1,reg%ntr ; if (reg%Tr(m)%id_tr_vardec > 0) then
716 m2 = reg%Tr(m)%ind_tr_squared
717 if (m2 < 1) call mom_error(fatal, "Bad value of Tr%ind_tr_squared for "//trim(reg%Tr(m)%name))
718 ! Update squared quantities
719 do k=1,nz ; do j=js,je ; do i=is,ie
720 work(i,j,k) = (reg%Tr(m2)%T(i,j,k) - reg%Tr(m)%T(i,j,k)**2) * idt
721 enddo ; enddo ; enddo
722 call post_data(reg%Tr(m)%id_tr_vardec, work, diag)
723 endif ; enddo
724
725end subroutine postale_tracer_diagnostics
726
727!> Post tracer diganostics when that should only be posted when MOM's state
728!! is self-consistent (also referred to as 'synchronized')
729subroutine post_tracer_diagnostics_at_sync(Reg, h, diag_prev, diag, G, GV, dt)
730 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
731 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
732 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
733 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
734 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
735 type(diag_grid_storage), intent(in) :: diag_prev !< Contains diagnostic grids from previous timestep
736 type(diag_ctrl), intent(inout) :: diag !< structure to regulate diagnostic output
737 real, intent(in) :: dt !< total time step for tracer updates [T ~> s]
738
739 real :: work3d(szi_(g),szj_(g),szk_(gv)) ! The time tendency of a diagnostic [CU T-1 ~> conc s-1]
740 real :: work2d(szi_(g),szj_(g)) ! The vertically integrated time tendency of a diagnostic
741 ! in [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1]
742 real :: idt ! The inverse of the time step [T-1 ~> s-1]
743 type(tracer_type), pointer :: tr=>null()
744 integer :: i, j, k, is, ie, js, je, nz, m
745 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
746
747 idt = 0. ; if (dt/=0.) idt = 1.0 / dt ! The "if" is in case the diagnostic is called for a zero length interval
748
749 ! Tendency diagnostics need to be posted on the grid from the last call to this routine
750 call diag_save_grids(diag)
751 call diag_copy_storage_to_diag(diag, diag_prev)
752 do m=1,reg%ntr ; if (reg%Tr(m)%registry_diags) then
753 tr => reg%Tr(m)
754 if (tr%id_tr > 0) call post_data(tr%id_tr, tr%t, diag)
755 if (tr%id_tendency > 0) then
756 work3d(:,:,:) = 0.0
757 do k=1,nz ; do j=js,je ; do i=is,ie
758 work3d(i,j,k) = (tr%t(i,j,k) - tr%t_prev(i,j,k))*idt
759 tr%t_prev(i,j,k) = tr%t(i,j,k)
760 enddo ; enddo ; enddo
761 call post_data(tr%id_tendency, work3d, diag, alt_h=diag_prev%h_state)
762 endif
763 if ((tr%id_trxh_tendency > 0) .or. (tr%id_trxh_tendency_2d > 0)) then
764 do k=1,nz ; do j=js,je ; do i=is,ie
765 work3d(i,j,k) = (tr%t(i,j,k)*h(i,j,k) - tr%Trxh_prev(i,j,k)) * idt
766 tr%Trxh_prev(i,j,k) = tr%t(i,j,k) * h(i,j,k)
767 enddo ; enddo ; enddo
768 if (tr%id_trxh_tendency > 0) call post_data(tr%id_trxh_tendency, work3d, diag, &
769 alt_h=diag_prev%h_state)
770 if (tr%id_trxh_tendency_2d > 0) then
771 work2d(:,:) = 0.0
772 do k=1,nz ; do j=js,je ; do i=is,ie
773 work2d(i,j) = work2d(i,j) + work3d(i,j,k)
774 enddo ; enddo ; enddo
775 call post_data(tr%id_trxh_tendency_2d, work2d, diag)
776 endif
777 endif
778 endif ; enddo
779 call diag_restore_grids(diag)
780
782
783!> Post the advective and diffusive tendencies
784subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag)
785 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
786 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
787 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
788 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
789 intent(in) :: h_diag !< Layer thicknesses on which to post fields [H ~> m or kg m-2]
790 type(diag_ctrl), intent(in) :: diag !< structure to regulate diagnostic output
791
792 integer :: i, j, k, is, ie, js, je, nz, m
793 real :: work2d(szi_(g),szj_(g)) ! The vertically integrated convergence of lateral advective
794 ! tracer fluxes [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1]
795 type(tracer_type), pointer :: tr=>null()
796
797 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
798
799 do m=1,reg%ntr ; if (reg%Tr(m)%registry_diags) then
800 tr => reg%Tr(m)
801 if (tr%id_tr_post_horzn> 0) call post_data(tr%id_tr_post_horzn, tr%t, diag)
802 if (tr%id_adx > 0) call post_data(tr%id_adx, tr%ad_x, diag, alt_h=h_diag)
803 if (tr%id_ady > 0) call post_data(tr%id_ady, tr%ad_y, diag, alt_h=h_diag)
804 if (tr%id_adx_resolved > 0) call post_data(tr%id_adx_resolved, tr%ad_x_resolved, diag, alt_h=h_diag)
805 if (tr%id_ady_resolved > 0) call post_data(tr%id_ady_resolved, tr%ad_y_resolved, diag, alt_h=h_diag)
806 if (tr%id_adx_param > 0) call post_data(tr%id_adx_param, tr%ad_x_param, diag, alt_h=h_diag)
807 if (tr%id_ady_param > 0) call post_data(tr%id_ady_param, tr%ad_y_param, diag, alt_h=h_diag)
808 if (tr%id_dfx > 0) call post_data(tr%id_dfx, tr%df_x, diag, alt_h=h_diag)
809 if (tr%id_dfy > 0) call post_data(tr%id_dfy, tr%df_y, diag, alt_h=h_diag)
810 if (tr%id_adx_2d > 0) call post_data(tr%id_adx_2d, tr%ad2d_x, diag)
811 if (tr%id_ady_2d > 0) call post_data(tr%id_ady_2d, tr%ad2d_y, diag)
812 if (tr%id_dfx_2d > 0) call post_data(tr%id_dfx_2d, tr%df2d_x, diag)
813 if (tr%id_dfy_2d > 0) call post_data(tr%id_dfy_2d, tr%df2d_y, diag)
814 if (tr%id_adv_xy > 0) call post_data(tr%id_adv_xy, tr%advection_xy, diag, alt_h=h_diag)
815 if (tr%id_adv_xy_2d > 0) then
816 work2d(:,:) = 0.0
817 do k=1,nz ; do j=js,je ; do i=is,ie
818 work2d(i,j) = work2d(i,j) + tr%advection_xy(i,j,k)
819 enddo ; enddo ; enddo
820 call post_data(tr%id_adv_xy_2d, work2d, diag)
821 endif
822 endif ; enddo
823
825
826!> Post diagnostics of vertically integrated tracer amouints
827subroutine post_tracer_integral_diagnostics(G, GV, US, Reg, h_diag, tv, diag)
828 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
829 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
830 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
831 type(tracer_registry_type), pointer :: reg !< pointer to the tracer registry
832 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
833 intent(in) :: h_diag !< Layer thicknesses on which to post fields [H ~> m or kg m-2]
834 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
835 !! thermodynamic variables.
836 type(diag_ctrl), intent(in) :: diag !< structure to regulate diagnostic output
837
838 integer :: i, j, k, is, ie, js, je, nz, m, khi
839 real :: work2d(szi_(g),szj_(g)) ! The vertically integrated tracer amounts [CU Z T-1 ~> conc m]
840 real :: dz(szi_(g),szj_(g),szk_(gv)) !< Geometric layer thicknesses in height units [Z ~> m]
841 real :: frac_under_100m(szi_(g),szj_(g),szk_(gv)) ! weights used to compute 100m vertical integrals [nondim]
842 real :: ztop(szi_(g),szj_(g)) ! position of the top interface [Z ~> m]
843 real :: zbot(szi_(g),szj_(g)) ! position of the bottom interface [Z ~> m]
844 real :: z_100 ! 100 m in depth units [Z ~> m]
845 logical :: dz_needed, dz100_used
846 type(tracer_type), pointer :: tr=>null()
847
848 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
849
850 dz_needed = .false.
851 dz100_used = .false.
852 do m=1,reg%ntr ; if (reg%Tr(m)%registry_diags) then
853 if (reg%Tr(m)%id_zint_100m > 0) dz100_used = .true.
854 if (reg%Tr(m)%id_zint > 0) dz_needed = .true.
855 endif ; enddo
856 if (dz100_used) dz_needed = .true.
857
858 if (dz_needed) then
859 ! Convert the layer thicknesses into geometric depths, using the pre-stored layer-mean specific
860 ! volumes when in non-Boussinesq mode.
861 call thickness_to_dz(h_diag, tv, dz, g, gv, us)
862 endif
863
864 if (dz100_used) then
865 ! If any tracers are posting 100m vertical integrals, compute weights
866 frac_under_100m(:,:,:) = 0.0
867 ! khi will be the largest layer index corresponding where ztop < 100m and ztop >= 100m
868 ! in any column (we can reduce computation of 100m integrals by only looping through khi
869 ! rather than GV%ke)
870 khi = 0
871
872 z_100 = 100.0*us%m_to_Z
873 zbot(:,:) = 0.0
874 do k=1,nz
875 do j=js,je ; do i=is,ie
876 ztop(i,j) = zbot(i,j)
877 zbot(i,j) = ztop(i,j) + dz(i,j,k)
878 if (zbot(i,j) <= z_100) then
879 frac_under_100m(i,j,k) = 1.0
880 elseif (ztop(i,j) < z_100) then
881 frac_under_100m(i,j,k) = (z_100 - ztop(i,j)) / (zbot(i,j) - ztop(i,j))
882 else
883 frac_under_100m(i,j,k) = 0.0
884 endif
885 ! frac_under_100m(i,j,k) = max(0, min(1.0, (Z_100 - ztop(i,j)) / (zbot(i,j) - ztop(i,j))))
886 enddo ; enddo
887 if (any(frac_under_100m(:,:,k) > 0)) khi = k
888 enddo
889 endif
890
891 do m=1,reg%ntr ; if (reg%Tr(m)%registry_diags) then
892 tr => reg%Tr(m)
893 ! A few diagnostics introduce with MARBL driver
894 ! Compute full-depth vertical integral
895 if (tr%id_zint > 0) then
896 work2d(:,:) = 0.0
897 do k=1,nz ; do j=js,je ; do i=is,ie
898 work2d(i,j) = work2d(i,j) + dz(i,j,k)*tr%t(i,j,k)
899 enddo ; enddo ; enddo
900 call post_data(tr%id_zint, work2d, diag)
901 endif
902
903 ! Compute 100m vertical integral
904 if (tr%id_zint_100m > 0) then
905 work2d(:,:) = 0.0
906 do k=1,khi ; do j=js,je ; do i=is,ie
907 work2d(i,j) = work2d(i,j) + frac_under_100m(i,j,k) * dz(i,j,k)*tr%t(i,j,k)
908 enddo ; enddo ; enddo
909 call post_data(tr%id_zint_100m, work2d, diag)
910 endif
911
912 ! Surface values of tracers
913 if (tr%id_SURF > 0) call post_data(tr%id_SURF, tr%t(:,:,1), diag)
914 endif ; enddo
915
917
918!> This subroutine writes out chksums for the first ntr registered tracers.
919subroutine tracer_array_chksum(mesg, Tr, ntr, G)
920 character(len=*), intent(in) :: mesg !< message that appears on the chksum lines
921 type(tracer_type), intent(in) :: Tr(:) !< array of all of registered tracers
922 integer, intent(in) :: ntr !< number of registered tracers
923 type(ocean_grid_type), intent(in) :: G !< ocean grid structure
924
925 integer :: m
926
927 do m=1,ntr
928 call hchksum(tr(m)%t, mesg//trim(tr(m)%name), g%HI, unscale=tr(m)%conc_scale)
929 enddo
930
931end subroutine tracer_array_chksum
932
933!> This subroutine writes out chksums for all the registered tracers.
934subroutine tracer_reg_chksum(mesg, Reg, G)
935 character(len=*), intent(in) :: mesg !< message that appears on the chksum lines
936 type(tracer_registry_type), pointer :: Reg !< pointer to the tracer registry
937 type(ocean_grid_type), intent(in) :: G !< ocean grid structure
938
939 integer :: m
940
941 if (.not.associated(reg)) return
942
943 do m=1,reg%ntr
944 call hchksum(reg%Tr(m)%t, mesg//trim(reg%Tr(m)%name), g%HI, unscale=reg%Tr(m)%conc_scale)
945 enddo
946
947end subroutine tracer_reg_chksum
948
949!> Calculates and prints the global inventory of the first ntr tracers in the registry.
950subroutine tracer_array_chkinv(mesg, G, GV, h, Tr, ntr)
951 character(len=*), intent(in) :: mesg !< message that appears on the chksum lines
952 type(ocean_grid_type), intent(in) :: G !< ocean grid structure
953 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
954 type(tracer_type), dimension(:), intent(in) :: Tr !< array of all of registered tracers
955 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
956 integer, intent(in) :: ntr !< number of registered tracers
957
958 ! Local variables
959 real :: vol_scale ! The dimensional scaling factor to convert volumes to m3 [m3 H-1 L-2 ~> 1] or cell
960 ! masses to kg [kg H-1 L-2 ~> 1], depending on whether the Boussinesq approximation is used
961 real :: tr_inv(SZI_(G),SZJ_(G),SZK_(GV)) ! Volumetric or mass-based tracer inventory in
962 ! each cell [conc m3] or [conc kg]
963 real :: total_inv ! The total amount of tracer [conc m3] or [conc kg]
964 integer :: is, ie, js, je, nz
965 integer :: i, j, k, m
966
967 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
968 vol_scale = gv%H_to_MKS*g%US%L_to_m**2
969 do m=1,ntr
970 do k=1,nz ; do j=js,je ; do i=is,ie
971 tr_inv(i,j,k) = tr(m)%conc_scale*tr(m)%t(i,j,k) * &
972 (vol_scale * h(i,j,k) * g%areaT(i,j)*g%mask2dT(i,j))
973 enddo ; enddo ; enddo
974 total_inv = reproducing_sum(tr_inv, is+(1-g%isd), ie+(1-g%isd), js+(1-g%jsd), je+(1-g%jsd))
975 if (is_root_pe()) write(0,'(A,1X,A5,1X,ES25.16,1X,A)') &
976 "h-point: inventory", tr(m)%name, total_inv, mesg
977 enddo
978
979end subroutine tracer_array_chkinv
980
981
982!> Calculates and prints the global inventory of all tracers in the registry.
983subroutine tracer_reg_chkinv(mesg, G, GV, h, Reg)
984 character(len=*), intent(in) :: mesg !< message that appears on the chksum lines
985 type(ocean_grid_type), intent(in) :: G !< ocean grid structure
986 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
987 type(tracer_registry_type), pointer :: Reg !< pointer to the tracer registry
988 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
989
990 ! Local variables
991 real :: vol_scale ! The dimensional scaling factor to convert volumes to m3 [m3 H-1 L-2 ~> 1] or cell
992 ! masses to kg [kg H-1 L-2 ~> 1], depending on whether the Boussinesq approximation is used
993 real :: tr_inv(SZI_(G),SZJ_(G),SZK_(GV)) ! Volumetric or mass-based tracer inventory in
994 ! each cell [conc m3] or [conc kg]
995 real :: total_inv ! The total amount of tracer [conc m3] or [conc kg]
996 integer :: is, ie, js, je, nz
997 integer :: i, j, k, m
998
999 if (.not.associated(reg)) return
1000
1001 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1002 vol_scale = gv%H_to_MKS*g%US%L_to_m**2
1003 do m=1,reg%ntr
1004 do k=1,nz ; do j=js,je ; do i=is,ie
1005 tr_inv(i,j,k) = reg%Tr(m)%conc_scale*reg%Tr(m)%t(i,j,k) * &
1006 (vol_scale * h(i,j,k) * g%areaT(i,j)*g%mask2dT(i,j))
1007 enddo ; enddo ; enddo
1008 total_inv = reproducing_sum(tr_inv, is+(1-g%isd), ie+(1-g%isd), js+(1-g%jsd), je+(1-g%jsd))
1009 if (is_root_pe()) write(0,'(A,1X,A5,1X,ES25.16,1X,A)') &
1010 "h-point: inventory", reg%Tr(m)%name, total_inv, mesg
1011 enddo
1012
1013end subroutine tracer_reg_chkinv
1014
1015
1016!> Find a tracer in the tracer registry by name.
1017subroutine tracer_name_lookup(Reg, n, tr_ptr, name)
1018 type(tracer_registry_type), pointer :: reg !< pointer to tracer registry
1019 type(tracer_type), pointer :: tr_ptr !< target or pointer to the tracer array
1020 character(len=32), intent(in) :: name !< tracer name
1021 integer, intent(out) :: n !< index to tracer registery
1022
1023 do n=1,reg%ntr
1024 if (lowercase(reg%Tr(n)%name) == lowercase(name)) then
1025 tr_ptr => reg%Tr(n)
1026 return
1027 endif
1028 enddo
1029
1030 call mom_error(fatal,"MOM cannot find registered tracer: "//name)
1031
1032end subroutine tracer_name_lookup
1033
1034!> Initialize the tracer registry.
1035subroutine tracer_registry_init(param_file, Reg)
1036 type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters
1037 type(tracer_registry_type), pointer :: reg !< pointer to tracer registry
1038
1039 integer, save :: init_calls = 0
1040
1041! This include declares and sets the variable "version".
1042#include "version_variable.h"
1043 character(len=40) :: mdl = "MOM_tracer_registry" ! This module's name.
1044 character(len=256) :: mesg ! Message for error messages.
1045
1046 if (.not.associated(reg)) then ; allocate(reg)
1047 else ; return ; endif
1048
1049 ! Read all relevant parameters and write them to the model log.
1050 call log_version(param_file, mdl, version, "", all_default=.true.)
1051
1052 init_calls = init_calls + 1
1053 if (init_calls > 1) then
1054 write(mesg,'("tracer_registry_init called ",I0, &
1055 &" times with different registry pointers.")') init_calls
1056 if (is_root_pe()) call mom_error(warning,"MOM_tracer "//mesg)
1057 endif
1058
1059end subroutine tracer_registry_init
1060
1061
1062!> This routine closes the tracer registry module.
1063subroutine tracer_registry_end(Reg)
1064 type(tracer_registry_type), pointer :: reg !< The tracer registry that will be deallocated
1065 if (associated(reg)) deallocate(reg)
1066end subroutine tracer_registry_end
1067
1068end module mom_tracer_registry