boundary_impulse_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!> Implements a boundary impulse response tracer to calculate Green's functions
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
18use mom_restart, only : register_restart_field, query_initialized, set_initialized, mom_restart_cs
21use mom_time_manager, only : time_type
22use mom_tracer_registry, only : register_tracer, tracer_registry_type
28
29implicit none ; private
30
31#include <MOM_memory.h>
32
36
37!> NTR_MAX is the maximum number of tracers in this module.
38integer, parameter :: ntr_max = 1
39
40!> The control structure for the boundary impulse tracer package
41type, public :: boundary_impulse_tracer_cs ; private
42 integer :: ntr=ntr_max !< The number of tracers that are actually used.
43 logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
44 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
45 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the tracer registry
46 real, pointer :: tr(:,:,:,:) => null() !< The array of tracers used in this subroutine, in [CU ~> conc] (g m-3)?
47 logical :: tracers_may_reinit !< If true, boundary_impulse can be initialized if not found in restart file
48 integer, dimension(NTR_MAX) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux if it is used and the
49 !! surface tracer concentrations are to be provided to the coupler.
50
51 integer :: nkml !< Number of layers in mixed layer
52 real, dimension(NTR_MAX) :: land_val = -1.0 !< A value to use to fill in tracers over land [CU ~> conc]
53 real :: remaining_source_time !< How much longer (same units as the timestep) to
54 !! inject the tracer at the surface [T ~> s]
55
56 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
57 !! regulate the timing of diagnostic output.
58 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the retart control structure
59
60 type(vardesc) :: tr_desc(ntr_max) !< Descriptions and metadata for the tracers
62
63contains
64
65!> Read in runtime options and add boundary impulse tracer to tracer registry
66function register_boundary_impulse_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
67 type(hor_index_type), intent(in ) :: hi !< A horizontal index type structure
68 type(verticalgrid_type), intent(in ) :: gv !< The ocean's vertical grid structure
69 type(unit_scale_type), intent(in ) :: us !< A dimensional unit scaling type
70 type(param_file_type), intent(in ) :: param_file !< A structure to parse for run-time parameters
71 type(boundary_impulse_tracer_cs), pointer :: cs !< The control structure returned by a previous
72 !! call to register_boundary_impulse_tracer.
73 type(tracer_registry_type), pointer :: tr_reg !< A pointer that is set to point to the control
74 !! structure for the tracer advection and
75 !! diffusion module
76 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control struct
77
78 ! Local variables
79 character(len=40) :: mdl = "boundary_impulse_tracer" ! This module's name.
80 character(len=48) :: var_name ! The variable's name.
81 character(len=48) :: flux_units ! The units for tracer fluxes, usually
82 ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
83 ! This include declares and sets the variable "version".
84# include "version_variable.h"
85 real, pointer :: tr_ptr(:,:,:) => null() ! The tracer concentration [CU ~> conc]
86 real, pointer :: rem_time_ptr => null() ! The ramaining injection time [T ~> s]
88 integer :: isd, ied, jsd, jed, nz, m
89 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
90
91 if (associated(cs)) then
92 call mom_error(fatal, "register_boundary_impulse_tracer called with an "// &
93 "associated control structure.")
94 endif
95 allocate(cs)
96
97 ! Read all relevant parameters and write them to the model log.
98 call log_version(param_file, mdl, version, "")
99 call get_param(param_file, mdl, "IMPULSE_SOURCE_TIME", cs%remaining_source_time, &
100 "Length of time for the boundary tracer to be injected "//&
101 "into the mixed layer. After this time has elapsed, the "//&
102 "surface becomes a sink for the boundary impulse tracer.", &
103 units="s", default=31536000.0, scale=us%s_to_T)
104 call get_param(param_file, mdl, "TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
105 "If true, tracers may go through the initialization code "//&
106 "if they are not found in the restart files. Otherwise "//&
107 "it is a fatal error if the tracers are not found in the "//&
108 "restart files of a restarted run.", default=.false.)
109 cs%ntr = ntr_max
110 allocate(cs%tr(isd:ied,jsd:jed,nz,cs%ntr), source=0.0)
111
112 cs%nkml = max(gv%nkml,1)
113
114 do m=1,cs%ntr
115 ! This is needed to force the compiler not to do a copy in the registration
116 ! calls. Curses on the designers and implementers of Fortran90.
117 cs%tr_desc(m) = var_desc(trim("boundary_impulse"), "kg kg-1", &
118 "Boundary impulse tracer", caller=mdl)
119 if (gv%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1"
120 else ; flux_units = "kg s-1" ; endif
121
122 tr_ptr => cs%tr(:,:,:,m)
123 call query_vardesc(cs%tr_desc(m), name=var_name, caller="register_boundary_impulse_tracer")
124 ! Register the tracer for horizontal advection, diffusion, and restarts.
125 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, tr_desc=cs%tr_desc(m), &
126 registry_diags=.true., flux_units=flux_units, &
127 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
128
129 ! Set coupled_tracers to be true (hard-coded above) to provide the surface
130 ! values to the coupler (if any). This is meta-code and its arguments will
131 ! currently (deliberately) give fatal errors if it is used.
132 if (cs%coupled_tracers) &
133 cs%ind_tr(m) = atmos_ocn_coupler_flux(trim(var_name)//'_flux', &
134 flux_type=' ', implementation=' ', caller="register_boundary_impulse_tracer")
135 enddo
136 ! Register remaining source time as a restart field
137 rem_time_ptr => cs%remaining_source_time
138 call register_restart_field(rem_time_ptr, "bir_remain_time", &
139 .not.cs%tracers_may_reinit, restart_cs, &
140 "Remaining time to apply BIR source", "s", conversion=us%T_to_s)
141
142 cs%tr_Reg => tr_reg
143 cs%restart_CSp => restart_cs
145
147
148!> Initialize tracer from restart or set to 1 at surface to initialize
149subroutine initialize_boundary_impulse_tracer(restart, day, G, GV, US, h, diag, OBC, CS, &
150 sponge_CSp, tv)
151 logical, intent(in) :: restart !< .true. if the fields have already
152 !! been read from a restart file.
153 type(time_type), target, intent(in) :: day !< Time of the start of the run.
154 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
155 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
156 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
157 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
158 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
159 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
160 !! diagnostic output.
161 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
162 !! whether, where, and what open boundary
163 !! conditions are used.
164 type(boundary_impulse_tracer_cs), pointer :: cs !< The control structure returned by a previous
165 !! call to register_boundary_impulse_tracer.
166 type(sponge_cs), pointer :: sponge_csp !< Pointer to the control structure for the sponges.
167 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
168 !! thermodynamic variables
169 ! Local variables
170 character(len=16) :: name ! A variable's name in a NetCDF file.
171 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
172 integer :: isdb, iedb, jsdb, jedb
173
174 if (.not.associated(cs)) return
175 if (cs%ntr < 1) return
176 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
177 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
178 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
179
180 cs%Time => day
181 cs%diag => diag
182 name = "boundary_impulse"
183
184 do m=1,cs%ntr
185 call query_vardesc(cs%tr_desc(m), name=name, caller="initialize_boundary_impulse_tracer")
186 if ((.not.restart) .or. (.not. query_initialized(cs%tr(:,:,:,m), name, cs%restart_CSp))) then
187 do k=1,cs%nkml ; do j=jsd,jed ; do i=isd,ied
188 cs%tr(i,j,k,m) = 1.0
189 enddo ; enddo ; enddo
190 call set_initialized(cs%tr(:,:,:,m), name, cs%restart_CSp)
191 endif
192 enddo ! Tracer loop
193
194 if (associated(obc)) then
195 ! Steal from updated DOME in the fullness of time.
196 endif
197
199
200!> Apply source or sink at boundary and do vertical diffusion
201subroutine boundary_impulse_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
202 tv, debug, evap_CFL_limit, minimum_forcing_depth)
203 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
204 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
205 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
206 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
207 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
208 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
209 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
210 intent(in) :: ea !< an array to which the amount of fluid entrained
211 !! from the layer above during this call will be
212 !! added [H ~> m or kg m-2].
213 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
214 intent(in) :: eb !< an array to which the amount of fluid entrained
215 !! from the layer below during this call will be
216 !! added [H ~> m or kg m-2].
217 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
218 !! and tracer forcing fields. Unused fields have NULL ptrs.
219 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
220 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
221 type(boundary_impulse_tracer_cs), pointer :: cs !< The control structure returned by a previous
222 !! call to register_boundary_impulse_tracer.
223 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
224 !! thermodynamic variables
225 logical, intent(in) :: debug !< If true calculate checksums
226 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
227 !! be fluxed out of the top layer in a timestep [nondim]
228 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
229 !! fluxes can be applied [H ~> m or kg m-2]
230
231! This subroutine applies diapycnal diffusion and any other column
232! tracer physics or chemistry to the tracers from this file.
233! This is a simple example of a set of advected passive tracers.
234
235! The arguments to this subroutine are redundant in that
236! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
237
238 ! Local variables
239 integer :: i, j, k, is, ie, js, je, nz, m
240 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
241
242 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
243
244 if (.not.associated(cs)) return
245 if (cs%ntr < 1) return
246
247 ! This uses applyTracerBoundaryFluxesInOut, usually in ALE mode
248 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
249 do k=1,nz ;do j=js,je ; do i=is,ie
250 h_work(i,j,k) = h_old(i,j,k)
251 enddo ; enddo ; enddo
252 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,1), dt, fluxes, h_work, &
253 evap_cfl_limit, minimum_forcing_depth)
254 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,1), g, gv)
255 else
256 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,1), g, gv)
257 endif
258
259 ! Set surface conditions
260 do m=1,1
261 if (cs%remaining_source_time > 0.0) then
262 do k=1,cs%nkml ; do j=js,je ; do i=is,ie
263 cs%tr(i,j,k,m) = 1.0
264 enddo ; enddo ; enddo
265 cs%remaining_source_time = cs%remaining_source_time-dt
266 else
267 do k=1,cs%nkml ; do j=js,je ; do i=is,ie
268 cs%tr(i,j,k,m) = 0.0
269 enddo ; enddo ; enddo
270 endif
271
272 enddo
273
275
276!> Calculate total inventory of tracer
277!> This function calculates the mass-weighted integral of the boundary impulse,
278!! tracer stocks returning the number of stocks it has calculated. If the stock_index
279!! is present, only the stock corresponding to that coded index is returned.
280function boundary_impulse_stock(h, stocks, G, GV, CS, names, units, stock_index)
281 type(ocean_grid_type), intent(in ) :: g !< The ocean's grid structure
282 type(verticalgrid_type), intent(in ) :: gv !< The ocean's vertical grid structure
283 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in ) :: h !< Layer thicknesses [H ~> m or kg m-2]
284 type(efp_type), dimension(:), intent( out) :: stocks !< The mass-weighted integrated amount of each
285 !! tracer, in kg times concentration units [kg conc]
286 type(boundary_impulse_tracer_cs), pointer :: cs !< The control structure returned by a previous
287 !! call to register_boundary_impulse_tracer.
288 character(len=*), dimension(:), intent( out) :: names !< The names of the stocks calculated.
289 character(len=*), dimension(:), intent( out) :: units !< The units of the stocks calculated.
290 integer, optional, intent(in ) :: stock_index !< The coded index of a specific stock
291 !! being sought.
292 integer :: boundary_impulse_stock !< Return value: the number of stocks calculated here.
293
294 ! Local variables
295 integer :: m
296
298 if (.not.associated(cs)) return
299 if (cs%ntr < 1) return
300
301 if (present(stock_index)) then ; if (stock_index > 0) then
302 ! Check whether this stock is available from this routine.
303
304 ! No stocks from this routine are being checked yet. Return 0.
305 return
306 endif ; endif
307
308 do m=1,1
309 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller="boundary_impulse_stock")
310 units(m) = trim(units(m))//" kg"
311 stocks(m) = global_mass_int_efp(h, g, gv, cs%tr(:,:,:,m), on_pe_only=.true.)
312 enddo
313
315
316end function boundary_impulse_stock
317
318!> This subroutine extracts the surface fields from this tracer package that
319!! are to be shared with the atmosphere in coupled configurations.
320!! This particular tracer package does not report anything back to the coupler.
321subroutine boundary_impulse_tracer_surface_state(sfc_state, h, G, GV, CS)
322 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
323 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
324 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
325 !! describe the surface state of the ocean.
326 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
327 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
328 type(boundary_impulse_tracer_cs), pointer :: cs !< The control structure returned by a previous
329 !! call to register_boundary_impulse_tracer.
330
331 ! This particular tracer package does not report anything back to the coupler.
332 ! The code that is here is just a rough guide for packages that would.
333
334 integer :: m, is, ie, js, je, isd, ied, jsd, jed
335 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
336 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
337
338 if (.not.associated(cs)) return
339
340 if (cs%coupled_tracers) then
341 do m=1,cs%ntr
342 ! This call loads the surface values into the appropriate array in the
343 ! coupler-type structure.
344 call set_coupler_type_data(cs%tr(:,:,1,m), cs%ind_tr(m), sfc_state%tr_fields, &
345 idim=(/isd, is, ie, ied/), jdim=(/jsd, js, je, jed/), turns=g%HI%turns)
346 enddo
347 endif
348
350
351!> Performs finalization of boundary impulse tracer
352subroutine boundary_impulse_tracer_end(CS)
353 type(boundary_impulse_tracer_cs), pointer :: cs !< The control structure returned by a previous
354 !! call to register_boundary_impulse_tracer.
355 if (associated(cs)) then
356 if (associated(cs%tr)) deallocate(cs%tr)
357 deallocate(cs)
358 endif
359end subroutine boundary_impulse_tracer_end
360
361!> \namespace boundary_impulse_tracer
362!!
363!! \section section_BIT_desc Boundary Impulse Response Tracer and Transit Time Distributions
364!! Transit time distributions (TTD) are the Green's function solution of the passive tracer equation between
365!! the oceanic surface and interior. The name derives from the idea that the 'age' (e.g. time since last
366!! contact with the atmosphere) of a water parcel is best characterized as a distribution of ages
367!! because water parcels leaving the surface arrive at a particular interior point at different times.
368!! The more commonly used ideal age tracer is the first moment of the TTD, equivalently referred to as the
369!! mean age.
370!!
371!! A boundary impulse response (BIR) is a passive tracer whose surface boundary condition is a rectangle
372!! function with width \f$\Delta t\f$. In the case of unsteady flow, multiple BIRs, initiated at different
373!! times in the model can be used to infer the transit time distribution or Green's function between
374!! the oceanic surface and interior. In the case of steady or cyclostationary flow, a single BIR is
375!! sufficient.
376!!
377!! In the References section, both the theoretical discussion of TTDs and BIRs are listed along with
378!! modeling studies which have this used framework in scientific investigations
379!!
380!! \section section_BIT_params Run-time parameters
381!! -DO_BOUNDARY_IMPULSE_TRACER: Enables the boundary impulse tracer model
382!! -IMPULSE_SOURCE_TIME: Length of time that the surface layer acts as a source of the BIR tracer
383!!
384!! \section section_BIT_refs References
385!! \subsection TTD and BIR Theory
386!! -Holzer, M., and T.M. Hall, 2000: Transit-time and tracer-age distributions in geophysical flows.
387!! J. Atmos. Sci., 57, 3539-3558, doi:10.1175/1520-0469(2000)057<3539:TTATAD>2.0.CO;2.
388!! -T.W.N. Haine, H. Zhang, D.W. Waugh, M. Holzer, On transit-time distributions in unsteady circulation
389!! models, Ocean Modelling, Volume 21, Issues 1–2, 2008, Pages 35-45, ISSN 1463-5003
390!! http://dx.doi.org/10.1016/j.ocemod.2007.11.004.
391!! \subsection section_BIT_apps Modelling applications
392!! -Peacock, S., and M. Maltrud (2006), Transit-time distributions in a global ocean model,
393!! J. Phys. Oceanogr., 36(3), 474–495, doi:10.1175/JPO2860.1.
394!! -Maltrud, M., Bryan, F. & Peacock, Boundary impulse response functions in a century-long eddying global
395!! ocean simulation, S. Environ Fluid Mech (2010) 10: 275. doi:10.1007/s10652-009-9154-3
396!!