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

14use mom_diag_mediator, only : diag_ctrl

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

20use mom_cvmix_kpp, only : kpp_nonlocaltransport, kpp_cs

21use mom_hor_index, only : hor_index_type

22use mom_interpolate, only : forcing_timeseries_dataset

23use mom_interpolate, only : forcing_timeseries_set_time_type_vars

24use mom_interpolate, only : map_model_time_to_forcing_time

25use mom_io, only : file_exists, mom_read_data, slasher, vardesc, var_desc, query_vardesc

26use mom_open_boundary, only : ocean_obc_type

27use mom_remapping, only : reintegrate_column

28use mom_remapping, only : remapping_cs, initialize_remapping, remapping_core_h

29use mom_restart, only : query_initialized, mom_restart_cs, register_restart_field

30use mom_spatial_means, only : global_mass_int_efp

31use mom_sponge, only : set_up_sponge_field, sponge_cs

32use mom_time_manager, only : time_type

33use mom_tracer_registry, only : register_tracer

34use mom_tracer_types, only : tracer_type, tracer_registry_type

35use mom_tracer_diabatic, only : tracer_vertdiff, applytracerboundaryfluxesinout

36use mom_tracer_initialization_from_z, only : mom_initialize_tracer_from_z

37use mom_tracer_z_init, only : read_z_edges

38use mom_unit_scaling, only : unit_scale_type

39use mom_variables, only : surface, thermo_var_ptrs

40use mom_verticalgrid, only : verticalgrid_type

41use mom_diag_mediator, only : register_diag_field, post_data!, safe_alloc_ptr

42

43use marbl_interface, only : marbl_interface_class

44use marbl_interface_public_types, only : marbl_diagnostics_type, marbl_saved_state_type

45

46use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux

47

48implicit none ; private

49

50#include <MOM_memory.h>

51

52public register_marbl_tracers, initialize_marbl_tracers

53public marbl_tracers_column_physics, marbl_tracers_surface_state

54public marbl_tracers_set_forcing

55public marbl_tracers_stock, marbl_tracers_get, marbl_tracers_end

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

92type :: saved_state_for_marbl_type

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

98end type saved_state_for_marbl_type

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) :: feventflux_file !< name of [netCDF] file containing iron vent flux

142 type(forcing_timeseries_dataset) :: d14c_dataset(3) !< File and time axis information for d14c forcing

143 real, dimension(3) :: d14c_bands !< forcing is organized into bands: [30 N, 90 N]; [30 S, 30 N]; [90 S, 30 S]

144 !! This variable contains D14C for each band [CU ~> conc]

145 integer :: d14c_id !< id for diagnostic field with d14c forcing

146 logical :: read_riv_fluxes !< If true, use river fluxes supplied from an input file.

147 !! This is temporary, we will always read river fluxes

148 type(forcing_timeseries_dataset) :: riv_flux_dataset !< File and time axis information for river fluxes

149 character(len=4) :: restoring_source !< location of tracer restoring data

150 !! valid values: file, none

151 integer :: restoring_nz !< number of levels in tracer restoring file

152 real, allocatable, dimension(:) :: &

153 restoring_z_edges !< The depths of the cell interfaces in the tracer restoring file [Z ~> m]

154 real, allocatable, dimension(:) :: &

155 restoring_dz !< The thickness of the cell layers in the tracer restoring file [H ~> m]

156 integer :: restoring_timescale_nz !< number of levels in tracer restoring timescale file

157 real, allocatable, dimension(:) :: &

158 restoring_timescale_z_edges !< The depths of the cell interfaces in the tracer restoring timescale file [Z ~> m]

159 real, allocatable, dimension(:) :: &

160 restoring_timescale_dz !< The thickness of the cell layers in the tracer restoring timescale file [H ~> m]

161 character(len=14) :: restoring_i_tau_source !< location of inverse restoring timescale data

162 !! valid values: file, grid_dependent

163 character(len=200) :: restoring_file !< name of [netCDF] file containing tracer restoring data

164 type(remapping_cs) :: restoring_remapcs !< Remapping parameters and work arrays for tracer restoring / timescale

165 character(len=200) :: restoring_i_tau_file !< name of [netCDF] file containing inverse restoring timescale

166 character(len=200) :: restoring_i_tau_var_name !< name of field containing inverse restoring timescale

167 character(len=35) :: marbl_settings_file !< name of [text] file containing MARBL settings

168

169 real :: bot_flux_mix_thickness !< for bottom flux -> tendency conversion, assume uniform mixing over

170 !! bottom layer of prescribed thickness [Z ~> m]

171 real :: ibfmt !< Reciprocal of bot_flux_mix_thickness [Z-1 ~> m-1]

172

173 type(temp_marbl_diag), allocatable :: surface_flux_diags(:) !< collect surface flux diagnostics from all columns

174 !! before posting

175 type(temp_marbl_diag), allocatable :: interior_tendency_diags(:) !< collect tendency diagnostics from all columns

176 !! before posting

177 type(saved_state_for_marbl_type), allocatable :: surface_flux_saved_state(:) !< surface_flux saved state

178 type(saved_state_for_marbl_type), allocatable :: interior_tendency_saved_state(:) !< interior_tendency saved state

179

180 ! TODO: If we can post data column by column, all we need are integer arrays for ids

181 ! integer, allocatable :: id_surface_flux_diags(:) !< array of indices for surface_flux diagnostics

182 ! integer, allocatable :: id_interior_tendency_diags(:) !< array of indices for interior_tendency diagnostics

183

184 type(tracer_ind_type) :: tracer_inds !< Indices to tracers that will have river fluxes added to STF

185

186 !> Need to store global output from both marbl_instance%surface_flux_compute() and

187 !! marbl_instance%interior_tendency_compute(). For the former, just need id to register

188 !! because we already copy data into CS%STF; latter requires copying data and indices

189 !! so currently using temp_MARBL_diag for that.

190 integer, allocatable :: id_surface_flux_out(:) !< register_diag indices for surface_flux output

191 integer, allocatable :: id_surface_flux_from_salt_flux(:) !< register_diag indices for surface_flux from salt_flux

192 type(temp_marbl_diag), allocatable :: interior_tendency_out(:) !< collect interior tendencies for diagnostic output

193 type(temp_marbl_diag), allocatable :: interior_tendency_out_zint(:) !< vertical integral of interior tendencies

194 !! (full column)

195 type(temp_marbl_diag), allocatable :: interior_tendency_out_zint_100m(:) !< vertical integral of interior tendencies

196 !! (top 100m)

197 integer :: bot_flux_to_tend_id !< register_diag index for BOT_FLUX_TO_TEND

198 integer, allocatable :: fracr_cat_id(:) !< register_diag index for per-category ice fraction

199 integer, allocatable :: qsw_cat_id(:) !< register_diag index for per-category shortwave

200

201 real :: dic_salt_ratio !< ratio to convert salt surface flux to DIC surface flux [conc ppt-1]

202 real :: alk_salt_ratio !< ratio to convert salt surface flux to ALK surface flux [conc ppt-1]

203

204 real, allocatable :: stf(:,:,:) !< surface fluxes returned from MARBL to use in tracer_vertdiff

205 !! (dims: i, j, tracer) [conc Z T-1 ~> conc m s-1]

206 real, pointer :: sfo(:,:,:) => null() !< surface flux output returned from MARBL for use in GCM

207 !! e.g. CO2 flux to pass to atmosphere (dims: i, j, num_sfo)

208 !! Units vary based on index of num_sfo dimension

209 real, pointer :: ito(:,:,:,:) => null() !< interior tendency output returned from MARBL for use in GCM

210 !! e.g. total chlorophyll to use in shortwave penetration

211 !! (dims: i, j, k, num_ito)

212 !! Units vary based on index of num_ito dimension

213

214 integer :: u10_sqr_ind !< index of MARBL forcing field array to copy 10-m wind (squared) into

215 integer :: sss_ind !< index of MARBL forcing field array to copy sea surface salinity into

216 integer :: sst_ind !< index of MARBL forcing field array to copy sea surface temperature into

217 integer :: ifrac_ind !< index of MARBL forcing field array to copy ice fraction into

218 integer :: dust_dep_ind !< index of MARBL forcing field array to copy dust flux into

219 integer :: fe_dep_ind !< index of MARBL forcing field array to copy iron flux into

220 integer :: nox_flux_ind !< index of MARBL forcing field array to copy NOx flux into

221 integer :: nhy_flux_ind !< index of MARBL forcing field array to copy NHy flux into

222 integer :: atmpress_ind !< index of MARBL forcing field array to copy atmospheric pressure into

223 integer :: xco2_ind !< index of MARBL forcing field array to copy CO2 flux into

224 integer :: xco2_alt_ind !< index of MARBL forcing field array to copy CO2 flux (alternate CO2) into

225 integer :: d14c_ind !< index of MARBL forcing field array to copy d14C into

226

227 !> external_field types for river fluxes (added to surface fluxes)

228 type(external_field) :: id_din_riv !< id for time_interp_external.

229 type(external_field) :: id_don_riv !< id for time_interp_external.

230 type(external_field) :: id_dip_riv !< id for time_interp_external.

231 type(external_field) :: id_dop_riv !< id for time_interp_external.

232 type(external_field) :: id_dsi_riv !< id for time_interp_external.

233 type(external_field) :: id_dfe_riv !< id for time_interp_external.

234 type(external_field) :: id_dic_riv !< id for time_interp_external.

235 type(external_field) :: id_alk_riv !< id for time_interp_external.

236 type(external_field) :: id_doc_riv !< id for time_interp_external.

237

238 !> external_field type for d14c (needed if abio_dic_on is True)

239 type(external_field) :: id_d14c(3) !< id for time_interp_external.

240

241 !> Indices for river fluxes (diagnostics)

242 integer :: no3_riv_flux !< NO3 riverine flux

243 integer :: po4_riv_flux !< PO4 riverine flux

244 integer :: don_riv_flux !< DON riverine flux

245 integer :: donr_riv_flux !< DONr riverine flux

246 integer :: dop_riv_flux !< DOP riverine flux

247 integer :: dopr_riv_flux !< DOPr riverine flux

248 integer :: sio3_riv_flux !< SiO3 riverine flux

249 integer :: fe_riv_flux !< Fe riverine flux

250 integer :: doc_riv_flux !< DOC riverine flux

251 integer :: docr_riv_flux !< DOCr riverine flux

252 integer :: alk_riv_flux !< ALK riverine flux

253 integer :: alk_alt_co2_riv_flux !< ALK (alternate CO2) riverine flux

254 integer :: dic_riv_flux !< DIC riverine flux

255 integer :: dic_alt_co2_riv_flux !< DIC (alternate CO2) riverine flux

256

257 !> Indices for forcing fields required to compute interior tendencies

258 integer :: dustflux_ind !< index of MARBL forcing field array to copy dust flux into

259 integer :: par_col_frac_ind !< index of MARBL forcing field array to copy PAR column fraction into

260 integer :: surf_shortwave_ind !< index of MARBL forcing field array to copy surface shortwave into

261 integer :: potemp_ind !< index of MARBL forcing field array to copy potential temperature into

262 integer :: salinity_ind !< index of MARBL forcing field array to copy salinity into

263 integer :: pressure_ind !< index of MARBL forcing field array to copy pressure into

264 integer :: fesedflux_ind !< index of MARBL forcing field array to copy iron sediment flux into

265 integer :: o2_scalef_ind !< index of MARBL forcing field array to copy O2 scale length into

266 integer :: remin_scalef_ind !< index of MARBL forcing field array to copy remin scale length into

267 type(external_field), allocatable :: id_tracer_restoring(:) !< id number for time_interp_external

268 integer, allocatable :: tracer_restoring_ind(:) !< index of MARBL forcing field to copy

269 !! per-tracer restoring field into

270 integer, allocatable :: tracer_i_tau_ind(:) !< index of MARBL forcing field to copy per-tracer

271 !! inverse restoring timescale into

272

273 !> Memory for storing river fluxes, tracer restoring fields, and abiotic forcing

274 real, allocatable :: d14c(:,:) !< d14c forcing for abiotic DIC and carbon isotope tracer modules

275 !! [mmol m-3 s-1]

276 real, allocatable :: riv_fluxes(:,:,:) !< river flux forcing for applyTracerBoundaryFluxesInOut

277 !! (needs to be time-integrated when passed to function!)

278 !! (dims: i, j, tracer) [conc m s-1]

279 character(len=15), allocatable :: tracer_restoring_varname(:) !< name of variable being restored

280 real, allocatable :: i_tau(:,:,:) !< inverse restoring timescale for marbl tracers (dims: i, j, k) [s-1]

281 real, allocatable, dimension(:,:,:,:) :: restoring_in !< Restoring fields read from file

282 !! (dims: i, j, restoring_nz, restoring_cnt) [tracer units]

283

284 !> Number of surface flux outputs as well as specific indices for each one

285 integer :: sfo_cnt !< number of surface flux outputs from MARBL

286 integer :: ito_cnt !< number of interior tendency outputs from MARBL

287 integer :: flux_co2_ind !< index to co2 flux surface flux output

288 integer :: total_chl_ind !< index to total chlorophyll interior tendency output

289

290 ! TODO: create generic 3D forcing input type to read z coordinate + values

291 real :: fesedflux_scale_factor !< scale factor for iron sediment flux [mmol umol-1 d s-1]

292 integer :: fesedflux_nz !< number of levels in iron sediment flux file

293 real, allocatable, dimension(:,:,:) :: fesedflux_in !< Field to read iron sediment flux into [conc m s-1]

294 real, allocatable, dimension(:,:,:) :: feventflux_in !< Field to read iron vent flux into [conc m s-1]

295 real, allocatable, dimension(:) :: &

296 fesedflux_z_edges !< The depths of the cell interfaces in the input data [Z ~> m]

297 ! TODO: this thickness does not need to be 3D, but it is easier to make thickness 0

298 ! below the surface on a per-column basis (could save memory by storing 1D

299 ! thickness from file and then computing a second 1D thickness array in (i,j) loop)

300 real, allocatable, dimension(:,:,:) :: &

301 fesedflux_dz !< The thickness of the cell layers in the input data [H ~> m]

302end type marbl_tracers_cs

303

304! Module parameters

305real, parameter :: atm_per_pa = 1./101325. !< convert from Pa -> atm [atm Pa-1]

306

307contains

308

309!> This subroutine is used to read marbl_in, configure MARBL accordingly, and then

310!! call MARBL's initialization routine

311subroutine configure_marbl_tracers(GV, US, param_file, CS)

312 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure

313 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

314 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters

315 type(marbl_tracers_cs), pointer :: CS !< A pointer that is set to point to the control

316 !! structure for this module

317

318# include "version_variable.h"

319 character(len=40) :: mdl = "MARBL_tracers" ! This module's name.

320 character(len=256) :: log_message

321 character(len=256) :: marbl_in_line(1)

322 character(len=256) :: forcing_sname, field_source

323 integer :: m, n, nz, marbl_settings_in, read_error, I_tau_count, fi

324 logical :: chl_from_file, forcing_processed

325 nz = gv%ke

326 marbl_settings_in = 615

327

328 ! (1) Read parameters necessary for general setup of MARBL

329 call log_version(param_file, mdl, version, "")

330 call get_param(param_file, mdl, "DEBUG", cs%debug, "If true, write out verbose debugging data.", &

331 default=.false., debuggingparam=.true.)

332 call get_param(param_file, mdl, "MARBL_IC_MIN_VAL", cs%IC_min, &

333 "Minimum value of tracer initial conditions (set to 1e-100 for dim scaling tests)", &

334 default=0., units="tracer units")

335 call get_param(param_file, mdl, "MARBL_SETTINGS_FILE", cs%marbl_settings_file, &

336 "The name of a file from which to read the run-time settings for MARBL.", default="marbl_in")

