DOME_tracer.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 that is used as a diagnostic in the DOME experiments
6module dome_tracer
7
8use mom_coupler_types, only : set_coupler_type_data, atmos_ocn_coupler_flux
10use mom_error_handler, only : mom_error, fatal, warning
11use mom_file_parser, only : get_param, log_param, log_version, param_file_type
12use mom_forcing_type, only : forcing
14use mom_grid, only : ocean_grid_type
15use mom_interface_heights, only : thickness_to_dz
16use mom_io, only : file_exists, mom_read_data, slasher, vardesc, var_desc, query_vardesc
21use mom_time_manager, only : time_type
22use mom_tracer_registry, only : register_tracer, tracer_registry_type
27
28implicit none ; private
29
30#include <MOM_memory.h>
31
34
35! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
36! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
37! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
38! vary with the Boussinesq approximation, the Boussinesq variant is given first.
39
40integer, parameter :: ntr = 11 !< The number of tracers in this module.
41
42!> The DOME_tracer control structure
43type, public :: dome_tracer_cs ; private
44 logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
45 character(len=200) :: tracer_ic_file !< The full path to the IC file, or " " to initialize internally.
46 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
47 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the tracer registry
48 real, pointer :: tr(:,:,:,:) => null() !< The array of tracers used in this package, perhaps in [g kg-1]
49 real :: land_val(ntr) = -1.0 !< The value of tr used where land is masked out, perhaps in [g kg-1]
50 logical :: use_sponge !< If true, sponges may be applied somewhere in the domain.
51
52 real :: stripe_width !< The meridional width of the vertical stripes in the initial condition
53 !! for some of the DOME tracers, in [km] or [degrees_N] or [m].
54 real :: stripe_s_lat !< The southern latitude of the first vertical stripe in the initial condition
55 !! for some of the DOME tracers, in [km] or [degrees_N] or [m].
56 real :: sheet_spacing !< The vertical spacing between successive horizontal sheets of tracer in the initial
57 !! conditions for some of the DOME tracers [Z ~> m], and twice the thickness of
58 !! these horizontal tracer sheets
59
60 integer, dimension(NTR) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux if it is used and the
61 !! surface tracer concentrations are to be provided to the coupler.
62
63 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
64 !! regulate the timing of diagnostic output.
65
66 type(vardesc) :: tr_desc(ntr) !< Descriptions and metadata for the tracers
67end type dome_tracer_cs
68
69contains
70
71!> Register tracer fields and subroutines to be used with MOM.
72function register_dome_tracer(G, GV, US, param_file, CS, tr_Reg, restart_CS)
73 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
74 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
75 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
76 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
77 type(dome_tracer_cs), pointer :: cs !< A pointer that is set to point to the
78 !! control structure for this module
79 type(tracer_registry_type), pointer :: tr_reg !< A pointer to the tracer registry.
80 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control struct
81
82 ! Local variables
83 character(len=80) :: name, longname
84 ! This include declares and sets the variable "version".
85# include "version_variable.h"
86 character(len=40) :: mdl = "DOME_tracer" ! This module's name.
87 character(len=48) :: flux_units ! The units for tracer fluxes, usually
88 ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
89 character(len=200) :: inputdir
90 real, pointer :: tr_ptr(:,:,:) => null() ! A pointer to one of the tracers, perhaps in [g kg-1]
91 logical :: register_dome_tracer
92 integer :: isd, ied, jsd, jed, nz, m
93 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
94
95 if (associated(cs)) then
96 call mom_error(fatal, "DOME_register_tracer called with an "// &
97 "associated control structure.")
98 endif
99 allocate(cs)
100
101 ! Read all relevant parameters and write them to the model log.
102 call log_version(param_file, mdl, version, "")
103 call get_param(param_file, mdl, "DOME_TRACER_IC_FILE", cs%tracer_IC_file, &
104 "The name of a file from which to read the initial "//&
105 "conditions for the DOME tracers, or blank to initialize "//&
106 "them internally.", default=" ")
107 if (len_trim(cs%tracer_IC_file) >= 1) then
108 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
109 inputdir = slasher(inputdir)
110 cs%tracer_IC_file = trim(inputdir)//trim(cs%tracer_IC_file)
111 call log_param(param_file, mdl, "INPUTDIR/DOME_TRACER_IC_FILE", &
112 cs%tracer_IC_file)
113 endif
114 call get_param(param_file, mdl, "DOME_TRACER_STRIPE_WIDTH", cs%stripe_width, &
115 "The meridional width of the vertical stripes in the initial condition "//&
116 "for the DOME tracers.", units=g%y_ax_unit_short, default=50.0)
117 call get_param(param_file, mdl, "DOME_TRACER_STRIPE_LAT", cs%stripe_s_lat, &
118 "The southern latitude of the first vertical stripe in the initial condition "//&
119 "for the DOME tracers.", units=g%y_ax_unit_short, default=350.0)
120 call get_param(param_file, mdl, "DOME_TRACER_SHEET_SPACING", cs%sheet_spacing, &
121 "The vertical spacing between successive horizontal sheets of tracer in the initial "//&
122 "conditions for the DOME tracers, and twice the thickness of these tracer sheets.", &
123 units="m", default=600.0, scale=us%m_to_Z)
124 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
125 "If true, sponges may be applied anywhere in the domain. "//&
126 "The exact location and properties of those sponges are "//&
127 "specified from MOM_initialization.F90.", default=.false.)
128
129 allocate(cs%tr(isd:ied,jsd:jed,nz,ntr), source=0.0)
130
131 do m=1,ntr
132 write(name,'("tr_D",I0)') m
133 write(longname,'("Concentration of DOME Tracer ",I2.2)') m
134 cs%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mdl)
135 if (gv%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1"
136 else ; flux_units = "kg s-1" ; endif
137
138 ! This is needed to force the compiler not to do a copy in the registration
139 ! calls. Curses on the designers and implementers of Fortran90.
140 tr_ptr => cs%tr(:,:,:,m)
141 ! Register the tracer for horizontal advection, diffusion, and restarts.
142 call register_tracer(tr_ptr, tr_reg, param_file, g%HI, gv, &
143 name=name, longname=longname, units="kg kg-1", &
144 registry_diags=.true., restart_cs=restart_cs, &
145 flux_units=trim(flux_units), flux_scale=gv%H_to_MKS)
146
147 ! Set coupled_tracers to be true (hard-coded above) to provide the surface
148 ! values to the coupler (if any). This is meta-code and its arguments will
149 ! currently (deliberately) give fatal errors if it is used.
150 if (cs%coupled_tracers) &
151 cs%ind_tr(m) = atmos_ocn_coupler_flux(trim(name)//'_flux', &
152 flux_type=' ', implementation=' ', caller="register_DOME_tracer")
153 enddo
154
155 cs%tr_Reg => tr_reg
156 register_dome_tracer = .true.
157end function register_dome_tracer
158
159!> Initializes the NTR tracer fields in tr(:,:,:,:) and sets up the tracer output.
160subroutine initialize_dome_tracer(restart, day, G, GV, US, h, diag, OBC, CS, &
161 sponge_CSp, tv)
162 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
163 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
164 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
165 logical, intent(in) :: restart !< .true. if the fields have already
166 !! been read from a restart file.
167 type(time_type), target, intent(in) :: day !< Time of the start of the run.
168 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
169 type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output.
170 type(ocean_obc_type), pointer :: obc !< Structure specifying open boundary options.
171 type(dome_tracer_cs), pointer :: cs !< The control structure returned by a previous
172 !! call to DOME_register_tracer.
173 type(sponge_cs), pointer :: sponge_csp !< A pointer to the control structure
174 !! for the sponges, if they are in use.
175 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
176
177 ! Local variables
178 real, allocatable :: temp(:,:,:) ! Target values for the tracers in the sponges, perhaps in [g kg-1]
179 character(len=16) :: name ! A variable's name in a NetCDF file.
180 real, pointer :: tr_ptr(:,:,:) => null() ! A pointer to one of the tracers, perhaps in [g kg-1]
181 real :: dz(szi_(g),szk_(gv)) ! Height change across layers [Z ~> m]
182 real :: tr_y ! Initial zonally uniform tracer concentrations, perhaps in [g kg-1]
183 real :: dz_neglect ! A thickness that is so small it is usually lost
184 ! in roundoff and can be neglected [Z ~> m]
185 real :: e(szk_(gv)+1) ! Interface heights relative to the sea surface (negative down) [Z ~> m]
186 real :: e_top ! Height of the top of the tracer band relative to the sea surface [Z ~> m]
187 real :: e_bot ! Height of the bottom of the tracer band relative to the sea surface [Z ~> m]
188 real :: d_tr ! A change in tracer concentrations, in tracer units, perhaps [g kg-1]
189 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
190
191 if (.not.associated(cs)) return
192 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
193 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
194
195 dz_neglect = gv%dz_subroundoff
196
197 cs%Time => day
198 cs%diag => diag
199
200 if (.not.restart) then
201 if (len_trim(cs%tracer_IC_file) >= 1) then
202 ! Read the tracer concentrations from a netcdf file.
203 if (.not.file_exists(cs%tracer_IC_file, g%Domain)) &
204 call mom_error(fatal, "DOME_initialize_tracer: Unable to open "// &
205 cs%tracer_IC_file)
206 do m=1,ntr
207 call query_vardesc(cs%tr_desc(m), name, caller="initialize_DOME_tracer")
208 call mom_read_data(cs%tracer_IC_file, trim(name), cs%tr(:,:,:,m), g%Domain)
209 enddo
210 else
211 do m=1,ntr
212 do k=1,nz ; do j=js,je ; do i=is,ie
213 cs%tr(i,j,k,m) = 1.0e-20 ! This could just as well be 0.
214 enddo ; enddo ; enddo
215 enddo
216
217! This sets a stripe of tracer across the basin.
218 do m=2,min(6,ntr) ; do j=js,je ; do i=is,ie
219 tr_y = 0.0
220 if ((g%geoLatT(i,j) > (cs%stripe_s_lat + cs%stripe_width*real(m-2))) .and. &
221 (g%geoLatT(i,j) < (cs%stripe_s_lat + cs%stripe_width*real(m-1)))) &
222 tr_y = 1.0
223 do k=1,nz
224! This adds the stripes of tracer to every layer.
225 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + tr_y
226 enddo
227 enddo ; enddo ; enddo
228
229 if (ntr >= 7) then
230 do j=js,je
231 call thickness_to_dz(h, tv, dz, j, g, gv)
232 do i=is,ie
233 e(1) = 0.0
234 do k=1,nz
235 e(k+1) = e(k) - dz(i,k)
236 do m=7,ntr
237 e_top = -cs%sheet_spacing * (real(m-6))
238 e_bot = -cs%sheet_spacing * (real(m-6) + 0.5)
239 if (e_top < e(k)) then
240 if (e_top < e(k+1)) then ; d_tr = 0.0
241 elseif (e_bot < e(k+1)) then
242 d_tr = 1.0 * (e_top-e(k+1)) / (dz(i,k)+dz_neglect)
243 else ; d_tr = 1.0 * (e_top-e_bot) / (dz(i,k)+dz_neglect)
244 endif
245 elseif (e_bot < e(k)) then
246 if (e_bot < e(k+1)) then ; d_tr = 1.0
247 else ; d_tr = 1.0 * (e(k)-e_bot) / (dz(i,k)+dz_neglect)
248 endif
249 else
250 d_tr = 0.0
251 endif
252 if (dz(i,k) < 2.0*gv%Angstrom_Z) d_tr=0.0
253 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + d_tr
254 enddo
255 enddo
256 enddo
257 enddo
258 endif
259
260 endif
261 endif ! restart
262
263 if ( cs%use_sponge ) then
264! If sponges are used, this example damps tracers in sponges in the
265! northern half of the domain to 1 and tracers in the southern half
266! to 0. For any tracers that are not damped in the sponge, the call
267! to set_up_sponge_field can simply be omitted.
268 if (.not.associated(sponge_csp)) &
269 call mom_error(fatal, "DOME_initialize_tracer: "// &
270 "The pointer to sponge_CSp must be associated if SPONGE is defined.")
271
272 allocate(temp(g%isd:g%ied,g%jsd:g%jed,nz))
273 do k=1,nz ; do j=js,je ; do i=is,ie
274 if (g%geoLatT(i,j) > 700.0 .and. (k > nz/2)) then
275 temp(i,j,k) = 1.0
276 else
277 temp(i,j,k) = 0.0
278 endif
279 enddo ; enddo ; enddo
280
281! do m=1,NTR
282 do m=1,1
283 ! This pointer is needed to force the compiler not to do a copy in the sponge calls.
284 tr_ptr => cs%tr(:,:,:,m)
285 call set_up_sponge_field(temp, tr_ptr, g, gv, nz, sponge_csp)
286 enddo
287 deallocate(temp)
288 endif
289
290end subroutine initialize_dome_tracer
291
292!> This subroutine applies diapycnal diffusion and any other column
293!! tracer physics or chemistry to the tracers from this file.
294!! This is a simple example of a set of advected passive tracers.
295!!
296!! The arguments to this subroutine are redundant in that
297!! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
298subroutine dome_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
299 evap_CFL_limit, minimum_forcing_depth)
300 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
301 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
302 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
303 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
304 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
305 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
306 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
307 intent(in) :: ea !< an array to which the amount of fluid entrained
308 !! from the layer above during this call will be
309 !! added [H ~> m or kg m-2].
310 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
311 intent(in) :: eb !< an array to which the amount of fluid entrained
312 !! from the layer below during this call will be
313 !! added [H ~> m or kg m-2].
314 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
315 !! and tracer forcing fields. Unused fields have NULL ptrs.
316 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
317 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
318 type(dome_tracer_cs), pointer :: cs !< The control structure returned by a previous
319 !! call to DOME_register_tracer.
320 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
321 !! be fluxed out of the top layer in a timestep [nondim]
322 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
323 !! fluxes can be applied [H ~> m or kg m-2]
324
325! Local variables
326 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
327 integer :: i, j, k, is, ie, js, je, nz, m
328 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
329
330 if (.not.associated(cs)) return
331
332 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
333 do m=1,ntr
334 do k=1,nz ;do j=js,je ; do i=is,ie
335 h_work(i,j,k) = h_old(i,j,k)
336 enddo ; enddo ; enddo
337 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,m), dt, fluxes, h_work, &
338 evap_cfl_limit, minimum_forcing_depth)
339 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
340 enddo
341 else
342 do m=1,ntr
343 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
344 enddo
345 endif
346
347end subroutine dome_tracer_column_physics
348
349!> This subroutine extracts the surface fields from this tracer package that
350!! are to be shared with the atmosphere in coupled configurations.
351!! This particular tracer package does not report anything back to the coupler.
352subroutine dome_tracer_surface_state(sfc_state, h, G, GV, CS)
353 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
354 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
355 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
356 !! describe the surface state of the ocean.
357 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
358 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
359 type(dome_tracer_cs), pointer :: cs !< The control structure returned by a previous
360 !! call to DOME_register_tracer.
361
362 ! This particular tracer package does not report anything back to the coupler.
363 ! The code that is here is just a rough guide for packages that would.
364
365 integer :: m, is, ie, js, je, isd, ied, jsd, jed
366 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
367 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
368
369 if (.not.associated(cs)) return
370
371 if (cs%coupled_tracers) then
372 do m=1,ntr
373 ! This call loads the surface values into the appropriate array in the
374 ! coupler-type structure.
375 call set_coupler_type_data(cs%tr(:,:,1,m), cs%ind_tr(m), sfc_state%tr_fields, &
376 idim=(/isd, is, ie, ied/), jdim=(/jsd, js, je, jed/), turns=g%HI%turns)
377 enddo
378 endif
379
380end subroutine dome_tracer_surface_state
381
382!> Clean up memory allocations, if any.
383subroutine dome_tracer_end(CS)
384 type(dome_tracer_cs), pointer :: cs !< The control structure returned by a previous
385 !! call to DOME_register_tracer.
386 if (associated(cs)) then
387 if (associated(cs%tr)) deallocate(cs%tr)
388 deallocate(cs)
389 endif
390end subroutine dome_tracer_end
391
392!> \namespace dome_tracer
393!!
394!! By Robert Hallberg, 2002
395!!
396!! This file contains an example of the code that is needed to set
397!! up and use a set (in this case eleven) of dynamically passive
398!! tracers. These tracers dye the inflowing water or water initially
399!! within a range of latitudes or water initially in a range of
400!! depths.
401!!
402!! A single subroutine is called from within each file to register
403!! each of the tracers for reinitialization and advection and to
404!! register the subroutine that initializes the tracers and set up
405!! their output and the subroutine that does any tracer physics or
406!! chemistry along with diapycnal mixing (included here because some
407!! tracers may float or swim vertically or dye diapycnal processes).
408
409end module dome_tracer