MARBL_tracers.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!> A tracer package for tracers computed in the MARBL library
6!!
7!! Currently configured for use with marbl0.36.0
8!! https://github.com/marbl-ecosys/MARBL/releases/tag/marbl0.36.0
9!! (clone entire repo into pkg/MARBL)
10module marbl_tracers
11
12use mom_coms, only : efp_type, root_pe, broadcast
13use mom_debugging, only : hchksum
15use mom_error_handler, only : is_root_pe, mom_error, fatal, warning, note
16use mom_file_parser, only : get_param, log_param, log_version, param_file_type
17use mom_forcing_type, only : forcing
18use mom_grid, only : ocean_grid_type
19use mom_interpolate, only : external_field, init_external_field, time_interp_external
25use mom_io, only : file_exists, mom_read_data, slasher, vardesc, var_desc, query_vardesc
29use mom_restart, only : query_initialized, mom_restart_cs, register_restart_field
32use mom_time_manager, only : time_type
39use mom_variables, only : surface
41use mom_diag_mediator, only : register_diag_field, post_data!, safe_alloc_ptr
42
45
47
48implicit none ; private
49
50#include <MOM_memory.h>
51
56
57! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
58! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
59! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
60! vary with the Boussinesq approximation, the Boussinesq variant is given first.
61
62!> Temporary type for diagnostic variables coming from MARBL
63!! Allocate exactly one of field_[23]d
64type :: temp_marbl_diag
65 integer :: id !< index into MOM diagnostic structure
66 real, allocatable :: field_2d(:,:) !< memory for 2D field
67 real, allocatable :: field_3d(:,:,:) !< memory for 3D field
68end type temp_marbl_diag
69
70!> MOM6 needs to know the index of some MARBL tracers to properly apply river fluxes
71type :: tracer_ind_type
72 integer :: no3_ind !< NO3 index
73 integer :: po4_ind !< PO4 index
74 integer :: don_ind !< DON index
75 integer :: donr_ind !< DONr index
76 integer :: dop_ind !< DOP index
77 integer :: dopr_ind !< DOPr index
78 integer :: sio3_ind !< SiO3 index
79 integer :: fe_ind !< Fe index
80 integer :: doc_ind !< DOC index
81 integer :: docr_ind !< DOCr index
82 integer :: alk_ind !< ALK index
83 integer :: alk_alt_co2_ind !< ALK_ALT_CO2 index
84 integer :: dic_ind !< DIC index
85 integer :: dic_alt_co2_ind !< DIC_ALT_CO2 index
86 integer :: abio_dic_ind !< ABIO_DIC index
87 integer :: abio_di14c_ind !< ABIO_DI14C index
88end type tracer_ind_type
89
90!> MOM needs to store some information about saved_state; besides providing these
91!! fields to MARBL, they are also written to restart files
93 character(len=200) :: short_name !< name of variable being saved
94 character(len=200) :: file_varname !< name of variable in restart file
95 character(len=200) :: units !< variable units
96 real, pointer :: field_2d(:,:) => null() !< memory for 2D field
97 real, pointer :: field_3d(:,:,:) => null() !< memory for 3D field
99
100!> All calls to MARBL are done via the interface class
101type(marbl_interface_class) :: marbl_instances
102
103!> Pointer to tracer concentration and to tracer_type in tracer registry
104type, private :: marbl_tracer_data
105 real, pointer :: tr(:,:,:) => null() !< Array of tracers used in this subroutine [CU ~> conc]
106 !! (ALK tracers use meq m-3 instead of mmol m-3)
107 type(tracer_type), pointer :: tr_ptr => null() !< pointer to tracer inside Tr_reg
108end type marbl_tracer_data
109
110!> The control structure for the MARBL tracer package
111type, public :: marbl_tracers_cs ; private
112 integer :: ntr !< The number of tracers that are actually used.
113 logical :: debug !< If true, write verbose checksums for debugging purposes.
114 logical :: base_bio_on !< Will MARBL use base biotic tracers?
115 logical :: abio_dic_on !< Will MARBL use abiotic DIC / DI14C tracers?
116 logical :: ciso_on !< Will MARBL use isotopic tracers?
117
118 integer :: restore_count !< The number of tracers MARBL is configured to restore
119 logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
120 logical :: use_ice_category_fields !< Forcing will include multiple ice categories for ice_frac and shortwave
121 logical :: request_chl_from_marbl !< MARBL can provide Chl to use in set_pen_shortwave()
122 integer :: ice_ncat !< Number of ice categories when use_ice_category_fields = True
123 real :: ic_min !< Minimum value for tracer initial conditions [CU ~> conc]
124 character(len=200) :: ic_file !< The file in which the age-tracer initial values cam be found.
125 logical :: ongrid !< True if IC_file is already interpolated to MOM grid
126 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the tracer registry
127 type(marbl_tracer_data), dimension(:), allocatable :: tracer_data !< type containing tracer data and pointer
128 !! into tracer registry
129
130 integer, allocatable, dimension(:) :: ind_tr !< Indices returned by aof_set_coupler_flux if it is used and the
131 !! surface tracer concentrations are to be provided to the coupler.
132
133 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
134 !! regulate the timing of diagnostic output.
135 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control structure
136
137 type(vardesc), allocatable :: tr_desc(:) !< Descriptions and metadata for the tracers
138 logical :: tracers_may_reinit !< If true the tracers may be initialized if not found in a restart file
139
140 character(len=200) :: fesedflux_file !< name of [netCDF] file containing iron sediment flux
141 character(len=200) :: fesedfluxred_file !< name of [netCDF] file containing reduced iron sediment flux
142 character(len=200) :: feventflux_file !< name of [netCDF] file containing iron vent flux
143 type(forcing_timeseries_dataset) :: d14c_dataset(3) !< File and time axis information for d14c forcing
144 real, dimension(3) :: d14c_bands !< forcing is organized into bands: [30 N, 90 N]; [30 S, 30 N]; [90 S, 30 S]
145 !! This variable contains D14C for each band [CU ~> conc]
146 integer :: d14c_id !< id for diagnostic field with d14c forcing
147 logical :: read_riv_fluxes !< If true, use river fluxes supplied from an input file.
148 !! This is temporary, we will always read river fluxes
149 type(forcing_timeseries_dataset) :: riv_flux_dataset !< File and time axis information for river fluxes
150 character(len=4) :: restoring_source !< location of tracer restoring data
151 !! valid values: file, none
152 integer :: restoring_nz !< number of levels in tracer restoring file
153 real, allocatable, dimension(:) :: &
154 restoring_z_edges !< The depths of the cell interfaces in the tracer restoring file [Z ~> m]
155 real, allocatable, dimension(:) :: &
156 restoring_dz !< The thickness of the cell layers in the tracer restoring file [H ~> m]
157 integer :: restoring_timescale_nz !< number of levels in tracer restoring timescale file
158 real, allocatable, dimension(:) :: &
159 restoring_timescale_z_edges !< The depths of the cell interfaces in the tracer restoring timescale file [Z ~> m]
160 real, allocatable, dimension(:) :: &
161 restoring_timescale_dz !< The thickness of the cell layers in the tracer restoring timescale file [H ~> m]
162 character(len=14) :: restoring_i_tau_source !< location of inverse restoring timescale data
163 !! valid values: file, grid_dependent
164 character(len=200) :: restoring_file !< name of [netCDF] file containing tracer restoring data
165 type(remapping_cs) :: restoring_remapcs !< Remapping parameters and work arrays for tracer restoring / timescale
166 character(len=200) :: restoring_i_tau_file !< name of [netCDF] file containing inverse restoring timescale
167 character(len=200) :: restoring_i_tau_var_name !< name of field containing inverse restoring timescale
168 character(len=35) :: marbl_settings_file !< name of [text] file containing MARBL settings
169
170 real :: bot_flux_mix_thickness !< for bottom flux -> tendency conversion, assume uniform mixing over
171 !! bottom layer of prescribed thickness [Z ~> m]
172 real :: ibfmt !< Reciprocal of bot_flux_mix_thickness [Z-1 ~> m-1]
173
174 type(temp_marbl_diag), allocatable :: surface_flux_diags(:) !< collect surface flux diagnostics from all columns
175 !! before posting
176 type(temp_marbl_diag), allocatable :: interior_tendency_diags(:) !< collect tendency diagnostics from all columns
177 !! before posting
178 type(saved_state_for_marbl_type), allocatable :: surface_flux_saved_state(:) !< surface_flux saved state
179 type(saved_state_for_marbl_type), allocatable :: interior_tendency_saved_state(:) !< interior_tendency saved state
180
181 ! TODO: If we can post data column by column, all we need are integer arrays for ids
182 ! integer, allocatable :: id_surface_flux_diags(:) !< array of indices for surface_flux diagnostics
183 ! integer, allocatable :: id_interior_tendency_diags(:) !< array of indices for interior_tendency diagnostics
184
185 type(tracer_ind_type) :: tracer_inds !< Indices to tracers that will have river fluxes added to STF
186
187 !> Need to store global output from both marbl_instance%surface_flux_compute() and
188 !! marbl_instance%interior_tendency_compute(). For the former, just need id to register
189 !! because we already copy data into CS%STF; latter requires copying data and indices
190 !! so currently using temp_MARBL_diag for that.
191 integer, allocatable :: id_surface_flux_out(:) !< register_diag indices for surface_flux output
192 integer, allocatable :: id_surface_flux_from_salt_flux(:) !< register_diag indices for surface_flux from salt_flux
193 type(temp_marbl_diag), allocatable :: interior_tendency_out(:) !< collect interior tendencies for diagnostic output
194 type(temp_marbl_diag), allocatable :: interior_tendency_out_zint(:) !< vertical integral of interior tendencies
195 !! (full column)
196 type(temp_marbl_diag), allocatable :: interior_tendency_out_zint_100m(:) !< vertical integral of interior tendencies
197 !! (top 100m)
198 integer :: bot_flux_to_tend_id !< register_diag index for BOT_FLUX_TO_TEND
199 integer, allocatable :: fracr_cat_id(:) !< register_diag index for per-category ice fraction
200 integer, allocatable :: qsw_cat_id(:) !< register_diag index for per-category shortwave
201
202 real :: dic_salt_ratio !< ratio to convert salt surface flux to DIC surface flux [conc ppt-1]
203 real :: alk_salt_ratio !< ratio to convert salt surface flux to ALK surface flux [conc ppt-1]
204
205 real, allocatable :: stf(:,:,:) !< surface fluxes returned from MARBL to use in tracer_vertdiff
206 !! (dims: i, j, tracer) [conc Z T-1 ~> conc m s-1]
207 real, pointer :: sfo(:,:,:) => null() !< surface flux output returned from MARBL for use in GCM
208 !! e.g. CO2 flux to pass to atmosphere (dims: i, j, num_sfo)
209 !! Units vary based on index of num_sfo dimension
210 real, pointer :: ito(:,:,:,:) => null() !< interior tendency output returned from MARBL for use in GCM
211 !! e.g. total chlorophyll to use in shortwave penetration
212 !! (dims: i, j, k, num_ito)
213 !! Units vary based on index of num_ito dimension
214
215 integer :: u10_sqr_ind !< index of MARBL forcing field array to copy 10-m wind (squared) into
216 integer :: sss_ind !< index of MARBL forcing field array to copy sea surface salinity into
217 integer :: sst_ind !< index of MARBL forcing field array to copy sea surface temperature into
218 integer :: ifrac_ind !< index of MARBL forcing field array to copy ice fraction into
219 integer :: dust_dep_ind !< index of MARBL forcing field array to copy dust flux into
220 integer :: fe_dep_ind !< index of MARBL forcing field array to copy iron flux into
221 integer :: nox_flux_ind !< index of MARBL forcing field array to copy NOx flux into
222 integer :: nhy_flux_ind !< index of MARBL forcing field array to copy NHy flux into
223 integer :: atmpress_ind !< index of MARBL forcing field array to copy atmospheric pressure into
224 integer :: xco2_ind !< index of MARBL forcing field array to copy CO2 flux into
225 integer :: xco2_alt_ind !< index of MARBL forcing field array to copy CO2 flux (alternate CO2) into
226 integer :: d14c_ind !< index of MARBL forcing field array to copy d14C into
227
228 !> external_field types for river fluxes (added to surface fluxes)
229 type(external_field) :: id_din_riv !< id for time_interp_external.
230 type(external_field) :: id_don_riv !< id for time_interp_external.
231 type(external_field) :: id_dip_riv !< id for time_interp_external.
232 type(external_field) :: id_dop_riv !< id for time_interp_external.
233 type(external_field) :: id_dsi_riv !< id for time_interp_external.
234 type(external_field) :: id_dfe_riv !< id for time_interp_external.
235 type(external_field) :: id_dic_riv !< id for time_interp_external.
236 type(external_field) :: id_alk_riv !< id for time_interp_external.
237 type(external_field) :: id_doc_riv !< id for time_interp_external.
238
239 !> external_field type for d14c (needed if abio_dic_on is True)
240 type(external_field) :: id_d14c(3) !< id for time_interp_external.
241
242 !> Indices for river fluxes (diagnostics)
243 integer :: no3_riv_flux !< NO3 riverine flux
244 integer :: po4_riv_flux !< PO4 riverine flux
245 integer :: don_riv_flux !< DON riverine flux
246 integer :: donr_riv_flux !< DONr riverine flux
247 integer :: dop_riv_flux !< DOP riverine flux
248 integer :: dopr_riv_flux !< DOPr riverine flux
249 integer :: sio3_riv_flux !< SiO3 riverine flux
250 integer :: fe_riv_flux !< Fe riverine flux
251 integer :: doc_riv_flux !< DOC riverine flux
252 integer :: docr_riv_flux !< DOCr riverine flux
253 integer :: alk_riv_flux !< ALK riverine flux
254 integer :: alk_alt_co2_riv_flux !< ALK (alternate CO2) riverine flux
255 integer :: dic_riv_flux !< DIC riverine flux
256 integer :: dic_alt_co2_riv_flux !< DIC (alternate CO2) riverine flux
257
258 !> Indices for forcing fields required to compute interior tendencies
259 integer :: dustflux_ind !< index of MARBL forcing field array to copy dust flux into
260 integer :: par_col_frac_ind !< index of MARBL forcing field array to copy PAR column fraction into
261 integer :: surf_shortwave_ind !< index of MARBL forcing field array to copy surface shortwave into
262 integer :: potemp_ind !< index of MARBL forcing field array to copy potential temperature into
263 integer :: salinity_ind !< index of MARBL forcing field array to copy salinity into
264 integer :: pressure_ind !< index of MARBL forcing field array to copy pressure into
265 integer :: fesedflux_ind !< index of MARBL forcing field array to copy iron sediment flux into
266 integer :: fesedfluxred_ind !< index of MARBL forcing field array to copy reduced iron sediment flux into
267 integer :: feventflux_ind !< index of MARBL forcing field array to copy iron vent flux into
268 integer :: o2_scalef_ind !< index of MARBL forcing field array to copy O2 scale length into
269 integer :: remin_scalef_ind !< index of MARBL forcing field array to copy remin scale length into
270 type(external_field), allocatable :: id_tracer_restoring(:) !< id number for time_interp_external
271 integer, allocatable :: tracer_restoring_ind(:) !< index of MARBL forcing field to copy
272 !! per-tracer restoring field into
273 integer, allocatable :: tracer_i_tau_ind(:) !< index of MARBL forcing field to copy per-tracer
274 !! inverse restoring timescale into
275
276 !> Memory for storing river fluxes, tracer restoring fields, and abiotic forcing
277 real, allocatable :: d14c(:,:) !< d14c forcing for abiotic DIC and carbon isotope tracer modules
278 !! [mmol m-3 s-1]
279 real, allocatable :: riv_fluxes(:,:,:) !< river flux forcing for applyTracerBoundaryFluxesInOut
280 !! (needs to be time-integrated when passed to function!)
281 !! (dims: i, j, tracer) [conc m s-1]
282 character(len=15), allocatable :: tracer_restoring_varname(:) !< name of variable being restored
283 real, allocatable :: i_tau(:,:,:) !< inverse restoring timescale for marbl tracers (dims: i, j, k) [s-1]
284 real, allocatable, dimension(:,:,:,:) :: restoring_in !< Restoring fields read from file
285 !! (dims: i, j, restoring_nz, restoring_cnt) [tracer units]
286
287 !> Number of surface flux outputs as well as specific indices for each one
288 integer :: sfo_cnt !< number of surface flux outputs from MARBL
289 integer :: ito_cnt !< number of interior tendency outputs from MARBL
290 integer :: flux_co2_ind !< index to co2 flux surface flux output
291 integer :: total_chl_ind !< index to total chlorophyll interior tendency output
292
293 ! TODO: create generic 3D forcing input type to read z coordinate + values
294 real :: fesedflux_scale_factor !< scale factor for iron sediment flux [mmol umol-1 d s-1]
295 integer :: fesedflux_nz !< number of levels in iron sediment flux file
296 real, allocatable, dimension(:,:,:) :: fesedflux_in !< Field to read iron sediment flux into [conc m s-1]
297 real, allocatable, dimension(:,:,:) :: fesedfluxred_in !< Field to read reduced iron sediment flux into [conc m s-1]
298 real, allocatable, dimension(:,:,:) :: feventflux_in !< Field to read iron vent flux into [conc m s-1]
299 real, allocatable, dimension(:) :: &
300 fesedflux_z_edges !< The depths of the cell interfaces in the input data [Z ~> m]
301 ! TODO: this thickness does not need to be 3D, but it is easier to make thickness 0
302 ! below the surface on a per-column basis (could save memory by storing 1D
303 ! thickness from file and then computing a second 1D thickness array in (i,j) loop)
304 real, allocatable, dimension(:,:,:) :: &
305 fesedflux_dz !< The thickness of the cell layers in the input data [H ~> m]
306end type marbl_tracers_cs
307
308! Module parameters
309real, parameter :: atm_per_pa = 1./101325. !< convert from Pa -> atm [atm Pa-1]
310
311contains
312
313!> This subroutine is used to read marbl_in, configure MARBL accordingly, and then
314!! call MARBL's initialization routine
315subroutine configure_marbl_tracers(GV, US, param_file, CS)
316 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
317 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
318 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
319 type(marbl_tracers_cs), pointer :: CS !< A pointer that is set to point to the control
320 !! structure for this module
321
322# include "version_variable.h"
323 character(len=40) :: mdl = "MARBL_tracers" ! This module's name.
324 character(len=256) :: log_message
325 character(len=256) :: marbl_in_line(1)
326 character(len=256) :: forcing_sname, field_source
327 integer :: m, n, nz, marbl_settings_in, read_error, I_tau_count, fi
328 logical :: chl_from_file, forcing_processed
329 nz = gv%ke
330 marbl_settings_in = 615
331
332 ! (1) Read parameters necessary for general setup of MARBL
333 call log_version(param_file, mdl, version, "")
334 call get_param(param_file, mdl, "DEBUG", cs%debug, "If true, write out verbose debugging data.", &
335 default=.false., debuggingparam=.true.)
336 call get_param(param_file, mdl, "MARBL_IC_MIN_VAL", cs%IC_min, &
337 "Minimum value of tracer initial conditions (set to 1e-100 for dim scaling tests)", &
338 default=0., units="tracer units")
339 call get_param(param_file, mdl, "MARBL_SETTINGS_FILE", cs%marbl_settings_file, &
340 "The name of a file from which to read the run-time settings for MARBL.", default="marbl_in")
341 call get_param(param_file, mdl, "BOT_FLUX_MIX_THICKNESS", cs%bot_flux_mix_thickness, &
342 "Bottom fluxes are uniformly mixed over layer of this thickness", default=1., units="m", &
343 scale=us%m_to_Z)
344 call get_param(param_file, mdl, "USE_ICE_CATEGORIES", cs%use_ice_category_fields, &
345 "If true, allocate memory for shortwave and ice fraction split by ice thickness category.", &
346 default=.false.)
347 call get_param(param_file, mdl, "ICE_NCAT", cs%ice_ncat, &
348 "Number of ice thickness categories in shortwave and ice fraction forcings.", default=0)
349 cs%Ibfmt = 1. / cs%bot_flux_mix_thickness
350
351 if (cs%use_ice_category_fields .and. (cs%ice_ncat == 0)) &
352 call mom_error(fatal, &
353 "Can not configure MARBL to use multiple ice categories without ice_ncat present")
354
355 ! (2) Read marbl settings file and call put_setting()
356 ! (2a) only master task opens file
357 if (is_root_pe()) then
358 ! read the marbl_in into buffer
359 open(unit=marbl_settings_in, file=cs%marbl_settings_file, iostat=read_error)
360 if (read_error .ne. 0) then
361 write(log_message, '(A, I0, 2A)') "IO ERROR ", read_error, " opening namelist file : ", &
362 trim(cs%marbl_settings_file)
363 call mom_error(fatal, log_message)
364 endif
365 endif
366
367 ! (2b) master task reads file and broadcasts line-by-line
368 marbl_in_line = ''
369 do
370 ! i. Read next line on master, iostat value out
371 ! (Exit loop if read is not successful; either read error or end of file)
372 if (is_root_pe()) read(marbl_settings_in, "(A)", iostat=read_error) marbl_in_line(1)
373 call broadcast(read_error, root_pe())
374 if (read_error .ne. 0) exit
375
376 ! ii. Broadcast line just read in on root PE to all tasks
377 call broadcast(marbl_in_line, 256, root_pe())
378
379 ! iii. All tasks call put_setting (TODO: openMP blocks?)
380 call marbl_instances%put_setting(marbl_in_line(1))
381 enddo
382
383 ! (2c) we should always reach the EOF to capture the entire file...
384 if (.not. is_iostat_end(read_error)) then
385 write(log_message, '(3A, I0)') "IO ERROR reading ", trim(cs%marbl_settings_file), ": ", &
386 read_error
387 call mom_error(fatal, log_message)
388 else
389 if (is_root_pe()) then
390 write(log_message, '(3A)') "Read '", trim(cs%marbl_settings_file), "' until EOF."
391 call mom_error(note, log_message)
392 endif
393 endif
394 if (is_root_pe()) close(marbl_settings_in)
395
396 ! (3) Initialize MARBL and configure MOM6 accordingly
397
398 ! (3a) call marbl%init()
399 ! TODO: We want to strip gcm_delta_z, gcm_zw, and gcm_zt values out of
400 ! init because MOM updates them every time step / every column
401 call marbl_instances%init(gcm_num_levels = nz, gcm_num_par_subcols = cs%ice_ncat + 1, &
402 gcm_num_elements_surface_flux = 1, & ! FIXME: change to number of grid cells on MPI task
403 gcm_delta_z = gv%sInterface(2:nz+1) - gv%sInterface(1:nz), gcm_zw = gv%sInterface(2:nz+1), &
404 gcm_zt = gv%sLayer, unit_system_opt = "mks", lgcm_has_global_ops = .false.) ! FIXME: add global ops
405 ! Regardless of vertical grid, MOM6 will always use GV%ke levels in all columns
406 marbl_instances%domain%kmt = gv%ke
407 if (marbl_instances%StatusLog%labort_marbl) &
408 call marbl_instances%StatusLog%log_error_trace("MARBL_instances%init", &
409 "configure_MARBL_tracers")
410 call print_marbl_log(marbl_instances%StatusLog)
411 call marbl_instances%StatusLog%erase()
412 cs%ntr = size(marbl_instances%tracer_metadata)
413 call marbl_instances%get_setting('base_bio_on', cs%base_bio_on)
414 call marbl_instances%get_setting('abio_dic_on', cs%abio_dic_on)
415 call marbl_instances%get_setting('ciso_on', cs%ciso_on)
416
417 ! (3b) Read parameters that depend on how MARBL is configured
418 if (cs%base_bio_on) then
419 call get_param(param_file, mdl, "CHL_FROM_FILE", chl_from_file, &
420 "If true, chl_a is read from a file.", default=.true.)
421 cs%request_Chl_from_MARBL = (.not. chl_from_file)
422 else
423 cs%request_Chl_from_MARBL = .false.
424 endif
425
426 ! (4) Request fields needed by MOM6
427 cs%sfo_cnt = 0
428 cs%ito_cnt = 0
429 cs%flux_co2_ind = -1
430 cs%total_Chl_ind = -1
431
432 if (cs%base_bio_on) then
433 ! CO2 Flux to the atmosphere
434 call marbl_instances%add_output_for_GCM(num_elements=1, field_name="flux_co2", &
435 output_id=cs%flux_co2_ind, field_source=field_source)
436 if (trim(field_source) == "surface_flux") then
437 cs%sfo_cnt = cs%sfo_cnt + 1
438 else if (trim(field_source) == "interior_tendency") then
439 cs%ito_cnt = cs%ito_cnt + 1
440 endif
441
442 ! Total 3D Chlorophyll
443 call marbl_instances%add_output_for_GCM(num_elements=1, num_levels=nz, field_name="total_Chl", &
444 output_id=cs%total_Chl_ind, field_source=field_source)
445 if (trim(field_source) == "surface_flux") then
446 cs%sfo_cnt = cs%sfo_cnt + 1
447 else if (trim(field_source) == "interior_tendency") then
448 cs%ito_cnt = cs%ito_cnt + 1
449 endif
450 endif
451
452 ! (5) Initialize forcing fields
453 ! i. store all surface forcing indices
454 cs%u10_sqr_ind = -1
455 cs%sss_ind = -1
456 cs%sst_ind = -1
457 cs%ifrac_ind = -1
458 cs%dust_dep_ind = -1
459 cs%fe_dep_ind = -1
460 cs%nox_flux_ind = -1
461 cs%nhy_flux_ind = -1
462 cs%atmpress_ind = -1
463 cs%xco2_ind = -1
464 cs%xco2_alt_ind = -1
465 do m=1,size(marbl_instances%surface_flux_forcings)
466 select case (trim(marbl_instances%surface_flux_forcings(m)%metadata%varname))
467 case('u10_sqr')
468 cs%u10_sqr_ind = m
469 case('sss')
470 cs%sss_ind = m
471 case('sst')
472 cs%sst_ind = m
473 case('Ice Fraction')
474 cs%ifrac_ind = m
475 case('Dust Flux')
476 cs%dust_dep_ind = m
477 case('Iron Flux')
478 cs%fe_dep_ind = m
479 case('NOx Flux')
480 cs%nox_flux_ind = m
481 case('NHy Flux')
482 cs%nhy_flux_ind = m
483 case('Atmospheric Pressure')
484 cs%atmpress_ind = m
485 case('xco2')
486 cs%xco2_ind = m
487 case('xco2_alt_co2')
488 cs%xco2_alt_ind = m
489 case('d14c')
490 cs%d14c_ind = m
491 case DEFAULT
492 write(log_message, "(A,1X,A)") &
493 trim(marbl_instances%surface_flux_forcings(m)%metadata%varname), &
494 'is not a valid surface flux forcing field name.'
495 call mom_error(fatal, log_message)
496 end select
497 enddo
498
499 ! ii. store all interior forcing indices
500 cs%dustflux_ind = -1
501 cs%PAR_col_frac_ind = -1
502 cs%surf_shortwave_ind = -1
503 cs%potemp_ind = -1
504 cs%salinity_ind = -1
505 cs%pressure_ind = -1
506 cs%fesedflux_ind = -1
507 cs%fesedfluxred_ind = -1
508 cs%feventflux_ind = -1
509 cs%o2_scalef_ind = -1
510 cs%remin_scalef_ind = -1
511 cs%d14c_ind = -1
512 allocate(cs%id_tracer_restoring(cs%ntr))
513 allocate(cs%tracer_restoring_varname(cs%ntr), source=' ') ! gfortran 13.2 bug?
514 ! source = '' does not blank out strings
515 allocate(cs%tracer_restoring_ind(cs%ntr), source=-1)
516 allocate(cs%tracer_I_tau_ind(cs%ntr), source=-1)
517 cs%restore_count = 0
518 i_tau_count = 0
519 do m=1,size(marbl_instances%interior_tendency_forcings)
520 select case (trim(marbl_instances%interior_tendency_forcings(m)%metadata%varname))
521 case('Dust Flux')
522 cs%dustflux_ind = m
523 case('PAR Column Fraction')
524 cs%PAR_col_frac_ind = m
525 case('Surface Shortwave')
526 cs%surf_shortwave_ind = m
527 case('Potential Temperature')
528 cs%potemp_ind = m
529 case('Salinity')
530 cs%salinity_ind = m
531 case('Pressure')
532 cs%pressure_ind = m
533 case('Iron Sediment Flux')
534 cs%fesedflux_ind = m
535 case('Iron Red Sediment Flux')
536 cs%fesedfluxred_ind = m
537 case('Iron Vent Flux')
538 cs%feventflux_ind = m
539 case('O2 Consumption Scale Factor')
540 cs%o2_scalef_ind = m
541 case('Particulate Remin Scale Factor')
542 cs%remin_scalef_ind = m
543 case DEFAULT
544 ! fi stands for forcing_index
545 fi = index(marbl_instances%interior_tendency_forcings(m)%metadata%varname, &
546 'Restoring Field')
547 if (fi > 0) then
548 cs%restore_count = cs%restore_count + 1
549 cs%tracer_restoring_ind(cs%restore_count) = m
550 cs%tracer_restoring_varname(cs%restore_count) = &
551 marbl_instances%interior_tendency_forcings(m)%metadata%varname(1:fi-2)
552 else
553 fi = index(marbl_instances%interior_tendency_forcings(m)%metadata%varname, &
554 'Restoring Inverse Timescale')
555 if (fi > 0) then
556 i_tau_count = i_tau_count + 1
557 cs%tracer_I_tau_ind(i_tau_count) = m
558 else
559 write(log_message, "(A,1X,A)") &
560 trim(marbl_instances%interior_tendency_forcings(m)%metadata%varname), &
561 'is not a valid interior tendency forcing field name.'
562 call mom_error(fatal, log_message)
563 endif
564 endif
565 end select
566 enddo
567end subroutine configure_marbl_tracers
568
569!> This subroutine is used to register tracer fields and subroutines
570!! to be used with MOM.
571function register_marbl_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS, MARBL_computes_chl)
572 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
573 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
574 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
575 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
576 type(marbl_tracers_cs), pointer :: cs !< A pointer that is set to point to the control
577 !! structure for this module
578 type(tracer_registry_type), pointer :: tr_reg !< A pointer that is set to point to the control
579 !! structure for the tracer advection and diffusion module.
580 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control struct
581 logical, intent(out) :: marbl_computes_chl !< If MARBL is computing chlorophyll, MOM
582 !! may use it to compute SW penetration
583
584! Local variables
585! This include declares and sets the variable "version".
586# include "version_variable.h"
587 character(len=40) :: mdl = "MARBL_tracers" ! This module's name.
588 character(len=256) :: log_message
589 character(len=200) :: inputdir ! The directory where the input files are.
590 character(len=48) :: var_name ! The variable's name.
591 character(len=128) :: desc_name ! The variable's descriptor.
592 character(len=48) :: units ! The variable's units.
593 character(len=96) :: file_name ! file name for d14c (looped over three bands)
594 real, pointer :: tr_ptr(:,:,:) => null() ! Pointer to 3D tracer array [CU ~> conc]
595 ! (ALK tracers use meq m-3 instead of mmol m-3)
596 integer :: forcing_file_start_year
597 integer :: forcing_file_end_year
598 integer :: forcing_file_data_ref_year
599 integer :: forcing_file_model_ref_year
600 integer :: forcing_file_forcing_year
601 logical :: register_marbl_tracers
602 ! read_Z_edges() has several mandatory arguments that we do not use given our expectation
603 ! of how the file being read in was created
604 logical :: z_edges_has_edges
605 logical :: z_edges_use_missing
606 real :: z_edges_missing ! required argument for read_Z_edges() [CU ~> conc]
607 integer :: isd, ied, jsd, jed, nz, m, k, kbot
608 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
609
610 if (associated(cs)) then
611 call mom_error(warning, "register_MARBL_tracers called with an associated control structure.")
612 return
613 endif
614 allocate(cs)
615
616 call configure_marbl_tracers(gv, us, param_file, cs)
617 marbl_computes_chl = cs%base_bio_on
618
619 ! Read all relevant parameters and write them to the model log.
620 call log_version(param_file, mdl, version, "")
621 ! ** Input directory
622 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
623 ! ** Tracer initial conditions
624 call get_param(param_file, mdl, "MARBL_TRACERS_IC_FILE", cs%IC_file, &
625 "The file in which the MARBL tracers initial values can be found.", &
626 default="ecosys_jan_IC_omip_latlon_1x1_180W_c230331.nc")
627 if (scan(cs%IC_file,'/') == 0) then
628 ! Add the directory if CS%IC_file is not already a complete path.
629 cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)
630 call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_IC_FILE", cs%IC_file)
631 endif
632 call get_param(param_file, mdl, "MARBL_TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
633 "If true, tracers may go through the initialization code if they are not found in the "//&
634 "restart files. Otherwise it is a fatal error if tracers are not found in the "//&
635 "restart files of a restarted run.", default=.false.)
636 call get_param(param_file, mdl, "MARBL_TRACERS_INIT_VERTICAL_REMAP_ONLY", cs%ongrid, &
637 "If true, initial conditions are on the model horizontal grid. Extrapolation over " //&
638 "missing ocean values is done using an ICE-9 procedure with vertical ALE remapping .", &
639 default=.false.)
640 if (cs%base_bio_on) then
641 ! ** FESEDFLUX
642 call get_param(param_file, mdl, "MARBL_FESEDFLUX_FILE", cs%fesedflux_file, &
643 "The file in which the iron sediment flux forcing field can be found.", &
644 default="fesedflux.nc")
645 if (scan(cs%fesedflux_file,'/') == 0) then
646 ! Add the directory if CS%fesedflux_file is not already a complete path.
647 cs%fesedflux_file = trim(slasher(inputdir))//trim(cs%fesedflux_file)
648 call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_FESEDFLUX_FILE", cs%fesedflux_file)
649 endif
650 ! ** FESEDFLUXRED
651 call get_param(param_file, mdl, "MARBL_FESEDFLUXRED_FILE", cs%fesedfluxred_file, &
652 "The file in which the iron sediment flux forcing field can be found.", &
653 default="fesedfluxred.nc")
654 if (scan(cs%fesedfluxred_file,'/') == 0) then
655 ! Add the directory if CS%fesedflux_file is not already a complete path.
656 cs%fesedfluxred_file = trim(slasher(inputdir))//trim(cs%fesedfluxred_file)
657 call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_FESEDFLUXRED_FILE", cs%fesedfluxred_file)
658 endif
659 ! ** FEVENTFLUX
660 call get_param(param_file, mdl, "MARBL_FEVENTFLUX_FILE", cs%feventflux_file, &
661 "The file in which the iron vent flux forcing field can be found.", &
662 default="feventflux.nc")
663 if (scan(cs%feventflux_file,'/') == 0) then
664 ! Add the directory if CS%feventflux_file is not already a complete path.
665 cs%feventflux_file = trim(slasher(inputdir))//trim(cs%feventflux_file)
666 call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_FEVENTFLUX_FILE", cs%feventflux_file)
667 endif
668 ! ** Scale factor for FESEDFLUX
669 call get_param(param_file, mdl, "MARBL_FESEDFLUX_SCALE_FACTOR", cs%fesedflux_scale_factor, &
670 "Conversion factor between FESEDFLUX file units and MARBL units", &
671 units="umol m-2 d-1 -> mmol m-2 s-1", default=0.001/86400.)
672
673 ! ** River fluxes
674 call get_param(param_file, mdl, "READ_RIV_FLUXES", cs%read_riv_fluxes, &
675 "If true, use river fluxes supplied from an input file", default=.true.)
676 if (cs%read_riv_fluxes) then
677 call get_param(param_file, mdl, "RIV_FLUX_FILE", cs%riv_flux_dataset%file_name, &
678 "The file in which the river fluxes can be found", &
679 default="riv_nut.gnews_gnm.JRA025m_to_tx0.66v1_nnsm_e333r100_190910.20210405.nc")
680 ! call get_param(param_file, mdl, "RIV_FLUX_OFFSET_YEAR", CS%riv)
681 if (scan(cs%riv_flux_dataset%file_name,'/') == 0) then
682 ! CS%riv_flux_dataset%file_name = trim(inputdir) // trim(CS%riv_flux_dataset%file_name)
683 cs%riv_flux_dataset%file_name = trim(slasher(inputdir)) //&
684 trim(cs%riv_flux_dataset%file_name)
685 call log_param(param_file, mdl, "INPUTDIR/RIV_FLUX_FILE", cs%riv_flux_dataset%file_name)
686 endif
687 call get_param(param_file, mdl, "RIV_FLUX_L_TIME_VARYING", &
688 cs%riv_flux_dataset%l_time_varying, &
689 ".true. for time-varying forcing, .false. for static forcing", default=.false.)
690 if (cs%riv_flux_dataset%l_time_varying) then
691 call get_param(param_file, mdl, "RIV_FLUX_FILE_START_YEAR", forcing_file_start_year, &
692 "First year of data to read in RIV_FLUX_FILE", default=1900)
693 call get_param(param_file, mdl, "RIV_FLUX_FILE_END_YEAR", forcing_file_end_year, &
694 "Last year of data to read in RIV_FLUX_FILE", default=2000)
695 call get_param(param_file, mdl, "RIV_FLUX_FILE_DATA_REF_YEAR", forcing_file_data_ref_year, &
696 "Align this year in RIV_FLUX_FILE with RIV_FLUX_FILE_MODEL_REF_YEAR in model", &
697 default=1900)
698 call get_param(param_file, mdl, "RIV_FLUX_FILE_MODEL_REF_YEAR", &
699 forcing_file_model_ref_year, &
700 "Align this year in model with RIV_FLUX_FILE_DATA_REF_YEAR in RIV_FLUX_FILE", &
701 default=1)
702 else
703 call get_param(param_file, mdl, "RIV_FLUX_FORCING_YEAR", forcing_file_forcing_year, &
704 "Year from RIV_FLUX_FILE to use for forcing", default=1900)
705 endif
706 call forcing_timeseries_set_time_type_vars(forcing_file_start_year, forcing_file_end_year, &
707 forcing_file_data_ref_year, forcing_file_model_ref_year, forcing_file_forcing_year, &
708 cs%riv_flux_dataset)
709 endif
710 endif
711
712 if (cs%abio_dic_on) then
713 call get_param(param_file, mdl, "D14C_L_TIME_VARYING", cs%d14c_dataset(1)%l_time_varying, &
714 ".true. for time-varying forcing, .false. for static forcing", default=.false.)
715 cs%d14c_dataset(2)%l_time_varying = cs%d14c_dataset(1)%l_time_varying
716 cs%d14c_dataset(3)%l_time_varying = cs%d14c_dataset(1)%l_time_varying
717 if (cs%d14c_dataset(1)%l_time_varying) then
718 call get_param(param_file, mdl, "D14C_FILE_START_YEAR", forcing_file_start_year, &
719 "First year of data to read in D14C_FILE", default=1850)
720 call get_param(param_file, mdl, "D14C_FILE_END_YEAR", forcing_file_end_year, &
721 "Last year of data to read in D14C_FILE", default=2015)
722 call get_param(param_file, mdl, "D14C_FILE_DATA_REF_YEAR", forcing_file_data_ref_year, &
723 "Align this year in D14C_FILE with D14C_FILE_MODEL_REF_YEAR in model", default=1850)
724 call get_param(param_file, mdl, "D14C_FILE_MODEL_REF_YEAR", forcing_file_model_ref_year, &
725 "Align this year in model with D14C_FILE_DATA_REF_YEAR in D14C_FILE", default=1)
726 else
727 call get_param(param_file, mdl, "D14C_FORCING_YEAR", forcing_file_forcing_year, &
728 "Year from D14C_FILE to use for forcing", default=1850)
729 endif
730 do m=1,3
731 write(var_name, "(A,I0)") "MARBL_D14C_FILE_", m
732 write(file_name, "(A,I0,A)") "atm_delta_C14_CMIP6_sector", m, &
733 "_global_1850-2015_yearly_v2.0_c240202.nc"
734 call get_param(param_file, mdl, var_name, cs%d14c_dataset(m)%file_name, &
735 "The file in which the d14c forcing field can be found.", default=file_name)
736 call forcing_timeseries_set_time_type_vars(forcing_file_start_year, forcing_file_end_year, &
737 forcing_file_data_ref_year, forcing_file_model_ref_year, forcing_file_forcing_year, &
738 cs%d14c_dataset(m))
739 if (scan(cs%d14c_dataset(m)%file_name,'/') == 0) then
740 ! Add the directory if CS%d14c_dataset%file_name is not already a complete path.
741 cs%d14c_dataset(m)%file_name = trim(slasher(inputdir))//trim(cs%d14c_dataset(m)%file_name)
742 call log_param(param_file, mdl, "INPUTDIR/D14C_FILE", cs%d14c_dataset(m)%file_name)
743 endif
744 enddo
745 endif
746
747 call get_param(param_file, mdl, "DIC_SALT_RATIO", cs%DIC_salt_ratio, &
748 "Ratio to convert salt surface flux to DIC surface flux", units="conc ppt-1", &
749 default=64.0)
750 call get_param(param_file, mdl, "ALK_SALT_RATIO", cs%ALK_salt_ratio, &
751 "Ratio to convert salt surface flux to ALK surface flux", units="conc ppt-1", &
752 default=70.0)
753
754 ! ** Tracer Restoring
755 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_SOURCE", cs%restoring_source, &
756 "Source of data for restoring MARBL tracers", default="none")
757 select case(cs%restoring_source)
758 case("none")
759 case("file")
760 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_FILE", cs%restoring_file, &
761 "File containing fields to restore MARBL tracers towards")
762 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_SOURCE", &
763 cs%restoring_I_tau_source, "Source of data for inverse timescale for restoring MARBL tracers")
764
765 ! Initialize remapping type
766 call initialize_remapping(cs%restoring_remapCS, 'PCM', boundary_extrapolation=.false., answer_date=99991231)
767
768 ! Set up array for thicknesses in restoring file
769 call read_z_edges(cs%restoring_file, "PO4", cs%restoring_z_edges, cs%restoring_nz, &
770 z_edges_has_edges, z_edges_use_missing, z_edges_missing, scale=us%m_to_Z, &
771 missing_scale=1.0)
772 allocate(cs%restoring_dz(cs%restoring_nz))
773 do k=cs%restoring_nz,1,-1
774 kbot = k + 1 ! level k is between z(k) and z(k+1)
775 cs%restoring_dz(k) = (cs%restoring_z_edges(k) - cs%restoring_z_edges(kbot)) * gv%Z_to_H
776 enddo
777
778 select case(cs%restoring_I_tau_source)
779 case("file")
780 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_FILE", &
781 cs%restoring_I_tau_file, &
782 "File containing the inverse timescale for restoring MARBL tracers")
783 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_VAR_NAME", &
784 cs%restoring_I_tau_var_name, &
785 "Field containing the inverse timescale for restoring MARBL tracers", &
786 default="I_TAU")
787 ! Set up array for thicknesses in restoring timescale file
788 call read_z_edges(cs%restoring_I_tau_file, cs%restoring_I_tau_var_name, cs%restoring_timescale_z_edges, &
789 cs%restoring_timescale_nz, z_edges_has_edges, z_edges_use_missing, z_edges_missing, scale=us%m_to_Z, &
790 missing_scale=1.0)
791 allocate(cs%restoring_timescale_dz(cs%restoring_timescale_nz))
792 do k=cs%restoring_timescale_nz,1,-1
793 kbot = k + 1 ! level k is between z(k) and z(k+1)
794 cs%restoring_timescale_dz(k) = (cs%restoring_timescale_z_edges(k) - &
795 cs%restoring_timescale_z_edges(kbot)) * gv%Z_to_H
796 enddo
797 case DEFAULT
798 write(log_message, "(3A)") "'", trim(cs%restoring_I_tau_source), &
799 "' is not a valid option for MARBL_TRACER_RESTORING_I_TAU_SOURCE"
800 call mom_error(fatal, log_message)
801 end select
802 case DEFAULT
803 write(log_message, "(3A)") "'", trim(cs%restoring_source), &
804 "' is not a valid option for MARBL_TRACER_RESTORING_SOURCE"
805 call mom_error(fatal, log_message)
806 end select
807
808 allocate(cs%ind_tr(cs%ntr))
809 allocate(cs%tr_desc(cs%ntr))
810 allocate(cs%tracer_data(cs%ntr))
811
812 do m=1,cs%ntr
813 allocate(cs%tracer_data(m)%tr(isd:ied,jsd:jed,nz), source=0.0)
814 write(var_name(:),'(A)') trim(marbl_instances%tracer_metadata(m)%short_name)
815 write(desc_name(:),'(A)') trim(marbl_instances%tracer_metadata(m)%long_name)
816 write(units(:),'(A)') trim(marbl_instances%tracer_metadata(m)%units)
817 cs%tr_desc(m) = var_desc(trim(var_name), trim(units), trim(desc_name), caller=mdl)
818
819 ! This is needed to force the compiler not to do a copy in the registration
820 ! calls. Curses on the designers and implementers of Fortran90.
821 tr_ptr => cs%tracer_data(m)%tr(:,:,:)
822 call query_vardesc(cs%tr_desc(m), name=var_name, &
823 caller="register_MARBL_tracers")
824 ! Register the tracer for horizontal advection, diffusion, and restarts.
825 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, units = units, &
826 tr_desc=cs%tr_desc(m), registry_diags=.true., &
827 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit, &
828 tr_out=cs%tracer_data(m)%tr_ptr)
829
830 ! Set coupled_tracers to be true (hard-coded above) to provide the surface
831 ! values to the coupler (if any). This is meta-code and its arguments will
832 ! currently (deliberately) give fatal errors if it is used.
833 if (cs%coupled_tracers) &
834 cs%ind_tr(m) = aof_set_coupler_flux(trim(var_name)//'_flux', &
835 flux_type=' ', implementation=' ', caller="register_MARBL_tracers")
836 enddo
837
838 ! Set up memory for saved state
839 call setup_saved_state(marbl_instances%surface_flux_saved_state, hi, gv, restart_cs, &
840 cs%tracers_may_reinit, cs%surface_flux_saved_state)
841 call setup_saved_state(marbl_instances%interior_tendency_saved_state, hi, gv, restart_cs, &
842 cs%tracers_may_reinit, cs%interior_tendency_saved_state)
843
844 ! Set up memory for additional output from MARBL and add to restart files
845 allocate(cs%SFO(szi_(hi), szj_(hi), cs%sfo_cnt), &
846 cs%ITO(szi_(hi), szj_(hi), szk_(gv), cs%ito_cnt), &
847 source=0.0)
848
849 do m=1,cs%sfo_cnt
850 write(var_name, "(2A)") 'MARBL_SFO_', &
851 trim(marbl_instances%surface_flux_output%outputs_for_GCM(m)%short_name)
852 call register_restart_field(cs%SFO(:,:,m), var_name, .false., restart_cs)
853 enddo
854
855 do m=1,cs%ito_cnt
856 write(var_name, "(2A)") 'MARBL_ITO_', &
857 trim(marbl_instances%interior_tendency_output%outputs_for_GCM(m)%short_name)
858 call register_restart_field(cs%ITO(:,:,:,m), var_name, .false., restart_cs)
859 enddo
860
861
862 cs%tr_Reg => tr_reg
863 cs%restart_CSp => restart_cs
864
867
868end function register_marbl_tracers
869
870!> This subroutine initializes the CS%ntr tracer fields in tr(:,:,:,:)
871!! and it sets up the tracer output.
872subroutine initialize_marbl_tracers(restart, day, G, GV, US, h, param_file, diag, OBC, CS, sponge_CSp)
873 logical, intent(in) :: restart !< .true. if the fields have already been
874 !! read from a restart file.
875 type(time_type), target, intent(in) :: day !< Time of the start of the run.
876 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
877 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
878 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
879 real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
880 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
881 type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output.
882 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
883 !! whether, where, and what open boundary
884 !! conditions are used.
885 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a previous
886 !! call to register_MARBL_tracers.
887 type(sponge_cs), pointer :: sponge_csp !< A pointer to the control structure
888 !! for the sponges, if they are in use.
889
890 ! Local variables
891 character(len=200) :: log_message
892 character(len=48) :: name ! A variable's name in a NetCDF file.
893 character(len=100) :: longname ! The long name of that variable.
894 character(len=48) :: units ! The units of the variable.
895 character(len=48) :: flux_units ! The units for age tracer fluxes, either
896 ! years m3 s-1 or years kg s-1.
897 character(len=48) :: tracer_name
898 logical :: fesedflux_has_edges, fesedflux_use_missing, tracer_init_from_z
899 real :: fesedflux_missing ! required argument for read_Z_edges() [CU ~> conc]
900 integer :: i, j, k, kbot, m, diag_size
901
902 if (.not.associated(cs)) return
903 if (cs%ntr < 1) return
904
905 cs%diag => diag
906
907 ! Allocate memory for surface tracer fluxes
908 allocate(cs%STF(szi_(g), szj_(g), cs%ntr), &
909 cs%RIV_FLUXES(szi_(g), szj_(g), cs%ntr), &
910 source=0.0)
911
912 ! Allocate memory for d14c forcing
913 if (cs%abio_dic_on) allocate(cs%d14c(szi_(g), szj_(g)))
914
915 ! Register diagnostics returned from MARBL (surface flux first, then interior tendency)
916 call register_marbl_diags(marbl_instances%surface_flux_diags, diag, day, g, cs%surface_flux_diags)
917 call register_marbl_diags(marbl_instances%interior_tendency_diags, diag, day, g, &
918 cs%interior_tendency_diags)
919
920 ! Register per-tracer diagnostics computed from MARBL surface flux / interior tendency values
921 allocate(cs%id_surface_flux_out(cs%ntr))
922 allocate(cs%id_surface_flux_from_salt_flux(cs%ntr))
923 allocate(cs%interior_tendency_out(cs%ntr))
924 allocate(cs%interior_tendency_out_zint(cs%ntr))
925 allocate(cs%interior_tendency_out_zint_100m(cs%ntr))
926 do m=1,cs%ntr
927 write(name, "(2A)") "STF_", trim(marbl_instances%tracer_metadata(m)%short_name)
928 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), " Surface Flux"
929 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), " m/s"
930 cs%id_surface_flux_out(m) = register_diag_field("ocean_model", trim(name), &
931 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid
932 day, trim(longname), trim(units), conversion=us%Z_to_m*us%s_to_T)
933
934 write(name, "(2A)") "STF_SALT_", trim(marbl_instances%tracer_metadata(m)%short_name)
935 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), " Surface Flux from Salt Flux"
936 cs%id_surface_flux_from_salt_flux(m) = register_diag_field("ocean_model", trim(name), &
937 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid
938 day, trim(longname), trim(units), conversion=us%Z_to_m*us%s_to_T)
939
940 write(name, "(2A)") "J_", trim(marbl_instances%tracer_metadata(m)%short_name)
941 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), " Source Sink Term"
942 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), "/s"
943 cs%interior_tendency_out(m)%id = register_diag_field("ocean_model", trim(name), &
944 diag%axesTL, & ! T=> tracer grid? L => layer center
945 day, trim(longname), trim(units))
946 if (cs%interior_tendency_out(m)%id > 0) &
947 allocate(cs%interior_tendency_out(m)%field_3d(szi_(g),szj_(g), szk_(g)), source=0.0)
948
949 write(name, "(2A)") "Jint_", trim(marbl_instances%tracer_metadata(m)%short_name)
950 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), &
951 " Source Sink Term Vertical Integral"
952 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), " m/s"
953 cs%interior_tendency_out_zint(m)%id = register_diag_field("ocean_model", trim(name), &
954 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
955 day, trim(longname), trim(units))
956 if (cs%interior_tendency_out_zint(m)%id > 0) &
957 allocate(cs%interior_tendency_out_zint(m)%field_2d(szi_(g),szj_(g)), source=0.0)
958
959 write(name, "(2A)") "Jint_100m_", trim(marbl_instances%tracer_metadata(m)%short_name)
960 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), &
961 " Source Sink Term Vertical Integral, 0-100m"
962 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), " m/s"
963 cs%interior_tendency_out_zint_100m(m)%id = register_diag_field("ocean_model", trim(name), &
964 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
965 day, trim(longname), trim(units))
966 if (cs%interior_tendency_out_zint_100m(m)%id > 0) &
967 allocate(cs%interior_tendency_out_zint_100m(m)%field_2d(szi_(g),szj_(g)), source=0.0)
968
969 enddo
970
971 ! Register diagnostics for MOM to report that are not tracer specific
972 cs%bot_flux_to_tend_id = register_diag_field("ocean_model", "BOT_FLUX_TO_TEND", &
973 diag%axesTL, & ! T=> tracer grid? L => layer center
974 day, "Conversion Factor for Bottom Flux -> Tend", "1/m")
975
976 ! Initialize tracers (if they weren't initialized from restart file)
977 tracer_init_from_z = .false.
978 do m=1,cs%ntr
979 call query_vardesc(cs%tr_desc(m), name=name, caller="initialize_MARBL_tracers")
980 if ((.not. restart) .or. &
981 (cs%tracers_may_reinit .and. &
982 .not. query_initialized(cs%tracer_data(m)%tr(:,:,:), name, cs%restart_CSp))) then
983 ! TODO: added the ongrid optional argument, but is there a good way to detect if the file is on grid?
984 call mom_initialize_tracer_from_z(h, cs%tracer_data(m)%tr, g, gv, us, param_file, &
985 cs%IC_file, name, ongrid=cs%ongrid)
986 tracer_init_from_z = .true.
987 do k=1,gv%ke ; do j=g%jsc, g%jec ; do i=g%isc, g%iec
988 ! Ensure tracer concentrations are at / above minimum value
989 if (cs%tracer_data(m)%tr(i,j,k) < cs%IC_min) cs%tracer_data(m)%tr(i,j,k) = cs%IC_min
990 enddo ; enddo ; enddo
991 endif
992 enddo
993 if (tracer_init_from_z) then
994 ! For each column, enforce consistency in MARBL tracers
995 ! (no negative concentrations; for a given autotroph, if one tracer is 0 they all are)
996 call mom_error(note, 'Enforcing consistency across autotroph tracer initial conditions')
997 do j=g%jsc, g%jec ; do i=g%isc, g%iec
998 ! Copy tracer data into flat array
999 do k=1,gv%ke ; do m=1, cs%ntr
1000 marbl_instances%tracers(m,k) = cs%tracer_data(m)%tr(i,j,k)
1001 enddo ; enddo
1002 ! call consistency enforcement
1003 call marbl_instances%autotroph_tracer_consistency_enforce()
1004 ! Copy tracer data out of flat array
1005 do k=1,gv%ke ; do m=1, cs%ntr
1006 cs%tracer_data(m)%tr(i,j,k) = marbl_instances%tracers(m,k)
1007 enddo ; enddo
1008 enddo ; enddo
1009 endif
1010
1011 ! Initialize total chlorophyll to get SW Pen correct (if it wasn't initialized from restart file)
1012 if ((cs%total_Chl_ind > 0) .and. &
1013 ((.not. restart) .or. &
1014 (.not. query_initialized(cs%ITO(:,:,:,cs%total_Chl_ind), "MARBL_ITO_total_Chl", cs%restart_CSp)))) then
1015 ! Three steps per column
1016 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1017 ! (i) Copy initial tracers into MARBL structure
1018 do k=1,gv%ke ; do m=1,cs%ntr
1019 marbl_instances%tracers(m,k) = max(cs%tracer_data(m)%tr(i,j,k), 0.)
1020 enddo ; enddo
1021 ! (ii) Compute total Chl for the column
1022 call marbl_instances%compute_totChl()
1023 ! (iii) Copy total Chl from MARBL data-structure into CS%ITO
1024 do k=1,gv%ke
1025 cs%ITO(i,j,k,cs%total_Chl_ind) = &
1026 marbl_instances%interior_tendency_output%outputs_for_GCM(cs%total_Chl_ind)%forcing_field_1d(1,k)
1027 enddo
1028 enddo ; enddo
1029 endif
1030
1031 ! Register diagnostics for river fluxes
1032 cs%no3_riv_flux = register_diag_field("ocean_model", "NO3_RIV_FLUX", &
1033 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1034 day, "Dissolved Inorganic Nitrate Riverine Flux", "mmol/m^3 m/s")
1035 cs%po4_riv_flux = register_diag_field("ocean_model", "PO4_RIV_FLUX", &
1036 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1037 day, "Dissolved Inorganic Phosphate Riverine Flux", "mmol/m^3 m/s")
1038 cs%don_riv_flux = register_diag_field("ocean_model", "DON_RIV_FLUX", &
1039 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1040 day, "Dissolved Organic Nitrogen Riverine Flux", "mmol/m^3 m/s")
1041 cs%donr_riv_flux = register_diag_field("ocean_model", "DONR_RIV_FLUX", &
1042 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1043 day, "Refractory DON Riverine Flux", "mmol/m^3 m/s")
1044 cs%dop_riv_flux = register_diag_field("ocean_model", "DOP_RIV_FLUX", &
1045 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1046 day, "Dissolved Organic Phosphorus Riverine Flux", "mmol/m^3 m/s")
1047 cs%dopr_riv_flux = register_diag_field("ocean_model", "DOPR_RIV_FLUX", &
1048 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1049 day, "Refractory DOP Riverine Flux", "mmol/m^3 m/s")
1050 cs%sio3_riv_flux = register_diag_field("ocean_model", "SiO3_RIV_FLUX", &
1051 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1052 day, "Dissolved Inorganic Silicate Riverine Flux", "mmol/m^3 m/s")
1053 cs%fe_riv_flux = register_diag_field("ocean_model", "Fe_RIV_FLUX", &
1054 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1055 day, "Dissolved Inorganic Iron Riverine Flux", "mmol/m^3 m/s")
1056 cs%doc_riv_flux = register_diag_field("ocean_model", "DOC_RIV_FLUX", &
1057 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1058 day, "Dissolved Organic Carbon Riverine Flux", "mmol/m^3 m/s")
1059 cs%docr_riv_flux = register_diag_field("ocean_model", "DOCR_RIV_FLUX", &
1060 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1061 day, "Refractory DOC Riverine Flux", "mmol/m^3 m/s")
1062 cs%alk_riv_flux = register_diag_field("ocean_model", "ALK_RIV_FLUX", &
1063 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1064 day, "Alkalinity Riverine Flux", "meq/m^3 m/s")
1065 cs%alk_alt_co2_riv_flux = register_diag_field("ocean_model", "ALK_ALT_CO2_RIV_FLUX", &
1066 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1067 day, "Alkalinity Riverine Flux, Alternative CO2", "meq/m^3 m/s")
1068 cs%dic_riv_flux = register_diag_field("ocean_model", "DIC_RIV_FLUX", &
1069 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1070 day, "Dissolved Inorganic Carbon Riverine Flux", "mmol/m^3 m/s")
1071 cs%dic_alt_co2_riv_flux = register_diag_field("ocean_model", "DIC_ALT_CO2_RIV_FLUX", &
1072 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1073 day, "Dissolved Inorganic Carbon Riverine Flux, Alternative CO2", "mmol/m^3 m/s")
1074
1075 ! Register diagnostics for d14c forcing
1076 if (cs%abio_dic_on) then
1077 cs%d14c_id = register_diag_field("ocean_model", "D14C_FORCING", &
1078 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
1079 day, "Delta-14C in atmospheric CO2", "per mil, relative to Modern")
1080 endif
1081
1082 ! Register diagnostics for per-category forcing fields
1083 if (cs%ice_ncat > 0) then
1084 allocate(cs%fracr_cat_id(cs%ice_ncat+1))
1085 allocate(cs%qsw_cat_id(cs%ice_ncat+1))
1086 do m=1,cs%ice_ncat+1
1087 write(name, "(A,I0)") "FRACR_CAT_", m
1088 write(longname, "(A,I0)") "Fraction of area in ice category ", m
1089 units = "fraction"
1090 cs%fracr_cat_id(m) = register_diag_field("ocean_model", trim(name), &
1091 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid
1092 day, trim(longname), trim(units))
1093 write(name, "(A,I0)") "QSW_CAT_", m
1094 write(longname, "(A,I0)") "Shortwave penetrating through ice category ", m
1095 units = "W m-2"
1096 cs%qsw_cat_id(m) = register_diag_field("ocean_model", trim(name), &
1097 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid
1098 day, trim(longname), trim(units))
1099 enddo
1100 endif
1101
1102 if (cs%base_bio_on) then
1103 ! Read initial fesedflux and feventflux fields
1104 ! (1) get vertical dimension
1105 ! -- comes from fesedflux_file, assume same dimension in feventflux
1106 ! (maybe these fields should be combined?)
1107 ! -- note: read_Z_edges treats depth as positive UP => 0 at surface, negative at depth
1108 fesedflux_use_missing = .false.
1109 call read_z_edges(cs%fesedflux_file, "FESEDFLUXIN", cs%fesedflux_z_edges, cs%fesedflux_nz, &
1110 fesedflux_has_edges, fesedflux_use_missing, fesedflux_missing, scale=us%m_to_Z, &
1111 missing_scale=1.0)
1112
1113 ! (2) Allocate memory for fesedflux and feventflux
1114 allocate(cs%fesedflux_in(szi_(g), szj_(g), cs%fesedflux_nz))
1115 allocate(cs%fesedfluxred_in(szi_(g), szj_(g), cs%fesedflux_nz))
1116 allocate(cs%feventflux_in(szi_(g), szj_(g), cs%fesedflux_nz))
1117 allocate(cs%fesedflux_dz(szi_(g), szj_(g), cs%fesedflux_nz))
1118
1119 ! (3) Read data
1120 ! TODO: Add US term to scale
1121 call mom_read_data(cs%fesedflux_file, "FESEDFLUXIN", cs%fesedflux_in(:,:,:), g%Domain, &
1122 scale=cs%fesedflux_scale_factor)
1123 call mom_read_data(cs%fesedfluxred_file, "FESEDFLUXIN", cs%fesedfluxred_in(:,:,:), g%Domain, &
1124 scale=cs%fesedflux_scale_factor)
1125 call mom_read_data(cs%feventflux_file, "FESEDFLUXIN", cs%feventflux_in(:,:,:), g%Domain, &
1126 scale=cs%fesedflux_scale_factor)
1127
1128 ! (4) Relocate values that are below ocean bottom to layer that intersects bathymetry
1129 ! Remember, fesedflux_z_edges = 0 at surface and is < 0 below surface
1130
1131 do k=cs%fesedflux_nz, 1, -1
1132 kbot = k + 1 ! level k is between z(k) and z(k+1)
1133 do j=g%jsc, g%jec
1134 do i=g%isc, g%iec
1135 if (g%mask2dT(i,j) == 0) cycle
1136 if (g%bathyT(i,j) + cs%fesedflux_z_edges(1) < 1e-8 * us%m_to_Z) then
1137 write(log_message, *) "Current implementation of fesedflux assumes G%bathyT >=", &
1138 " first edge;first edge = ", -cs%fesedflux_z_edges(1), "bathyT = ", g%bathyT(i,j)
1139 call mom_error(fatal, log_message)
1140 endif
1141 ! Also figure out layer thickness while we're here
1142 cs%fesedflux_dz(i,j,k) = (cs%fesedflux_z_edges(k) - cs%fesedflux_z_edges(kbot)) * gv%Z_to_H
1143 ! If top interface is at or below ocean bottom, move flux in current layer up one
1144 ! and set thickness of current level to 0
1145 if (g%bathyT(i,j) + cs%fesedflux_z_edges(k) < 1e-8 * us%m_to_Z) then
1146 cs%fesedflux_in(i,j,k-1) = cs%fesedflux_in(i,j,k-1) + cs%fesedflux_in(i,j,k)
1147 cs%fesedflux_in(i,j,k) = 0.
1148 cs%fesedfluxred_in(i,j,k-1) = cs%fesedfluxred_in(i,j,k-1) + cs%fesedfluxred_in(i,j,k)
1149 cs%fesedfluxred_in(i,j,k) = 0.
1150 cs%feventflux_in(i,j,k-1) = cs%feventflux_in(i,j,k-1) + cs%feventflux_in(i,j,k)
1151 cs%feventflux_in(i,j,k) = 0.
1152 cs%fesedflux_dz(i,j,k) = 0.
1153 elseif (g%bathyT(i,j) + cs%fesedflux_z_edges(kbot) < 1e-8 * us%m_to_Z) then
1154 ! Otherwise, if lower interface is below bathymetry move interface to ocean bottom
1155 cs%fesedflux_dz(i,j,k) = (g%bathyT(i,j) + cs%fesedflux_z_edges(k)) * gv%Z_to_H
1156 endif
1157 enddo
1158 enddo
1159 enddo
1160
1161 ! Initialize external field for river fluxes
1162 if (cs%read_riv_fluxes) then
1163 cs%id_din_riv = init_external_field(cs%riv_flux_dataset%file_name, 'din_riv_flux', &
1164 domain=g%Domain%mpp_domain)
1165 cs%id_don_riv = init_external_field(cs%riv_flux_dataset%file_name, 'don_riv_flux', &
1166 domain=g%Domain%mpp_domain)
1167 cs%id_dip_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dip_riv_flux', &
1168 domain=g%Domain%mpp_domain)
1169 cs%id_dop_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dop_riv_flux', &
1170 domain=g%Domain%mpp_domain)
1171 cs%id_dsi_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dsi_riv_flux', &
1172 domain=g%Domain%mpp_domain)
1173 cs%id_dfe_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dfe_riv_flux', &
1174 domain=g%Domain%mpp_domain)
1175 cs%id_dic_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dic_riv_flux', &
1176 domain=g%Domain%mpp_domain)
1177 cs%id_alk_riv = init_external_field(cs%riv_flux_dataset%file_name, 'alk_riv_flux', &
1178 domain=g%Domain%mpp_domain)
1179 cs%id_doc_riv = init_external_field(cs%riv_flux_dataset%file_name, 'doc_riv_flux', &
1180 domain=g%Domain%mpp_domain)
1181 endif
1182 endif
1183
1184 if (cs%abio_dic_on) then
1185 ! Initialize external field for d14c forcing
1186 do m=1,3
1187 cs%id_d14c(m) = init_external_field(cs%d14c_dataset(m)%file_name, "Delta14co2_in_air", &
1188 ignore_axis_atts=.true.)
1189 enddo
1190 endif
1191
1192 ! Initialize external field for restoring
1193 if (cs%restoring_I_tau_source == "file") then
1194 select case(cs%restoring_source)
1195 case("file")
1196 ! Set up array for reading in raw restoring data
1197 allocate(cs%restoring_in(szi_(g), szj_(g), cs%restoring_nz, cs%restore_count), source=0.)
1198 do m=1,cs%restore_count
1199 cs%id_tracer_restoring(m) = init_external_field(cs%restoring_file, &
1200 trim(cs%tracer_restoring_varname(m)), domain=g%Domain%mpp_domain)
1201 enddo
1202 end select
1203 select case(cs%restoring_I_tau_source)
1204 case("file")
1205 allocate(cs%I_tau(szi_(g), szj_(g), cs%restoring_timescale_nz), source=0.)
1206 call mom_read_data(cs%restoring_I_tau_file, "RTAU", cs%I_tau(:,:,:), g%Domain)
1207 end select
1208 endif
1209
1210end subroutine initialize_marbl_tracers
1211
1212!> This subroutine is used to register tracer fields and subroutines
1213!! to be used with MOM.
1214subroutine register_marbl_diags(MARBL_diags, diag, day, G, id_diags)
1215
1216 type(marbl_diagnostics_type), intent(in) :: MARBL_diags !< MARBL diagnostics from MARBL_instances
1217 type(time_type), target, intent(in) :: day !< Time of the start of the run.
1218 type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output.
1219 !integer, allocatable, intent(inout) :: id_diags(:) !< allocatable array storing diagnostic index number
1220 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1221 type(temp_marbl_diag), allocatable, intent(inout) :: id_diags(:) !< allocatable array storing diagnostic index
1222 !! number and buffer space for collecting diags
1223 !! from all columns
1224
1225 integer :: m, diag_size
1226
1227 diag_size = size(marbl_diags%diags)
1228 allocate(id_diags(diag_size))
1229 do m = 1, diag_size
1230 id_diags(m)%id = -1
1231 if (trim(marbl_diags%diags(m)%vertical_grid) == "none") then ! 2D field
1232 id_diags(m)%id = register_diag_field("ocean_model", &
1233 trim(marbl_diags%diags(m)%short_name), &
1234 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid
1235 day, &
1236 trim(marbl_diags%diags(m)%long_name), &
1237 trim(marbl_diags%diags(m)%units))
1238 if (id_diags(m)%id > 0) allocate(id_diags(m)%field_2d(szi_(g),szj_(g)), source=0.0)
1239 else ! 3D field
1240 ! TODO: MARBL should provide v_extensive through MARBL_diags
1241 ! (for now, FESEDFLUX is the only one that should be true)
1242 ! Also, known issue where passing v_extensive=.false. isn't
1243 ! treated the same as not passing v_extensive
1244 if ((trim(marbl_diags%diags(m)%short_name) == "FESEDFLUX") .or. &
1245 (trim(marbl_diags%diags(m)%short_name) == "FEREDSEDFLUX") .or. &
1246 (trim(marbl_diags%diags(m)%short_name) == "FEVENTFLUX")) then
1247 id_diags(m)%id = register_diag_field("ocean_model", &
1248 trim(marbl_diags%diags(m)%short_name), &
1249 diag%axesTL, & ! T=> tracer grid? L => layer center
1250 day, &
1251 trim(marbl_diags%diags(m)%long_name), &
1252 trim(marbl_diags%diags(m)%units), &
1253 v_extensive=.true.)
1254 else
1255 id_diags(m)%id = register_diag_field("ocean_model", &
1256 trim(marbl_diags%diags(m)%short_name), &
1257 diag%axesTL, & ! T=> tracer grid? L => layer center
1258 day, &
1259 trim(marbl_diags%diags(m)%long_name), &
1260 trim(marbl_diags%diags(m)%units))
1261 endif
1262 if (id_diags(m)%id > 0) allocate(id_diags(m)%field_3d(szi_(g),szj_(g), szk_(g)), source=0.0)
1263 endif
1264 enddo
1265
1266end subroutine register_marbl_diags
1267
1268!> This subroutine allocates memory for saved state fields and registers them in the restart files
1269subroutine setup_saved_state(MARBL_saved_state, HI, GV, restart_CS, tracers_may_reinit, &
1270 local_saved_state)
1271
1272 type(marbl_saved_state_type), intent(in) :: MARBL_saved_state !< MARBL saved state from
1273 !! MARBL_instances
1274 type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
1275 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1276 type(mom_restart_cs), pointer, intent(in) :: restart_CS !< control structure to add saved state
1277 !! to restarts
1278 logical, intent(in) :: tracers_may_reinit !< used to determine mandatory
1279 !! flag in restart
1280 type(saved_state_for_marbl_type), allocatable, intent(inout) :: local_saved_state(:) !< allocatable array for local
1281 !! saved state
1282
1283 integer :: num_fields, m
1284 character(len=200) :: log_message, varname
1285
1286 num_fields = marbl_saved_state%saved_state_cnt
1287 allocate(local_saved_state(num_fields))
1288
1289 do m=1,num_fields
1290 write(varname, "(2A)") "MARBL_", trim(marbl_saved_state%state(m)%short_name)
1291 select case (marbl_saved_state%state(m)%rank)
1292 case (2)
1293 allocate(local_saved_state(m)%field_2d(szi_(hi),szj_(hi)), source=0.0)
1294 call register_restart_field(local_saved_state(m)%field_2d, varname, &
1295 .not.tracers_may_reinit, restart_cs)
1296 case (3)
1297 if (trim(marbl_saved_state%state(m)%vertical_grid).eq."layer_avg") then
1298 allocate(local_saved_state(m)%field_3d(szi_(hi),szj_(hi), szk_(gv)), source=0.0)
1299 call register_restart_field(local_saved_state(m)%field_3d, varname, &
1300 .not.tracers_may_reinit, restart_cs)
1301 else
1302 write(log_message, "(3A, I0, A)") "'", trim(marbl_saved_state%state(m)%vertical_grid), &
1303 "' is an invalid vertical grid for saved state (ind = ", m, ")"
1304 call mom_error(fatal, log_message)
1305 endif
1306 case DEFAULT
1307 write(log_message, "(I0, A, I0, A)") marbl_saved_state%state(m)%rank, &
1308 " is an invalid rank for saved state (ind = ", m, ")"
1309 call mom_error(fatal, log_message)
1310 end select
1311 local_saved_state(m)%short_name = trim(marbl_saved_state%state(m)%short_name)
1312 write(local_saved_state(m)%file_varname, "(2A)") "MARBL_", trim(local_saved_state(m)%short_name)
1313 local_saved_state(m)%units = trim(marbl_saved_state%state(m)%units)
1314 enddo
1315
1316end subroutine setup_saved_state
1317
1318!> This subroutine applies diapycnal diffusion and any other column
1319!! tracer physics or chemistry to the tracers from this file.
1320subroutine marbl_tracers_column_physics(h_old, ea, eb, fluxes, dt, G, GV, US, CS, &
1321 prediabatic_T, prediabatic_S, KPP_CSp, nonLocalTrans, evap_CFL_limit, minimum_forcing_depth)
1322
1323 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1324 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1325 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1326 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
1327 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1328 intent(in) :: ea !< an array to which the amount of fluid entrained
1329 !! from the layer above during this call will be
1330 !! added [H ~> m or kg m-2].
1331 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1332 intent(in) :: eb !< an array to which the amount of fluid entrained
1333 !! from the layer below during this call will be
1334 !! added [H ~> m or kg m-2].
1335 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
1336 !! and tracer forcing fields. Unused fields have NULL ptrs.
1337 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
1338 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1339 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a previous
1340 !! call to register_MARBL_tracers.
1341 real, dimension(:,:,:), intent(in) :: prediabatic_t !< Temperature prior to calling diabatic driver [C ~> degC]
1342 real, dimension(:,:,:), intent(in) :: prediabatic_s !< Salinity prior to calling diabatic driver [S ~> ppt]
1343
1344 type(kpp_cs), optional, pointer :: kpp_csp !< KPP control structure
1345 real, optional, intent(in) :: nonlocaltrans(:,:,:) !< Non-local transport [1]
1346 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
1347 !! be fluxed out of the top layer in a timestep [1]
1348 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
1349 !! fluxes can be applied [m]
1350
1351 ! Local variables
1352 character(len=256) :: log_message
1353 real, dimension(SZI_(G),SZJ_(G)) :: net_salt_rate ! Surface salt flux into the ocean
1354 ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1].
1355 real, dimension(SZI_(G),SZJ_(G)) :: flux_from_salt_flux ! Surface tracer flux from salt flux
1356 ! [conc Z T-1 ~> conc m s-1].
1357 real, dimension(SZI_(G),SZJ_(G)) :: ref_mask ! Mask for 2D MARBL diags using ref_depth [1]
1358 real, dimension(SZI_(G),SZJ_(G)) :: riv_flux_loc ! Local copy of CS%RIV_FLUXES*dt [conc H ~> mmol m-2]
1359 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
1360 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: bot_flux_to_tend ! Conversion factor for bottom tlux -> tend
1361 ! [Z-1 ~> m-1]
1362 real :: cum_bftt_dz ! sum of bot_flux_to_tend * dz from the bottom layer to current layer [1]
1363 real, dimension(0:GV%ke) :: zi ! z-coordinate interface depth [Z ~> m]
1364 real, dimension(GV%ke) :: zc ! z-coordinate layer center depth [Z ~> m]
1365 real, dimension(GV%ke) :: dz ! z-coordinate cell thickness [H ~> m]
1366 integer :: i, j, k, is, ie, js, je, nz, m
1367
1368 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1369
1370 if (.not.associated(cs)) return
1371
1372 ! (1) Compute surface fluxes and interior tendencies
1373 ! FIXME: MARBL can handle computing surface fluxes for all columns simultaneously
1374 ! I was just thinking going column-by-column at first might be easier
1375 bot_flux_to_tend(:, :, :) = 0.
1376 do j=js,je
1377 do i=is,ie
1378 ! Surface fluxes
1379 ! i. only want ocean points in this loop
1380 if (g%mask2dT(i,j) == 0) cycle
1381
1382 ! ii. Load proper column data
1383 ! * surface flux forcings
1384 ! These fields are getting the correct data
1385 ! TODO: if top layer is vanishly thin, do we actually want (e.g.) top 5m average temp / salinity?
1386 ! How does MOM pass SST and SSS to GFDL coupler? (look in core.F90?)
1387 if (cs%sss_ind > 0) &
1388 marbl_instances%surface_flux_forcings(cs%sss_ind)%field_0d(1) = prediabatic_s(i,j,1) * us%S_to_ppt
1389 if (cs%sst_ind > 0) &
1390 marbl_instances%surface_flux_forcings(cs%sst_ind)%field_0d(1) = prediabatic_t(i,j,1) * us%C_to_degC
1391 if (cs%ifrac_ind > 0) &
1392 marbl_instances%surface_flux_forcings(cs%ifrac_ind)%field_0d(1) = fluxes%ice_fraction(i,j)
1393
1394 ! MARBL wants u10_sqr in (m/s)^2
1395 if (cs%u10_sqr_ind > 0) &
1396 marbl_instances%surface_flux_forcings(cs%u10_sqr_ind)%field_0d(1) = fluxes%u10_sqr(i,j) * &
1397 ((us%L_T_to_m_s)**2)
1398
1399 ! mct_driver/ocn_cap_methods:93 -- ice_ocean_boundary%p(i,j) comes from coupler
1400 ! We may need a new ice_ocean_boundary%p_atm because %p includes ice in GFDL driver
1401 if (cs%atmpress_ind > 0) then
1402 if (associated(fluxes%p_surf_full)) then
1403 marbl_instances%surface_flux_forcings(cs%atmpress_ind)%field_0d(1) = &
1404 fluxes%p_surf_full(i,j) * ((us%R_to_kg_m3 * (us%L_T_to_m_s**2)) * atm_per_pa)
1405 else
1406 ! hardcode value of 1 atm (can't figure out how to get this from solo_driver)
1407 marbl_instances%surface_flux_forcings(cs%atmpress_ind)%field_0d(1) = 1.
1408 endif
1409 endif
1410
1411 ! These are okay, but need option to come in from coupler
1412 if (cs%xco2_ind > 0) &
1413 marbl_instances%surface_flux_forcings(cs%xco2_ind)%field_0d(1) = fluxes%atm_co2(i,j)
1414 if (cs%xco2_alt_ind > 0) &
1415 marbl_instances%surface_flux_forcings(cs%xco2_alt_ind)%field_0d(1) = fluxes%atm_alt_co2(i,j)
1416
1417 ! These are okay, but need option to read in from file
1418 if (cs%dust_dep_ind > 0) &
1419 marbl_instances%surface_flux_forcings(cs%dust_dep_ind)%field_0d(1) = &
1420 fluxes%dust_flux(i,j) * us%RZ_T_to_kg_m2s
1421
1422 if (cs%fe_dep_ind > 0) &
1423 marbl_instances%surface_flux_forcings(cs%fe_dep_ind)%field_0d(1) = &
1424 fluxes%iron_flux(i,j) * (us%Z_to_m * us%s_to_T)
1425
1426 ! MARBL wants ndep in (mmol/m^2/s)
1427 if (cs%nox_flux_ind > 0) &
1428 marbl_instances%surface_flux_forcings(cs%nox_flux_ind)%field_0d(1) = fluxes%noy_dep(i,j) * &
1429 (us%Z_to_m * us%s_to_T)
1430 if (cs%nhy_flux_ind > 0) &
1431 marbl_instances%surface_flux_forcings(cs%nhy_flux_ind)%field_0d(1) = fluxes%nhx_dep(i,j) * &
1432 (us%Z_to_m * us%s_to_T)
1433
1434 if (cs%d14c_ind > 0) &
1435 marbl_instances%surface_flux_forcings(cs%d14c_ind)%field_0d(1) = cs%d14c(i,j)
1436
1437 ! * tracers at surface
1438 ! TODO: average over some shallow depth (e.g. 5m)
1439 do m=1,cs%ntr
1440 marbl_instances%tracers_at_surface(1,m) = cs%tracer_data(m)%tr(i,j,1)
1441 enddo
1442
1443 ! * surface flux saved state
1444 do m=1,size(marbl_instances%surface_flux_saved_state%state)
1445 ! (currently only 2D fields are saved from surface_flux_compute())
1446 marbl_instances%surface_flux_saved_state%state(m)%field_2d(1) = &
1447 cs%surface_flux_saved_state(m)%field_2d(i,j)
1448 enddo
1449
1450 ! iii. Compute surface fluxes in MARBL
1451 call marbl_instances%surface_flux_compute()
1452 if (marbl_instances%StatusLog%labort_marbl) then
1453 call marbl_instances%StatusLog%log_error_trace("MARBL_instances%surface_flux_compute()", &
1454 "MARBL_tracers_column_physics")
1455 endif
1456 call print_marbl_log(marbl_instances%StatusLog, g, i, j)
1457 call marbl_instances%StatusLog%erase()
1458
1459 ! iv. Copy output that MOM6 needs to hold on to
1460 ! * saved state
1461 do m=1,size(marbl_instances%surface_flux_saved_state%state)
1462 cs%surface_flux_saved_state(m)%field_2d(i,j) = &
1463 marbl_instances%surface_flux_saved_state%state(m)%field_2d(1)
1464 enddo
1465
1466 ! * diagnostics
1467 do m=1,size(marbl_instances%surface_flux_diags%diags)
1468 ! All diags are 2D coming from surface
1469 if (cs%surface_flux_diags(m)%id > 0) &
1470 cs%surface_flux_diags(m)%field_2d(i,j) = &
1471 real(marbl_instances%surface_flux_diags%diags(m)%field_2d(1))
1472 enddo
1473
1474 ! * Surface tracer flux
1475 cs%STF(i,j,:) = marbl_instances%surface_fluxes(1,:) * (us%m_to_Z * us%T_to_s)
1476
1477 ! * Surface flux output
1478 do m=1,cs%sfo_cnt
1479 cs%SFO(i,j,m) = marbl_instances%surface_flux_output%outputs_for_GCM(m)%forcing_field_0d(1)
1480 enddo
1481
1482 ! interior tendencies
1483 ! i. Set up vertical domain and bot_flux_to_tend
1484 ! Calculate depth of interface by building up thicknesses from the bottom (top interface is always 0)
1485 ! MARBL wants this to be positive-down
1486 zi(gv%ke) = g%bathyT(i,j)
1487 marbl_instances%bot_flux_to_tend(:) = 0.
1488 cum_bftt_dz = 0.
1489 do k = gv%ke, 1, -1
1490 dz(k) = h_old(i,j,k) ! cell thickness
1491 zc(k) = zi(k) - 0.5 * (dz(k)*gv%H_to_Z)
1492 zi(k-1) = zi(k) - (dz(k)*gv%H_to_Z)
1493 if (g%bathyT(i,j) - zi(k-1) <= cs%bot_flux_mix_thickness) then
1494 marbl_instances%bot_flux_to_tend(k) = us%m_to_Z * cs%Ibfmt
1495 cum_bftt_dz = cum_bftt_dz + marbl_instances%bot_flux_to_tend(k) * (gv%H_to_m * dz(k))
1496 elseif (g%bathyT(i,j) - zi(k) < cs%bot_flux_mix_thickness) then
1497 ! MARBL_instances%bot_flux_to_tend(k) = (1. - (G%bathyT(i,j) - zi(k)) * CS%Ibfmt) / dz(k)
1498 marbl_instances%bot_flux_to_tend(k) = (1. - cum_bftt_dz) / (gv%H_to_m * dz(k))
1499 endif
1500 enddo
1501 if (g%bathyT(i,j) - zi(0) < cs%bot_flux_mix_thickness) &
1502 marbl_instances%bot_flux_to_tend(:) = marbl_instances%bot_flux_to_tend(:) * &
1503 cs%bot_flux_mix_thickness / (g%bathyT(i,j) - zi(0))
1504 if (cs%bot_flux_to_tend_id > 0) &
1505 bot_flux_to_tend(i, j, :) = marbl_instances%bot_flux_to_tend(:)
1506
1507 ! zw(1:nz) is bottom cell depth so no element of zw = 0, it is assumed to be top layer depth
1508 marbl_instances%domain%zw(:) = us%Z_to_m * zi(1:gv%ke)
1509 marbl_instances%domain%zt(:) = us%Z_to_m * zc(:)
1510 marbl_instances%domain%delta_z(:) = gv%H_to_m * dz(:)
1511
1512 ! ii. Load proper column data
1513 ! * Forcing Fields
1514 ! These fields are getting the correct data
1515 if (cs%potemp_ind > 0) &
1516 marbl_instances%interior_tendency_forcings(cs%potemp_ind)%field_1d(1,:) = prediabatic_t(i,j,:) * us%C_to_degC
1517 if (cs%salinity_ind > 0) &
1518 marbl_instances%interior_tendency_forcings(cs%salinity_ind)%field_1d(1,:) = prediabatic_s(i,j,:) * us%S_to_ppt
1519
1520 ! This is okay, but need option to read in from file
1521 ! (Same as dust_dep_ind for surface_flux_forcings)
1522 if (cs%dustflux_ind > 0) &
1523 marbl_instances%interior_tendency_forcings(cs%dustflux_ind)%field_0d(1) = &
1524 fluxes%dust_flux(i,j) * us%RZ_T_to_kg_m2s
1525
1526 ! TODO: Support PAR (currently just using single subcolumn)
1527 ! (Look for Pen_sw_bnd?)
1528 if (cs%PAR_col_frac_ind > 0) then
1529 ! second index is num_subcols, not depth
1530 !MARBL_instances%interior_tendency_forcings(CS%PAR_col_frac_ind)%field_1d(1,:) = fluxes%fracr_cat(i,j,:)
1531 if (cs%use_ice_category_fields) then
1532 marbl_instances%interior_tendency_forcings(cs%PAR_col_frac_ind)%field_1d(1,:) = &
1533 fluxes%fracr_cat(i,j,:)
1534 else
1535 marbl_instances%interior_tendency_forcings(cs%PAR_col_frac_ind)%field_1d(1,1) = 1.
1536 endif
1537 endif
1538
1539 if (cs%surf_shortwave_ind > 0) then
1540 ! second index is num_subcols, not depth
1541 if (cs%use_ice_category_fields) then
1542 marbl_instances%interior_tendency_forcings(cs%surf_shortwave_ind)%field_1d(1,:) = &
1543 fluxes%qsw_cat(i,j,:) * us%QRZ_T_to_W_m2
1544 else
1545 marbl_instances%interior_tendency_forcings(cs%surf_shortwave_ind)%field_1d(1,1) = &
1546 fluxes%sw(i,j) * us%QRZ_T_to_W_m2
1547 endif
1548 endif
1549 ! Tracer restoring
1550 do m=1,cs%restore_count
1551 marbl_instances%interior_tendency_forcings(cs%tracer_restoring_ind(m))%field_1d(1,:) = 0.
1552 call remapping_core_h(cs%restoring_remapCS, cs%restoring_nz, cs%restoring_dz(:), &
1553 cs%restoring_in(i,j,:,m), gv%ke, dz(:), &
1554 marbl_instances%interior_tendency_forcings(cs%tracer_restoring_ind(m))%field_1d(1,:))
1555 if (m==1) then
1556 call remapping_core_h(cs%restoring_remapCS, cs%restoring_timescale_nz, &
1557 cs%restoring_timescale_dz(:), cs%I_tau(i,j,:), gv%ke, dz(:), &
1558 marbl_instances%interior_tendency_forcings(cs%tracer_I_tau_ind(m))%field_1d(1,:))
1559 else
1560 marbl_instances%interior_tendency_forcings(cs%tracer_I_tau_ind(m))%field_1d(1,:) = &
1561 marbl_instances%interior_tendency_forcings(cs%tracer_I_tau_ind(1))%field_1d(1,:)
1562 endif
1563 enddo
1564
1565 ! TODO: In POP, pressure comes from a function in state_mod.F90; I don't see a similar function here
1566 ! This formulation is from Levitus 1994, and I think it belongs in MOM_EOS.F90?
1567 ! Converts depth [m] -> pressure [bars]
1568 ! NOTE: Andrew recommends using GV%H_to_Pa
1569 if (cs%pressure_ind > 0) &
1570 marbl_instances%interior_tendency_forcings(cs%pressure_ind)%field_1d(1,:) = &
1571 (0.0598088 * (exp(-0.025*us%Z_to_m * zc(:)) - 1.)) + &
1572 (0.100766 * us%Z_to_m * zc(:)) + (2.28405e-7*((us%Z_to_m * zc(:))**2))
1573
1574 if (cs%fesedflux_ind > 0) then
1575 marbl_instances%interior_tendency_forcings(cs%fesedflux_ind)%field_1d(1,:) = 0.
1576 call reintegrate_column(cs%fesedflux_nz, &
1577 cs%fesedflux_dz(i,j,:) * (sum(dz(:) * gv%H_to_Z) / g%bathyT(i,j)), &
1578 cs%fesedflux_in(i,j,:), gv%ke, dz(:), &
1579 marbl_instances%interior_tendency_forcings(cs%fesedflux_ind)%field_1d(1,:))
1580 endif
1581 if (cs%fesedfluxred_ind > 0) then
1582 marbl_instances%interior_tendency_forcings(cs%fesedfluxred_ind)%field_1d(1,:) = 0.
1583 call reintegrate_column(cs%fesedflux_nz, &
1584 cs%fesedflux_dz(i,j,:) * (sum(dz(:) * gv%H_to_Z) / g%bathyT(i,j)), &
1585 cs%fesedfluxred_in(i,j,:), gv%ke, dz(:), &
1586 marbl_instances%interior_tendency_forcings(cs%fesedfluxred_ind)%field_1d(1,:))
1587 endif
1588 if (cs%feventflux_ind > 0) then
1589 marbl_instances%interior_tendency_forcings(cs%feventflux_ind)%field_1d(1,:) = 0.
1590 call reintegrate_column(cs%fesedflux_nz, &
1591 cs%fesedflux_dz(i,j,:) * (sum(dz(:) * gv%H_to_Z) / g%bathyT(i,j)), &
1592 cs%feventflux_in(i,j,:), gv%ke, dz(:), &
1593 marbl_instances%interior_tendency_forcings(cs%feventflux_ind)%field_1d(1,:))
1594 endif
1595
1596 ! TODO: add ability to read these fields from file
1597 ! also, add constant values to CS
1598 if (cs%o2_scalef_ind > 0) &
1599 marbl_instances%interior_tendency_forcings(cs%o2_scalef_ind)%field_1d(1,:) = 1.
1600 if (cs%remin_scalef_ind > 0) &
1601 marbl_instances%interior_tendency_forcings(cs%remin_scalef_ind)%field_1d(1,:) = 1.
1602
1603 ! * Column Tracers
1604 do m=1,cs%ntr
1605 marbl_instances%tracers(m, :) = cs%tracer_data(m)%tr(i,j,:)
1606 enddo
1607
1608 ! * interior tendency saved state
1609 ! (currently only 3D fields are saved from interior_tendency_compute())
1610 do m=1,size(marbl_instances%interior_tendency_saved_state%state)
1611 marbl_instances%interior_tendency_saved_state%state(m)%field_3d(:,1) = &
1612 cs%interior_tendency_saved_state(m)%field_3d(i,j,:)
1613 enddo
1614
1615 ! iii. Compute interior tendencies in MARBL
1616 call marbl_instances%interior_tendency_compute()
1617 if (marbl_instances%StatusLog%labort_marbl) then
1618 call marbl_instances%StatusLog%log_error_trace(&
1619 "MARBL_instances%interior_tendency_compute()", "MARBL_tracers_column_physics")
1620 endif
1621 call print_marbl_log(marbl_instances%StatusLog, g, i, j)
1622 call marbl_instances%StatusLog%erase()
1623
1624 ! iv. Apply tendencies immediately
1625 ! First pass - Euler step; if stability issues, we can do something different (subcycle?)
1626 do m=1,cs%ntr
1627 cs%tracer_data(m)%tr(i,j,:) = cs%tracer_data(m)%tr(i,j,:) + (dt * us%T_to_s) * &
1628 marbl_instances%interior_tendencies(m,:)
1629 enddo
1630
1631 ! v. Copy output that MOM6 needs to hold on to
1632 ! * saved state
1633 do m=1,size(marbl_instances%interior_tendency_saved_state%state)
1634 cs%interior_tendency_saved_state(m)%field_3d(i,j,:) = &
1635 marbl_instances%interior_tendency_saved_state%state(m)%field_3d(:,1)
1636 enddo
1637
1638 ! * diagnostics
1639 do m=1,size(marbl_instances%interior_tendency_diags%diags)
1640 if (cs%interior_tendency_diags(m)%id > 0) then
1641 if (allocated(cs%interior_tendency_diags(m)%field_2d)) then
1642 ! Only copy values if ref_depth < bathyT
1643 if (g%bathyT(i,j) > real(marbl_instances%interior_tendency_diags%diags(m)%ref_depth)) then
1644 cs%interior_tendency_diags(m)%field_2d(i,j) = &
1645 real(marbl_instances%interior_tendency_diags%diags(m)%field_2d(1))
1646 endif
1647 else ! not a 2D diagnostic
1648 cs%interior_tendency_diags(m)%field_3d(i,j,:) = &
1649 real(marbl_instances%interior_tendency_diags%diags(m)%field_3d(:,1))
1650 endif
1651 endif
1652 enddo
1653
1654 ! * tendency values themselves (and vertical integrals of them)
1655 do m=1,cs%ntr
1656 if (allocated(cs%interior_tendency_out(m)%field_3d)) &
1657 cs%interior_tendency_out(m)%field_3d(i,j,:) = marbl_instances%interior_tendencies(m,:)
1658
1659 if (allocated(cs%interior_tendency_out_zint(m)%field_2d)) &
1660 cs%interior_tendency_out_zint(m)%field_2d(i,j) = (sum(dz(:) * &
1661 marbl_instances%interior_tendencies(m,:)))
1662
1663 if (allocated(cs%interior_tendency_out_zint_100m(m)%field_2d)) then
1664 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) = 0.
1665 do k=1,gv%ke
1666 if (zi(k) < us%m_to_Z * 100.) then
1667 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) = &
1668 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) + gv%H_to_m * dz(k) * &
1669 marbl_instances%interior_tendencies(m,k)
1670 elseif (zi(k-1) < us%m_to_Z * 100.) then
1671 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) = &
1672 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) + gv%H_to_m * dz(k) * &
1673 ((us%m_to_Z * 100. - zi(k-1)) / (zi(k) - zi(k-1))) * &
1674 marbl_instances%interior_tendencies(m,k)
1675 else
1676 exit
1677 endif
1678 enddo
1679 endif
1680 enddo
1681
1682 ! * Interior tendency output
1683 do m=1,cs%ito_cnt
1684 cs%ITO(i,j,:,m) = &
1685 marbl_instances%interior_tendency_output%outputs_for_GCM(m)%forcing_field_1d(1,:)
1686 enddo
1687
1688 enddo
1689 enddo
1690
1691 if (cs%debug) then
1692 do m=1,cs%ntr
1693 call hchksum(cs%tracer_data(m)%tr(:,:,m), &
1694 trim(marbl_instances%tracer_metadata(m)%short_name)//' post source-sink', g%HI)
1695 enddo
1696 endif
1697
1698 if (associated(fluxes%salt_flux)) then
1699 ! convert salt flux to tracer fluxes and add to STF
1700 do j=js,je ; do i=is,ie
1701 net_salt_rate(i,j) = (1000.0*us%ppt_to_S * fluxes%salt_flux(i,j)) * gv%RZ_to_H
1702 enddo ; enddo
1703
1704 ! DIC related tracers
1705 do j=js,je ; do i=is,ie
1706 flux_from_salt_flux(i,j) = (cs%DIC_salt_ratio * gv%H_to_Z) * net_salt_rate(i,j)
1707 enddo ; enddo
1708 m = cs%tracer_inds%dic_ind
1709 if (m > 0) then
1710 do j=js,je ; do i=is,ie
1711 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)
1712 enddo ; enddo
1713 if (cs%id_surface_flux_from_salt_flux(m) > 0) &
1714 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)
1715 endif
1716 m = cs%tracer_inds%dic_alt_co2_ind
1717 if (m > 0) then
1718 do j=js,je ; do i=is,ie
1719 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)
1720 enddo ; enddo
1721 if (cs%id_surface_flux_from_salt_flux(m) > 0) &
1722 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)
1723 endif
1724 m = cs%tracer_inds%abio_dic_ind
1725 if (m > 0) then
1726 do j=js,je ; do i=is,ie
1727 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)
1728 enddo ; enddo
1729 if (cs%id_surface_flux_from_salt_flux(m) > 0) &
1730 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)
1731 endif
1732 m = cs%tracer_inds%abio_di14c_ind
1733 if (m > 0) then
1734 do j=js,je ; do i=is,ie
1735 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)
1736 enddo ; enddo
1737 if (cs%id_surface_flux_from_salt_flux(m) > 0) &
1738 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)
1739 endif
1740
1741 ! ALK related tracers
1742 do j=js,je ; do i=is,ie
1743 flux_from_salt_flux(i,j) = (cs%ALK_salt_ratio * gv%H_to_Z) * net_salt_rate(i,j)
1744 enddo ; enddo
1745 m = cs%tracer_inds%alk_ind
1746 if (m > 0) then
1747 do j=js,je ; do i=is,ie
1748 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)
1749 enddo ; enddo
1750 if (cs%id_surface_flux_from_salt_flux(m) > 0) &
1751 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)
1752 endif
1753 m = cs%tracer_inds%alk_alt_co2_ind
1754 if (m > 0) then
1755 do j=js,je ; do i=is,ie
1756 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)
1757 enddo ; enddo
1758 if (cs%id_surface_flux_from_salt_flux(m) > 0) &
1759 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)
1760 endif
1761 endif
1762
1763 if (cs%debug) then
1764 do m=1,cs%ntr
1765 call hchksum(cs%STF(:,:,m), &
1766 trim(marbl_instances%tracer_metadata(m)%short_name)//" sfc_flux", g%HI, &
1767 unscale=us%Z_to_m*us%s_to_T)
1768 enddo
1769 endif
1770
1771 ! (2) Apply surface fluxes via vertical diffusion
1772 ! Compute KPP nonlocal term if necessary
1773 if (present(kpp_csp)) then
1774 if (associated(kpp_csp) .and. present(nonlocaltrans)) then
1775 do m=1,cs%ntr
1776 call kpp_nonlocaltransport(kpp_csp, g, gv, h_old, nonlocaltrans, cs%STF(:,:,m), dt, &
1777 cs%diag, cs%tracer_data(m)%tr_ptr, cs%tracer_data(m)%tr(:,:,:), &
1778 flux_scale=gv%Z_to_H)
1779 enddo
1780 endif
1781 if (cs%debug) then
1782 do m=1,cs%ntr
1783 call hchksum(cs%tracer_data(m)%tr(:,:,m), &
1784 trim(marbl_instances%tracer_metadata(m)%short_name)//' post KPP', g%HI)
1785 enddo
1786 endif
1787 endif
1788
1789 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
1790 do m=1,cs%ntr
1791 do k=1,nz ;do j=js,je ; do i=is,ie
1792 h_work(i,j,k) = h_old(i,j,k)
1793 enddo ; enddo ; enddo
1794 ! CS%RIV_FLUXES is conc m/s, in_flux_optional expects time-integrated flux (conc H)
1795 do j=js,je ; do i=is,ie
1796 riv_flux_loc(i,j) = (cs%RIV_FLUXES(i,j,m) * (dt*us%T_to_s)) * gv%m_to_H
1797 enddo ; enddo
1798 if (cs%debug) &
1799 call hchksum(riv_flux_loc(:,:), &
1800 trim(marbl_instances%tracer_metadata(m)%short_name)//' riv flux', g%HI, unscale=gv%H_to_m)
1801 call applytracerboundaryfluxesinout(g, gv, cs%tracer_data(m)%tr(:,:,:) , dt, fluxes, h_work, &
1802 evap_cfl_limit, minimum_forcing_depth, in_flux_optional=riv_flux_loc)
1803 call tracer_vertdiff(h_work, ea, eb, dt, cs%tracer_data(m)%tr(:,:,:), g, gv, &
1804 sfc_flux=gv%Rho0 * cs%STF(:,:,m))
1805 enddo
1806 else
1807 ! TODO: do we want to support these options? does not apply river fluxes!
1808 ! an alternative would be to require evap_CFL_limit and minimum_forcing_depth.
1809 ! Much like we now require prediabatic_T and prediabatic_S, we can abort
1810 ! in tracer flow control if they are not present.
1811 do m=1,cs%ntr
1812 call tracer_vertdiff(h_old, ea, eb, dt, cs%tracer_data(m)%tr(:,:,:), g, gv, &
1813 sfc_flux=gv%Rho0 * cs%STF(:,:,m))
1814 enddo
1815 endif
1816
1817 if (cs%debug) then
1818 do m=1,cs%ntr
1819 call hchksum(cs%tracer_data(m)%tr(:,:,m), &
1820 trim(marbl_instances%tracer_metadata(m)%short_name)//' post tracer_vertdiff', g%HI)
1821 enddo
1822 endif
1823
1824 ! (3) Post diagnostics from our buffer
1825 ! i. surface fluxes and their diagnostics (currently all 2D)
1826 ! ii. Interior tendency diagnostics (mix of 2D and 3D)
1827 ! iii. Interior tendencies themselves
1828 ! iv. Forcing fields
1829 do m=1,cs%ntr
1830 if (cs%id_surface_flux_out(m) > 0) &
1831 call post_data(cs%id_surface_flux_out(m), cs%STF(:,:,m), cs%diag)
1832 enddo
1833
1834 do m=1,size(cs%surface_flux_diags)
1835 if (cs%surface_flux_diags(m)%id > 0) &
1836 call post_data(cs%surface_flux_diags(m)%id, cs%surface_flux_diags(m)%field_2d(:,:), cs%diag)
1837 enddo
1838
1839 if (cs%bot_flux_to_tend_id > 0) &
1840 call post_data(cs%bot_flux_to_tend_id, bot_flux_to_tend(:, :, :), cs%diag)
1841
1842 do m=1,size(cs%interior_tendency_diags)
1843 if (cs%interior_tendency_diags(m)%id > 0) then
1844 if (allocated(cs%interior_tendency_diags(m)%field_2d)) then
1845 if (real(marbl_instances%interior_tendency_diags%diags(m)%ref_depth) == 0.) then
1846 call post_data(cs%interior_tendency_diags(m)%id, &
1847 cs%interior_tendency_diags(m)%field_2d(:,:), cs%diag)
1848 else ! non-zero ref-depth
1849 ref_mask(:, :) = 0.
1850 do j=js,je ; do i=is,ie
1851 if (g%bathyT(i,j) > real(marbl_instances%interior_tendency_diags%diags(m)%ref_depth)) &
1852 ref_mask(i,j) = 1.
1853 enddo ; enddo
1854 call post_data(cs%interior_tendency_diags(m)%id, &
1855 cs%interior_tendency_diags(m)%field_2d(:,:), cs%diag, mask=ref_mask(:,:))
1856 endif
1857 elseif (allocated(cs%interior_tendency_diags(m)%field_3d)) then
1858 call post_data(cs%interior_tendency_diags(m)%id, &
1859 cs%interior_tendency_diags(m)%field_3d(:,:,:), cs%diag)
1860 else
1861 write(log_message, "(A, I0, A, I0, A)") "Diagnostic number ", m, " post id ", &
1862 cs%interior_tendency_diags(m)%id," did not allocate 2D or 3D array"
1863 call mom_error(fatal, log_message)
1864 endif
1865 endif
1866 enddo
1867
1868 do m=1,cs%ntr
1869 if (allocated(cs%interior_tendency_out(m)%field_3d)) &
1870 call post_data(cs%interior_tendency_out(m)%id, &
1871 cs%interior_tendency_out(m)%field_3d(:,:,:), cs%diag)
1872 if (allocated(cs%interior_tendency_out_zint(m)%field_2d)) &
1873 call post_data(cs%interior_tendency_out_zint(m)%id, &
1874 cs%interior_tendency_out_zint(m)%field_2d(:,:), cs%diag)
1875 if (allocated(cs%interior_tendency_out_zint_100m(m)%field_2d)) &
1876 call post_data(cs%interior_tendency_out_zint_100m(m)%id, &
1877 cs%interior_tendency_out_zint_100m(m)%field_2d(:,:), cs%diag)
1878 enddo
1879
1880 if (cs%ice_ncat > 0) then
1881 do m=1,cs%ice_ncat+1
1882 if (cs%fracr_cat_id(m) > 0) &
1883 call post_data(cs%fracr_cat_id(m), fluxes%fracr_cat(:,:,m), cs%diag)
1884 if (cs%qsw_cat_id(m) > 0) &
1885 call post_data(cs%qsw_cat_id(m), fluxes%qsw_cat(:,:,m), cs%diag)
1886 enddo
1887 endif
1888
1889
1890end subroutine marbl_tracers_column_physics
1891
1892!> This subroutine reads time-varying forcing from files
1893subroutine marbl_tracers_set_forcing(day_start, G, CS)
1894
1895 type(time_type), intent(in) :: day_start !< Start time of the fluxes.
1896 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1897 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a
1898
1899 real, parameter :: donriv_refract = 0.1 ! Fraction of DON river nutrients in refractory pools [1]
1900 real, parameter :: docriv_refract = 0.2 ! Fraction of DOC river nutrients in refractory pools [1]
1901 real, parameter :: dopriv_refract = 0.025 ! Fraction of DOP river nutrients in refractory pools [1]
1902
1903 real, dimension(SZI_(G),SZJ_(G)) :: riv_flux_in !< The field read in from forcing file with time dimension
1904 !! [mmol m-2 s-1]
1905 type(time_type) :: time_forcing !< For reading river flux fields, we use a modified version of Time
1906 integer :: i, j, k, is, ie, js, je, m
1907
1908 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1909
1910 ! Abiotic DIC forcing
1911 if (cs%abio_dic_on) then
1912 ! Read d14c bands
1913 do m=1,3
1914 time_forcing = map_model_time_to_forcing_time(day_start, cs%d14c_dataset(m))
1915 call time_interp_external(cs%id_d14c(m), time_forcing, cs%d14c_bands(m))
1916 enddo
1917
1918 ! Set d14c according to the bands
1919 do j=js,je ; do i=is,ie
1920 if (g%geoLatT(i,j) > 30.) then
1921 cs%d14c(i,j) = cs%d14c_bands(1)
1922 elseif (g%geoLatT(i,j) > -30.) then
1923 cs%d14c(i,j) = cs%d14c_bands(2)
1924 else
1925 cs%d14c(i,j) = cs%d14c_bands(3)
1926 endif
1927 enddo ; enddo
1928 endif
1929
1930 ! River fluxes
1931 if (cs%read_riv_fluxes) then
1932 cs%RIV_FLUXES(:,:,:) = 0.
1933 time_forcing = map_model_time_to_forcing_time(day_start, cs%riv_flux_dataset)
1934
1935 ! DIN river flux affects NO3, ALK, and ALK_ALT_CO2
1936 call time_interp_external(cs%id_din_riv, time_forcing, riv_flux_in)
1937
1938 if (cs%tracer_inds%no3_ind > 0) then
1939 do j=js,je ; do i=is,ie
1940 cs%RIV_FLUXES(i,j,cs%tracer_inds%no3_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)
1941 enddo ; enddo
1942 endif
1943 if (cs%tracer_inds%alk_ind > 0) then
1944 do j=js,je ; do i=is,ie
1945 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) = cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) - &
1946 g%mask2dT(i,j) *riv_flux_in(i,j)
1947 enddo ; enddo
1948 endif
1949 if (cs%tracer_inds%alk_alt_co2_ind > 0) then
1950 do j=js,je ; do i=is,ie
1951 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) = &
1952 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) - g%mask2dT(i,j) *riv_flux_in(i,j)
1953 enddo ; enddo
1954 endif
1955
1956 call time_interp_external(cs%id_dip_riv, time_forcing, riv_flux_in)
1957 if (cs%tracer_inds%po4_ind > 0) then
1958 do j=js,je ; do i=is,ie
1959 cs%RIV_FLUXES(i,j,cs%tracer_inds%po4_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)
1960 enddo ; enddo
1961 endif
1962
1963 call time_interp_external(cs%id_don_riv, time_forcing, riv_flux_in)
1964 if (cs%tracer_inds%don_ind > 0) then
1965 do j=js,je ; do i=is,ie
1966 cs%RIV_FLUXES(i,j,cs%tracer_inds%don_ind) = g%mask2dT(i,j) * (1. - donriv_refract) * &
1967 riv_flux_in(i,j)
1968 enddo ; enddo
1969 endif
1970 if (cs%tracer_inds%donr_ind > 0) then
1971 do j=js,je ; do i=is,ie
1972 cs%RIV_FLUXES(i,j,cs%tracer_inds%donr_ind) = g%mask2dT(i,j) * donriv_refract * &
1973 riv_flux_in(i,j)
1974 enddo ; enddo
1975 endif
1976
1977 call time_interp_external(cs%id_dop_riv, time_forcing, riv_flux_in)
1978 if (cs%tracer_inds%dop_ind > 0) then
1979 do j=js,je ; do i=is,ie
1980 cs%RIV_FLUXES(i,j,cs%tracer_inds%dop_ind) = g%mask2dT(i,j) * (1. - dopriv_refract) * &
1981 riv_flux_in(i,j)
1982 enddo ; enddo
1983 endif
1984 if (cs%tracer_inds%dopr_ind > 0) then
1985 do j=js,je ; do i=is,ie
1986 cs%RIV_FLUXES(i,j,cs%tracer_inds%dopr_ind) = g%mask2dT(i,j) * dopriv_refract * &
1987 riv_flux_in(i,j)
1988 enddo ; enddo
1989 endif
1990
1991 call time_interp_external(cs%id_dsi_riv, time_forcing, riv_flux_in)
1992 if (cs%tracer_inds%sio3_ind > 0) then
1993 do j=js,je ; do i=is,ie
1994 cs%RIV_FLUXES(i,j,cs%tracer_inds%sio3_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)
1995 enddo ; enddo
1996 endif
1997
1998 call time_interp_external(cs%id_dfe_riv, time_forcing, riv_flux_in)
1999 if (cs%tracer_inds%fe_ind > 0) then
2000 do j=js,je ; do i=is,ie
2001 cs%RIV_FLUXES(i,j,cs%tracer_inds%fe_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)
2002 enddo ; enddo
2003 endif
2004
2005 call time_interp_external(cs%id_dic_riv, time_forcing, riv_flux_in)
2006 if (cs%tracer_inds%dic_ind > 0) then
2007 do j=js,je ; do i=is,ie
2008 cs%RIV_FLUXES(i,j,cs%tracer_inds%dic_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)
2009 enddo ; enddo
2010 endif
2011 if (cs%tracer_inds%dic_alt_co2_ind > 0) then
2012 do j=js,je ; do i=is,ie
2013 cs%RIV_FLUXES(i,j,cs%tracer_inds%dic_alt_co2_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)
2014 enddo ; enddo
2015 endif
2016
2017 call time_interp_external(cs%id_alk_riv, time_forcing, riv_flux_in)
2018 if (cs%tracer_inds%alk_ind > 0) then
2019 do j=js,je ; do i=is,ie
2020 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) = cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) + &
2021 g%mask2dT(i,j) *riv_flux_in(i,j)
2022 enddo ; enddo
2023 endif
2024 if (cs%tracer_inds%alk_alt_co2_ind > 0) then
2025 do j=js,je ; do i=is,ie
2026 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) = &
2027 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) + g%mask2dT(i,j) * riv_flux_in(i,j)
2028 enddo ; enddo
2029 endif
2030
2031 call time_interp_external(cs%id_doc_riv, time_forcing, riv_flux_in)
2032 if (cs%tracer_inds%doc_ind > 0) then
2033 do j=js,je ; do i=is,ie
2034 cs%RIV_FLUXES(i,j,cs%tracer_inds%doc_ind) = g%mask2dT(i,j) * (1. - docriv_refract) * &
2035 riv_flux_in(i,j)
2036 enddo ; enddo
2037 endif
2038 if (cs%tracer_inds%docr_ind > 0) then
2039 do j=js,je ; do i=is,ie
2040 cs%RIV_FLUXES(i,j,cs%tracer_inds%docr_ind) = g%mask2dT(i,j) * docriv_refract * &
2041 riv_flux_in(i,j)
2042 enddo ; enddo
2043 endif
2044 endif
2045
2046 ! Tracer restoring
2047 do m=1,cs%restore_count
2048 call time_interp_external(cs%id_tracer_restoring(m),day_start,cs%restoring_in(:,:,:,m))
2049 do k=1,cs%restoring_nz ; do j=js,je ; do i=is,ie
2050 cs%restoring_in(i,j,k,m) = g%mask2dT(i,j) * cs%restoring_in(i,j,k,m)
2051 enddo ; enddo ; enddo
2052 enddo
2053
2054 ! Post Forcing to Diagnostics
2055 if (cs%read_riv_fluxes) then
2056 if (cs%no3_riv_flux > 0 .and. cs%tracer_inds%no3_ind > 0) &
2057 call post_data(cs%no3_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%no3_ind), cs%diag)
2058 if (cs%po4_riv_flux > 0 .and. cs%tracer_inds%po4_ind > 0) &
2059 call post_data(cs%po4_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%po4_ind), cs%diag)
2060 if (cs%don_riv_flux > 0 .and. cs%tracer_inds%don_ind > 0) &
2061 call post_data(cs%don_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%don_ind), cs%diag)
2062 if (cs%donr_riv_flux > 0 .and. cs%tracer_inds%donr_ind > 0) &
2063 call post_data(cs%donr_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%donr_ind), cs%diag)
2064 if (cs%dop_riv_flux > 0 .and. cs%tracer_inds%dop_ind > 0) &
2065 call post_data(cs%dop_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dop_ind), cs%diag)
2066 if (cs%dopr_riv_flux > 0 .and. cs%tracer_inds%dopr_ind > 0) &
2067 call post_data(cs%dopr_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dopr_ind), cs%diag)
2068 if (cs%sio3_riv_flux > 0 .and. cs%tracer_inds%sio3_ind > 0) &
2069 call post_data(cs%sio3_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%sio3_ind), cs%diag)
2070 if (cs%fe_riv_flux > 0 .and. cs%tracer_inds%fe_ind > 0) &
2071 call post_data(cs%fe_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%fe_ind), cs%diag)
2072 if (cs%doc_riv_flux > 0 .and. cs%tracer_inds%doc_ind > 0) &
2073 call post_data(cs%doc_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%doc_ind), cs%diag)
2074 if (cs%docr_riv_flux > 0 .and. cs%tracer_inds%docr_ind > 0) &
2075 call post_data(cs%docr_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%docr_ind), cs%diag)
2076 if (cs%alk_riv_flux > 0 .and. cs%tracer_inds%alk_ind > 0) &
2077 call post_data(cs%alk_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%alk_ind), cs%diag)
2078 if (cs%alk_alt_co2_riv_flux > 0 .and. cs%tracer_inds%alk_alt_co2_ind > 0) &
2079 call post_data(cs%alk_alt_co2_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%alk_alt_co2_ind), &
2080 cs%diag)
2081 if (cs%dic_riv_flux > 0 .and. cs%tracer_inds%dic_ind > 0) &
2082 call post_data(cs%dic_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dic_ind), cs%diag)
2083 if (cs%dic_alt_co2_riv_flux > 0 .and. cs%tracer_inds%dic_alt_co2_ind > 0) &
2084 call post_data(cs%dic_alt_co2_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dic_alt_co2_ind), &
2085 cs%diag)
2086 endif
2087 if (cs%abio_dic_on) then
2088 if (cs%d14c_id > 0) &
2089 call post_data(cs%d14c_id, cs%d14c, cs%diag)
2090 endif
2091
2092end subroutine marbl_tracers_set_forcing
2093
2094!> This function calculates the mass-weighted integral of all tracer stocks,
2095!! returning the number of stocks it has calculated. If the stock_index
2096!! is present, only the stock corresponding to that coded index is returned.
2097function marbl_tracers_stock(h, stocks, G, GV, CS, names, units, stock_index)
2098 real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
2099 type(efp_type), dimension(:), intent(out) :: stocks !< the mass-weighted integrated amount of
2100 !! each tracer, in kg times concentration units
2101 !! [kg conc].
2102 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
2103 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
2104 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a
2105 !! previous call to register_MARBL_tracers.
2106 character(len=*), dimension(:), intent(out) :: names !< the names of the stocks calculated.
2107 character(len=*), dimension(:), intent(out) :: units !< the units of the stocks calculated.
2108 integer, optional, intent(in) :: stock_index !< the coded index of a specific stock
2109 !! being sought.
2110 integer :: marbl_tracers_stock !< Return value: the number of stocks
2111 !! calculated here.
2112
2113 ! Local variables
2114 integer :: i, j, k, is, ie, js, je, nz, m
2115 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2116
2118 if (.not.associated(cs)) return
2119 if (cs%ntr < 1) return
2120
2121 if (present(stock_index)) then ; if (stock_index > 0) then
2122 ! Check whether this stock is available from this routine.
2123
2124 ! No stocks from this routine are being checked yet. Return 0.
2125 return
2126 endif ; endif
2127
2128 do m=1,cs%ntr
2129 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller="MARBL_tracers_stock")
2130 units(m) = trim(units(m))//" kg"
2131 stocks(m) = global_mass_int_efp(h, g, gv, cs%tracer_data(m)%tr(:,:,:), on_pe_only=.true.)
2132 enddo
2133 marbl_tracers_stock = cs%ntr
2134
2135end function marbl_tracers_stock
2136
2137!> This subroutine extracts the surface fields from this tracer package that
2138!! are to be shared with the atmosphere in coupled configurations.
2139subroutine marbl_tracers_surface_state(sfc_state, G, US, CS)
2140 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
2141 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
2142 !! describe the surface state of the ocean.
2143 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2144 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a previous
2145 !! call to register_MARBL_tracers.
2146
2147 ! Local variables
2148 integer :: i, j, is, ie, js, je
2149
2150 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2151
2152 if (.not.associated(cs)) return
2153
2154 if (allocated(sfc_state%fco2)) then
2155 do j=js,je ; do i=is,ie
2156 ! 44e-6 converts mmol/m^2/s (positive down) to kg CO2/m^2/s (positive down)
2157 sfc_state%fco2(i,j) = us%kg_m2s_to_RZ_T * (44.0e-6 * cs%SFO(i,j,cs%flux_co2_ind))
2158 enddo ; enddo
2159 endif
2160
2161end subroutine marbl_tracers_surface_state
2162
2163!> Copy the requested interior tendency output field into an array.
2164subroutine marbl_tracers_get(name, G, GV, array, CS)
2165
2166 character(len=*), intent(in) :: name !< Name of requested tracer.
2167 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
2168 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
2169 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2170 intent(inout) :: array !< Array filled by this routine.
2171 type(marbl_tracers_cs), pointer :: cs !< Pointer to the control structure for this module.
2172
2173 character(len=128), parameter :: sub_name = 'MARBL_tracers_get'
2174 character(len=128) :: log_message
2175
2176 array(:,:,:) = 0.0
2177 select case(trim(name))
2178 case ('Chl')
2179 array(:,:,:) = cs%ITO(:,:,:,cs%total_Chl_ind)
2180 case DEFAULT
2181 write(log_message, "(3A)") "'", trim(name), &
2182 "' is not a valid interior tendency output field name"
2183 call mom_error(fatal, log_message)
2184 end select
2185
2186end subroutine marbl_tracers_get
2187
2188!> Clean up any allocated memory after the run.
2189subroutine marbl_tracers_end(CS)
2190 type(marbl_tracers_cs), pointer, intent(inout) :: cs !< The control structure returned by a previous
2191 !! call to register_MARBL_tracers.
2192
2193 integer :: m
2194
2195 call print_marbl_log(marbl_instances%StatusLog)
2196 call marbl_instances%StatusLog%erase()
2197 call marbl_instances%shutdown()
2198 ! TODO: print MARBL timers to stdout as well
2199
2200 if (associated(cs)) then
2201 if (allocated(cs%tracer_data)) then
2202 do m=1,cs%ntr
2203 if (associated(cs%tracer_data(m)%tr)) deallocate(cs%tracer_data(m)%tr)
2204 enddo
2205 deallocate(cs%tracer_data)
2206 endif
2207 if (allocated(cs%ind_tr)) deallocate(cs%ind_tr)
2208 if (allocated(cs%id_surface_flux_out)) deallocate(cs%id_surface_flux_out)
2209 if (allocated(cs%interior_tendency_out)) deallocate(cs%interior_tendency_out)
2210 if (allocated(cs%interior_tendency_out_zint)) deallocate(cs%interior_tendency_out_zint)
2211 if (allocated(cs%interior_tendency_out_zint_100m)) &
2212 deallocate(cs%interior_tendency_out_zint_100m)
2213 if (allocated(cs%fracr_cat_id)) deallocate(cs%fracr_cat_id)
2214 if (allocated(cs%qsw_cat_id)) deallocate(cs%qsw_cat_id)
2215 if (allocated(cs%STF)) deallocate(cs%STF)
2216 if (allocated(cs%RIV_FLUXES)) deallocate(cs%RIV_FLUXES)
2217 if (associated(cs%SFO)) then
2218 deallocate(cs%SFO)
2219 nullify(cs%SFO)
2220 endif
2221 if (associated(cs%ITO)) then
2222 deallocate(cs%ITO)
2223 nullify(cs%ITO)
2224 endif
2225 if (allocated(cs%tracer_restoring_ind)) deallocate(cs%tracer_restoring_ind)
2226 if (allocated(cs%tracer_I_tau_ind)) deallocate(cs%tracer_I_tau_ind)
2227 if (allocated(cs%fesedflux_in)) deallocate(cs%fesedflux_in)
2228 if (allocated(cs%fesedfluxred_in)) deallocate(cs%fesedfluxred_in)
2229 if (allocated(cs%feventflux_in)) deallocate(cs%feventflux_in)
2230 if (allocated(cs%I_tau)) deallocate(cs%I_tau)
2231 deallocate(cs)
2232 endif
2233end subroutine marbl_tracers_end
2234
2235subroutine set_riv_flux_tracer_inds(CS)
2236
2237 type(marbl_tracers_cs), pointer, intent(inout) :: CS !< The MARBL tracers control structure
2238
2239 character(len=256) :: log_message
2240 character(len=48) :: name ! A variable's name in a NetCDF file.
2241 integer :: m
2242
2243 ! Initialize tracers from file (unless they were initialized by restart file)
2244 ! Also save indices of tracers that have river fluxes
2245 cs%tracer_inds%no3_ind = 0
2246 cs%tracer_inds%po4_ind = 0
2247 cs%tracer_inds%don_ind = 0
2248 cs%tracer_inds%donr_ind = 0
2249 cs%tracer_inds%dop_ind = 0
2250 cs%tracer_inds%dopr_ind = 0
2251 cs%tracer_inds%sio3_ind = 0
2252 cs%tracer_inds%fe_ind = 0
2253 cs%tracer_inds%doc_ind = 0
2254 cs%tracer_inds%docr_ind = 0
2255 cs%tracer_inds%alk_ind = 0
2256 cs%tracer_inds%alk_alt_co2_ind = 0
2257 cs%tracer_inds%dic_ind = 0
2258 cs%tracer_inds%dic_alt_co2_ind = 0
2259 cs%tracer_inds%abio_dic_ind = 0
2260 cs%tracer_inds%abio_di14c_ind = 0
2261 do m=1,cs%ntr
2262 name = marbl_instances%tracer_metadata(m)%short_name
2263 if (trim(name) == "NO3") then
2264 cs%tracer_inds%no3_ind = m
2265 elseif (trim(name) == "PO4") then
2266 cs%tracer_inds%po4_ind = m
2267 elseif (trim(name) == "DON") then
2268 cs%tracer_inds%don_ind = m
2269 elseif (trim(name) == "DONr") then
2270 cs%tracer_inds%donr_ind = m
2271 elseif (trim(name) == "DOP") then
2272 cs%tracer_inds%dop_ind = m
2273 elseif (trim(name) == "DOPr") then
2274 cs%tracer_inds%dopr_ind = m
2275 elseif (trim(name) == "SiO3") then
2276 cs%tracer_inds%sio3_ind = m
2277 elseif (trim(name) == "Fe") then
2278 cs%tracer_inds%fe_ind = m
2279 elseif (trim(name) == "DOC") then
2280 cs%tracer_inds%doc_ind = m
2281 elseif (trim(name) == "DOCr") then
2282 cs%tracer_inds%docr_ind = m
2283 elseif (trim(name) == "ALK") then
2284 cs%tracer_inds%alk_ind = m
2285 elseif (trim(name) == "ALK_ALT_CO2") then
2286 cs%tracer_inds%alk_alt_co2_ind = m
2287 elseif (trim(name) == "DIC") then
2288 cs%tracer_inds%dic_ind = m
2289 elseif (trim(name) == "DIC_ALT_CO2") then
2290 cs%tracer_inds%dic_alt_co2_ind = m
2291 elseif (trim(name) == "ABIO_DIC") then
2292 cs%tracer_inds%abio_dic_ind = m
2293 elseif (trim(name) == "ABIO_DI14C") then
2294 cs%tracer_inds%abio_di14c_ind = m
2295 endif
2296 enddo
2297
2298 ! Log indices for each tracer to ensure we set them all correctly
2299 write(log_message, "(A,I0)") "NO3 index: ", cs%tracer_inds%no3_ind
2300 call mom_error(note, log_message)
2301 write(log_message, "(A,I0)") "PO4 index: ", cs%tracer_inds%po4_ind
2302 call mom_error(note, log_message)
2303 write(log_message, "(A,I0)") "DON index: ", cs%tracer_inds%don_ind
2304 call mom_error(note, log_message)
2305 write(log_message, "(A,I0)") "DONr index: ", cs%tracer_inds%donr_ind
2306 call mom_error(note, log_message)
2307 write(log_message, "(A,I0)") "DOP index: ", cs%tracer_inds%dop_ind
2308 call mom_error(note, log_message)
2309 write(log_message, "(A,I0)") "DOPr index: ", cs%tracer_inds%dopr_ind
2310 call mom_error(note, log_message)
2311 write(log_message, "(A,I0)") "SiO3 index: ", cs%tracer_inds%sio3_ind
2312 call mom_error(note, log_message)
2313 write(log_message, "(A,I0)") "Fe index: ", cs%tracer_inds%fe_ind
2314 call mom_error(note, log_message)
2315 write(log_message, "(A,I0)") "DOC index: ", cs%tracer_inds%doc_ind
2316 call mom_error(note, log_message)
2317 write(log_message, "(A,I0)") "DOCr index: ", cs%tracer_inds%docr_ind
2318 call mom_error(note, log_message)
2319 write(log_message, "(A,I0)") "ALK index: ", cs%tracer_inds%alk_ind
2320 call mom_error(note, log_message)
2321 write(log_message, "(A,I0)") "ALK_ALT_CO2 index: ", cs%tracer_inds%alk_alt_co2_ind
2322 call mom_error(note, log_message)
2323 write(log_message, "(A,I0)") "DIC index: ", cs%tracer_inds%dic_ind
2324 call mom_error(note, log_message)
2325 write(log_message, "(A,I0)") "DIC_ALT_CO2 index: ", cs%tracer_inds%dic_alt_co2_ind
2326 call mom_error(note, log_message)
2327
2328end subroutine set_riv_flux_tracer_inds
2329
2330! TODO: some log messages come from a specific grid point, and this routine
2331! needs to include the location in the preamble
2332!> This subroutine writes the contents of the MARBL log using MOM_error(NOTE, ...).
2333subroutine print_marbl_log(log_to_print, G, i, j)
2334
2336 use marbl_logging, only : marbl_log_type
2337 use mom_coms, only : pe_here
2338
2339 class(marbl_log_type), intent(in) :: log_to_print !< MARBL log to include in MOM6 logfile
2340 type(ocean_grid_type), optional, intent(in) :: G !< The ocean's grid structure
2341 integer, optional, intent(in) :: i !< i of (i,j) index of column providing the log
2342 integer, optional, intent(in) :: j !< j of (i,j) index of column providing the log
2343
2344 character(len=*), parameter :: subname = 'MARBL_tracers:print_marbl_log'
2345 character(len=256) :: message_prefix, message_location, log_message
2346 type(marbl_status_log_entry_type), pointer :: tmp
2347 integer :: msg_lev, elem_old
2348
2349 ! elem_old is used to keep track of whether all messages are coming from the same point
2350 elem_old = -1
2351 write(message_prefix, "(A,I0,A)") '(Task ', pe_here(), ')'
2352
2353 tmp => log_to_print%FullLog
2354 do while (associated(tmp))
2355 ! 1) Do I need to write this message? Yes, if all tasks should write this
2356 ! or if I am master_task
2357 if ((.not. tmp%lonly_master_writes) .or. is_root_pe()) then
2358 ! 2) Print message location? (only if ElementInd changed and is positive; requires G)
2359 if ((present(g)) .and. (tmp%ElementInd .ne. elem_old)) then
2360 if (tmp%ElementInd .gt. 0) then
2361 if (present(i) .and. present(j)) then
2362 write(message_location, "(A,F8.3,A,F7.3,A,I0,A,I0,A,I0)") &
2363 'Message from (lon, lat) (', g%geoLonT(i,j), ', ', g%geoLatT(i,j), &
2364 '), which is global (i,j) (', i + g%HI%idg_offset, ', ', j + g%HI%jdg_offset, &
2365 '). Level: ', tmp%ElementInd
2366 else
2367 write(message_location, "(A)") "Grid cell responsible for message is unknown"
2368 endif ! i,j present
2369 ! master task does not need prefix
2370 if (is_root_pe()) then
2371 write(log_message, "(A)") trim(message_location)
2372 msg_lev = note
2373 else
2374 write(log_message, "(A,1X,A)") trim(message_prefix), trim(message_location)
2375 msg_lev = warning
2376 endif ! print message prefix?
2377 call mom_error(msg_lev, log_message, all_print=.true.)
2378 endif ! ElementInd > 0
2379 elem_old = tmp%ElementInd
2380 endif ! ElementInd /= elem_old
2381
2382 ! 3) Write message from the log
2383 ! master task does not need prefix
2384 if (is_root_pe()) then
2385 write(log_message, "(A)") trim(tmp%LogMessage)
2386 msg_lev = note
2387 else
2388 write(log_message, "(A,1X,A)") trim(message_prefix), trim(tmp%LogMessage)
2389 msg_lev = warning
2390 endif ! print message prefix?
2391 call mom_error(msg_lev, log_message, all_print=.true.)
2392 endif ! write the message?
2393 tmp => tmp%next
2394 enddo
2395
2396 if (log_to_print%labort_marbl) then
2397 call mom_error(warning, 'ERROR reported from MARBL library', all_print=.true.)
2398 call mom_error(fatal, 'Stopping in ' // subname)
2399 endif
2400
2401end subroutine print_marbl_log
2402
2403!> \namespace MARBL_tracers
2404!!
2405!! This module contains the code that is needed to provide
2406!! the MARBL BGC tracer library with necessary forcings and
2407!! apply the resulting surface fluxes and tendencies to the
2408!! requested tracers.
2409
2410end module marbl_tracers