337 call get_param(param_file, mdl, "BOT_FLUX_MIX_THICKNESS", cs%bot_flux_mix_thickness, &

338 "Bottom fluxes are uniformly mixed over layer of this thickness", default=1., units="m", &

339 scale=us%m_to_Z)

340 call get_param(param_file, mdl, "USE_ICE_CATEGORIES", cs%use_ice_category_fields, &

341 "If true, allocate memory for shortwave and ice fraction split by ice thickness category.", &

342 default=.false.)

343 call get_param(param_file, mdl, "ICE_NCAT", cs%ice_ncat, &

344 "Number of ice thickness categories in shortwave and ice fraction forcings.", default=0)

345 cs%Ibfmt = 1. / cs%bot_flux_mix_thickness

346

347 if (cs%use_ice_category_fields .and. (cs%ice_ncat == 0)) &

348 call mom_error(fatal, &

349 "Can not configure MARBL to use multiple ice categories without ice_ncat present")

350

351 ! (2) Read marbl settings file and call put_setting()

352 ! (2a) only master task opens file

353 if (is_root_pe()) then

354 ! read the marbl_in into buffer

355 open(unit=marbl_settings_in, file=cs%marbl_settings_file, iostat=read_error)

356 if (read_error .ne. 0) then

357 write(log_message, '(A, I0, 2A)') "IO ERROR ", read_error, " opening namelist file : ", &

358 trim(cs%marbl_settings_file)

359 call mom_error(fatal, log_message)

360 endif

361 endif

362

363 ! (2b) master task reads file and broadcasts line-by-line

364 marbl_in_line = ''

365 do

366 ! i. Read next line on master, iostat value out

367 ! (Exit loop if read is not successful; either read error or end of file)

368 if (is_root_pe()) read(marbl_settings_in, "(A)", iostat=read_error) marbl_in_line(1)

369 call broadcast(read_error, root_pe())

370 if (read_error .ne. 0) exit

371

372 ! ii. Broadcast line just read in on root PE to all tasks

373 call broadcast(marbl_in_line, 256, root_pe())

374

375 ! iii. All tasks call put_setting (TODO: openMP blocks?)

376 call marbl_instances%put_setting(marbl_in_line(1))

377 enddo

378

379 ! (2c) we should always reach the EOF to capture the entire file...

380 if (.not. is_iostat_end(read_error)) then

381 write(log_message, '(3A, I0)') "IO ERROR reading ", trim(cs%marbl_settings_file), ": ", &

382 read_error

383 call mom_error(fatal, log_message)

384 else

385 if (is_root_pe()) then

386 write(log_message, '(3A)') "Read '", trim(cs%marbl_settings_file), "' until EOF."

387 call mom_error(note, log_message)

388 endif

389 endif

390 if (is_root_pe()) close(marbl_settings_in)

391

392 ! (3) Initialize MARBL and configure MOM6 accordingly

393

394 ! (3a) call marbl%init()

395 ! TODO: We want to strip gcm_delta_z, gcm_zw, and gcm_zt values out of

396 ! init because MOM updates them every time step / every column

397 call marbl_instances%init(gcm_num_levels = nz, gcm_num_par_subcols = cs%ice_ncat + 1, &

398 gcm_num_elements_surface_flux = 1, & ! FIXME: change to number of grid cells on MPI task

399 gcm_delta_z = gv%sInterface(2:nz+1) - gv%sInterface(1:nz), gcm_zw = gv%sInterface(2:nz+1), &

400 gcm_zt = gv%sLayer, unit_system_opt = "mks", lgcm_has_global_ops = .false.) ! FIXME: add global ops

401 ! Regardless of vertical grid, MOM6 will always use GV%ke levels in all columns

402 marbl_instances%domain%kmt = gv%ke

403 if (marbl_instances%StatusLog%labort_marbl) &

404 call marbl_instances%StatusLog%log_error_trace("MARBL_instances%init", &

405 "configure_MARBL_tracers")

406 call print_marbl_log(marbl_instances%StatusLog)

407 call marbl_instances%StatusLog%erase()

408 cs%ntr = size(marbl_instances%tracer_metadata)

409 call marbl_instances%get_setting('base_bio_on', cs%base_bio_on)

410 call marbl_instances%get_setting('abio_dic_on', cs%abio_dic_on)

411 call marbl_instances%get_setting('ciso_on', cs%ciso_on)

412

413 ! (3b) Read parameters that depend on how MARBL is configured

414 if (cs%base_bio_on) then

415 call get_param(param_file, mdl, "CHL_FROM_FILE", chl_from_file, &

416 "If true, chl_a is read from a file.", default=.true.)

417 cs%request_Chl_from_MARBL = (.not. chl_from_file)

418 else

419 cs%request_Chl_from_MARBL = .false.

420 endif

421

422 ! (4) Request fields needed by MOM6

423 cs%sfo_cnt = 0

424 cs%ito_cnt = 0

425 cs%flux_co2_ind = -1

426 cs%total_Chl_ind = -1

427

428 if (cs%base_bio_on) then

429 ! CO2 Flux to the atmosphere

430 call marbl_instances%add_output_for_GCM(num_elements=1, field_name="flux_co2", &

431 output_id=cs%flux_co2_ind, field_source=field_source)

432 if (trim(field_source) == "surface_flux") then

433 cs%sfo_cnt = cs%sfo_cnt + 1

434 else if (trim(field_source) == "interior_tendency") then

435 cs%ito_cnt = cs%ito_cnt + 1

436 endif

437

438 ! Total 3D Chlorophyll

439 call marbl_instances%add_output_for_GCM(num_elements=1, num_levels=nz, field_name="total_Chl", &

440 output_id=cs%total_Chl_ind, field_source=field_source)

441 if (trim(field_source) == "surface_flux") then

442 cs%sfo_cnt = cs%sfo_cnt + 1

443 else if (trim(field_source) == "interior_tendency") then

444 cs%ito_cnt = cs%ito_cnt + 1

445 endif

446 endif

447

448 ! (5) Initialize forcing fields

449 ! i. store all surface forcing indices

450 cs%u10_sqr_ind = -1

451 cs%sss_ind = -1

452 cs%sst_ind = -1

453 cs%ifrac_ind = -1

454 cs%dust_dep_ind = -1

455 cs%fe_dep_ind = -1

456 cs%nox_flux_ind = -1

457 cs%nhy_flux_ind = -1

458 cs%atmpress_ind = -1

459 cs%xco2_ind = -1

460 cs%xco2_alt_ind = -1

461 do m=1,size(marbl_instances%surface_flux_forcings)

462 select case (trim(marbl_instances%surface_flux_forcings(m)%metadata%varname))

463 case('u10_sqr')

464 cs%u10_sqr_ind = m

465 case('sss')

466 cs%sss_ind = m

467 case('sst')

468 cs%sst_ind = m

469 case('Ice Fraction')

470 cs%ifrac_ind = m

471 case('Dust Flux')

472 cs%dust_dep_ind = m

473 case('Iron Flux')

474 cs%fe_dep_ind = m

475 case('NOx Flux')

476 cs%nox_flux_ind = m

477 case('NHy Flux')

478 cs%nhy_flux_ind = m

479 case('Atmospheric Pressure')

480 cs%atmpress_ind = m

481 case('xco2')

482 cs%xco2_ind = m

483 case('xco2_alt_co2')

484 cs%xco2_alt_ind = m

485 case('d14c')

486 cs%d14c_ind = m

487 case DEFAULT

488 write(log_message, "(A,1X,A)") &

489 trim(marbl_instances%surface_flux_forcings(m)%metadata%varname), &

490 'is not a valid surface flux forcing field name.'

491 call mom_error(fatal, log_message)

492 end select

493 enddo

494

495 ! ii. store all interior forcing indices

496 cs%dustflux_ind = -1

497 cs%PAR_col_frac_ind = -1

498 cs%surf_shortwave_ind = -1

499 cs%potemp_ind = -1

500 cs%salinity_ind = -1

501 cs%pressure_ind = -1

502 cs%fesedflux_ind = -1

503 cs%o2_scalef_ind = -1

504 cs%remin_scalef_ind = -1

505 cs%d14c_ind = -1

506 allocate(cs%id_tracer_restoring(cs%ntr))

507 allocate(cs%tracer_restoring_varname(cs%ntr), source=' ') ! gfortran 13.2 bug?

508 ! source = '' does not blank out strings

509 allocate(cs%tracer_restoring_ind(cs%ntr), source=-1)

510 allocate(cs%tracer_I_tau_ind(cs%ntr), source=-1)

511 cs%restore_count = 0

512 i_tau_count = 0

513 do m=1,size(marbl_instances%interior_tendency_forcings)

514 select case (trim(marbl_instances%interior_tendency_forcings(m)%metadata%varname))

515 case('Dust Flux')

516 cs%dustflux_ind = m

517 case('PAR Column Fraction')

518 cs%PAR_col_frac_ind = m

519 case('Surface Shortwave')

520 cs%surf_shortwave_ind = m

521 case('Potential Temperature')

522 cs%potemp_ind = m

523 case('Salinity')

524 cs%salinity_ind = m

525 case('Pressure')

526 cs%pressure_ind = m

527 case('Iron Sediment Flux')

528 cs%fesedflux_ind = m

529 case('O2 Consumption Scale Factor')

530 cs%o2_scalef_ind = m

531 case('Particulate Remin Scale Factor')

532 cs%remin_scalef_ind = m

533 case DEFAULT

534 ! fi stands for forcing_index

535 fi = index(marbl_instances%interior_tendency_forcings(m)%metadata%varname, &

536 'Restoring Field')

537 if (fi > 0) then

538 cs%restore_count = cs%restore_count + 1

539 cs%tracer_restoring_ind(cs%restore_count) = m

540 cs%tracer_restoring_varname(cs%restore_count) = &

541 marbl_instances%interior_tendency_forcings(m)%metadata%varname(1:fi-2)

542 else

543 fi = index(marbl_instances%interior_tendency_forcings(m)%metadata%varname, &

544 'Restoring Inverse Timescale')

545 if (fi > 0) then

546 i_tau_count = i_tau_count + 1

547 cs%tracer_I_tau_ind(i_tau_count) = m

548 else

549 write(log_message, "(A,1X,A)") &

550 trim(marbl_instances%interior_tendency_forcings(m)%metadata%varname), &

551 'is not a valid interior tendency forcing field name.'

552 call mom_error(fatal, log_message)

553 endif

554 endif

555 end select

556 enddo

557end subroutine configure_marbl_tracers

558

559!> This subroutine is used to register tracer fields and subroutines

560!! to be used with MOM.

561function register_marbl_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS, MARBL_computes_chl)

562 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.

563 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure

564 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

565 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters

566 type(marbl_tracers_cs), pointer :: cs !< A pointer that is set to point to the control

567 !! structure for this module

568 type(tracer_registry_type), pointer :: tr_reg !< A pointer that is set to point to the control

569 !! structure for the tracer advection and diffusion module.

570 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control struct

571 logical, intent(out) :: marbl_computes_chl !< If MARBL is computing chlorophyll, MOM

572 !! may use it to compute SW penetration

573

574! Local variables

575! This include declares and sets the variable "version".

576# include "version_variable.h"

577 character(len=40) :: mdl = "MARBL_tracers" ! This module's name.

578 character(len=256) :: log_message

579 character(len=200) :: inputdir ! The directory where the input files are.

580 character(len=48) :: var_name ! The variable's name.

581 character(len=128) :: desc_name ! The variable's descriptor.

582 character(len=48) :: units ! The variable's units.

583 character(len=96) :: file_name ! file name for d14c (looped over three bands)

584 real, pointer :: tr_ptr(:,:,:) => null() ! Pointer to 3D tracer array [CU ~> conc]

585 ! (ALK tracers use meq m-3 instead of mmol m-3)

586 integer :: forcing_file_start_year

587 integer :: forcing_file_end_year

588 integer :: forcing_file_data_ref_year

589 integer :: forcing_file_model_ref_year

590 integer :: forcing_file_forcing_year

591 logical :: register_marbl_tracers

592 ! read_Z_edges() has several mandatory arguments that we do not use given our expectation

593 ! of how the file being read in was created

594 logical :: z_edges_has_edges

595 logical :: z_edges_use_missing

596 real :: z_edges_missing ! required argument for read_Z_edges() [CU ~> conc]

597 integer :: isd, ied, jsd, jed, nz, m, k, kbot

598 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke

599

600 if (associated(cs)) then

601 call mom_error(warning, "register_MARBL_tracers called with an associated control structure.")

602 return

603 endif

604 allocate(cs)

605

606 call configure_marbl_tracers(gv, us, param_file, cs)

607 marbl_computes_chl = cs%base_bio_on

608

609 ! Read all relevant parameters and write them to the model log.

610 call log_version(param_file, mdl, version, "")

611 ! ** Input directory

612 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")

613 ! ** Tracer initial conditions

614 call get_param(param_file, mdl, "MARBL_TRACERS_IC_FILE", cs%IC_file, &

615 "The file in which the MARBL tracers initial values can be found.", &

616 default="ecosys_jan_IC_omip_latlon_1x1_180W_c230331.nc")

617 if (scan(cs%IC_file,'/') == 0) then

618 ! Add the directory if CS%IC_file is not already a complete path.

619 cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)

620 call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_IC_FILE", cs%IC_file)

621 endif

622 call get_param(param_file, mdl, "MARBL_TRACERS_MAY_REINIT", cs%tracers_may_reinit, &

623 "If true, tracers may go through the initialization code if they are not found in the "//&

624 "restart files. Otherwise it is a fatal error if tracers are not found in the "//&

625 "restart files of a restarted run.", default=.false.)

626 call get_param(param_file, mdl, "MARBL_TRACERS_INIT_VERTICAL_REMAP_ONLY", cs%ongrid, &

627 "If true, initial conditions are on the model horizontal grid. Extrapolation over " //&

628 "missing ocean values is done using an ICE-9 procedure with vertical ALE remapping .", &

629 default=.false.)

630 if (cs%base_bio_on) then

631 ! ** FESEDFLUX

632 call get_param(param_file, mdl, "MARBL_FESEDFLUX_FILE", cs%fesedflux_file, &

633 "The file in which the iron sediment flux forcing field can be found.", &

634 default="fesedflux_total_reduce_oxic_tx0.66v1.c230817.nc")

635 if (scan(cs%fesedflux_file,'/') == 0) then

636 ! Add the directory if CS%fesedflux_file is not already a complete path.

637 cs%fesedflux_file = trim(slasher(inputdir))//trim(cs%fesedflux_file)

638 call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_FESEDFLUX_FILE", cs%fesedflux_file)

639 endif

640 ! ** FEVENTFLUX

641 call get_param(param_file, mdl, "MARBL_FEVENTFLUX_FILE", cs%feventflux_file, &

642 "The file in which the iron vent flux forcing field can be found.", &

643 default="feventflux_5gmol_tx0.66v1.c230817.nc")

644 if (scan(cs%feventflux_file,'/') == 0) then

645 ! Add the directory if CS%feventflux_file is not already a complete path.

646 cs%feventflux_file = trim(slasher(inputdir))//trim(cs%feventflux_file)

647 call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_FEVENTFLUX_FILE", cs%feventflux_file)

648 endif

649 ! ** Scale factor for FESEDFLUX

650 call get_param(param_file, mdl, "MARBL_FESEDFLUX_SCALE_FACTOR", cs%fesedflux_scale_factor, &

651 "Conversion factor between FESEDFLUX file units and MARBL units", &

652 units="umol m-2 d-1 -> mmol m-2 s-1", default=0.001/86400.)

653

654 ! ** River fluxes

655 call get_param(param_file, mdl, "READ_RIV_FLUXES", cs%read_riv_fluxes, &

656 "If true, use river fluxes supplied from an input file", default=.true.)

657 if (cs%read_riv_fluxes) then

