nw2_tracers.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!> Ideal tracers designed to help diagnose a tracer diffusivity tensor in NeverWorld2
6module nw2_tracers
7
9use mom_error_handler, only : mom_error, fatal, warning
10use mom_file_parser, only : get_param, log_param, log_version, param_file_type
11use mom_forcing_type, only : forcing
12use mom_grid, only : ocean_grid_type
14use mom_interface_heights, only : thickness_to_dz
15use mom_io, only : file_exists, mom_read_data, slasher, vardesc, var_desc
16use mom_restart, only : query_initialized, set_initialized, mom_restart_cs
17use mom_time_manager, only : time_type
18use mom_tracer_registry, only : register_tracer, tracer_registry_type
23
24implicit none ; private
25
26#include <MOM_memory.h>
27
31public nw2_tracers_end
32
33!> The control structure for the nw2_tracers package
34type, public :: nw2_tracers_cs ; private
35 integer :: ntr = 0 !< The number of tracers that are actually used.
36 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
37 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the tracer registry
38 real, pointer :: tr(:,:,:,:) => null() !< The array of tracers used in this package, in [conc] (g m-3)?
39 real, allocatable , dimension(:) :: restore_rate !< The rate at which the tracer is damped toward
40 !! its target profile [T-1 ~> s-1]
41 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
42 !! regulate the timing of diagnostic output.
43 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart controls structure
44end type nw2_tracers_cs
45
46contains
47
48!> Register the NW2 tracer fields to be used with MOM.
49logical function register_nw2_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
50 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure
51 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
52 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
53 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
54 type(nw2_tracers_cs), pointer :: cs !< The control structure returned by a previous
55 !! call to register_nw2_tracer.
56 type(tracer_registry_type), pointer :: tr_reg !< A pointer that is set to point to the control
57 !! structure for the tracer advection and
58 !! diffusion module
59 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control struct
60
61 ! This include declares and sets the variable "version".
62# include "version_variable.h"
63 character(len=40) :: mdl = "nw2_tracers" ! This module's name.
64 character(len=8) :: var_name ! The variable's name.
65 real, pointer :: tr_ptr(:,:,:) => null() ! The tracer concentration [conc]
66 integer :: isd, ied, jsd, jed, nz, m, ig
67 integer :: n_groups ! Number of groups of three tracers (i.e. # tracers/3)
68 real, allocatable, dimension(:) :: timescale_in_days ! Damping timescale [days]
69 type(vardesc) :: tr_desc ! Descriptions and metadata for the tracers
70 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
71
72 if (associated(cs)) then
73 call mom_error(fatal, "register_nw2_tracer called with an "// &
74 "associated control structure.")
75 endif
76 allocate(cs)
77
78 ! Read all relevant parameters and write them to the model log.
79 call log_version(param_file, mdl, version, "")
80 call get_param(param_file, mdl, "NW2_TRACER_GROUPS", n_groups, &
81 "The number of tracer groups where a group is of three tracers "//&
82 "initialized and restored to sin(x), y and z, respectively. Each "//&
83 "group is restored with an independent restoration rate.", &
84 default=3)
85 allocate(timescale_in_days(n_groups))
86 timescale_in_days = (/365., 730., 1460./)
87 call get_param(param_file, mdl, "NW2_TRACER_RESTORE_TIMESCALE", timescale_in_days, &
88 "A list of timescales, one for each tracer group.", &
89 units="days")
90
91 cs%ntr = 3 * n_groups
92 allocate(cs%tr(isd:ied,jsd:jed,nz,cs%ntr), source=0.0)
93 allocate(cs%restore_rate(cs%ntr))
94
95 do m=1,cs%ntr
96 write(var_name(1:8),'(a6,i2.2)') 'tracer',m
97 tr_desc = var_desc(var_name, "1", "Ideal Tracer", caller=mdl)
98 ! This is needed to force the compiler not to do a copy in the registration
99 ! calls. Curses on the designers and implementers of Fortran90.
100 tr_ptr => cs%tr(:,:,:,m)
101 ! Register the tracer for horizontal advection, diffusion, and restarts.
102 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, tr_desc=tr_desc, &
103 registry_diags=.true., restart_cs=restart_cs, mandatory=.false.)
104 ig = int( (m+2)/3 ) ! maps (1,2,3)->1, (4,5,6)->2, ...
105 cs%restore_rate(m) = 1.0 / ( timescale_in_days(ig) * 86400.0*us%s_to_T )
106 enddo
107
108 cs%tr_Reg => tr_reg
109 cs%restart_CSp => restart_cs
110 register_nw2_tracers = .true.
111end function register_nw2_tracers
112
113!> Sets the NW2 traces to their initial values and sets up the tracer output
114subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS)
115 logical, intent(in) :: restart !< .true. if the fields have already
116 !! been read from a restart file.
117 type(time_type), target, intent(in) :: day !< Time of the start of the run.
118 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
119 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
120 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
121 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
122 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
123 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
124 !! thermodynamic variables
125 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
126 !! diagnostic output.
127 type(nw2_tracers_cs), pointer :: cs !< The control structure returned by a previous
128 !! call to register_nw2_tracer.
129 ! Local variables
130 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Interface heights [Z ~> m]
131 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Vertical extent of layers [Z ~> m]
132 real :: rscl ! z* scaling factor [nondim]
133 character(len=8) :: var_name ! The variable's name.
134 integer :: i, j, k, m
135
136 if (.not.associated(cs)) return
137
138 cs%Time => day
139 cs%diag => diag
140
141 ! Calculate z* interface positions
142 call thickness_to_dz(h, tv, dz, g, gv, us)
143
144 if (gv%Boussinesq) then
145 ! First calculate interface positions in z-space (m)
146 do j=g%jsc,g%jec ; do i=g%isc,g%iec
147 eta(i,j,gv%ke+1) = - g%mask2dT(i,j) * g%bathyT(i,j)
148 enddo ; enddo
149 do k=gv%ke,1,-1 ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
150 eta(i,j,k) = eta(i,j,k+1) + g%mask2dT(i,j) * dz(i,j,k)
151 enddo ; enddo ; enddo
152 ! Re-calculate for interface positions in z*-space (m)
153 do j=g%jsc,g%jec ; do i=g%isc,g%iec
154 if (g%bathyT(i,j)>0.) then
155 rscl = g%bathyT(i,j) / ( eta(i,j,1) + g%bathyT(i,j) )
156 do k=gv%ke, 1, -1
157 eta(i,j,k) = eta(i,j,k+1) + g%mask2dT(i,j) * dz(i,j,k) * rscl
158 enddo
159 endif
160 enddo ; enddo
161 else
162 call mom_error(fatal, "NW2 tracers assume Boussinesq mode")
163 endif
164
165 do m=1,cs%ntr
166 ! Initialize only if this is not a restart or we are using a restart
167 ! in which the tracers were not present
168 write(var_name(1:8),'(a6,i2.2)') 'tracer',m
169 if ((.not.restart) .or. &
170 (.not. query_initialized(cs%tr(:,:,:,m), var_name, cs%restart_CSp))) then
171 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
172 cs%tr(i,j,k,m) = nw2_tracer_dist(m, g, gv, eta, i, j, k)
173 enddo ; enddo ; enddo
174 call set_initialized(cs%tr(:,:,:,m), var_name, cs%restart_CSp)
175 endif ! restart
176 enddo ! Tracer loop
177
178end subroutine initialize_nw2_tracers
179
180!> Applies diapycnal diffusion, aging and regeneration at the surface to the NW2 tracers
181subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, tv, CS, &
182 evap_CFL_limit, minimum_forcing_depth)
183 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
184 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
185 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
186 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
187 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
188 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
189 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
190 intent(in) :: ea !< an array to which the amount of fluid entrained
191 !! from the layer above during this call will be
192 !! added [H ~> m or kg m-2].
193 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
194 intent(in) :: eb !< an array to which the amount of fluid entrained
195 !! from the layer below during this call will be
196 !! added [H ~> m or kg m-2].
197 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
198 !! and tracer forcing fields. Unused fields have NULL ptrs.
199 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
200 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
201 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
202 type(nw2_tracers_cs), pointer :: cs !< The control structure returned by a previous
203 !! call to register_nw2_tracer.
204 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
205 !! be fluxed out of the top layer in a timestep [nondim]
206 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
207 !! fluxes can be applied [H ~> m or kg m-2]
208! This subroutine applies diapycnal diffusion and any other column
209! tracer physics or chemistry to the tracers from this file.
210! This is a simple example of a set of advected passive tracers.
211
212! The arguments to this subroutine are redundant in that
213! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
214 ! Local variables
215 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
216 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Interface heights [Z ~> m]
217 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Vertical extent of layers [Z ~> m]
218 integer :: i, j, k, m
219 real :: dt_x_rate ! dt * restoring rate [nondim]
220 real :: rscl ! z* scaling factor [nondim]
221 real :: target_value ! tracer target value for damping [conc]
222
223! if (.not.associated(CS)) return
224
225 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
226 do m=1,cs%ntr
227 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
228 h_work(i,j,k) = h_old(i,j,k)
229 enddo ; enddo ; enddo
230 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,m), dt, fluxes, h_work, &
231 evap_cfl_limit, minimum_forcing_depth)
232 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
233 enddo
234 else
235 do m=1,cs%ntr
236 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
237 enddo
238 endif
239
240 ! Calculate z* interface positions
241 call thickness_to_dz(h_new, tv, dz, g, gv, us)
242
243 if (gv%Boussinesq) then
244 ! First calculate interface positions in z-space [Z ~> m]
245 do j=g%jsc,g%jec ; do i=g%isc,g%iec
246 eta(i,j,gv%ke+1) = - g%mask2dT(i,j) * g%bathyT(i,j)
247 enddo ; enddo
248 do k=gv%ke,1,-1 ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
249 eta(i,j,k) = eta(i,j,k+1) + g%mask2dT(i,j) * dz(i,j,k)
250 enddo ; enddo ; enddo
251 ! Re-calculate for interface positions in z*-space [Z ~> m]
252 do j=g%jsc,g%jec ; do i=g%isc,g%iec
253 if (g%bathyT(i,j)>0.) then
254 rscl = g%bathyT(i,j) / ( eta(i,j,1) + g%bathyT(i,j) )
255 do k=gv%ke, 1, -1
256 eta(i,j,k) = eta(i,j,k+1) + g%mask2dT(i,j) * dz(i,j,k) * rscl
257 enddo
258 endif
259 enddo ; enddo
260 else
261 call mom_error(fatal, "NW2 tracers assume Boussinesq mode")
262 endif
263
264 do m=1,cs%ntr
265 dt_x_rate = dt * cs%restore_rate(m)
266 !$OMP parallel do default(shared) private(target_value)
267 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
268 target_value = nw2_tracer_dist(m, g, gv, eta, i, j, k)
269 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + g%mask2dT(i,j) * dt_x_rate * ( target_value - cs%tr(i,j,k,m) )
270 enddo ; enddo ; enddo
271 enddo
272
273end subroutine nw2_tracer_column_physics
274
275!> The target value of a NeverWorld2 tracer label m [conc] at non-dimensional
276!! position x=lon/Lx, y=lat/Ly, z=eta/H
277real function nw2_tracer_dist(m, G, GV, eta, i, j, k)
278 integer, intent(in) :: m !< Indicates the NW2 tracer
279 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
280 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
281 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
282 intent(in) :: eta !< Interface position [Z ~> m]
283 integer, intent(in) :: i !< Cell index i
284 integer, intent(in) :: j !< Cell index j
285 integer, intent(in) :: k !< Layer index k
286 ! Local variables
287 real :: pi ! 3.1415... [nondim]
288 real :: x, y, z ! non-dimensional relative positions [nondim]
289 pi = 2.*acos(0.)
290 x = ( g%geolonT(i,j) - g%west_lon ) / g%len_lon ! 0 ... 1
291 y = -g%geolatT(i,j) / g%south_lat ! -1 ... 1
292 z = - 0.5 * ( eta(i,j,k) + eta(i,j,k+1) ) / gv%max_depth ! 0 ... 1
293 select case ( mod(m-1,3) )
294 case (0) ! sin(2 pi x/L)
295 nw2_tracer_dist = sin( 2.0 * pi * x )
296 case (1) ! y/L
298 case (2) ! -z/L
299 nw2_tracer_dist = -z
300 case default
301 stop 'This should not happen. Died in nw2_tracer_dist()!'
302 end select
303 nw2_tracer_dist = nw2_tracer_dist * g%mask2dT(i,j)
304end function nw2_tracer_dist
305
306!> Deallocate any memory associated with this tracer package
307subroutine nw2_tracers_end(CS)
308 type(nw2_tracers_cs), pointer :: cs !< The control structure returned by a previous
309 !! call to register_nw2_tracers.
310
311 if (associated(cs)) then
312 if (associated(cs%tr)) deallocate(cs%tr)
313 deallocate(cs)
314 endif
315end subroutine nw2_tracers_end
316
317!> \namespace nw2_tracers
318!!
319!! TBD
320
321end module nw2_tracers