advection_test_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!> This tracer package is used to test advection schemes
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 : slasher, vardesc, var_desc, query_vardesc
18use mom_restart, only : query_initialized, set_initialized, mom_restart_cs
21use mom_time_manager, only : time_type
22use 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
35
36integer, parameter :: ntr = 11 !< The number of tracers in this module.
37
38!> The control structure for the advect_test_tracer module
39type, public :: advection_test_tracer_cs ; private
40 integer :: ntr = ntr !< Number of tracers in this module
41 logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
42 character(len=200) :: tracer_ic_file !< The full path to the IC file, or " " to initialize internally.
43 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
44 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the MOM tracer registry
45 real, pointer :: tr(:,:,:,:) => null() !< The array of tracers used in this subroutine [conc]
46 real :: land_val(ntr) = -1.0 !< The value of tr used where land is masked out [conc]
47 logical :: use_sponge !< If true, sponges may be applied somewhere in the domain.
48 logical :: tracers_may_reinit !< If true, the tracers may be set up via the initialization code if
49 !! they are not found in the restart files. Otherwise it is a fatal error
50 !! if the tracers are not found in the restart files of a restarted run.
51 real :: x_origin !< Starting x-position of the tracer [m] or [km] or [degrees_E]
52 real :: x_width !< Initial size in the x-direction of the tracer patch [m] or [km] or [degrees_E]
53 real :: y_origin !< Starting y-position of the tracer [m] or [km] or [degrees_N]
54 real :: y_width !< Initial size in the y-direction of the tracer patch [m] or [km] or [degrees_N]
55
56 integer, dimension(NTR) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux if it is used and
57 !! the surface tracer concentrations are to be provided to the coupler.
58
59 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
60 !! regulate the timing of diagnostic output.
61 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control structure.
62
63 type(vardesc) :: tr_desc(ntr) !< Descriptions and metadata for the tracers
65
66contains
67
68!> Register tracer fields and subroutines to be used with MOM.
69function register_advection_test_tracer(G, GV, param_file, CS, tr_Reg, restart_CS)
70 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
71 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
72 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
73 type(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
74 !! call to register_advection_test_tracer.
75 type(tracer_registry_type), pointer :: tr_reg !< A pointer that is set to point to the control
76 !! structure for the tracer advection and
77 !! diffusion module
78 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control struct
79
80 ! Local variables
81 character(len=80) :: name, longname
82 ! This include declares and sets the variable "version".
83# include "version_variable.h"
84 character(len=40) :: mdl = "advection_test_tracer" ! This module's name.
85 character(len=200) :: inputdir ! The directory where the input file can be found
86 character(len=48) :: flux_units ! The units for tracer fluxes, usually
87 ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
88 real, pointer :: tr_ptr(:,:,:) => null() ! A pointer to a tracer array [conc]
90 integer :: isd, ied, jsd, jed, nz, m
91 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
92
93 if (associated(cs)) then
94 call mom_error(fatal, "register_advection_test_tracer called with an "// &
95 "associated control structure.")
96 endif
97 allocate(cs)
98
99 ! Read all relevant parameters and write them to the model log.
100 call log_version(param_file, mdl, version, "")
101
102 call get_param(param_file, mdl, "ADVECTION_TEST_X_ORIGIN", cs%x_origin, &
103 "The x-coordinate of the center of the test-functions.", units=g%x_ax_unit_short, default=0.)
104 call get_param(param_file, mdl, "ADVECTION_TEST_Y_ORIGIN", cs%y_origin, &
105 "The y-coordinate of the center of the test-functions.", units=g%y_ax_unit_short, default=0.)
106 call get_param(param_file, mdl, "ADVECTION_TEST_X_WIDTH", cs%x_width, &
107 "The x-width of the test-functions.", units=g%x_ax_unit_short, default=0.)
108 call get_param(param_file, mdl, "ADVECTION_TEST_Y_WIDTH", cs%y_width, &
109 "The y-width of the test-functions.", units=g%y_ax_unit_short, default=0.)
110 call get_param(param_file, mdl, "ADVECTION_TEST_TRACER_IC_FILE", cs%tracer_IC_file, &
111 "The name of a file from which to read the initial "//&
112 "conditions for the tracers, or blank to initialize "//&
113 "them internally.", default=" ")
114
115 if (len_trim(cs%tracer_IC_file) >= 1) then
116 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
117 cs%tracer_IC_file = trim(slasher(inputdir))//trim(cs%tracer_IC_file)
118 call log_param(param_file, mdl, "INPUTDIR/ADVECTION_TEST_TRACER_IC_FILE", &
119 cs%tracer_IC_file)
120 endif
121 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
122 "If true, sponges may be applied anywhere in the domain. "//&
123 "The exact location and properties of those sponges are "//&
124 "specified from MOM_initialization.F90.", default=.false.)
125
126 call get_param(param_file, mdl, "TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
127 "If true, tracers may go through the initialization code "//&
128 "if they are not found in the restart files. Otherwise "//&
129 "it is a fatal error if the tracers are not found in the "//&
130 "restart files of a restarted run.", default=.false.)
131
132
133 allocate(cs%tr(isd:ied,jsd:jed,nz,ntr), source=0.0)
134
135 do m=1,ntr
136 write(name,'("tr",I0)') m
137 write(longname,'("Concentration of Tracer ",I2.2)') m
138 cs%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mdl)
139 if (gv%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1"
140 else ; flux_units = "kg s-1" ; endif
141
142
143 ! This is needed to force the compiler not to do a copy in the registration
144 ! calls. Curses on the designers and implementers of Fortran90.
145 tr_ptr => cs%tr(:,:,:,m)
146 ! Register the tracer for horizontal advection, diffusion, and restarts.
147 call register_tracer(tr_ptr, tr_reg, param_file, g%HI, gv, &
148 name=name, longname=longname, units="kg kg-1", &
149 registry_diags=.true., flux_units=flux_units, &
150 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
151
152 ! Set coupled_tracers to be true (hard-coded above) to provide the surface
153 ! values to the coupler (if any). This is meta-code and its arguments will
154 ! currently (deliberately) give fatal errors if it is used.
155 if (cs%coupled_tracers) &
156 cs%ind_tr(m) = atmos_ocn_coupler_flux(trim(name)//'_flux', &
157 flux_type=' ', implementation=' ', caller="register_advection_test_tracer")
158 enddo
159
160 cs%tr_Reg => tr_reg
161 cs%restart_CSp => restart_cs
164
165!> Initializes the NTR tracer fields in tr(:,:,:,:) and it sets up the tracer output.
166subroutine initialize_advection_test_tracer(restart, day, G, GV, h,diag, OBC, CS, &
167 sponge_CSp)
168 logical, intent(in) :: restart !< .true. if the fields have already
169 !! been read from a restart file.
170 type(time_type), target, intent(in) :: day !< Time of the start of the run.
171 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
172 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
173 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
174 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
175 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
176 !! diagnostic output.
177 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
178 !! whether, where, and what open boundary
179 !! conditions are used.
180 type(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
181 !! call to register_advection_test_tracer.
182 type(sponge_cs), pointer :: sponge_csp !< Pointer to the control structure for the sponges.
183
184 ! Local variables
185 character(len=16) :: name ! A variable's name in a NetCDF file.
186 real :: locx, locy ! x- and y- positions relative to the center of the tracer patch
187 ! normalized by its size [nondim]
188 real :: h_neglect ! A thickness that is so small it is usually lost
189 ! in roundoff and can be neglected [H ~> m or kg m-2].
190 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
191 integer :: isdb, iedb, jsdb, jedb
192
193 if (.not.associated(cs)) return
194 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
195 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
196 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
197 h_neglect = gv%H_subroundoff
198
199 cs%diag => diag
200 cs%ntr = ntr
201 do m=1,ntr
202 call query_vardesc(cs%tr_desc(m), name=name, &
203 caller="initialize_advection_test_tracer")
204 if ((.not.restart) .or. (cs%tracers_may_reinit .and. .not. &
205 query_initialized(cs%tr(:,:,:,m), name, cs%restart_CSp))) then
206 do k=1,nz ; do j=js,je ; do i=is,ie
207 cs%tr(i,j,k,m) = 0.0
208 enddo ; enddo ; enddo
209 k=1 ! Square wave
210 do j=js,je ; do i=is,ie
211 if (abs(g%geoLonT(i,j)-cs%x_origin)<0.5*cs%x_width .and. &
212 abs(g%geoLatT(i,j)-cs%y_origin)<0.5*cs%y_width) cs%tr(i,j,k,m) = 1.0
213 enddo ; enddo
214 k=2 ! Triangle wave
215 do j=js,je ; do i=is,ie
216 locx = abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
217 locy = abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
218 cs%tr(i,j,k,m) = max(0.0, 1.0-locx)*max(0.0, 1.0-locy)
219 enddo ; enddo
220 k=3 ! Cosine bell
221 do j=js,je ; do i=is,ie
222 locx = min(1.0, abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width) * (acos(0.0)*2.)
223 locy = min(1.0, abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width) * (acos(0.0)*2.)
224 cs%tr(i,j,k,m) = (1.0+cos(locx))*(1.0+cos(locy))*0.25
225 enddo ; enddo
226 k=4 ! Cylinder
227 do j=js,je ; do i=is,ie
228 locx = abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
229 locy = abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
230 if ((locx**2) + (locy**2) <= 1.0) cs%tr(i,j,k,m) = 1.0
231 enddo ; enddo
232 k=5 ! Cut cylinder
233 do j=js,je ; do i=is,ie
234 locx = (g%geoLonT(i,j)-cs%x_origin)/cs%x_width
235 locy = (g%geoLatT(i,j)-cs%y_origin)/cs%y_width
236 if ((locx**2) + (locy**2) <= 1.0) cs%tr(i,j,k,m) = 1.0
237 if (locx>0.0 .and. abs(locy)<0.2) cs%tr(i,j,k,m) = 0.0
238 enddo ; enddo
239
240 call set_initialized(cs%tr(:,:,:,m), name, cs%restart_CSp)
241 endif ! restart
242 enddo
243
244
246
247
248!> Applies diapycnal diffusion and any other column tracer physics or chemistry to the tracers
249!! from this package. This is a simple example of a set of advected passive tracers.
250subroutine advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
251 evap_CFL_limit, minimum_forcing_depth)
252 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
253 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
254 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
255 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
256 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
257 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
258 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
259 intent(in) :: ea !< an array to which the amount of fluid entrained
260 !! from the layer above during this call will be
261 !! added [H ~> m or kg m-2].
262 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
263 intent(in) :: eb !< an array to which the amount of fluid entrained
264 !! from the layer below during this call will be
265 !! added [H ~> m or kg m-2].
266 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
267 !! and tracer forcing fields. Unused fields have NULL ptrs.
268 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
269 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
270 type(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
271 !! call to register_advection_test_tracer.
272 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
273 !! be fluxed out of the top layer in a timestep [nondim]
274 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
275 !! fluxes can be applied [H ~> m or kg m-2]
276! This subroutine applies diapycnal diffusion and any other column
277! tracer physics or chemistry to the tracers from this file.
278! This is a simple example of a set of advected passive tracers.
279
280! The arguments to this subroutine are redundant in that
281! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
282 ! Local variables
283 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
284 integer :: i, j, k, is, ie, js, je, nz, m
285 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
286
287 if (.not.associated(cs)) return
288
289 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
290 do m=1,cs%ntr
291 do k=1,nz ;do j=js,je ; do i=is,ie
292 h_work(i,j,k) = h_old(i,j,k)
293 enddo ; enddo ; enddo
294 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,m), dt, fluxes, h_work, &
295 evap_cfl_limit, minimum_forcing_depth)
296 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
297 enddo
298 else
299 do m=1,ntr
300 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
301 enddo
302 endif
303
305
306!> This subroutine extracts the surface fields from this tracer package that
307!! are to be shared with the atmosphere in coupled configurations.
308!! This particular tracer package does not report anything back to the coupler.
309subroutine advection_test_tracer_surface_state(sfc_state, h, G, GV, CS)
310 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
311 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
312 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
313 !! describe the surface state of the ocean.
314 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
315 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
316 type(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
317 !! call to register_advection_test_tracer.
318
319 ! This particular tracer package does not report anything back to the coupler.
320 ! The code that is here is just a rough guide for packages that would.
321
322 integer :: m, is, ie, js, je, isd, ied, jsd, jed
323 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
324 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
325
326 if (.not.associated(cs)) return
327
328 if (cs%coupled_tracers) then
329 do m=1,cs%ntr
330 ! This call loads the surface values into the appropriate array in the
331 ! coupler-type structure.
332 call set_coupler_type_data(cs%tr(:,:,1,m), cs%ind_tr(m), sfc_state%tr_fields, &
333 idim=(/isd, is, ie, ied/), jdim=(/jsd, js, je, jed/), turns=g%HI%turns)
334 enddo
335 endif
336
338
339!> Calculate the mass-weighted integral of all tracer stocks, returning the number of stocks it has calculated.
340!! If the stock_index is present, only the stock corresponding to that coded index is returned.
341function advection_test_stock(h, stocks, G, GV, CS, names, units, stock_index)
342 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
343 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
344 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
345 type(efp_type), dimension(:), intent(out) :: stocks !< the mass-weighted integrated amount of each
346 !! tracer, in kg times concentration units [kg conc].
347 type(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
348 !! call to register_advection_test_tracer.
349 character(len=*), dimension(:), intent(out) :: names !< the names of the stocks calculated.
350 character(len=*), dimension(:), intent(out) :: units !< the units of the stocks calculated.
351 integer, optional, intent(in) :: stock_index !< the coded index of a specific stock being sought.
352 integer :: advection_test_stock !< the number of stocks calculated here.
353
354 ! Local variables
355 integer :: is, ie, js, je, nz, m
356 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
357
359 if (.not.associated(cs)) return
360 if (cs%ntr < 1) return
361
362 if (present(stock_index)) then ; if (stock_index > 0) then
363 ! Check whether this stock is available from this routine.
364
365 ! No stocks from this routine are being checked yet. Return 0.
366 return
367 endif ; endif
368
369 do m=1,cs%ntr
370 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller="advection_test_stock")
371 stocks(m) = global_mass_int_efp(h, g, gv, cs%tr(:,:,:,m), on_pe_only=.true.)
372 enddo
373 advection_test_stock = cs%ntr
374
375end function advection_test_stock
376
377!> Deallocate memory associated with this module
378subroutine advection_test_tracer_end(CS)
379 type(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
380 !! call to register_advection_test_tracer.
381 if (associated(cs)) then
382 if (associated(cs%tr)) deallocate(cs%tr)
383 deallocate(cs)
384 endif
385end subroutine advection_test_tracer_end
386
387end module advection_test_tracer