658 call get_param(param_file, mdl, "RIV_FLUX_FILE", cs%riv_flux_dataset%file_name, &

659 "The file in which the river fluxes can be found", &

660 default="riv_nut.gnews_gnm.JRA025m_to_tx0.66v1_nnsm_e333r100_190910.20210405.nc")

661 ! call get_param(param_file, mdl, "RIV_FLUX_OFFSET_YEAR", CS%riv)

662 if (scan(cs%riv_flux_dataset%file_name,'/') == 0) then

663 ! CS%riv_flux_dataset%file_name = trim(inputdir) // trim(CS%riv_flux_dataset%file_name)

664 cs%riv_flux_dataset%file_name = trim(slasher(inputdir)) //&

665 trim(cs%riv_flux_dataset%file_name)

666 call log_param(param_file, mdl, "INPUTDIR/RIV_FLUX_FILE", cs%riv_flux_dataset%file_name)

667 endif

668 call get_param(param_file, mdl, "RIV_FLUX_L_TIME_VARYING", &

669 cs%riv_flux_dataset%l_time_varying, &

670 ".true. for time-varying forcing, .false. for static forcing", default=.false.)

671 if (cs%riv_flux_dataset%l_time_varying) then

672 call get_param(param_file, mdl, "RIV_FLUX_FILE_START_YEAR", forcing_file_start_year, &

673 "First year of data to read in RIV_FLUX_FILE", default=1900)

674 call get_param(param_file, mdl, "RIV_FLUX_FILE_END_YEAR", forcing_file_end_year, &

675 "Last year of data to read in RIV_FLUX_FILE", default=2000)

676 call get_param(param_file, mdl, "RIV_FLUX_FILE_DATA_REF_YEAR", forcing_file_data_ref_year, &

677 "Align this year in RIV_FLUX_FILE with RIV_FLUX_FILE_MODEL_REF_YEAR in model", &

678 default=1900)

679 call get_param(param_file, mdl, "RIV_FLUX_FILE_MODEL_REF_YEAR", &

680 forcing_file_model_ref_year, &

681 "Align this year in model with RIV_FLUX_FILE_DATA_REF_YEAR in RIV_FLUX_FILE", &

682 default=1)

683 else

684 call get_param(param_file, mdl, "RIV_FLUX_FORCING_YEAR", forcing_file_forcing_year, &

685 "Year from RIV_FLUX_FILE to use for forcing", default=1900)

686 endif

687 call forcing_timeseries_set_time_type_vars(forcing_file_start_year, forcing_file_end_year, &

688 forcing_file_data_ref_year, forcing_file_model_ref_year, forcing_file_forcing_year, &

689 cs%riv_flux_dataset)

690 endif

691 endif

692

693 if (cs%abio_dic_on) then

694 call get_param(param_file, mdl, "D14C_L_TIME_VARYING", cs%d14c_dataset(1)%l_time_varying, &

695 ".true. for time-varying forcing, .false. for static forcing", default=.false.)

696 cs%d14c_dataset(2)%l_time_varying = cs%d14c_dataset(1)%l_time_varying

697 cs%d14c_dataset(3)%l_time_varying = cs%d14c_dataset(1)%l_time_varying

698 if (cs%d14c_dataset(1)%l_time_varying) then

699 call get_param(param_file, mdl, "D14C_FILE_START_YEAR", forcing_file_start_year, &

700 "First year of data to read in D14C_FILE", default=1850)

701 call get_param(param_file, mdl, "D14C_FILE_END_YEAR", forcing_file_end_year, &

702 "Last year of data to read in D14C_FILE", default=2015)

703 call get_param(param_file, mdl, "D14C_FILE_DATA_REF_YEAR", forcing_file_data_ref_year, &

704 "Align this year in D14C_FILE with D14C_FILE_MODEL_REF_YEAR in model", default=1850)

705 call get_param(param_file, mdl, "D14C_FILE_MODEL_REF_YEAR", forcing_file_model_ref_year, &

706 "Align this year in model with D14C_FILE_DATA_REF_YEAR in D14C_FILE", default=1)

707 else

708 call get_param(param_file, mdl, "D14C_FORCING_YEAR", forcing_file_forcing_year, &

709 "Year from D14C_FILE to use for forcing", default=1850)

710 endif

711 do m=1,3

712 write(var_name, "(A,I0)") "MARBL_D14C_FILE_", m

713 write(file_name, "(A,I0,A)") "atm_delta_C14_CMIP6_sector", m, &

714 "_global_1850-2015_yearly_v2.0_c240202.nc"

715 call get_param(param_file, mdl, var_name, cs%d14c_dataset(m)%file_name, &

716 "The file in which the d14c forcing field can be found.", default=file_name)

717 call forcing_timeseries_set_time_type_vars(forcing_file_start_year, forcing_file_end_year, &

718 forcing_file_data_ref_year, forcing_file_model_ref_year, forcing_file_forcing_year, &

719 cs%d14c_dataset(m))

720 if (scan(cs%d14c_dataset(m)%file_name,'/') == 0) then

721 ! Add the directory if CS%d14c_dataset%file_name is not already a complete path.

722 cs%d14c_dataset(m)%file_name = trim(slasher(inputdir))//trim(cs%d14c_dataset(m)%file_name)

723 call log_param(param_file, mdl, "INPUTDIR/D14C_FILE", cs%d14c_dataset(m)%file_name)

724 endif

725 enddo

726 endif

727

728 call get_param(param_file, mdl, "DIC_SALT_RATIO", cs%DIC_salt_ratio, &

729 "Ratio to convert salt surface flux to DIC surface flux", units="conc ppt-1", &

730 default=64.0)

731 call get_param(param_file, mdl, "ALK_SALT_RATIO", cs%ALK_salt_ratio, &

732 "Ratio to convert salt surface flux to ALK surface flux", units="conc ppt-1", &

733 default=70.0)

734

735 ! ** Tracer Restoring

736 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_SOURCE", cs%restoring_source, &

737 "Source of data for restoring MARBL tracers", default="none")

738 select case(cs%restoring_source)

739 case("none")

740 case("file")

741 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_FILE", cs%restoring_file, &

742 "File containing fields to restore MARBL tracers towards")

743 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_SOURCE", &

744 cs%restoring_I_tau_source, "Source of data for inverse timescale for restoring MARBL tracers")

745

746 ! Initialize remapping type

747 call initialize_remapping(cs%restoring_remapCS, 'PCM', boundary_extrapolation=.false., answer_date=99991231)

748

749 ! Set up array for thicknesses in restoring file

750 call read_z_edges(cs%restoring_file, "PO4", cs%restoring_z_edges, cs%restoring_nz, &

751 z_edges_has_edges, z_edges_use_missing, z_edges_missing, scale=us%m_to_Z, &

752 missing_scale=1.0)

753 allocate(cs%restoring_dz(cs%restoring_nz))

754 do k=cs%restoring_nz,1,-1

755 kbot = k + 1 ! level k is between z(k) and z(k+1)

756 cs%restoring_dz(k) = (cs%restoring_z_edges(k) - cs%restoring_z_edges(kbot)) * gv%Z_to_H

757 enddo

758

759 select case(cs%restoring_I_tau_source)

760 case("file")

761 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_FILE", &

762 cs%restoring_I_tau_file, &

763 "File containing the inverse timescale for restoring MARBL tracers")

764 call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_VAR_NAME", &

765 cs%restoring_I_tau_var_name, &

766 "Field containing the inverse timescale for restoring MARBL tracers", &

767 default="I_TAU")

768 ! Set up array for thicknesses in restoring timescale file

769 call read_z_edges(cs%restoring_I_tau_file, cs%restoring_I_tau_var_name, cs%restoring_timescale_z_edges, &

770 cs%restoring_timescale_nz, 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_timescale_dz(cs%restoring_timescale_nz))

773 do k=cs%restoring_timescale_nz,1,-1

774 kbot = k + 1 ! level k is between z(k) and z(k+1)

775 cs%restoring_timescale_dz(k) = (cs%restoring_timescale_z_edges(k) - &

776 cs%restoring_timescale_z_edges(kbot)) * gv%Z_to_H

777 enddo

778 case DEFAULT

779 write(log_message, "(3A)") "'", trim(cs%restoring_I_tau_source), &

780 "' is not a valid option for MARBL_TRACER_RESTORING_I_TAU_SOURCE"

781 call mom_error(fatal, log_message)

782 end select

783 case DEFAULT

784 write(log_message, "(3A)") "'", trim(cs%restoring_source), &

785 "' is not a valid option for MARBL_TRACER_RESTORING_SOURCE"

786 call mom_error(fatal, log_message)

787 end select

788

789 allocate(cs%ind_tr(cs%ntr))

790 allocate(cs%tr_desc(cs%ntr))

791 allocate(cs%tracer_data(cs%ntr))

792

793 do m=1,cs%ntr

794 allocate(cs%tracer_data(m)%tr(isd:ied,jsd:jed,nz), source=0.0)

795 write(var_name(:),'(A)') trim(marbl_instances%tracer_metadata(m)%short_name)

796 write(desc_name(:),'(A)') trim(marbl_instances%tracer_metadata(m)%long_name)

797 write(units(:),'(A)') trim(marbl_instances%tracer_metadata(m)%units)

798 cs%tr_desc(m) = var_desc(trim(var_name), trim(units), trim(desc_name), caller=mdl)

799

800 ! This is needed to force the compiler not to do a copy in the registration

801 ! calls. Curses on the designers and implementers of Fortran90.

802 tr_ptr => cs%tracer_data(m)%tr(:,:,:)

803 call query_vardesc(cs%tr_desc(m), name=var_name, &

804 caller="register_MARBL_tracers")

805 ! Register the tracer for horizontal advection, diffusion, and restarts.

806 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, units = units, &

807 tr_desc=cs%tr_desc(m), registry_diags=.true., &

808 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit, &

809 tr_out=cs%tracer_data(m)%tr_ptr)

810

811 ! Set coupled_tracers to be true (hard-coded above) to provide the surface

812 ! values to the coupler (if any). This is meta-code and its arguments will

813 ! currently (deliberately) give fatal errors if it is used.

814 if (cs%coupled_tracers) &

815 cs%ind_tr(m) = aof_set_coupler_flux(trim(var_name)//'_flux', &

816 flux_type=' ', implementation=' ', caller="register_MARBL_tracers")

817 enddo

818

819 ! Set up memory for saved state

820 call setup_saved_state(marbl_instances%surface_flux_saved_state, hi, gv, restart_cs, &

821 cs%tracers_may_reinit, cs%surface_flux_saved_state)

822 call setup_saved_state(marbl_instances%interior_tendency_saved_state, hi, gv, restart_cs, &

823 cs%tracers_may_reinit, cs%interior_tendency_saved_state)

824

825 ! Set up memory for additional output from MARBL and add to restart files

826 allocate(cs%SFO(szi_(hi), szj_(hi), cs%sfo_cnt), &

827 cs%ITO(szi_(hi), szj_(hi), szk_(gv), cs%ito_cnt), &

828 source=0.0)

829

830 do m=1,cs%sfo_cnt

831 write(var_name, "(2A)") 'MARBL_SFO_', &

832 trim(marbl_instances%surface_flux_output%outputs_for_GCM(m)%short_name)

833 call register_restart_field(cs%SFO(:,:,m), var_name, .false., restart_cs)

834 enddo

835

836 do m=1,cs%ito_cnt

837 write(var_name, "(2A)") 'MARBL_ITO_', &

838 trim(marbl_instances%interior_tendency_output%outputs_for_GCM(m)%short_name)

839 call register_restart_field(cs%ITO(:,:,:,m), var_name, .false., restart_cs)

840 enddo

841

842

843 cs%tr_Reg => tr_reg

844 cs%restart_CSp => restart_cs

845

846 call set_riv_flux_tracer_inds(cs)

847 register_marbl_tracers = .true.

848

849end function register_marbl_tracers

850

851!> This subroutine initializes the CS%ntr tracer fields in tr(:,:,:,:)

852!! and it sets up the tracer output.

853subroutine initialize_marbl_tracers(restart, day, G, GV, US, h, param_file, diag, OBC, CS, sponge_CSp)

854 logical, intent(in) :: restart !< .true. if the fields have already been

855 !! read from a restart file.

856 type(time_type), target, intent(in) :: day !< Time of the start of the run.

857 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure

858 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure

859 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

860 real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]

861 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters

862 type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output.

863 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies

864 !! whether, where, and what open boundary

865 !! conditions are used.

866 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a previous

867 !! call to register_MARBL_tracers.

868 type(sponge_cs), pointer :: sponge_csp !< A pointer to the control structure

869 !! for the sponges, if they are in use.

870

871 ! Local variables

872 character(len=200) :: log_message

873 character(len=48) :: name ! A variable's name in a NetCDF file.

874 character(len=100) :: longname ! The long name of that variable.

875 character(len=48) :: units ! The units of the variable.

876 character(len=48) :: flux_units ! The units for age tracer fluxes, either

877 ! years m3 s-1 or years kg s-1.

878 character(len=48) :: tracer_name

879 logical :: fesedflux_has_edges, fesedflux_use_missing

880 real :: fesedflux_missing ! required argument for read_Z_edges() [CU ~> conc]

881 integer :: i, j, k, kbot, m, diag_size

882

883 if (.not.associated(cs)) return

884 if (cs%ntr < 1) return

885

886 cs%diag => diag

887

888 ! Allocate memory for surface tracer fluxes

889 allocate(cs%STF(szi_(g), szj_(g), cs%ntr), &

890 cs%RIV_FLUXES(szi_(g), szj_(g), cs%ntr), &

891 source=0.0)

892

893 ! Allocate memory for d14c forcing

894 if (cs%abio_dic_on) allocate(cs%d14c(szi_(g), szj_(g)))

895

896 ! Register diagnostics returned from MARBL (surface flux first, then interior tendency)

897 call register_marbl_diags(marbl_instances%surface_flux_diags, diag, day, g, cs%surface_flux_diags)

898 call register_marbl_diags(marbl_instances%interior_tendency_diags, diag, day, g, &

899 cs%interior_tendency_diags)

900

901 ! Register per-tracer diagnostics computed from MARBL surface flux / interior tendency values

902 allocate(cs%id_surface_flux_out(cs%ntr))

903 allocate(cs%id_surface_flux_from_salt_flux(cs%ntr))

904 allocate(cs%interior_tendency_out(cs%ntr))

905 allocate(cs%interior_tendency_out_zint(cs%ntr))

906 allocate(cs%interior_tendency_out_zint_100m(cs%ntr))

907 do m=1,cs%ntr

908 write(name, "(2A)") "STF_", trim(marbl_instances%tracer_metadata(m)%short_name)

909 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), " Surface Flux"

910 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), " m/s"

911 cs%id_surface_flux_out(m) = register_diag_field("ocean_model", trim(name), &

912 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid

913 day, trim(longname), trim(units), conversion=us%Z_to_m*us%s_to_T)

914

915 write(name, "(2A)") "STF_SALT_", trim(marbl_instances%tracer_metadata(m)%short_name)

916 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), " Surface Flux from Salt Flux"

917 cs%id_surface_flux_from_salt_flux(m) = register_diag_field("ocean_model", trim(name), &

918 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid

919 day, trim(longname), trim(units), conversion=us%Z_to_m*us%s_to_T)

920

921 write(name, "(2A)") "J_", trim(marbl_instances%tracer_metadata(m)%short_name)

922 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), " Source Sink Term"

923 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), "/s"

924 cs%interior_tendency_out(m)%id = register_diag_field("ocean_model", trim(name), &

925 diag%axesTL, & ! T=> tracer grid? L => layer center

926 day, trim(longname), trim(units))

927 if (cs%interior_tendency_out(m)%id > 0) &

928 allocate(cs%interior_tendency_out(m)%field_3d(szi_(g),szj_(g), szk_(g)), source=0.0)

