pseudo_salt_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 mimics salinity
7
8use mom_coms, only : efp_type
9use mom_debugging, only : hchksum
10use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
12use mom_error_handler, only : mom_error, fatal, warning
13use mom_file_parser, only : get_param, log_param, log_version, param_file_type
14use mom_forcing_type, only : forcing
15use mom_grid, only : ocean_grid_type
20use mom_restart, only : query_initialized, set_initialized, mom_restart_cs
23use mom_time_manager, only : time_type
24use mom_tracer_registry, only : register_tracer, tracer_registry_type, tracer_type
30
31implicit none ; private
32
33#include <MOM_memory.h>
34
38
39!> The control structure for the pseudo-salt tracer
40type, public :: pseudo_salt_tracer_cs ; private
41 type(tracer_type), pointer :: tr_ptr !< pointer to tracer inside Tr_reg
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 MOM tracer registry
44 real, pointer :: ps(:,:,:) => null() !< The array of pseudo-salt tracer used in this
45 !! subroutine [ppt]
46 real, allocatable :: diff(:,:,:) !< The difference between the pseudo-salt
47 !! tracer and the real salt [ppt].
48 logical :: pseudo_salt_may_reinit = .true. !< Hard coding since this should not matter
49
50 integer :: id_psd = -1 !< A diagnostic ID
51
52 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate
53 !! the timing of diagnostic output.
54 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control structure
55
56 type(vardesc) :: tr_desc !< A description and metadata for the pseudo-salt tracer
58
59contains
60
61!> Register the pseudo-salt tracer with MOM6, and return .true. if the tracer is to be used.
62function register_pseudo_salt_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
63 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure
64 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
65 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
66 type(pseudo_salt_tracer_cs), pointer :: cs !< The control structure returned by a previous
67 !! call to register_pseudo_salt_tracer.
68 type(tracer_registry_type), pointer :: tr_reg !< A pointer that is set to point to the control
69 !! structure for the tracer advection and
70 !! diffusion module
71 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control structure
72
73 ! Local variables
74 character(len=40) :: mdl = "pseudo_salt_tracer" ! This module's name.
75 character(len=48) :: var_name ! The variable's name.
76 ! This include declares and sets the variable "version".
77# include "version_variable.h"
78 real, pointer :: tr_ptr(:,:,:) => null() ! The tracer concentration [ppt]
80 integer :: isd, ied, jsd, jed, nz
81 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
82
83 if (associated(cs)) then
84 call mom_error(fatal, "register_pseudo_salt_tracer called with an "// &
85 "associated control structure.")
86 endif
87 allocate(cs)
88
89 ! Read all relevant parameters and write them to the model log.
90 call log_version(param_file, mdl, version, "")
91
92 allocate(cs%ps(isd:ied,jsd:jed,nz), source=0.0)
93
94 cs%tr_desc = var_desc(trim("pseudo_salt"), "psu", &
95 "Pseudo salt passive tracer", caller=mdl)
96
97 tr_ptr => cs%ps(:,:,:)
98 call query_vardesc(cs%tr_desc, name=var_name, caller="register_pseudo_salt_tracer")
99 ! Register the tracer for horizontal advection, diffusion, and restarts.
100 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, name="pseudo_salt", &
101 longname="Pseudo salt passive tracer", units="psu", &
102 registry_diags=.true., restart_cs=restart_cs, &
103 mandatory=.not.cs%pseudo_salt_may_reinit, tr_out=cs%tr_ptr)
104
105 cs%tr_Reg => tr_reg
106 cs%restart_CSp => restart_cs
108
110
111!> Initialize the pseudo-salt tracer
112subroutine initialize_pseudo_salt_tracer(restart, day, G, GV, US, h, diag, OBC, CS, &
113 sponge_CSp, tv)
114 logical, intent(in) :: restart !< .true. if the fields have already
115 !! been read from a restart file.
116 type(time_type), target, intent(in) :: day !< Time of the start of the run
117 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
118 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
119 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
120 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
121 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
122 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
123 !! diagnostic output
124 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
125 !! whether, where, and what open boundary
126 !! conditions are used.
127 type(pseudo_salt_tracer_cs), pointer :: cs !< The control structure returned by a previous
128 !! call to register_pseudo_salt_tracer
129 type(sponge_cs), pointer :: sponge_csp !< Pointer to the control structure for the sponges
130 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing various thermodynamic variables
131
132 ! This subroutine initializes the tracer fields in CS%ps(:,:,:).
133
134 ! Local variables
135 character(len=16) :: name ! A variable's name in a NetCDF file
136 integer :: i, j, k, isd, ied, jsd, jed, nz
137
138 if (.not.associated(cs)) return
139
140 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
141
142 cs%Time => day
143 cs%diag => diag
144 name = "pseudo_salt"
145
146 call query_vardesc(cs%tr_desc, name=name, caller="initialize_pseudo_salt_tracer")
147 if ((.not.restart) .or. (.not.query_initialized(cs%ps, name, cs%restart_CSp))) then
148 do k=1,nz ; do j=jsd,jed ; do i=isd,ied
149 cs%ps(i,j,k) = us%S_to_ppt*tv%S(i,j,k)
150 enddo ; enddo ; enddo
151 call set_initialized(cs%ps, name, cs%restart_CSp)
152 endif
153
154 if (associated(obc)) then
155 ! Steal from updated DOME in the fullness of time.
156 endif
157
158 cs%id_psd = register_diag_field("ocean_model", "pseudo_salt_diff", cs%diag%axesTL, &
159 day, "Difference between pseudo salt passive tracer and salt tracer", "psu")
160 if (.not.allocated(cs%diff)) allocate(cs%diff(isd:ied,jsd:jed,nz), source=0.0)
161
163
164!> Apply sources, sinks and diapycnal diffusion to the tracers in this package.
165subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, tv, debug, &
166 KPP_CSp, nonLocalTrans, evap_CFL_limit, minimum_forcing_depth)
167 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
168 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
169 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
170 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]
171 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
172 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]
173 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
174 intent(in) :: ea !< The amount of fluid entrained from the layer above
175 !! during this call [H ~> m or kg m-2]
176 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
177 intent(in) :: eb !< The amount of fluid entrained from the layer below
178 !! during this call [H ~> m or kg m-2]
179 type(forcing), intent(in) :: fluxes !< A structure containing thermodynamic and
180 !! tracer forcing fields
181 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
182 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
183 type(pseudo_salt_tracer_cs), pointer :: cs !< The control structure returned by a previous
184 !! call to register_pseudo_salt_tracer
185 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
186 logical, intent(in) :: debug !< If true calculate checksums
187 type(kpp_cs), optional, pointer :: kpp_csp !< KPP control structure
188 real, optional, intent(in) :: nonlocaltrans(:,:,:) !< Non-local transport [nondim]
189 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
190 !! be fluxed out of the top layer in a timestep [nondim]
191 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
192 !! fluxes can be applied [H ~> m or kg m-2]
193
194 ! This subroutine applies diapycnal diffusion and any other column
195 ! tracer physics or chemistry to the tracers from this file.
196
197 ! The arguments to this subroutine are redundant in that
198 ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
199
200 ! Local variables
201 real :: net_salt_rate(szi_(g),szj_(g)) ! Net salt flux into the ocean
202 ! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
203 real :: net_salt(szi_(g),szj_(g)) ! Net salt flux into the ocean integrated over
204 ! a timestep [ppt H ~> ppt m or ppt kg m-2]
205 real :: htot(szi_(g)) ! Total ocean depth [H ~> m or kg m-2]
206 real :: fluxrescaledepth ! Minimum total ocean depth at which fluxes start to be scaled
207 ! away [H ~> m or kg m-2]
208 real :: ih_limit ! Inverse of FluxRescaleDepth or 0 for no limiting [H-1 ~> m-1 or m2 kg-1]
209 real :: scale ! Scale scales away fluxes if depth < FluxRescaleDepth [nondim]
210 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
211 integer :: i, j, k, is, ie, js, je, nz
212
213 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
214
215 if (.not.associated(cs)) return
216 if (.not.associated(cs%ps)) return
217
218 if (debug) then
219 call hchksum(tv%S,"salt pre pseudo-salt vertdiff", g%HI, unscale=us%S_to_ppt)
220 call hchksum(cs%ps,"pseudo_salt pre pseudo-salt vertdiff", g%HI)
221 endif
222
223 fluxrescaledepth = max( gv%Angstrom_H, 1.e-30*gv%m_to_H )
224 ih_limit = 0.0 ; if (fluxrescaledepth > 0.0) ih_limit = 1.0 / fluxrescaledepth
225
226 ! Compute KPP nonlocal term if necessary
227 if (present(kpp_csp)) then
228 if (associated(kpp_csp) .and. present(nonlocaltrans)) then
229 ! Determine the salt flux, including limiting for small total ocean depths.
230 net_salt_rate(:,:) = 0.0
231 if (associated(fluxes%salt_flux)) then
232 do j=js,je
233 do i=is,ie ; htot(i) = h_old(i,j,1) ; enddo
234 do k=2,nz ; do i=is,ie ; htot(i) = htot(i) + h_old(i,j,k) ; enddo ; enddo
235 do i=is,ie
236 scale = 1.0 ; if ((ih_limit > 0.0) .and. (htot(i)*ih_limit < 1.0)) scale = htot(i)*ih_limit
237 net_salt_rate(i,j) = (scale * (1000.0 * fluxes%salt_flux(i,j))) * gv%RZ_to_H
238 enddo
239 enddo
240 endif
241 call kpp_nonlocaltransport(kpp_csp, g, gv, h_old, nonlocaltrans, net_salt_rate, &
242 dt, cs%diag, cs%tr_ptr, cs%ps(:,:,:))
243 endif
244 endif
245
246 ! This uses applyTracerBoundaryFluxesInOut, usually in ALE mode
247 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
248 ! This option uses applyTracerBoundaryFluxesInOut, usually in ALE mode
249
250 ! Determine the time-integrated salt flux, including limiting for small total ocean depths.
251 net_salt(:,:) = 0.0
252 do j=js,je
253 do i=is,ie ; htot(i) = h_old(i,j,1) ; enddo
254 do k=2,nz ; do i=is,ie ; htot(i) = htot(i) + h_old(i,j,k) ; enddo ; enddo
255 do i=is,ie
256 scale = 1.0 ; if ((ih_limit > 0.0) .and. (htot(i)*ih_limit < 1.0)) scale = htot(i)*ih_limit
257 net_salt(i,j) = (scale * dt * (1000.0 * fluxes%salt_flux(i,j))) * gv%RZ_to_H
258 enddo
259 enddo
260
261 do k=1,nz ; do j=js,je ; do i=is,ie
262 h_work(i,j,k) = h_old(i,j,k)
263 enddo ; enddo ; enddo
264 call applytracerboundaryfluxesinout(g, gv, cs%ps, dt, fluxes, h_work, evap_cfl_limit, &
265 minimum_forcing_depth, out_flux_optional=net_salt)
266 call tracer_vertdiff(h_work, ea, eb, dt, cs%ps, g, gv)
267 else
268 call tracer_vertdiff(h_old, ea, eb, dt, cs%ps, g, gv)
269 endif
270
271 if (debug) then
272 call hchksum(tv%S, "salt post pseudo-salt vertdiff", g%HI, unscale=us%S_to_ppt)
273 call hchksum(cs%ps, "pseudo_salt post pseudo-salt vertdiff", g%HI)
274 endif
275
276 if (allocated(cs%diff)) then
277 do k=1,nz ; do j=js,je ; do i=is,ie
278 cs%diff(i,j,k) = cs%ps(i,j,k) - us%S_to_ppt*tv%S(i,j,k)
279 enddo ; enddo ; enddo
280 if (cs%id_psd>0) call post_data(cs%id_psd, cs%diff, cs%diag)
281 endif
282
284
285
286!> Calculates the mass-weighted integral of all tracer stocks, returning the number of stocks it has
287!! calculated. If the stock_index is present, only the stock corresponding to that coded index is returned.
288function pseudo_salt_stock(h, stocks, G, GV, CS, names, units, stock_index)
289 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
290 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
291 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
292 type(efp_type), dimension(:), intent(out) :: stocks !< The mass-weighted integrated amount of each
293 !! tracer, in kg times concentration units [kg conc]
294 type(pseudo_salt_tracer_cs), pointer :: cs !< The control structure returned by a previous
295 !! call to register_pseudo_salt_tracer
296 character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated
297 character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated
298 integer, optional, intent(in) :: stock_index !< The coded index of a specific stock
299 !! being sought
300 integer :: pseudo_salt_stock !< Return value: the number of
301 !! stocks calculated here
302
303
305 if (.not.associated(cs)) return
306 if (.not.allocated(cs%diff)) return
307
308 if (present(stock_index)) then ; if (stock_index > 0) then
309 ! Check whether this stock is available from this routine.
310
311 ! No stocks from this routine are being checked yet. Return 0.
312 return
313 endif ; endif
314
315 call query_vardesc(cs%tr_desc, name=names(1), units=units(1), caller="pseudo_salt_stock")
316 units(1) = trim(units(1))//" kg"
317 stocks(1) = global_mass_int_efp(h, g, gv, cs%diff, on_pe_only=.true.)
318
320
321end function pseudo_salt_stock
322
323!> This subroutine extracts the surface fields from this tracer package that
324!! are to be shared with the atmosphere in coupled configurations.
325!! This particular tracer package does not report anything back to the coupler.
326subroutine pseudo_salt_tracer_surface_state(sfc_state, h, G, GV, CS)
327 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
328 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
329 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
330 !! describe the surface state of the ocean
331 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
332 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
333 type(pseudo_salt_tracer_cs), pointer :: cs !< The control structure returned by a previous
334 !! call to register_pseudo_salt_tracer
335
336 ! This particular tracer package does not report anything back to the coupler.
337 ! The code that is here is just a rough guide for packages that would.
338
339 if (.not.associated(cs)) return
340
341 ! By design, this tracer package does not return any surface states.
342
344
345!> Deallocate memory associated with this tracer package
346subroutine pseudo_salt_tracer_end(CS)
347 type(pseudo_salt_tracer_cs), pointer :: cs !< The control structure returned by a previous
348 !! call to register_pseudo_salt_tracer
349
350 if (associated(cs)) then
351 if (associated(cs%ps)) deallocate(cs%ps)
352 if (allocated(cs%diff)) deallocate(cs%diff)
353 deallocate(cs)
354 endif
355end subroutine pseudo_salt_tracer_end
356
357!> \namespace pseudo_salt_tracer
358!!
359!! By Andrew Shao, 2016
360!!
361!! This file contains the routines necessary to model a passive
362!! tracer that uses the same boundary fluxes as salinity. At the
363!! beginning of the run, salt is set to the same as tv%S. Any
364!! deviations between this salt-like tracer and tv%S signifies a
365!! difference between how active and passive tracers are treated.
366
367end module pseudo_salt_tracer