tracer_example.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 sample tracer package that has striped initial conditions
7
8use mom_coms, only : efp_type
9use mom_coupler_types, only : set_coupler_type_data, atmos_ocn_coupler_flux
11use mom_error_handler, only : mom_error, fatal, warning
12use mom_file_parser, only : get_param, log_param, log_version, param_file_type
13use mom_forcing_type, only : forcing
14use mom_grid, only : ocean_grid_type
16use mom_io, only : file_exists, mom_read_data, slasher
22use mom_time_manager, only : time_type
23use mom_tracer_registry, only : register_tracer, tracer_registry_type
25use mom_variables, only : surface
27
28implicit none ; private
29
30#include <MOM_memory.h>
31
34
35integer, parameter :: ntr = 1 !< The number of tracers in this module.
36
37!> The control structure for the USER_tracer_example module
38type, public :: user_tracer_example_cs ; private
39 logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
40 character(len=200) :: tracer_ic_file !< The full path to the IC file, or " "
41 !! to initialize internally.
42 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
43 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the tracer registry
44 real, pointer :: tr(:,:,:,:) => null() !< The array of tracers used in this subroutine, perhaps in [g kg-1]?
45 real :: land_val(ntr) = -1.0 !< The value of tr that is used where land is masked out, perhaps in [g kg-1]?
46
47 real :: stripe_width !< The Gaussian width of the stripe in the initial condition
48 !! for the tracer_example tracers [L ~> m]
49 real :: stripe_lat !< The central latitude of the stripe in the initial condition
50 !! for the tracer_example tracers, in [degrees_N] or [km] or [m].
51 logical :: use_sponge !< If true, sponges may be applied somewhere in the domain.
52
53 integer, dimension(NTR) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux if it is used and the
54 !! surface tracer concentrations are to be provided to the coupler.
55
56 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the timing of diagnostic output.
57
58 type(vardesc) :: tr_desc(ntr) !< Descriptions of each of the tracers.
60
61contains
62
63!> This subroutine is used to register tracer fields and subroutines to be used with MOM.
64function user_register_tracer_example(G, GV, US, param_file, CS, tr_Reg, restart_CS)
65 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
66 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
67 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
68 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
69 type(user_tracer_example_cs), pointer :: cs !< A pointer that is set to point to the control
70 !! structure for this module
71 type(tracer_registry_type), pointer :: tr_reg !< A pointer that is set to point to the control
72 !! structure for the tracer advection and
73 !! diffusion module
74 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control struct
75
76! Local variables
77 character(len=80) :: name, longname
78 ! This include declares and sets the variable "version".
79# include "version_variable.h"
80 character(len=40) :: mdl = "tracer_example" ! This module's name.
81 character(len=200) :: inputdir
82 character(len=48) :: flux_units ! The units for tracer fluxes, usually
83 ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
84 real, pointer :: tr_ptr(:,:,:) => null() ! A pointer to one of the tracers, perhaps in [g kg-1]
86 integer :: isd, ied, jsd, jed, nz, m
87 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
88
89 if (associated(cs)) then
90 call mom_error(fatal, "USER_register_tracer_example called with an "// &
91 "associated control structure.")
92 endif
93 allocate(cs)
94
95 ! Read all relevant parameters and write them to the model log.
96 call log_version(param_file, mdl, version, "")
97 call get_param(param_file, mdl, "TRACER_EXAMPLE_IC_FILE", cs%tracer_IC_file, &
98 "The name of a file from which to read the initial conditions for "//&
99 "the tracer_example tracers, or blank to initialize them internally.", &
100 default=" ")
101 if (len_trim(cs%tracer_IC_file) >= 1) then
102 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
103 cs%tracer_IC_file = trim(slasher(inputdir))//trim(cs%tracer_IC_file)
104 call log_param(param_file, mdl, "INPUTDIR/TRACER_EXAMPLE_IC_FILE", &
105 cs%tracer_IC_file)
106 endif
107 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
108 "If true, sponges may be applied anywhere in the domain. "//&
109 "The exact location and properties of those sponges are "//&
110 "specified from MOM_initialization.F90.", default=.false.)
111 call get_param(param_file, mdl, "TRACER_EXAMPLE_STRIPE_WIDTH", cs%stripe_width, &
112 "The Gaussian width of the stripe in the initial condition for the "//&
113 "tracer_example tracers.", units="m", default=1.0e5, scale=us%m_to_L)
114 call get_param(param_file, mdl, "TRACER_EXAMPLE_STRIPE_LAT", cs%stripe_lat, &
115 "The central latitude of the stripe in the initial condition for the "//&
116 "tracer_example tracers.", units=g%y_ax_unit_short, default=40.0)
117
118 allocate(cs%tr(isd:ied,jsd:jed,nz,ntr), source=0.0)
119
120 do m=1,ntr
121 write(name,'("tr",I0)') m
122 write(longname,'("Concentration of Tracer ",I2.2)') m
123 cs%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mdl)
124
125 ! This needs to be changed if the units of tracer are changed above.
126 if (gv%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1"
127 else ; flux_units = "kg s-1" ; endif
128
129 ! This pointer is needed to force the compiler not to do a copy in the registration calls.
130 tr_ptr => cs%tr(:,:,:,m)
131 ! Register the tracer for horizontal advection, diffusion, and restarts.
132 call register_tracer(tr_ptr, tr_reg, param_file, g%HI, gv, &
133 name=name, longname=longname, units="kg kg-1", &
134 registry_diags=.true., flux_units=flux_units, &
135 restart_cs=restart_cs)
136
137 ! Set coupled_tracers to be true (hard-coded above) to provide the surface
138 ! values to the coupler (if any). This is meta-code and its arguments will
139 ! currently (deliberately) give fatal errors if it is used.
140 if (cs%coupled_tracers) &
141 cs%ind_tr(m) = atmos_ocn_coupler_flux(trim(name)//'_flux', &
142 flux_type=' ', implementation=' ', caller="USER_register_tracer_example")
143 enddo
144
145 cs%tr_Reg => tr_reg
148
149!> This subroutine initializes the NTR tracer fields in tr(:,:,:,:)
150!! and it sets up the tracer output.
151subroutine user_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS, &
152 sponge_CSp)
153 logical, intent(in) :: restart !< .true. if the fields have already
154 !! been read from a restart file.
155 type(time_type), target, intent(in) :: day !< Time of the start of the run.
156 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
157 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
158 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
159 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
160 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
161 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
162 !! diagnostic output.
163 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
164 !! whether, where, and what open boundary
165 !! conditions are used.
166 type(user_tracer_example_cs), pointer :: cs !< The control structure returned by a previous
167 !! call to USER_register_tracer_example.
168 type(sponge_cs), pointer :: sponge_csp !< A pointer to the control structure
169 !! for the sponges, if they are in use.
170
171! Local variables
172 real, allocatable :: temp(:,:,:) ! Target values for the tracers in the sponges, perhaps in [g kg-1]
173 character(len=32) :: name ! A variable's name in a NetCDF file.
174 real, pointer :: tr_ptr(:,:,:) => null() ! A pointer to one of the tracers, perhaps in [g kg-1]
175 real :: pi ! 3.1415926... calculated as 4*atan(1) [nondim]
176 real :: tr_y ! Initial zonally uniform tracer concentrations, perhaps in [g kg-1]
177 real :: dist2 ! The distance squared from a line [L2 ~> m2].
178 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
179 integer :: isdb, iedb, jsdb, jedb, lntr
180
181 if (.not.associated(cs)) return
182 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
183 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
184 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
185
186 lntr = ntr ! Avoids compile-time warning when NTR<2
187 cs%Time => day
188 cs%diag => diag
189
190 if (.not.restart) then
191 if (len_trim(cs%tracer_IC_file) >= 1) then
192! Read the tracer concentrations from a netcdf file.
193 if (.not.file_exists(cs%tracer_IC_file, g%Domain)) &
194 call mom_error(fatal, "USER_initialize_tracer: Unable to open "// &
195 cs%tracer_IC_file)
196 do m=1,ntr
197 call query_vardesc(cs%tr_desc(m), name, caller="USER_initialize_tracer")
198 call mom_read_data(cs%tracer_IC_file, trim(name), cs%tr(:,:,:,m), g%Domain)
199 enddo
200 else
201 do m=1,ntr
202 do k=1,nz ; do j=js,je ; do i=is,ie
203 cs%tr(i,j,k,m) = 1.0e-20 ! This could just as well be 0.
204 enddo ; enddo ; enddo
205 enddo
206
207! This sets a stripe of tracer across the basin.
208 pi = 4.0*atan(1.0)
209 do j=js,je
210 dist2 = (g%Rad_Earth_L * pi / 180.0)**2 * (g%geoLatT(i,j) - cs%stripe_lat)**2
211 tr_y = 0.5 * exp( -dist2 / cs%stripe_width**2 )
212
213 do k=1,nz ; do i=is,ie
214! This adds the stripes of tracer to every layer.
215 cs%tr(i,j,k,1) = cs%tr(i,j,k,1) + tr_y
216 enddo ; enddo
217 enddo
218 endif
219 endif ! restart
220
221 if ( cs%use_sponge ) then
222! If sponges are used, this example damps tracers in sponges in the
223! northern half of the domain to 1 and tracers in the southern half
224! to 0. For any tracers that are not damped in the sponge, the call
225! to set_up_sponge_field can simply be omitted.
226 if (.not.associated(sponge_csp)) &
227 call mom_error(fatal, "USER_initialize_tracer: "// &
228 "The pointer to sponge_CSp must be associated if SPONGE is defined.")
229
230 allocate(temp(g%isd:g%ied,g%jsd:g%jed,nz))
231 do k=1,nz ; do j=js,je ; do i=is,ie
232 if ((g%geoLatT(i,j) > 0.5*g%len_lat + g%south_lat) .and. (k > nz/2)) then
233 temp(i,j,k) = 1.0
234 else
235 temp(i,j,k) = 0.0
236 endif
237 enddo ; enddo ; enddo
238
239! do m=1,NTR
240 do m=1,1
241 ! This pointer is needed to force the compiler not to do a copy in the sponge calls.
242 tr_ptr => cs%tr(:,:,:,m)
243 call set_up_sponge_field(temp, tr_ptr, g, gv, nz, sponge_csp)
244 enddo
245 deallocate(temp)
246 endif
247
248 if (associated(obc)) then
249 call query_vardesc(cs%tr_desc(1), name, caller="USER_initialize_tracer")
250 if (obc%specified_v_BCs_exist_globally) then
251 ! Steal from updated DOME in the fullness of time.
252 else
253 ! Steal from updated DOME in the fullness of time.
254 endif
255 ! All tracers but the first have 0 concentration in their inflows. As this
256 ! is the default value, the following calls are unnecessary.
257 !do m=2,lntr
258 do m=2,ntr
259 call query_vardesc(cs%tr_desc(m), name, caller="USER_initialize_tracer")
260 ! Steal from updated DOME in the fullness of time.
261 enddo
262 endif
263
264end subroutine user_initialize_tracer
265
266!> This subroutine applies diapycnal diffusion and any other column
267!! tracer physics or chemistry to the tracers from this file.
268!! This is a simple example of a set of advected passive tracers.
269!! The arguments to this subroutine are redundant in that
270!! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
271subroutine tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS)
272 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
273 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
274 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
275 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
276 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
277 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
278 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
279 intent(in) :: ea !< an array to which the amount of fluid entrained
280 !! from the layer above during this call will be
281 !! added [H ~> m or kg m-2].
282 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
283 intent(in) :: eb !< an array to which the amount of fluid entrained
284 !! from the layer below during this call will be
285 !! added [H ~> m or kg m-2].
286 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
287 !! and tracer forcing fields. Unused fields have NULL ptrs.
288 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
289 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
290 type(user_tracer_example_cs), pointer :: cs !< The control structure returned by a previous
291 !! call to USER_register_tracer_example.
292
293! Local variables
294 real :: hold0(szi_(g)) ! The original topmost layer thickness,
295 ! with surface mass fluxes added back [H ~> m or kg m-2].
296 real :: b1(szi_(g)) ! b1 is a variable used by the tridiagonal solver [H ~> m or kg m-2].
297 real :: c1(szi_(g),szk_(gv)) ! c1 is a variable used by the tridiagonal solver [nondim].
298 real :: d1(szi_(g)) ! d1=1-c1 is used by the tridiagonal solver [nondim].
299 real :: h_neglect ! A thickness that is so small it is usually lost
300 ! in roundoff and can be neglected [H ~> m or kg m-2].
301 real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
302 real :: diapyc_filt ! A multiplicative filter that can be set to 0 to disable diapycnal
303 ! advection of the tracer [nondim]
304 real :: dye_up ! The tracer concentration of upwelled water, perhaps in [g kg-1]?
305 real :: dye_down ! The tracer concentration of downwelled water, perhaps in [g kg-1]?
306 integer :: i, j, k, is, ie, js, je, nz, m
307
308 ! These are the settings for most "physical" tracers, which
309 ! are advected diapycnally in the usual manner.
310 diapyc_filt = 1.0 ; dye_down = 0.0 ; dye_down = 0.0
311
312 ! Uncomment the following line to dye downwelling.
313! diapyc_filt = 0.0 ; dye_down = 1.0
314 ! Uncomment the following line to dye upwelling.
315! diapyc_filt = 0.0 ; dye_up = 1.0
316 ! Uncomment the following line for tracer concentrations to be set
317 ! to zero in any diapycnal motions.
318! diapyc_filt = 0.0 ; dye_down = 0.0 ; dye_down = 0.0
319
320 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
321
322 if (.not.associated(cs)) return
323 h_neglect = gv%H_subroundoff
324
325 do j=js,je
326 do i=is,ie
327! The following line is appropriate for quantities like salinity
328! that are left behind by evaporation, and any surface fluxes would
329! be explicitly included in the flux structure.
330 hold0(i) = h_old(i,j,1)
331! The following line is appropriate for quantities like temperature
332! that can be assumed to have the same concentration in evaporation
333! as they had in the water. The explicit surface fluxes here would
334! reflect differences in concentration from the ambient water, not
335! the absolute fluxes.
336 ! hold0(i) = h_old(i,j,1) + ea(i,j,1)
337 b_denom_1 = h_old(i,j,1) + ea(i,j,1) + h_neglect
338 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
339! d1(i) = b_denom_1 * b1(i)
340 d1(i) = diapyc_filt * (b_denom_1 * b1(i)) + (1.0 - diapyc_filt)
341 do m=1,ntr
342 cs%tr(i,j,1,m) = b1(i)*(hold0(i)*cs%tr(i,j,1,m) + dye_up*eb(i,j,1))
343 ! Add any surface tracer fluxes to the preceding line.
344 enddo
345 enddo
346 do k=2,nz ; do i=is,ie
347 c1(i,k) = diapyc_filt * eb(i,j,k-1) * b1(i)
348 b_denom_1 = h_old(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
349 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
350 d1(i) = diapyc_filt * (b_denom_1 * b1(i)) + (1.0 - diapyc_filt)
351 do m=1,ntr
352 cs%tr(i,j,k,m) = b1(i) * (h_old(i,j,k)*cs%tr(i,j,k,m) + &
353 ea(i,j,k)*(diapyc_filt*cs%tr(i,j,k-1,m) + dye_down) + &
354 eb(i,j,k)*dye_up)
355 enddo
356 enddo ; enddo
357 do m=1,ntr ; do k=nz-1,1,-1 ; do i=is,ie
358 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + c1(i,k+1)*cs%tr(i,j,k+1,m)
359 enddo ; enddo ; enddo
360 enddo
361
362end subroutine tracer_column_physics
363
364!> This function calculates the mass-weighted integral of all tracer stocks,
365!! returning the number of stocks it has calculated. If the stock_index
366!! is present, only the stock corresponding to that coded index is returned.
367function user_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index)
368 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
369 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
370 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
371 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
372 type(efp_type), dimension(:), intent(out) :: stocks !< The mass-weighted integrated amount of each
373 !! tracer, in kg times concentration units [kg conc]
374 type(user_tracer_example_cs), pointer :: cs !< The control structure returned by a
375 !! previous call to register_USER_tracer.
376 character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated.
377 character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated.
378 integer, optional, intent(in) :: stock_index !< The coded index of a specific stock
379 !! being sought.
380 integer :: user_tracer_stock !< Return value: the number of
381 !! stocks calculated here.
382
383 ! Local variables
384 integer :: m
385
387 if (.not.associated(cs)) return
388
389 if (present(stock_index)) then ; if (stock_index > 0) then
390 ! Check whether this stock is available from this routine.
391
392 ! No stocks from this routine are being checked yet. Return 0.
393 return
394 endif ; endif
395
396 do m=1,ntr
397 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller="USER_tracer_stock")
398 units(m) = trim(units(m))//" kg"
399 stocks(m) = global_mass_int_efp(h, g, gv, cs%tr(:,:,:,m), on_pe_only=.true.)
400 enddo
402
403end function user_tracer_stock
404
405!> This subroutine extracts the surface fields from this tracer package that
406!! are to be shared with the atmosphere in coupled configurations.
407subroutine user_tracer_surface_state(sfc_state, h, G, GV, CS)
408 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
409 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
410 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
411 !! describe the surface state of the ocean.
412 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
413 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
414 type(user_tracer_example_cs), pointer :: cs !< The control structure returned by a previous
415 !! call to register_USER_tracer.
416
417 ! This particular tracer package does not report anything back to the coupler.
418 ! The code that is here is just a rough guide for packages that would.
419
420 integer :: m, is, ie, js, je, isd, ied, jsd, jed
421 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
422 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
423
424 if (.not.associated(cs)) return
425
426 if (cs%coupled_tracers) then
427 do m=1,ntr
428 ! This call loads the surface values into the appropriate array in the
429 ! coupler-type structure.
430 call set_coupler_type_data(cs%tr(:,:,1,m), cs%ind_tr(m), sfc_state%tr_fields, &
431 idim=(/isd, is, ie, ied/), jdim=(/jsd, js, je, jed/), turns=g%HI%turns)
432 enddo
433 endif
434
435end subroutine user_tracer_surface_state
436
437!> Clean up allocated memory at the end.
438subroutine user_tracer_example_end(CS)
439 type(user_tracer_example_cs), pointer :: cs !< The control structure returned by a previous
440 !! call to register_USER_tracer.
441
442 if (associated(cs)) then
443 if (associated(cs%tr)) deallocate(cs%tr)
444 deallocate(cs)
445 endif
446end subroutine user_tracer_example_end
447
448!> \namespace user_tracer_example
449!!
450!! Original by Robert Hallberg, 2002
451!!
452!! This file contains an example of the code that is needed to set
453!! up and use a set (in this case one) of dynamically passive tracers.
454!!
455!! A single subroutine is called from within each file to register
456!! each of the tracers for reinitialization and advection and to
457!! register the subroutine that initializes the tracers and set up
458!! their output and the subroutine that does any tracer physics or
459!! chemistry along with diapycnal mixing (included here because some
460!! tracers may float or swim vertically or dye diapycnal processes).
461end module user_tracer_example