929

930 write(name, "(2A)") "Jint_", trim(marbl_instances%tracer_metadata(m)%short_name)

931 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), &

932 " Source Sink Term Vertical Integral"

933 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), " m/s"

934 cs%interior_tendency_out_zint(m)%id = register_diag_field("ocean_model", trim(name), &

935 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

936 day, trim(longname), trim(units))

937 if (cs%interior_tendency_out_zint(m)%id > 0) &

938 allocate(cs%interior_tendency_out_zint(m)%field_2d(szi_(g),szj_(g)), source=0.0)

939

940 write(name, "(2A)") "Jint_100m_", trim(marbl_instances%tracer_metadata(m)%short_name)

941 write(longname, "(2A)") trim(marbl_instances%tracer_metadata(m)%long_name), &

942 " Source Sink Term Vertical Integral, 0-100m"

943 write(units, "(2A)") trim(marbl_instances%tracer_metadata(m)%units), " m/s"

944 cs%interior_tendency_out_zint_100m(m)%id = register_diag_field("ocean_model", trim(name), &

945 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

946 day, trim(longname), trim(units))

947 if (cs%interior_tendency_out_zint_100m(m)%id > 0) &

948 allocate(cs%interior_tendency_out_zint_100m(m)%field_2d(szi_(g),szj_(g)), source=0.0)

949

950 enddo

951

952 ! Register diagnostics for MOM to report that are not tracer specific

953 cs%bot_flux_to_tend_id = register_diag_field("ocean_model", "BOT_FLUX_TO_TEND", &

954 diag%axesTL, & ! T=> tracer grid? L => layer center

955 day, "Conversion Factor for Bottom Flux -> Tend", "1/m")

956

957 ! Initialize tracers (if they weren't initialized from restart file)

958 do m=1,cs%ntr

959 call query_vardesc(cs%tr_desc(m), name=name, caller="initialize_MARBL_tracers")

960 if ((.not. restart) .or. &

961 (cs%tracers_may_reinit .and. &

962 .not. query_initialized(cs%tracer_data(m)%tr(:,:,:), name, cs%restart_CSp))) then

963 ! TODO: added the ongrid optional argument, but is there a good way to detect if the file is on grid?

964 call mom_initialize_tracer_from_z(h, cs%tracer_data(m)%tr, g, gv, us, param_file, &

965 cs%IC_file, name, ongrid=cs%ongrid)

966 do k=1,gv%ke ; do j=g%jsc, g%jec ; do i=g%isc, g%iec

967 ! Ensure tracer concentrations are at / above minimum value

968 if (cs%tracer_data(m)%tr(i,j,k) < cs%IC_min) cs%tracer_data(m)%tr(i,j,k) = cs%IC_min

969 enddo ; enddo ; enddo

970 endif

971 enddo

972

973 ! Initialize total chlorophyll to get SW Pen correct (if it wasn't initialized from restart file)

974 if ((cs%total_Chl_ind > 0) .and. &

975 ((.not. restart) .or. &

976 (.not. query_initialized(cs%ITO(:,:,:,cs%total_Chl_ind), "MARBL_ITO_total_Chl", cs%restart_CSp)))) then

977 ! Three steps per column

978 do j=g%jsc, g%jec ; do i=g%isc, g%iec

979 ! (i) Copy initial tracers into MARBL structure

980 do k=1,gv%ke ; do m=1,cs%ntr

981 marbl_instances%tracers(m,k) = max(cs%tracer_data(m)%tr(i,j,k), 0.)

982 enddo ; enddo

983 ! (ii) Compute total Chl for the column

984 call marbl_instances%compute_totChl()

985 ! (iii) Copy total Chl from MARBL data-structure into CS%ITO

986 do k=1,gv%ke

987 cs%ITO(i,j,k,cs%total_Chl_ind) = &

988 marbl_instances%interior_tendency_output%outputs_for_GCM(cs%total_Chl_ind)%forcing_field_1d(1,k)

989 enddo

990 enddo ; enddo

991 endif

992

993 ! Register diagnostics for river fluxes

994 cs%no3_riv_flux = register_diag_field("ocean_model", "NO3_RIV_FLUX", &

995 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

996 day, "Dissolved Inorganic Nitrate Riverine Flux", "mmol/m^3 m/s")

997 cs%po4_riv_flux = register_diag_field("ocean_model", "PO4_RIV_FLUX", &

998 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

999 day, "Dissolved Inorganic Phosphate Riverine Flux", "mmol/m^3 m/s")

1000 cs%don_riv_flux = register_diag_field("ocean_model", "DON_RIV_FLUX", &

1001 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1002 day, "Dissolved Organic Nitrogen Riverine Flux", "mmol/m^3 m/s")

1003 cs%donr_riv_flux = register_diag_field("ocean_model", "DONR_RIV_FLUX", &

1004 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1005 day, "Refractory DON Riverine Flux", "mmol/m^3 m/s")

1006 cs%dop_riv_flux = register_diag_field("ocean_model", "DOP_RIV_FLUX", &

1007 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1008 day, "Dissolved Organic Phosphorus Riverine Flux", "mmol/m^3 m/s")

1009 cs%dopr_riv_flux = register_diag_field("ocean_model", "DOPR_RIV_FLUX", &

1010 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1011 day, "Refractory DOP Riverine Flux", "mmol/m^3 m/s")

1012 cs%sio3_riv_flux = register_diag_field("ocean_model", "SiO3_RIV_FLUX", &

1013 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1014 day, "Dissolved Inorganic Silicate Riverine Flux", "mmol/m^3 m/s")

1015 cs%fe_riv_flux = register_diag_field("ocean_model", "Fe_RIV_FLUX", &

1016 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1017 day, "Dissolved Inorganic Iron Riverine Flux", "mmol/m^3 m/s")

1018 cs%doc_riv_flux = register_diag_field("ocean_model", "DOC_RIV_FLUX", &

1019 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1020 day, "Dissolved Organic Carbon Riverine Flux", "mmol/m^3 m/s")

1021 cs%docr_riv_flux = register_diag_field("ocean_model", "DOCR_RIV_FLUX", &

1022 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1023 day, "Refractory DOC Riverine Flux", "mmol/m^3 m/s")

1024 cs%alk_riv_flux = register_diag_field("ocean_model", "ALK_RIV_FLUX", &

1025 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1026 day, "Alkalinity Riverine Flux", "meq/m^3 m/s")

1027 cs%alk_alt_co2_riv_flux = register_diag_field("ocean_model", "ALK_ALT_CO2_RIV_FLUX", &

1028 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1029 day, "Alkalinity Riverine Flux, Alternative CO2", "meq/m^3 m/s")

1030 cs%dic_riv_flux = register_diag_field("ocean_model", "DIC_RIV_FLUX", &

1031 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1032 day, "Dissolved Inorganic Carbon Riverine Flux", "mmol/m^3 m/s")

1033 cs%dic_alt_co2_riv_flux = register_diag_field("ocean_model", "DIC_ALT_CO2_RIV_FLUX", &

1034 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1035 day, "Dissolved Inorganic Carbon Riverine Flux, Alternative CO2", "mmol/m^3 m/s")

1036

1037 ! Register diagnostics for d14c forcing

1038 if (cs%abio_dic_on) then

1039 cs%d14c_id = register_diag_field("ocean_model", "D14C_FORCING", &

1040 diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid

1041 day, "Delta-14C in atmospheric CO2", "per mil, relative to Modern")

1042 endif

1043

1044 ! Register diagnostics for per-category forcing fields

1045 if (cs%ice_ncat > 0) then

1046 allocate(cs%fracr_cat_id(cs%ice_ncat+1))

1047 allocate(cs%qsw_cat_id(cs%ice_ncat+1))

1048 do m=1,cs%ice_ncat+1

1049 write(name, "(A,I0)") "FRACR_CAT_", m

1050 write(longname, "(A,I0)") "Fraction of area in ice category ", m

1051 units = "fraction"

1052 cs%fracr_cat_id(m) = register_diag_field("ocean_model", trim(name), &

1053 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid

1054 day, trim(longname), trim(units))

1055 write(name, "(A,I0)") "QSW_CAT_", m

1056 write(longname, "(A,I0)") "Shortwave penetrating through ice category ", m

1057 units = "W m-2"

1058 cs%qsw_cat_id(m) = register_diag_field("ocean_model", trim(name), &

1059 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid

1060 day, trim(longname), trim(units))

1061 enddo

1062 endif

1063

1064 if (cs%base_bio_on) then

1065 ! Read initial fesedflux and feventflux fields

1066 ! (1) get vertical dimension

1067 ! -- comes from fesedflux_file, assume same dimension in feventflux

1068 ! (maybe these fields should be combined?)

1069 ! -- note: read_Z_edges treats depth as positive UP => 0 at surface, negative at depth

1070 fesedflux_use_missing = .false.

1071 call read_z_edges(cs%fesedflux_file, "FESEDFLUXIN", cs%fesedflux_z_edges, cs%fesedflux_nz, &

1072 fesedflux_has_edges, fesedflux_use_missing, fesedflux_missing, scale=us%m_to_Z, &

1073 missing_scale=1.0)

1074

1075 ! (2) Allocate memory for fesedflux and feventflux

1076 allocate(cs%fesedflux_in(szi_(g), szj_(g), cs%fesedflux_nz))

1077 allocate(cs%feventflux_in(szi_(g), szj_(g), cs%fesedflux_nz))

1078 allocate(cs%fesedflux_dz(szi_(g), szj_(g), cs%fesedflux_nz))

1079

1080 ! (3) Read data

1081 ! TODO: Add US term to scale

1082 call mom_read_data(cs%fesedflux_file, "FESEDFLUXIN", cs%fesedflux_in(:,:,:), g%Domain, &

1083 scale=cs%fesedflux_scale_factor)

1084 call mom_read_data(cs%feventflux_file, "FESEDFLUXIN", cs%feventflux_in(:,:,:), g%Domain, &

1085 scale=cs%fesedflux_scale_factor)

1086

1087 ! (4) Relocate values that are below ocean bottom to layer that intersects bathymetry

1088 ! Remember, fesedflux_z_edges = 0 at surface and is < 0 below surface

1089

1090 do k=cs%fesedflux_nz, 1, -1

1091 kbot = k + 1 ! level k is between z(k) and z(k+1)

1092 do j=g%jsc, g%jec

1093 do i=g%isc, g%iec

1094 if (g%mask2dT(i,j) == 0) cycle

1095 if (g%bathyT(i,j) + cs%fesedflux_z_edges(1) < 1e-8 * us%m_to_Z) then

1096 write(log_message, *) "Current implementation of fesedflux assumes G%bathyT >=", &

1097 " first edge;first edge = ", -cs%fesedflux_z_edges(1), "bathyT = ", g%bathyT(i,j)

1098 call mom_error(fatal, log_message)

1099 endif

1100 ! Also figure out layer thickness while we're here

1101 cs%fesedflux_dz(i,j,k) = (cs%fesedflux_z_edges(k) - cs%fesedflux_z_edges(kbot)) * gv%Z_to_H

1102 ! If top interface is at or below ocean bottom, move flux in current layer up one

1103 ! and set thickness of current level to 0

1104 if (g%bathyT(i,j) + cs%fesedflux_z_edges(k) < 1e-8 * us%m_to_Z) then

1105 cs%fesedflux_in(i,j,k-1) = cs%fesedflux_in(i,j,k-1) + cs%fesedflux_in(i,j,k)

1106 cs%fesedflux_in(i,j,k) = 0.

1107 cs%feventflux_in(i,j,k-1) = cs%feventflux_in(i,j,k-1) + cs%feventflux_in(i,j,k)

1108 cs%feventflux_in(i,j,k) = 0.

1109 cs%fesedflux_dz(i,j,k) = 0.

1110 elseif (g%bathyT(i,j) + cs%fesedflux_z_edges(kbot) < 1e-8 * us%m_to_Z) then

1111 ! Otherwise, if lower interface is below bathymetry move interface to ocean bottom

1112 cs%fesedflux_dz(i,j,k) = (g%bathyT(i,j) + cs%fesedflux_z_edges(k)) * gv%Z_to_H

1113 endif

1114 enddo

1115 enddo

1116 enddo

1117

1118 ! Initialize external field for river fluxes

1119 if (cs%read_riv_fluxes) then

1120 cs%id_din_riv = init_external_field(cs%riv_flux_dataset%file_name, 'din_riv_flux', &

1121 domain=g%Domain%mpp_domain)

1122 cs%id_don_riv = init_external_field(cs%riv_flux_dataset%file_name, 'don_riv_flux', &

1123 domain=g%Domain%mpp_domain)

1124 cs%id_dip_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dip_riv_flux', &

1125 domain=g%Domain%mpp_domain)

1126 cs%id_dop_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dop_riv_flux', &

1127 domain=g%Domain%mpp_domain)

1128 cs%id_dsi_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dsi_riv_flux', &

1129 domain=g%Domain%mpp_domain)

1130 cs%id_dfe_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dfe_riv_flux', &

1131 domain=g%Domain%mpp_domain)

1132 cs%id_dic_riv = init_external_field(cs%riv_flux_dataset%file_name, 'dic_riv_flux', &

1133 domain=g%Domain%mpp_domain)

1134 cs%id_alk_riv = init_external_field(cs%riv_flux_dataset%file_name, 'alk_riv_flux', &

1135 domain=g%Domain%mpp_domain)

1136 cs%id_doc_riv = init_external_field(cs%riv_flux_dataset%file_name, 'doc_riv_flux', &

1137 domain=g%Domain%mpp_domain)

1138 endif

1139 endif

1140

1141 if (cs%abio_dic_on) then

1142 ! Initialize external field for d14c forcing

1143 do m=1,3

1144 cs%id_d14c(m) = init_external_field(cs%d14c_dataset(m)%file_name, "Delta14co2_in_air", &

1145 ignore_axis_atts=.true.)

1146 enddo

1147 endif

1148

1149 ! Initialize external field for restoring

1150 if (cs%restoring_I_tau_source == "file") then

1151 select case(cs%restoring_source)

1152 case("file")

1153 ! Set up array for reading in raw restoring data

1154 allocate(cs%restoring_in(szi_(g), szj_(g), cs%restoring_nz, cs%restore_count), source=0.)

1155 do m=1,cs%restore_count

1156 cs%id_tracer_restoring(m) = init_external_field(cs%restoring_file, &

1157 trim(cs%tracer_restoring_varname(m)), domain=g%Domain%mpp_domain)

1158 enddo

1159 end select

1160 select case(cs%restoring_I_tau_source)

1161 case("file")

1162 allocate(cs%I_tau(szi_(g), szj_(g), cs%restoring_timescale_nz), source=0.)

1163 call mom_read_data(cs%restoring_I_tau_file, "RTAU", cs%I_tau(:,:,:), g%Domain)

1164 end select

1165 endif

1166

1167end subroutine initialize_marbl_tracers

1168

1169!> This subroutine is used to register tracer fields and subroutines

1170!! to be used with MOM.

1171subroutine register_marbl_diags(MARBL_diags, diag, day, G, id_diags)

1172

1173 type(marbl_diagnostics_type), intent(in) :: MARBL_diags !< MARBL diagnostics from MARBL_instances

1174 type(time_type), target, intent(in) :: day !< Time of the start of the run.

1175 type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output.

1176 !integer, allocatable, intent(inout) :: id_diags(:) !< allocatable array storing diagnostic index number

1177 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure

1178 type(temp_marbl_diag), allocatable, intent(inout) :: id_diags(:) !< allocatable array storing diagnostic index

1179 !! number and buffer space for collecting diags

1180 !! from all columns

1181

1182 integer :: m, diag_size

1183

1184 diag_size = size(marbl_diags%diags)

1185 allocate(id_diags(diag_size))

1186 do m = 1, diag_size

1187 id_diags(m)%id = -1

1188 if (trim(marbl_diags%diags(m)%vertical_grid) .eq. "none") then ! 2D field

1189 id_diags(m)%id = register_diag_field("ocean_model", &

1190 trim(marbl_diags%diags(m)%short_name), &

1191 diag%axesT1, & ! T => tracer grid? 1 => no vertical grid

1192 day, &

1193 trim(marbl_diags%diags(m)%long_name), &

1194 trim(marbl_diags%diags(m)%units))

1195 if (id_diags(m)%id > 0) allocate(id_diags(m)%field_2d(szi_(g),szj_(g)), source=0.0)

1196 else ! 3D field

1197 ! TODO: MARBL should provide v_extensive through MARBL_diags

1198 ! (for now, FESEDFLUX is the only one that should be true)

1199 ! Also, known issue where passing v_extensive=.false. isn't

1200 ! treated the same as not passing v_extensive

1201 if (trim(marbl_diags%diags(m)%short_name).eq."FESEDFLUX") then

1202 id_diags(m)%id = register_diag_field("ocean_model", &

1203 trim(marbl_diags%diags(m)%short_name), &

1204 diag%axesTL, & ! T=> tracer grid? L => layer center

1205 day, &

1206 trim(marbl_diags%diags(m)%long_name), &

1207 trim(marbl_diags%diags(m)%units), &

1208 v_extensive=.true.)

1209 else

1210 id_diags(m)%id = register_diag_field("ocean_model", &

1211 trim(marbl_diags%diags(m)%short_name), &

1212 diag%axesTL, & ! T=> tracer grid? L => layer center

1213 day, &

1214 trim(marbl_diags%diags(m)%long_name), &

1215 trim(marbl_diags%diags(m)%units))

1216 endif

1217 if (id_diags(m)%id > 0) allocate(id_diags(m)%field_3d(szi_(g),szj_(g), szk_(g)), source=0.0)

1218 endif

1219 enddo

1220

1221end subroutine register_marbl_diags

1222

1223!> This subroutine allocates memory for saved state fields and registers them in the restart files

1224subroutine setup_saved_state(MARBL_saved_state, HI, GV, restart_CS, tracers_may_reinit, &

1225 local_saved_state)

1226

1227 type(marbl_saved_state_type), intent(in) :: MARBL_saved_state !< MARBL saved state from

1228 !! MARBL_instances

1229 type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.

1230 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure

1231 type(mom_restart_cs), pointer, intent(in) :: restart_CS !< control structure to add saved state

1232 !! to restarts

1233 logical, intent(in) :: tracers_may_reinit !< used to determine mandatory

1234 !! flag in restart

1235 type(saved_state_for_marbl_type), allocatable, intent(inout) :: local_saved_state(:) !< allocatable array for local

1236 !! saved state

1237

1238 integer :: num_fields, m

1239 character(len=200) :: log_message, varname

1240

1241 num_fields = marbl_saved_state%saved_state_cnt

1242 allocate(local_saved_state(num_fields))

1243

1244 do m=1,num_fields

1245 write(varname, "(2A)") "MARBL_", trim(marbl_saved_state%state(m)%short_name)

1246 select case (marbl_saved_state%state(m)%rank)

1247 case (2)

1248 allocate(local_saved_state(m)%field_2d(szi_(hi),szj_(hi)), source=0.0)

1249 call register_restart_field(local_saved_state(m)%field_2d, varname, &

1250 .not.tracers_may_reinit, restart_cs)

1251 case (3)

1252 if (trim(marbl_saved_state%state(m)%vertical_grid).eq."layer_avg") then

1253 allocate(local_saved_state(m)%field_3d(szi_(hi),szj_(hi), szk_(gv)), source=0.0)

1254 call register_restart_field(local_saved_state(m)%field_3d, varname, &

1255 .not.tracers_may_reinit, restart_cs)

1256 else

1257 write(log_message, "(3A, I0, A)") "'", trim(marbl_saved_state%state(m)%vertical_grid), &

1258 "' is an invalid vertical grid for saved state (ind = ", m, ")"

1259 call mom_error(fatal, log_message)

1260 endif

1261 case DEFAULT

1262 write(log_message, "(I0, A, I0, A)") marbl_saved_state%state(m)%rank, &

1263 " is an invalid rank for saved state (ind = ", m, ")"

1264 call mom_error(fatal, log_message)

1265 end select

1266 local_saved_state(m)%short_name = trim(marbl_saved_state%state(m)%short_name)

1267 write(local_saved_state(m)%file_varname, "(2A)") "MARBL_", trim(local_saved_state(m)%short_name)

1268 local_saved_state(m)%units = trim(marbl_saved_state%state(m)%units)

1269 enddo

1270

1271end subroutine setup_saved_state

1272

1273!> This subroutine applies diapycnal diffusion and any other column

1274!! tracer physics or chemistry to the tracers from this file.

1275subroutine marbl_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, tv, &

1276 KPP_CSp, nonLocalTrans, evap_CFL_limit, minimum_forcing_depth)

1277

1278 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure

1279 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure

1280 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &

1281 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].

1282 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &

1283 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].

1284 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &

1285 intent(in) :: ea !< an array to which the amount of fluid entrained

1286 !! from the layer above during this call will be

1287 !! added [H ~> m or kg m-2].

1288 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &

1289 intent(in) :: eb !< an array to which the amount of fluid entrained

1290 !! from the layer below during this call will be

1291 !! added [H ~> m or kg m-2].

1292 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic

1293 !! and tracer forcing fields. Unused fields have NULL ptrs.

1294 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]

1295 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

1296 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a previous

1297 !! call to register_MARBL_tracers.

1298 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables

1299 type(kpp_cs), optional, pointer :: kpp_csp !< KPP control structure

1300 real, optional, intent(in) :: nonlocaltrans(:,:,:) !< Non-local transport [1]

1301 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can

1302 !! be fluxed out of the top layer in a timestep [1]

1303 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which

1304 !! fluxes can be applied [m]

1305

1306 ! Local variables

1307 character(len=256) :: log_message

1308 real, dimension(SZI_(G),SZJ_(G)) :: net_salt_rate ! Surface salt flux into the ocean

1309 ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1].

1310 real, dimension(SZI_(G),SZJ_(G)) :: flux_from_salt_flux ! Surface tracer flux from salt flux

1311 ! [conc Z T-1 ~> conc m s-1].

1312 real, dimension(SZI_(G),SZJ_(G)) :: ref_mask ! Mask for 2D MARBL diags using ref_depth [1]

1313 real, dimension(SZI_(G),SZJ_(G)) :: riv_flux_loc ! Local copy of CS%RIV_FLUXES*dt [conc H ~> mmol m-2]

1314 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]

1315 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: bot_flux_to_tend ! Conversion factor for bottom tlux -> tend

1316 ! [Z-1 ~> m-1]

1317 real :: cum_bftt_dz ! sum of bot_flux_to_tend * dz from the bottom layer to current layer [1]

1318 real, dimension(0:GV%ke) :: zi ! z-coordinate interface depth [Z ~> m]

1319 real, dimension(GV%ke) :: zc ! z-coordinate layer center depth [Z ~> m]

1320 real, dimension(GV%ke) :: dz ! z-coordinate cell thickness [H ~> m]

1321 integer :: i, j, k, is, ie, js, je, nz, m

1322

1323 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

1324

1325 if (.not.associated(cs)) return

1326

1327 ! (1) Compute surface fluxes

1328 ! FIXME: MARBL can handle computing surface fluxes for all columns simultaneously

1329 ! I was just thinking going column-by-column at first might be easier

1330 do j=js,je

1331 do i=is,ie

1332 ! i. only want ocean points in this loop

1333 if (g%mask2dT(i,j) == 0) cycle

1334

1335 ! ii. Load proper column data

1336 ! * surface flux forcings

1337 ! These fields are getting the correct data

1338 ! TODO: if top layer is vanishly thin, do we actually want (e.g.) top 5m average temp / salinity?

1339 ! How does MOM pass SST and SSS to GFDL coupler? (look in core.F90?)

1340 if (cs%sss_ind > 0) &

1341 marbl_instances%surface_flux_forcings(cs%sss_ind)%field_0d(1) = tv%S(i,j,1) * us%S_to_ppt

1342 if (cs%sst_ind > 0) &

1343 marbl_instances%surface_flux_forcings(cs%sst_ind)%field_0d(1) = tv%T(i,j,1) * us%C_to_degC

1344 if (cs%ifrac_ind > 0) &

1345 marbl_instances%surface_flux_forcings(cs%ifrac_ind)%field_0d(1) = fluxes%ice_fraction(i,j)

1346

1347 ! MARBL wants u10_sqr in (m/s)^2

1348 if (cs%u10_sqr_ind > 0) &

1349 marbl_instances%surface_flux_forcings(cs%u10_sqr_ind)%field_0d(1) = fluxes%u10_sqr(i,j) * &

1350 ((us%L_T_to_m_s)**2)

1351

1352 ! mct_driver/ocn_cap_methods:93 -- ice_ocean_boundary%p(i,j) comes from coupler

1353 ! We may need a new ice_ocean_boundary%p_atm because %p includes ice in GFDL driver

1354 if (cs%atmpress_ind > 0) then

1355 if (associated(fluxes%p_surf_full)) then

1356 marbl_instances%surface_flux_forcings(cs%atmpress_ind)%field_0d(1) = &

1357 fluxes%p_surf_full(i,j) * ((us%R_to_kg_m3 * (us%L_T_to_m_s**2)) * atm_per_pa)

1358 else

1359 ! hardcode value of 1 atm (can't figure out how to get this from solo_driver)

1360 marbl_instances%surface_flux_forcings(cs%atmpress_ind)%field_0d(1) = 1.

1361 endif

1362 endif

1363

1364 ! These are okay, but need option to come in from coupler

1365 if (cs%xco2_ind > 0) &

1366 marbl_instances%surface_flux_forcings(cs%xco2_ind)%field_0d(1) = fluxes%atm_co2(i,j)

1367 if (cs%xco2_alt_ind > 0) &

1368 marbl_instances%surface_flux_forcings(cs%xco2_alt_ind)%field_0d(1) = fluxes%atm_alt_co2(i,j)

1369

1370 ! These are okay, but need option to read in from file

1371 if (cs%dust_dep_ind > 0) &

1372 marbl_instances%surface_flux_forcings(cs%dust_dep_ind)%field_0d(1) = &

1373 fluxes%dust_flux(i,j) * us%RZ_T_to_kg_m2s

1374

1375 if (cs%fe_dep_ind > 0) &

1376 marbl_instances%surface_flux_forcings(cs%fe_dep_ind)%field_0d(1) = &

1377 fluxes%iron_flux(i,j) * (us%Z_to_m * us%s_to_T)

1378

1379 ! MARBL wants ndep in (mmol/m^2/s)

1380 if (cs%nox_flux_ind > 0) &

1381 marbl_instances%surface_flux_forcings(cs%nox_flux_ind)%field_0d(1) = fluxes%noy_dep(i,j) * &

1382 (us%Z_to_m * us%s_to_T)

1383 if (cs%nhy_flux_ind > 0) &

1384 marbl_instances%surface_flux_forcings(cs%nhy_flux_ind)%field_0d(1) = fluxes%nhx_dep(i,j) * &

1385 (us%Z_to_m * us%s_to_T)

1386

1387 if (cs%d14c_ind > 0) &

1388 marbl_instances%surface_flux_forcings(cs%d14c_ind)%field_0d(1) = cs%d14c(i,j)

1389

1390 ! * tracers at surface

1391 ! TODO: average over some shallow depth (e.g. 5m)

1392 do m=1,cs%ntr

1393 marbl_instances%tracers_at_surface(1,m) = cs%tracer_data(m)%tr(i,j,1)

1394 enddo

1395

1396 ! * surface flux saved state

1397 do m=1,size(marbl_instances%surface_flux_saved_state%state)

1398 ! (currently only 2D fields are saved from surface_flux_compute())

1399 marbl_instances%surface_flux_saved_state%state(m)%field_2d(1) = &

1400 cs%surface_flux_saved_state(m)%field_2d(i,j)

1401 enddo

1402

1403 ! iii. Compute surface fluxes in MARBL

1404 call marbl_instances%surface_flux_compute()

1405 if (marbl_instances%StatusLog%labort_marbl) then

1406 call marbl_instances%StatusLog%log_error_trace("MARBL_instances%surface_flux_compute()", &

1407 "MARBL_tracers_column_physics")

1408 endif

1409 call print_marbl_log(marbl_instances%StatusLog)

1410 call marbl_instances%StatusLog%erase()

1411

1412 ! iv. Copy output that MOM6 needs to hold on to

1413 ! * saved state

1414 do m=1,size(marbl_instances%surface_flux_saved_state%state)

1415 cs%surface_flux_saved_state(m)%field_2d(i,j) = &

1416 marbl_instances%surface_flux_saved_state%state(m)%field_2d(1)

1417 enddo

1418

1419 ! * diagnostics

1420 do m=1,size(marbl_instances%surface_flux_diags%diags)

1421 ! All diags are 2D coming from surface

1422 if (cs%surface_flux_diags(m)%id > 0) &

1423 cs%surface_flux_diags(m)%field_2d(i,j) = &

1424 real(marbl_instances%surface_flux_diags%diags(m)%field_2d(1))

1425 enddo

1426

1427 ! * Surface tracer flux

1428 cs%STF(i,j,:) = marbl_instances%surface_fluxes(1,:) * (us%m_to_Z * us%T_to_s)

1429

1430 ! * Surface flux output

1431 do m=1,cs%sfo_cnt

1432 cs%SFO(i,j,m) = marbl_instances%surface_flux_output%outputs_for_GCM(m)%forcing_field_0d(1)

1433 enddo

1434

1435 enddo

1436 enddo

1437

1438 if (associated(fluxes%salt_flux)) then

1439 ! convert salt flux to tracer fluxes and add to STF

1440 do j=js,je ; do i=is,ie

1441 net_salt_rate(i,j) = (1000.0*us%ppt_to_S * fluxes%salt_flux(i,j)) * gv%RZ_to_H

1442 enddo ; enddo

1443

1444 ! DIC related tracers

1445 do j=js,je ; do i=is,ie

1446 flux_from_salt_flux(i,j) = (cs%DIC_salt_ratio * gv%H_to_Z) * net_salt_rate(i,j)

1447 enddo ; enddo

1448 m = cs%tracer_inds%dic_ind

1449 if (m > 0) then

1450 do j=js,je ; do i=is,ie

1451 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)

1452 enddo ; enddo

1453 if (cs%id_surface_flux_from_salt_flux(m) > 0) &

1454 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)

1455 endif

1456 m = cs%tracer_inds%dic_alt_co2_ind

1457 if (m > 0) then

1458 do j=js,je ; do i=is,ie

1459 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)

1460 enddo ; enddo

1461 if (cs%id_surface_flux_from_salt_flux(m) > 0) &

1462 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)

1463 endif

1464 m = cs%tracer_inds%abio_dic_ind

1465 if (m > 0) then

1466 do j=js,je ; do i=is,ie

1467 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)

1468 enddo ; enddo

1469 if (cs%id_surface_flux_from_salt_flux(m) > 0) &

1470 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)

1471 endif

1472 m = cs%tracer_inds%abio_di14c_ind

1473 if (m > 0) then

1474 do j=js,je ; do i=is,ie

1475 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)

1476 enddo ; enddo

1477 if (cs%id_surface_flux_from_salt_flux(m) > 0) &

1478 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)

1479 endif

1480

1481 ! ALK related tracers

1482 do j=js,je ; do i=is,ie

1483 flux_from_salt_flux(i,j) = (cs%ALK_salt_ratio * gv%H_to_Z) * net_salt_rate(i,j)

1484 enddo ; enddo

1485 m = cs%tracer_inds%alk_ind

1486 if (m > 0) then

1487 do j=js,je ; do i=is,ie

1488 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)

1489 enddo ; enddo

1490 if (cs%id_surface_flux_from_salt_flux(m) > 0) &

1491 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)

1492 endif

1493 m = cs%tracer_inds%alk_alt_co2_ind

1494 if (m > 0) then

1495 do j=js,je ; do i=is,ie

1496 cs%STF(i,j,m) = cs%STF(i,j,m) + flux_from_salt_flux(i,j)

1497 enddo ; enddo

1498 if (cs%id_surface_flux_from_salt_flux(m) > 0) &

1499 call post_data(cs%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, cs%diag)

1500 endif

1501 endif

1502

1503 if (cs%debug) then

1504 do m=1,cs%ntr

1505 call hchksum(cs%STF(:,:,m), &

1506 trim(marbl_instances%tracer_metadata(m)%short_name)//" sfc_flux", g%HI, &

1507 unscale=us%Z_to_m*us%s_to_T)

1508 enddo

1509 endif

1510

1511 ! (2) Post surface fluxes and their diagnostics (currently all 2D)

1512 do m=1,cs%ntr

1513 if (cs%id_surface_flux_out(m) > 0) &

1514 call post_data(cs%id_surface_flux_out(m), cs%STF(:,:,m), cs%diag)

1515 enddo

1516 do m=1,size(cs%surface_flux_diags)

1517 if (cs%surface_flux_diags(m)%id > 0) &

1518 call post_data(cs%surface_flux_diags(m)%id, cs%surface_flux_diags(m)%field_2d(:,:), cs%diag)

1519 enddo

1520

1521 ! (3) Apply surface fluxes via vertical diffusion

1522 ! Compute KPP nonlocal term if necessary

1523 if (present(kpp_csp)) then

1524 if (associated(kpp_csp) .and. present(nonlocaltrans)) then

1525 do m=1,cs%ntr

1526 call kpp_nonlocaltransport(kpp_csp, g, gv, h_old, nonlocaltrans, cs%STF(:,:,m), dt, &

1527 cs%diag, cs%tracer_data(m)%tr_ptr, cs%tracer_data(m)%tr(:,:,:), &

1528 flux_scale=gv%Z_to_H)

1529 enddo

1530 endif

1531 if (cs%debug) then

1532 do m=1,cs%ntr

1533 call hchksum(cs%tracer_data(m)%tr(:,:,m), &

1534 trim(marbl_instances%tracer_metadata(m)%short_name)//' post KPP', g%HI)

1535 enddo

1536 endif

1537 endif

1538

1539 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then

1540 do m=1,cs%ntr

1541 do k=1,nz ;do j=js,je ; do i=is,ie

1542 h_work(i,j,k) = h_old(i,j,k)

1543 enddo ; enddo ; enddo

1544 ! CS%RIV_FLUXES is conc m/s, in_flux_optional expects time-integrated flux (conc H)

1545 do j=js,je ; do i=is,ie

1546 riv_flux_loc(i,j) = (cs%RIV_FLUXES(i,j,m) * (dt*us%T_to_s)) * gv%m_to_H

1547 enddo ; enddo

1548 if (cs%debug) &

1549 call hchksum(riv_flux_loc(:,:), &

1550 trim(marbl_instances%tracer_metadata(m)%short_name)//' riv flux', g%HI, unscale=gv%H_to_m)

1551 call applytracerboundaryfluxesinout(g, gv, cs%tracer_data(m)%tr(:,:,:) , dt, fluxes, h_work, &

1552 evap_cfl_limit, minimum_forcing_depth, in_flux_optional=riv_flux_loc)

1553 call tracer_vertdiff(h_work, ea, eb, dt, cs%tracer_data(m)%tr(:,:,:), g, gv, &

1554 sfc_flux=gv%Rho0 * cs%STF(:,:,m))

1555 enddo

1556 else

1557 do m=1,cs%ntr

1558 call tracer_vertdiff(h_old, ea, eb, dt, cs%tracer_data(m)%tr(:,:,:), g, gv, &

1559 sfc_flux=gv%Rho0 * cs%STF(:,:,m))

1560 enddo

1561 endif

1562

1563 if (cs%debug) then

1564 do m=1,cs%ntr

1565 call hchksum(cs%tracer_data(m)%tr(:,:,m), &

1566 trim(marbl_instances%tracer_metadata(m)%short_name)//' post tracer_vertdiff', g%HI)

1567 enddo

1568 endif

1569

1570 ! (4) Compute interior tendencies

1571

1572 bot_flux_to_tend(:, :, :) = 0.

1573 do j=js,je

1574 do i=is,ie

1575 ! i. only want ocean points in this loop

1576 if (g%mask2dT(i,j) == 0) cycle

1577

1578 ! ii. Set up vertical domain and bot_flux_to_tend

1579 ! Calculate depth of interface by building up thicknesses from the bottom (top interface is always 0)

1580 ! MARBL wants this to be positive-down

1581 zi(gv%ke) = g%bathyT(i,j)

1582 marbl_instances%bot_flux_to_tend(:) = 0.

1583 cum_bftt_dz = 0.

1584 do k = gv%ke, 1, -1

1585 ! TODO: if we move this above vertical mixing, use h_old

1586 dz(k) = h_new(i,j,k) ! cell thickness

1587 zc(k) = zi(k) - 0.5 * (dz(k)*gv%H_to_Z)

1588 zi(k-1) = zi(k) - (dz(k)*gv%H_to_Z)

1589 if (g%bathyT(i,j) - zi(k-1) <= cs%bot_flux_mix_thickness) then

1590 marbl_instances%bot_flux_to_tend(k) = us%m_to_Z * cs%Ibfmt

1591 cum_bftt_dz = cum_bftt_dz + marbl_instances%bot_flux_to_tend(k) * (gv%H_to_m * dz(k))

1592 elseif (g%bathyT(i,j) - zi(k) < cs%bot_flux_mix_thickness) then

1593 ! MARBL_instances%bot_flux_to_tend(k) = (1. - (G%bathyT(i,j) - zi(k)) * CS%Ibfmt) / dz(k)

1594 marbl_instances%bot_flux_to_tend(k) = (1. - cum_bftt_dz) / (gv%H_to_m * dz(k))

1595 endif

1596 enddo

1597 if (g%bathyT(i,j) - zi(0) < cs%bot_flux_mix_thickness) &

1598 marbl_instances%bot_flux_to_tend(:) = marbl_instances%bot_flux_to_tend(:) * &

1599 cs%bot_flux_mix_thickness / (g%bathyT(i,j) - zi(0))

1600 if (cs%bot_flux_to_tend_id > 0) &

1601 bot_flux_to_tend(i, j, :) = marbl_instances%bot_flux_to_tend(:)

1602

1603 ! zw(1:nz) is bottom cell depth so no element of zw = 0, it is assumed to be top layer depth

1604 marbl_instances%domain%zw(:) = us%Z_to_m * zi(1:gv%ke)

1605 marbl_instances%domain%zt(:) = us%Z_to_m * zc(:)

1606 marbl_instances%domain%delta_z(:) = gv%H_to_m * dz(:)

1607

1608 ! iii. Load proper column data

1609 ! * Forcing Fields

1610 ! These fields are getting the correct data

1611 if (cs%potemp_ind > 0) &

1612 marbl_instances%interior_tendency_forcings(cs%potemp_ind)%field_1d(1,:) = tv%T(i,j,:) * us%C_to_degC

1613 if (cs%salinity_ind > 0) &

1614 marbl_instances%interior_tendency_forcings(cs%salinity_ind)%field_1d(1,:) = tv%S(i,j,:) * us%S_to_ppt

1615

1616 ! This are okay, but need option to read in from file

1617 ! (Same as dust_dep_ind for surface_flux_forcings)

1618 if (cs%dustflux_ind > 0) &

1619 marbl_instances%interior_tendency_forcings(cs%dustflux_ind)%field_0d(1) = &

1620 fluxes%dust_flux(i,j) * us%RZ_T_to_kg_m2s

1621

1622 ! TODO: Support PAR (currently just using single subcolumn)

1623 ! (Look for Pen_sw_bnd?)

1624 if (cs%PAR_col_frac_ind > 0) then

1625 ! second index is num_subcols, not depth

1626 !MARBL_instances%interior_tendency_forcings(CS%PAR_col_frac_ind)%field_1d(1,:) = fluxes%fracr_cat(i,j,:)

1627 if (cs%use_ice_category_fields) then

1628 marbl_instances%interior_tendency_forcings(cs%PAR_col_frac_ind)%field_1d(1,:) = &

1629 fluxes%fracr_cat(i,j,:)

1630 else

1631 marbl_instances%interior_tendency_forcings(cs%PAR_col_frac_ind)%field_1d(1,1) = 1.

1632 endif

1633 endif

1634

1635 if (cs%surf_shortwave_ind > 0) then

1636 ! second index is num_subcols, not depth

1637 if (cs%use_ice_category_fields) then

1638 marbl_instances%interior_tendency_forcings(cs%surf_shortwave_ind)%field_1d(1,:) = &

1639 fluxes%qsw_cat(i,j,:) * us%QRZ_T_to_W_m2

1640 else

1641 marbl_instances%interior_tendency_forcings(cs%surf_shortwave_ind)%field_1d(1,1) = &

1642 fluxes%sw(i,j) * us%QRZ_T_to_W_m2

1643 endif

1644 endif

1645 ! Tracer restoring

1646 do m=1,cs%restore_count

1647 marbl_instances%interior_tendency_forcings(cs%tracer_restoring_ind(m))%field_1d(1,:) = 0.

1648 call remapping_core_h(cs%restoring_remapCS, cs%restoring_nz, cs%restoring_dz(:), &

1649 cs%restoring_in(i,j,:,m), gv%ke, dz(:), &

1650 marbl_instances%interior_tendency_forcings(cs%tracer_restoring_ind(m))%field_1d(1,:))

1651 if (m==1) then

1652 call remapping_core_h(cs%restoring_remapCS, cs%restoring_timescale_nz, &

1653 cs%restoring_timescale_dz(:), cs%I_tau(i,j,:), gv%ke, dz(:), &

1654 marbl_instances%interior_tendency_forcings(cs%tracer_I_tau_ind(m))%field_1d(1,:))

1655 else

1656 marbl_instances%interior_tendency_forcings(cs%tracer_I_tau_ind(m))%field_1d(1,:) = &

1657 marbl_instances%interior_tendency_forcings(cs%tracer_I_tau_ind(1))%field_1d(1,:)

1658 endif

1659 enddo

1660

1661 ! TODO: In POP, pressure comes from a function in state_mod.F90; I don't see a similar function here

1662 ! This formulation is from Levitus 1994, and I think it belongs in MOM_EOS.F90?

1663 ! Converts depth [m] -> pressure [bars]

1664 ! NOTE: Andrew recommends using GV%H_to_Pa

1665 if (cs%pressure_ind > 0) &

1666 marbl_instances%interior_tendency_forcings(cs%pressure_ind)%field_1d(1,:) = &

1667 (0.0598088 * (exp(-0.025*us%Z_to_m * zc(:)) - 1.)) + &

1668 (0.100766 * us%Z_to_m * zc(:)) + (2.28405e-7*((us%Z_to_m * zc(:))**2))

1669

1670 if (cs%fesedflux_ind > 0) then

1671 marbl_instances%interior_tendency_forcings(cs%fesedflux_ind)%field_1d(1,:) = 0.

1672 call reintegrate_column(cs%fesedflux_nz, &

1673 cs%fesedflux_dz(i,j,:) * (sum(dz(:) * gv%H_to_Z) / g%bathyT(i,j)), &

1674 cs%fesedflux_in(i,j,:) + cs%feventflux_in(i,j,:), gv%ke, dz(:), &

1675 marbl_instances%interior_tendency_forcings(cs%fesedflux_ind)%field_1d(1,:))

1676 endif

1677

1678 ! TODO: add ability to read these fields from file

1679 ! also, add constant values to CS

1680 if (cs%o2_scalef_ind > 0) &

1681 marbl_instances%interior_tendency_forcings(cs%o2_scalef_ind)%field_1d(1,:) = 1.

1682 if (cs%remin_scalef_ind > 0) &

1683 marbl_instances%interior_tendency_forcings(cs%remin_scalef_ind)%field_1d(1,:) = 1.

1684

1685 ! * Column Tracers

1686 do m=1,cs%ntr

1687 marbl_instances%tracers(m, :) = cs%tracer_data(m)%tr(i,j,:)

1688 enddo

1689

1690 ! * interior tendency saved state

1691 ! (currently only 3D fields are saved from interior_tendency_compute())

1692 do m=1,size(marbl_instances%interior_tendency_saved_state%state)

1693 marbl_instances%interior_tendency_saved_state%state(m)%field_3d(:,1) = &

1694 cs%interior_tendency_saved_state(m)%field_3d(i,j,:)

1695 enddo

1696

1697 ! iv. Compute interior tendencies in MARBL

1698 call marbl_instances%interior_tendency_compute()

1699 if (marbl_instances%StatusLog%labort_marbl) then

1700 call marbl_instances%StatusLog%log_error_trace(&

1701 "MARBL_instances%interior_tendency_compute()", "MARBL_tracers_column_physics")

1702 endif

1703 call print_marbl_log(marbl_instances%StatusLog, g, i, j)

1704 call marbl_instances%StatusLog%erase()

1705

1706 ! v. Apply tendencies immediately

1707 ! First pass - Euler step; if stability issues, we can do something different (subcycle?)

1708 do m=1,cs%ntr

1709 cs%tracer_data(m)%tr(i,j,:) = cs%tracer_data(m)%tr(i,j,:) + (dt * us%T_to_s) * &

1710 marbl_instances%interior_tendencies(m,:)

1711 enddo

1712

1713 ! vi. Copy output that MOM6 needs to hold on to

1714 ! * saved state

1715 do m=1,size(marbl_instances%interior_tendency_saved_state%state)

1716 cs%interior_tendency_saved_state(m)%field_3d(i,j,:) = &

1717 marbl_instances%interior_tendency_saved_state%state(m)%field_3d(:,1)

1718 enddo

1719

1720 ! * diagnostics

1721 do m=1,size(marbl_instances%interior_tendency_diags%diags)

1722 if (cs%interior_tendency_diags(m)%id > 0) then

1723 if (allocated(cs%interior_tendency_diags(m)%field_2d)) then

1724 ! Only copy values if ref_depth < bathyT

1725 if (g%bathyT(i,j) > real(marbl_instances%interior_tendency_diags%diags(m)%ref_depth)) then

1726 cs%interior_tendency_diags(m)%field_2d(i,j) = &

1727 real(marbl_instances%interior_tendency_diags%diags(m)%field_2d(1))

1728 endif

1729 else ! not a 2D diagnostic

1730 cs%interior_tendency_diags(m)%field_3d(i,j,:) = &

1731 real(marbl_instances%interior_tendency_diags%diags(m)%field_3d(:,1))

1732 endif

1733 endif

1734 enddo

1735

1736 ! * tendency values themselves (and vertical integrals of them)

1737 do m=1,cs%ntr

1738 if (allocated(cs%interior_tendency_out(m)%field_3d)) &

1739 cs%interior_tendency_out(m)%field_3d(i,j,:) = marbl_instances%interior_tendencies(m,:)

1740

1741 if (allocated(cs%interior_tendency_out_zint(m)%field_2d)) &

1742 cs%interior_tendency_out_zint(m)%field_2d(i,j) = (sum(dz(:) * &

1743 marbl_instances%interior_tendencies(m,:)))

1744

1745 if (allocated(cs%interior_tendency_out_zint_100m(m)%field_2d)) then

1746 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) = 0.

1747 do k=1,gv%ke

1748 if (zi(k) < us%m_to_Z * 100.) then

1749 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) = &

1750 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) + gv%H_to_m * dz(k) * &

1751 marbl_instances%interior_tendencies(m,k)

1752 elseif (zi(k-1) < us%m_to_Z * 100.) then

1753 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) = &

1754 cs%interior_tendency_out_zint_100m(m)%field_2d(i,j) + gv%H_to_m * dz(k) * &

1755 ((us%m_to_Z * 100. - zi(k-1)) / (zi(k) - zi(k-1))) * &

1756 marbl_instances%interior_tendencies(m,k)

1757 else

1758 exit

1759 endif

1760 enddo

1761 endif

1762 enddo

1763

1764 ! * Interior tendency output

1765 do m=1,cs%ito_cnt

1766 cs%ITO(i,j,:,m) = &

1767 marbl_instances%interior_tendency_output%outputs_for_GCM(m)%forcing_field_1d(1,:)

1768 enddo

1769

1770 enddo

1771 enddo

1772

1773 if (cs%debug) then

1774 do m=1,cs%ntr

1775 call hchksum(cs%tracer_data(m)%tr(:,:,m), &

1776 trim(marbl_instances%tracer_metadata(m)%short_name)//' post source-sink', g%HI)

1777 enddo

1778 endif

1779

1780 ! (5) Post diagnostics from our buffer

1781 ! i. Interior tendency diagnostics (mix of 2D and 3D)

1782 ! ii. Interior tendencies themselves

1783 ! iii. Forcing fields

1784 if (cs%bot_flux_to_tend_id > 0) &

1785 call post_data(cs%bot_flux_to_tend_id, bot_flux_to_tend(:, :, :), cs%diag)

1786

1787 do m=1,size(cs%interior_tendency_diags)

1788 if (cs%interior_tendency_diags(m)%id > 0) then

1789 if (allocated(cs%interior_tendency_diags(m)%field_2d)) then

1790 if (real(marbl_instances%interior_tendency_diags%diags(m)%ref_depth) == 0.) then

1791 call post_data(cs%interior_tendency_diags(m)%id, &

1792 cs%interior_tendency_diags(m)%field_2d(:,:), cs%diag)

1793 else ! non-zero ref-depth

1794 ref_mask(:, :) = 0.

1795 do j=js,je ; do i=is,ie

1796 if (g%bathyT(i,j) > real(marbl_instances%interior_tendency_diags%diags(m)%ref_depth)) &

1797 ref_mask(i,j) = 1.

1798 enddo ; enddo

1799 call post_data(cs%interior_tendency_diags(m)%id, &

1800 cs%interior_tendency_diags(m)%field_2d(:,:), cs%diag, mask=ref_mask(:,:))

1801 endif

1802 elseif (allocated(cs%interior_tendency_diags(m)%field_3d)) then

1803 call post_data(cs%interior_tendency_diags(m)%id, &

1804 cs%interior_tendency_diags(m)%field_3d(:,:,:), cs%diag)

1805 else

1806 write(log_message, "(A, I0, A, I0, A)") "Diagnostic number ", m, " post id ", &

1807 cs%interior_tendency_diags(m)%id," did not allocate 2D or 3D array"

1808 call mom_error(fatal, log_message)

1809 endif

1810 endif

1811 enddo

1812

1813 do m=1,cs%ntr

1814 if (allocated(cs%interior_tendency_out(m)%field_3d)) &

1815 call post_data(cs%interior_tendency_out(m)%id, &

1816 cs%interior_tendency_out(m)%field_3d(:,:,:), cs%diag)

1817 if (allocated(cs%interior_tendency_out_zint(m)%field_2d)) &

1818 call post_data(cs%interior_tendency_out_zint(m)%id, &

1819 cs%interior_tendency_out_zint(m)%field_2d(:,:), cs%diag)

1820 if (allocated(cs%interior_tendency_out_zint_100m(m)%field_2d)) &

1821 call post_data(cs%interior_tendency_out_zint_100m(m)%id, &

1822 cs%interior_tendency_out_zint_100m(m)%field_2d(:,:), cs%diag)

1823 enddo

1824

1825 if (cs%ice_ncat > 0) then

1826 do m=1,cs%ice_ncat+1

1827 if (cs%fracr_cat_id(m) > 0) &

1828 call post_data(cs%fracr_cat_id(m), fluxes%fracr_cat(:,:,m), cs%diag)

1829 if (cs%qsw_cat_id(m) > 0) &

1830 call post_data(cs%qsw_cat_id(m), fluxes%qsw_cat(:,:,m), cs%diag)

1831 enddo

1832 endif

1833

1834

1835end subroutine marbl_tracers_column_physics

1836

1837!> This subroutine reads time-varying forcing from files

1838subroutine marbl_tracers_set_forcing(day_start, G, CS)

1839

1840 type(time_type), intent(in) :: day_start !< Start time of the fluxes.

1841 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

1842 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a

1843

1844 real, parameter :: donriv_refract = 0.1 ! Fraction of DON river nutrients in refractory pools [1]

1845 real, parameter :: docriv_refract = 0.2 ! Fraction of DOC river nutrients in refractory pools [1]

1846 real, parameter :: dopriv_refract = 0.025 ! Fraction of DOP river nutrients in refractory pools [1]

1847

1848 real, dimension(SZI_(G),SZJ_(G)) :: riv_flux_in !< The field read in from forcing file with time dimension

1849 !! [mmol m-2 s-1]

1850 type(time_type) :: time_forcing !< For reading river flux fields, we use a modified version of Time

1851 integer :: i, j, k, is, ie, js, je, m

1852

1853 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

1854

1855 ! Abiotic DIC forcing

1856 if (cs%abio_dic_on) then

1857 ! Read d14c bands

1858 do m=1,3

1859 time_forcing = map_model_time_to_forcing_time(day_start, cs%d14c_dataset(m))

1860 call time_interp_external(cs%id_d14c(m),time_forcing,cs%d14c_bands(m))

1861 enddo

1862

1863 ! Set d14c according to the bands

1864 do j=js,je ; do i=is,ie

1865 if (g%geoLatT(i,j) > 30.) then

1866 cs%d14c(i,j) = cs%d14c_bands(1)

1867 elseif (g%geoLatT(i,j) > -30.) then

1868 cs%d14c(i,j) = cs%d14c_bands(2)

1869 else

1870 cs%d14c(i,j) = cs%d14c_bands(3)

1871 endif

1872 enddo ; enddo

1873 endif

1874

1875 ! River fluxes

1876 if (cs%read_riv_fluxes) then

1877 cs%RIV_FLUXES(:,:,:) = 0.

1878 time_forcing = map_model_time_to_forcing_time(day_start, cs%riv_flux_dataset)

1879

1880 ! DIN river flux affects NO3, ALK, and ALK_ALT_CO2

1881 call time_interp_external(cs%id_din_riv,time_forcing,riv_flux_in)

1882

1883 if (cs%tracer_inds%no3_ind > 0) then

1884 do j=js,je ; do i=is,ie

1885 cs%RIV_FLUXES(i,j,cs%tracer_inds%no3_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)

1886 enddo ; enddo

1887 endif

1888 if (cs%tracer_inds%alk_ind > 0) then

1889 do j=js,je ; do i=is,ie

1890 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) = cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) - &

1891 g%mask2dT(i,j) *riv_flux_in(i,j)

1892 enddo ; enddo

1893 endif

1894 if (cs%tracer_inds%alk_alt_co2_ind > 0) then

1895 do j=js,je ; do i=is,ie

1896 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) = &

1897 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) - g%mask2dT(i,j) *riv_flux_in(i,j)

1898 enddo ; enddo

1899 endif

1900

1901 call time_interp_external(cs%id_dip_riv,time_forcing,riv_flux_in)

1902 if (cs%tracer_inds%po4_ind > 0) then

1903 do j=js,je ; do i=is,ie

1904 cs%RIV_FLUXES(i,j,cs%tracer_inds%po4_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)

1905 enddo ; enddo

1906 endif

1907

1908 call time_interp_external(cs%id_don_riv,time_forcing,riv_flux_in)

1909 if (cs%tracer_inds%don_ind > 0) then

1910 do j=js,je ; do i=is,ie

1911 cs%RIV_FLUXES(i,j,cs%tracer_inds%don_ind) = g%mask2dT(i,j) * (1. - donriv_refract) * &

1912 riv_flux_in(i,j)

1913 enddo ; enddo

1914 endif

1915 if (cs%tracer_inds%donr_ind > 0) then

1916 do j=js,je ; do i=is,ie

1917 cs%RIV_FLUXES(i,j,cs%tracer_inds%donr_ind) = g%mask2dT(i,j) * donriv_refract * &

1918 riv_flux_in(i,j)

1919 enddo ; enddo

1920 endif

1921

1922 call time_interp_external(cs%id_dop_riv,time_forcing,riv_flux_in)

1923 if (cs%tracer_inds%dop_ind > 0) then

1924 do j=js,je ; do i=is,ie

1925 cs%RIV_FLUXES(i,j,cs%tracer_inds%dop_ind) = g%mask2dT(i,j) * (1. - dopriv_refract) * &

1926 riv_flux_in(i,j)

1927 enddo ; enddo

1928 endif

1929 if (cs%tracer_inds%dopr_ind > 0) then

1930 do j=js,je ; do i=is,ie

1931 cs%RIV_FLUXES(i,j,cs%tracer_inds%dopr_ind) = g%mask2dT(i,j) * dopriv_refract * &

1932 riv_flux_in(i,j)

1933 enddo ; enddo

1934 endif

1935

1936 call time_interp_external(cs%id_dsi_riv,time_forcing,riv_flux_in)

1937 if (cs%tracer_inds%sio3_ind > 0) then

1938 do j=js,je ; do i=is,ie

1939 cs%RIV_FLUXES(i,j,cs%tracer_inds%sio3_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)

1940 enddo ; enddo

1941 endif

1942

1943 call time_interp_external(cs%id_dfe_riv,time_forcing,riv_flux_in)

1944 if (cs%tracer_inds%fe_ind > 0) then

1945 do j=js,je ; do i=is,ie

1946 cs%RIV_FLUXES(i,j,cs%tracer_inds%fe_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)

1947 enddo ; enddo

1948 endif

1949

1950 call time_interp_external(cs%id_dic_riv,time_forcing,riv_flux_in)

1951 if (cs%tracer_inds%dic_ind > 0) then

1952 do j=js,je ; do i=is,ie

1953 cs%RIV_FLUXES(i,j,cs%tracer_inds%dic_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)

1954 enddo ; enddo

1955 endif

1956 if (cs%tracer_inds%dic_alt_co2_ind > 0) then

1957 do j=js,je ; do i=is,ie

1958 cs%RIV_FLUXES(i,j,cs%tracer_inds%dic_alt_co2_ind) = g%mask2dT(i,j) * riv_flux_in(i,j)

1959 enddo ; enddo

1960 endif

1961

1962 call time_interp_external(cs%id_alk_riv,time_forcing,riv_flux_in)

1963 if (cs%tracer_inds%alk_ind > 0) then

1964 do j=js,je ; do i=is,ie

1965 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) = cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_ind) + &

1966 g%mask2dT(i,j) *riv_flux_in(i,j)

1967 enddo ; enddo

1968 endif

1969 if (cs%tracer_inds%alk_alt_co2_ind > 0) then

1970 do j=js,je ; do i=is,ie

1971 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) = &

1972 cs%RIV_FLUXES(i,j,cs%tracer_inds%alk_alt_co2_ind) + g%mask2dT(i,j) * riv_flux_in(i,j)

1973 enddo ; enddo

1974 endif

1975

1976 call time_interp_external(cs%id_doc_riv,time_forcing,riv_flux_in)

1977 if (cs%tracer_inds%doc_ind > 0) then

1978 do j=js,je ; do i=is,ie

1979 cs%RIV_FLUXES(i,j,cs%tracer_inds%doc_ind) = g%mask2dT(i,j) * (1. - docriv_refract) * &

1980 riv_flux_in(i,j)

1981 enddo ; enddo

1982 endif

1983 if (cs%tracer_inds%docr_ind > 0) then

1984 do j=js,je ; do i=is,ie

1985 cs%RIV_FLUXES(i,j,cs%tracer_inds%docr_ind) = g%mask2dT(i,j) * docriv_refract * &

1986 riv_flux_in(i,j)

1987 enddo ; enddo

1988 endif

1989 endif

1990

1991 ! Tracer restoring

1992 do m=1,cs%restore_count

1993 call time_interp_external(cs%id_tracer_restoring(m),day_start,cs%restoring_in(:,:,:,m))

1994 do k=1,cs%restoring_nz ; do j=js,je ; do i=is,ie

1995 cs%restoring_in(i,j,k,m) = g%mask2dT(i,j) * cs%restoring_in(i,j,k,m)

1996 enddo ; enddo ; enddo

1997 enddo

1998

1999 ! Post Forcing to Diagnostics

2000 if (cs%read_riv_fluxes) then

2001 if (cs%no3_riv_flux > 0 .and. cs%tracer_inds%no3_ind > 0) &

2002 call post_data(cs%no3_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%no3_ind), cs%diag)

2003 if (cs%po4_riv_flux > 0 .and. cs%tracer_inds%po4_ind > 0) &

2004 call post_data(cs%po4_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%po4_ind), cs%diag)

2005 if (cs%don_riv_flux > 0 .and. cs%tracer_inds%don_ind > 0) &

2006 call post_data(cs%don_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%don_ind), cs%diag)

2007 if (cs%donr_riv_flux > 0 .and. cs%tracer_inds%donr_ind > 0) &

2008 call post_data(cs%donr_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%donr_ind), cs%diag)

2009 if (cs%dop_riv_flux > 0 .and. cs%tracer_inds%dop_ind > 0) &

2010 call post_data(cs%dop_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dop_ind), cs%diag)

2011 if (cs%dopr_riv_flux > 0 .and. cs%tracer_inds%dopr_ind > 0) &

2012 call post_data(cs%dopr_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dopr_ind), cs%diag)

2013 if (cs%sio3_riv_flux > 0 .and. cs%tracer_inds%sio3_ind > 0) &

2014 call post_data(cs%sio3_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%sio3_ind), cs%diag)

2015 if (cs%fe_riv_flux > 0 .and. cs%tracer_inds%fe_ind > 0) &

2016 call post_data(cs%fe_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%fe_ind), cs%diag)

2017 if (cs%doc_riv_flux > 0 .and. cs%tracer_inds%doc_ind > 0) &

2018 call post_data(cs%doc_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%doc_ind), cs%diag)

2019 if (cs%docr_riv_flux > 0 .and. cs%tracer_inds%docr_ind > 0) &

2020 call post_data(cs%docr_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%docr_ind), cs%diag)

2021 if (cs%alk_riv_flux > 0 .and. cs%tracer_inds%alk_ind > 0) &

2022 call post_data(cs%alk_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%alk_ind), cs%diag)

2023 if (cs%alk_alt_co2_riv_flux > 0 .and. cs%tracer_inds%alk_alt_co2_ind > 0) &

2024 call post_data(cs%alk_alt_co2_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%alk_alt_co2_ind), &

2025 cs%diag)

2026 if (cs%dic_riv_flux > 0 .and. cs%tracer_inds%dic_ind > 0) &

2027 call post_data(cs%dic_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dic_ind), cs%diag)

2028 if (cs%dic_alt_co2_riv_flux > 0 .and. cs%tracer_inds%dic_alt_co2_ind > 0) &

2029 call post_data(cs%dic_alt_co2_riv_flux, cs%RIV_FLUXES(:,:,cs%tracer_inds%dic_alt_co2_ind), &

2030 cs%diag)

2031 endif

2032 if (cs%abio_dic_on) then

2033 if (cs%d14c_id > 0) &

2034 call post_data(cs%d14c_id, cs%d14c, cs%diag)

2035 endif

2036

2037end subroutine marbl_tracers_set_forcing

2038

2039!> This function calculates the mass-weighted integral of all tracer stocks,

2040!! returning the number of stocks it has calculated. If the stock_index

2041!! is present, only the stock corresponding to that coded index is returned.

2042function marbl_tracers_stock(h, stocks, G, GV, CS, names, units, stock_index)

2043 real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]

2044 type(efp_type), dimension(:), intent(out) :: stocks !< the mass-weighted integrated amount of

2045 !! each tracer, in kg times concentration units

2046 !! [kg conc].

2047 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure

2048 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure

2049 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a

2050 !! previous call to register_MARBL_tracers.

2051 character(len=*), dimension(:), intent(out) :: names !< the names of the stocks calculated.

2052 character(len=*), dimension(:), intent(out) :: units !< the units of the stocks calculated.

2053 integer, optional, intent(in) :: stock_index !< the coded index of a specific stock

2054 !! being sought.

2055 integer :: marbl_tracers_stock !< Return value: the number of stocks

2056 !! calculated here.

2057

2058 ! Local variables

2059 integer :: i, j, k, is, ie, js, je, nz, m

2060 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

2061

2062 marbl_tracers_stock = 0

2063 if (.not.associated(cs)) return

2064 if (cs%ntr < 1) return

2065

2066 if (present(stock_index)) then ; if (stock_index > 0) then

2067 ! Check whether this stock is available from this routine.

2068

2069 ! No stocks from this routine are being checked yet. Return 0.

2070 return

2071 endif ; endif

2072

2073 do m=1,cs%ntr

2074 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller="MARBL_tracers_stock")

2075 units(m) = trim(units(m))//" kg"

2076 stocks(m) = global_mass_int_efp(h, g, gv, cs%tracer_data(m)%tr(:,:,:), on_pe_only=.true.)

2077 enddo

2078 marbl_tracers_stock = cs%ntr

2079

2080end function marbl_tracers_stock

2081

2082!> This subroutine extracts the surface fields from this tracer package that

2083!! are to be shared with the atmosphere in coupled configurations.

2084subroutine marbl_tracers_surface_state(sfc_state, G, US, CS)

2085 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

2086 type(surface), intent(inout) :: sfc_state !< A structure containing fields that

2087 !! describe the surface state of the ocean.

2088 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2089 type(marbl_tracers_cs), pointer :: cs !< The control structure returned by a previous

2090 !! call to register_MARBL_tracers.

2091

2092 ! Local variables

2093 integer :: i, j, is, ie, js, je

2094

2095 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

2096

2097 if (.not.associated(cs)) return

2098

2099 if (allocated(sfc_state%fco2)) then

2100 do j=js,je ; do i=is,ie

2101 ! 44e-6 converts mmol/m^2/s (positive down) to kg CO2/m^2/s (positive down)

2102 sfc_state%fco2(i,j) = us%kg_m2s_to_RZ_T * (44.0e-6 * cs%SFO(i,j,cs%flux_co2_ind))

2103 enddo ; enddo

2104 endif

2105

2106end subroutine marbl_tracers_surface_state

2107

2108!> Copy the requested interior tendency output field into an array.

2109subroutine marbl_tracers_get(name, G, GV, array, CS)

2110

2111 character(len=*), intent(in) :: name !< Name of requested tracer.

2112 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

2113 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.

2114 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

2115 intent(inout) :: array !< Array filled by this routine.

2116 type(marbl_tracers_cs), pointer :: cs !< Pointer to the control structure for this module.

2117

2118 character(len=128), parameter :: sub_name = 'MARBL_tracers_get'

2119 character(len=128) :: log_message

2120

2121 array(:,:,:) = 0.0

2122 select case(trim(name))

2123 case ('Chl')

2124 array(:,:,:) = cs%ITO(:,:,:,cs%total_Chl_ind)

2125 case DEFAULT

2126 write(log_message, "(3A)") "'", trim(name), &

2127 "' is not a valid interior tendency output field name"

2128 call mom_error(fatal, log_message)

2129 end select

2130

2131end subroutine marbl_tracers_get

2132

2133!> Clean up any allocated memory after the run.

2134subroutine marbl_tracers_end(CS)

2135 type(marbl_tracers_cs), pointer, intent(inout) :: cs !< The control structure returned by a previous

2136 !! call to register_MARBL_tracers.

2137

2138 integer :: m

2139

2140 call print_marbl_log(marbl_instances%StatusLog)

2141 call marbl_instances%StatusLog%erase()

2142 call marbl_instances%shutdown()

2143 ! TODO: print MARBL timers to stdout as well

2144

2145 if (associated(cs)) then

2146 if (allocated(cs%tracer_data)) then

2147 do m=1,cs%ntr

2148 if (associated(cs%tracer_data(m)%tr)) deallocate(cs%tracer_data(m)%tr)

2149 enddo

2150 deallocate(cs%tracer_data)

2151 endif

2152 if (allocated(cs%ind_tr)) deallocate(cs%ind_tr)

2153 if (allocated(cs%id_surface_flux_out)) deallocate(cs%id_surface_flux_out)

2154 if (allocated(cs%interior_tendency_out)) deallocate(cs%interior_tendency_out)

2155 if (allocated(cs%interior_tendency_out_zint)) deallocate(cs%interior_tendency_out_zint)

2156 if (allocated(cs%interior_tendency_out_zint_100m)) &

2157 deallocate(cs%interior_tendency_out_zint_100m)

2158 if (allocated(cs%fracr_cat_id)) deallocate(cs%fracr_cat_id)

2159 if (allocated(cs%qsw_cat_id)) deallocate(cs%qsw_cat_id)

2160 if (allocated(cs%STF)) deallocate(cs%STF)

2161 if (allocated(cs%RIV_FLUXES)) deallocate(cs%RIV_FLUXES)

2162 if (associated(cs%SFO)) then

2163 deallocate(cs%SFO)

2164 nullify(cs%SFO)

2165 endif

2166 if (associated(cs%ITO)) then

2167 deallocate(cs%ITO)

2168 nullify(cs%ITO)

2169 endif

2170 if (allocated(cs%tracer_restoring_ind)) deallocate(cs%tracer_restoring_ind)

2171 if (allocated(cs%tracer_I_tau_ind)) deallocate(cs%tracer_I_tau_ind)

2172 if (allocated(cs%fesedflux_in)) deallocate(cs%fesedflux_in)

2173 if (allocated(cs%feventflux_in)) deallocate(cs%feventflux_in)

2174 if (allocated(cs%I_tau)) deallocate(cs%I_tau)

2175 deallocate(cs)

2176 endif

2177end subroutine marbl_tracers_end

2178

2179subroutine set_riv_flux_tracer_inds(CS)

2180

2181 type(marbl_tracers_cs), pointer, intent(inout) :: CS !< The MARBL tracers control structure

2182

2183 character(len=256) :: log_message

2184 character(len=48) :: name ! A variable's name in a NetCDF file.

2185 integer :: m

2186

2187 ! Initialize tracers from file (unless they were initialized by restart file)

2188 ! Also save indices of tracers that have river fluxes

2189 cs%tracer_inds%no3_ind = 0

2190 cs%tracer_inds%po4_ind = 0

2191 cs%tracer_inds%don_ind = 0

2192 cs%tracer_inds%donr_ind = 0

2193 cs%tracer_inds%dop_ind = 0

2194 cs%tracer_inds%dopr_ind = 0

2195 cs%tracer_inds%sio3_ind = 0

2196 cs%tracer_inds%fe_ind = 0

2197 cs%tracer_inds%doc_ind = 0

2198 cs%tracer_inds%docr_ind = 0

2199 cs%tracer_inds%alk_ind = 0

2200 cs%tracer_inds%alk_alt_co2_ind = 0

2201 cs%tracer_inds%dic_ind = 0

2202 cs%tracer_inds%dic_alt_co2_ind = 0

2203 cs%tracer_inds%abio_dic_ind = 0

2204 cs%tracer_inds%abio_di14c_ind = 0

2205 do m=1,cs%ntr

2206 name = marbl_instances%tracer_metadata(m)%short_name

2207 if (trim(name) == "NO3") then

2208 cs%tracer_inds%no3_ind = m

2209 elseif (trim(name) == "PO4") then

2210 cs%tracer_inds%po4_ind = m

2211 elseif (trim(name) == "DON") then

2212 cs%tracer_inds%don_ind = m

2213 elseif (trim(name) == "DONr") then

2214 cs%tracer_inds%donr_ind = m

2215 elseif (trim(name) == "DOP") then

2216 cs%tracer_inds%dop_ind = m

2217 elseif (trim(name) == "DOPr") then

2218 cs%tracer_inds%dopr_ind = m

2219 elseif (trim(name) == "SiO3") then

2220 cs%tracer_inds%sio3_ind = m

2221 elseif (trim(name) == "Fe") then

2222 cs%tracer_inds%fe_ind = m

2223 elseif (trim(name) == "DOC") then

2224 cs%tracer_inds%doc_ind = m

2225 elseif (trim(name) == "DOCr") then

2226 cs%tracer_inds%docr_ind = m

2227 elseif (trim(name) == "ALK") then

2228 cs%tracer_inds%alk_ind = m

2229 elseif (trim(name) == "ALK_ALT_CO2") then

2230 cs%tracer_inds%alk_alt_co2_ind = m

2231 elseif (trim(name) == "DIC") then

2232 cs%tracer_inds%dic_ind = m

2233 elseif (trim(name) == "DIC_ALT_CO2") then

2234 cs%tracer_inds%dic_alt_co2_ind = m

2235 elseif (trim(name) == "ABIO_DIC") then

2236 cs%tracer_inds%abio_dic_ind = m

2237 elseif (trim(name) == "ABIO_DI14C") then

2238 cs%tracer_inds%abio_di14c_ind = m

2239 endif

2240 enddo

2241

2242 ! Log indices for each tracer to ensure we set them all correctly

2243 write(log_message, "(A,I0)") "NO3 index: ", cs%tracer_inds%no3_ind

2244 call mom_error(note, log_message)

2245 write(log_message, "(A,I0)") "PO4 index: ", cs%tracer_inds%po4_ind

2246 call mom_error(note, log_message)

2247 write(log_message, "(A,I0)") "DON index: ", cs%tracer_inds%don_ind

2248 call mom_error(note, log_message)

2249 write(log_message, "(A,I0)") "DONr index: ", cs%tracer_inds%donr_ind

2250 call mom_error(note, log_message)

2251 write(log_message, "(A,I0)") "DOP index: ", cs%tracer_inds%dop_ind

2252 call mom_error(note, log_message)

2253 write(log_message, "(A,I0)") "DOPr index: ", cs%tracer_inds%dopr_ind

2254 call mom_error(note, log_message)

2255 write(log_message, "(A,I0)") "SiO3 index: ", cs%tracer_inds%sio3_ind

2256 call mom_error(note, log_message)

2257 write(log_message, "(A,I0)") "Fe index: ", cs%tracer_inds%fe_ind

2258 call mom_error(note, log_message)

2259 write(log_message, "(A,I0)") "DOC index: ", cs%tracer_inds%doc_ind

2260 call mom_error(note, log_message)

2261 write(log_message, "(A,I0)") "DOCr index: ", cs%tracer_inds%docr_ind

2262 call mom_error(note, log_message)

2263 write(log_message, "(A,I0)") "ALK index: ", cs%tracer_inds%alk_ind

2264 call mom_error(note, log_message)

2265 write(log_message, "(A,I0)") "ALK_ALT_CO2 index: ", cs%tracer_inds%alk_alt_co2_ind

2266 call mom_error(note, log_message)

2267 write(log_message, "(A,I0)") "DIC index: ", cs%tracer_inds%dic_ind

2268 call mom_error(note, log_message)

2269 write(log_message, "(A,I0)") "DIC_ALT_CO2 index: ", cs%tracer_inds%dic_alt_co2_ind

2270 call mom_error(note, log_message)

2271

2272end subroutine set_riv_flux_tracer_inds

2273

2274! TODO: some log messages come from a specific grid point, and this routine

2275! needs to include the location in the preamble

2276!> This subroutine writes the contents of the MARBL log using MOM_error(NOTE, ...).

2277subroutine print_marbl_log(log_to_print, G, i, j)

2278

2279 use marbl_logging, only : marbl_status_log_entry_type

2280 use marbl_logging, only : marbl_log_type

2281 use mom_coms, only : pe_here

2282

2283 class(marbl_log_type), intent(in) :: log_to_print !< MARBL log to include in MOM6 logfile

2284 type(ocean_grid_type), optional, intent(in) :: G !< The ocean's grid structure

2285 integer, optional, intent(in) :: i !< i of (i,j) index of column providing the log

2286 integer, optional, intent(in) :: j !< j of (i,j) index of column providing the log

2287

2288 character(len=*), parameter :: subname = 'MARBL_tracers:print_marbl_log'

2289 character(len=256) :: message_prefix, message_location, log_message

2290 type(marbl_status_log_entry_type), pointer :: tmp

2291 integer :: msg_lev, elem_old

2292

2293 ! elem_old is used to keep track of whether all messages are coming from the same point

2294 elem_old = -1

2295 write(message_prefix, "(A,I0,A)") '(Task ', pe_here(), ')'

2296

2297 tmp => log_to_print%FullLog

2298 do while (associated(tmp))

2299 ! 1) Do I need to write this message? Yes, if all tasks should write this

2300 ! or if I am master_task

2301 if ((.not. tmp%lonly_master_writes) .or. is_root_pe()) then

2302 ! 2) Print message location? (only if ElementInd changed and is positive; requires G)

2303 if ((present(g)) .and. (tmp%ElementInd .ne. elem_old)) then

2304 if (tmp%ElementInd .gt. 0) then

2305 if (present(i) .and. present(j)) then

2306 write(message_location, "(A,F8.3,A,F7.3,A,I0,A,I0,A,I0)") &

2307 'Message from (lon, lat) (', g%geoLonT(i,j), ', ', g%geoLatT(i,j), &

2308 '), which is global (i,j) (', i + g%HI%idg_offset, ', ', j + g%HI%jdg_offset, &

2309 '). Level: ', tmp%ElementInd

2310 else

2311 write(message_location, "(A)") "Grid cell responsible for message is unknown"

2312 endif ! i,j present

2313 ! master task does not need prefix

2314 if (is_root_pe()) then

2315 write(log_message, "(A)") trim(message_location)

2316 msg_lev = note

2317 else

2318 write(log_message, "(A,1X,A)") trim(message_prefix), trim(message_location)

2319 msg_lev = warning

2320 endif ! print message prefix?

2321 call mom_error(msg_lev, log_message, all_print=.true.)

2322 endif ! ElementInd > 0

2323 elem_old = tmp%ElementInd

2324 endif ! ElementInd /= elem_old

2325

2326 ! 3) Write message from the log

2327 ! master task does not need prefix

2328 if (is_root_pe()) then

2329 write(log_message, "(A)") trim(tmp%LogMessage)

2330 msg_lev = note

2331 else

2332 write(log_message, "(A,1X,A)") trim(message_prefix), trim(tmp%LogMessage)

2333 msg_lev = warning

2334 endif ! print message prefix?

2335 call mom_error(msg_lev, log_message, all_print=.true.)

2336 endif ! write the message?

2337 tmp => tmp%next

2338 enddo

2339

2340 if (log_to_print%labort_marbl) then

2341 call mom_error(warning, 'ERROR reported from MARBL library', all_print=.true.)

2342 call mom_error(fatal, 'Stopping in ' // subname)

2343 endif

2344

2345end subroutine print_marbl_log

2346

2347!> \namespace MARBL_tracers

2348!!

2349!! This module contains the code that is needed to provide

2350!! the MARBL BGC tracer library with necessary forcings and

2351!! apply the resulting surface fluxes and tendencies to the

2352!! requested tracers.

2353

2354end module marbl_tracers