RGC_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 module contains the routines used to set up a
6!! dynamically passive tracer.
7!! Set up and use passive tracers requires the following:
8!! (1) register_RGC_tracer
9!! (2) apply diffusion, physics/chemistry and advect the tracer
10
11!********+*********+*********+*********+*********+*********+*********+**
12!* *
13!* By Elizabeth Yankovsky, June 2019 *
14!*********+*********+*********+*********+*********+*********+***********
15
16module rgc_tracer
17
19use mom_error_handler, only : mom_error, fatal, warning
20use mom_file_parser, only : get_param, log_param, log_version, param_file_type
21use mom_forcing_type, only : forcing
23use mom_grid, only : ocean_grid_type
24use mom_io, only : file_exists, mom_read_data, slasher, vardesc, var_desc, query_vardesc
26use mom_ale_sponge, only : set_up_ale_sponge_field, ale_sponge_cs, get_ale_sponge_nz_data
28use mom_time_manager, only : time_type
29use mom_tracer_registry, only : register_tracer, tracer_registry_type
32use mom_variables, only : surface
35
36implicit none ; private
37
38#include <MOM_memory.h>
39
40!< Publicly available functions
43
44integer, parameter :: ntr = 1 !< The number of tracers in this module.
45
46!> tracer control structure
47type, public :: rgc_tracer_cs ; private
48 logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
49 character(len = 200) :: tracer_ic_file !< The full path to the IC file, or " " to initialize internally.
50 type(time_type), pointer :: time !< A pointer to the ocean model's clock.
51 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the tracer registry.
52 real, pointer :: tr(:,:,:,:) => null() !< The array of tracers used in this package [kg kg-1]
53 real, pointer :: tr_aux(:,:,:,:) => null() !< The masked tracer concentration [kg kg-1]
54 real :: land_val(ntr) = -1.0 !< The value of tr used where land is masked out [kg kg-1]
55 real :: csl !< The length of the continental shelf (x direction) [km]
56 real :: lensponge !< the length of the sponge layer [km]
57 logical :: mask_tracers !< If true, tracers are masked out in massless layers.
58 logical :: use_sponge !< If true, sponges may be applied somewhere in the domain.
59 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the timing of diagnostic output.
60 type(vardesc) :: tr_desc(ntr) !< Descriptions and metadata for the tracers.
61end type rgc_tracer_cs
62
63contains
64
65!> This subroutine is used to register tracer fields
66function register_rgc_tracer(G, GV, param_file, CS, tr_Reg, restart_CS)
67 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
68 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
69 type(param_file_type), intent(in) :: param_file !<A structure indicating the open file to parse
70 !! for model parameter values.
71 type(rgc_tracer_cs), pointer :: cs !< A pointer that is set to point to the control
72 !! structure for this module (in/out).
73 type(tracer_registry_type), pointer :: tr_reg !< A pointer to the tracer registry.
74 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
75
76 character(len=80) :: name, longname
77 ! This include declares and sets the variable "version".
78# include "version_variable.h"
79 character(len=40) :: mdl = "RGC_tracer" ! This module's name.
80 character(len=200) :: inputdir
81 real, pointer :: tr_ptr(:,:,:) => null() ! A pointer to one of the tracers in this module [kg kg-1]
82 logical :: register_rgc_tracer
83 integer :: isd, ied, jsd, jed, nz, m
84 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
85
86 if (associated(cs)) then
87 call mom_error(fatal, "RGC_register_tracer called with an "// &
88 "associated control structure.")
89 endif
90 allocate(cs)
91
92 ! Read all relevant parameters and write them to the model log.
93 call log_version(param_file, mdl, version, "")
94 call get_param(param_file, mdl, "RGC_TRACER_IC_FILE", cs%tracer_IC_file, &
95 "The name of a file from which to read the initial \n"//&
96 "conditions for the RGC tracers, or blank to initialize \n"//&
97 "them internally.", default=" ")
98 if (len_trim(cs%tracer_IC_file) >= 1) then
99 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
100 inputdir = slasher(inputdir)
101 cs%tracer_IC_file = trim(inputdir)//trim(cs%tracer_IC_file)
102 call log_param(param_file, mdl, "INPUTDIR/RGC_TRACER_IC_FILE", &
103 cs%tracer_IC_file)
104 endif
105 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
106 "If true, sponges may be applied anywhere in the domain. \n"//&
107 "The exact location and properties of those sponges are \n"//&
108 "specified from MOM_initialization.F90.", default=.false.)
109
110 call get_param(param_file, mdl, "CONT_SHELF_LENGTH", cs%CSL, &
111 "The length of the continental shelf (x dir, km).", &
112 units=g%x_ax_unit_short, default=15.0)
113
114 call get_param(param_file, mdl, "LENSPONGE", cs%lensponge, &
115 "The length of the sponge layer (km).", &
116 units=g%x_ax_unit_short, default=10.0)
117
118 allocate(cs%tr(isd:ied,jsd:jed,nz,ntr), source=0.0)
119 if (cs%mask_tracers) then
120 allocate(cs%tr_aux(isd:ied,jsd:jed,nz,ntr), source=0.0)
121 endif
122
123 do m=1,ntr
124 write(name,'("tr_RGC",I0)') m
125 write(longname,'("Concentration of RGC Tracer ",I2.2)') m
126 cs%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mdl)
127
128 ! This is needed to force the compiler not to do a copy in the registration calls.
129 tr_ptr => cs%tr(:,:,:,m)
130 ! Register the tracer for horizontal advection & diffusion.
131 call register_tracer(tr_ptr, tr_reg, param_file, g%HI, gv, &
132 name=name, longname=longname, units="kg kg-1", &
133 registry_diags=.true., flux_units="kg/s", &
134 restart_cs=restart_cs)
135 enddo
136
137 cs%tr_Reg => tr_reg
138 register_rgc_tracer = .true.
139end function register_rgc_tracer
140
141!> Initializes the NTR tracer fields in tr(:,:,:,:)
142!! and it sets up the tracer output.
143subroutine initialize_rgc_tracer(restart, day, G, GV, h, diag, OBC, CS, &
144 layer_CSp, sponge_CSp)
145
146 type(ocean_grid_type), intent(in) :: g !< Grid structure.
147 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
148 logical, intent(in) :: restart !< .true. if the fields have already
149 !! been read from a restart file.
150 type(time_type), target, intent(in) :: day !< Time of the start of the run.
151 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
152 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
153 type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output.
154 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
155 !! whether, where, and what open boundary
156 !! conditions are used. This is not being used for now.
157 type(rgc_tracer_cs), pointer :: cs !< The control structure returned by a previous
158 !! call to RGC_register_tracer.
159 type(sponge_cs), pointer :: layer_csp !< A pointer to the control structure
160 type(ale_sponge_cs), pointer :: sponge_csp !< A pointer to the control structure for the
161 !! sponges, if they are in use. Otherwise this may be unassociated.
162
163 real, allocatable :: temp(:,:,:) ! A temporary array used for several sponge target values [various]
164 character(len=16) :: name ! A variable's name in a NetCDF file.
165 real, pointer :: tr_ptr(:,:,:) => null() ! A pointer to one of the tracers in this module [kg kg-1]
166 real :: h_neglect ! A thickness that is so small it is usually lost
167 ! in roundoff and can be neglected [H ~> m or kg m-2].
168 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
169 integer :: isdb, iedb, jsdb, jedb
170 integer :: nzdata
171
172 if (.not.associated(cs)) return
173 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
174 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
175 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
176 h_neglect = gv%H_subroundoff
177
178 cs%Time => day
179 cs%diag => diag
180
181 if (.not.restart) then
182 if (len_trim(cs%tracer_IC_file) >= 1) then
183 ! Read the tracer concentrations from a netcdf file.
184 if (.not.file_exists(cs%tracer_IC_file, g%Domain)) &
185 call mom_error(fatal, "RGC_initialize_tracer: Unable to open "// &
186 cs%tracer_IC_file)
187 do m=1,ntr
188 call query_vardesc(cs%tr_desc(m), name, caller="initialize_RGC_tracer")
189 call mom_read_data(cs%tracer_IC_file, trim(name), cs%tr(:,:,:,m), g%Domain)
190 enddo
191 else
192 do m=1,ntr
193 do k=1,nz ; do j=js,je ; do i=is,ie
194 cs%tr(i,j,k,m) = 0.0
195 enddo ; enddo ; enddo
196 enddo
197 m=1
198 do j=js,je ; do i=is,ie
199 !set tracer to 1.0 in the surface of the continental shelf
200 if (g%geoLonT(i,j) <= (cs%CSL)) then
201 cs%tr(i,j,1,m) = 1.0 !first layer
202 endif
203 enddo ; enddo
204
205 endif
206 endif ! restart
207
208 if ( cs%use_sponge ) then
209! If sponges are used, this damps values to zero in the offshore boundary.
210! For any tracers that are not damped in the sponge, the call
211! to set_up_sponge_field can simply be omitted.
212 if (associated(sponge_csp)) then !ALE mode
213 nzdata = get_ale_sponge_nz_data(sponge_csp)
214 if (nzdata>0) then
215 allocate(temp(g%isd:g%ied,g%jsd:g%jed,nzdata))
216 do k=1,nzdata ; do j=js,je ; do i=is,ie
217 if (g%geoLonT(i,j) >= (g%len_lon - cs%lensponge) .AND. g%geoLonT(i,j) <= g%len_lon) then
218 temp(i,j,k) = 0.0
219 endif
220 enddo ; enddo ; enddo
221 do m=1,1
222 ! This is needed to force the compiler not to do a copy in the sponge calls.
223 tr_ptr => cs%tr(:,:,:,m)
224 call set_up_ale_sponge_field(temp, g, gv, tr_ptr, sponge_csp, 'RGC_tracer')
225 enddo
226 deallocate(temp)
227 endif
228
229 elseif (associated(layer_csp)) then !layer mode
230 if (nz>0) then
231 allocate(temp(g%isd:g%ied,g%jsd:g%jed,nz))
232 do k=1,nz ; do j=js,je ; do i=is,ie
233 if (g%geoLonT(i,j) >= (g%len_lon - cs%lensponge) .AND. g%geoLonT(i,j) <= g%len_lon) then
234 temp(i,j,k) = 0.0
235 endif
236 enddo ; enddo ; enddo
237 do m=1,1
238 tr_ptr => cs%tr(:,:,:,m)
239 call set_up_sponge_field(temp, tr_ptr, g, gv, nz, layer_csp)
240 enddo
241 deallocate(temp)
242 endif
243 else
244 call mom_error(fatal, "RGC_initialize_tracer: "// &
245 "The pointer to sponge_CSp must be associated if SPONGE is defined.")
246 endif !selecting mode/calling error if no pointer
247 endif !using sponge
248
249end subroutine initialize_rgc_tracer
250
251!> This subroutine applies diapycnal diffusion and any other column
252!! tracer physics or chemistry to the tracers from this file.
253!! This is a simple example of a set of advected passive tracers.
254subroutine rgc_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
255 evap_CFL_limit, minimum_forcing_depth)
256 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
257 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
258 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
259 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
260 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
261 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
262 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
263 intent(in) :: ea !< an array to which the amount of fluid entrained
264 !! from the layer above during this call will be
265 !! added [H ~> m or kg m-2].
266 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
267 intent(in) :: eb !< an array to which the amount of fluid entrained
268 !! from the layer below during this call will be
269 !! added [H ~> m or kg m-2].
270 type(forcing), intent(in) :: fluxes !< A structure containing pointers to any possible
271 !! forcing fields. Unused fields have NULL ptrs.
272 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
273 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
274 type(rgc_tracer_cs), pointer :: cs !< The control structure returned by a previous call.
275 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can be
276 !! fluxed out of the top layer in a timestep [nondim].
277 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which fluxes
278 !! can be applied [H ~> m or kg m-2].
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
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
285 integer :: i, j, k, is, ie, js, je, nz, m
286 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
287
288 if (.not.associated(cs)) return
289
290 m=1
291 do j=js,je ; do i=is,ie
292 ! set tracer to 1.0 in the surface of the continental shelf
293 if (g%geoLonT(i,j) <= (cs%CSL)) then
294 cs%tr(i,j,1,m) = 1.0 !first layer
295 endif
296 enddo ; enddo
297
298 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
299 do m=1,ntr
300 do k=1,nz ; do j=js,je ; do i=is,ie
301 h_work(i,j,k) = h_old(i,j,k)
302 enddo ; enddo ; enddo
303 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,m) , dt, fluxes, h_work, &
304 evap_cfl_limit, minimum_forcing_depth)
305
306 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
307 enddo
308 else
309 do m=1,ntr
310 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
311 enddo
312 endif
313
314end subroutine rgc_tracer_column_physics
315
316subroutine rgc_tracer_end(CS)
317 type(rgc_tracer_cs), pointer :: cs !< The control structure returned by a previous call to RGC_register_tracer.
318
319 if (associated(cs)) then
320 if (associated(cs%tr)) deallocate(cs%tr)
321 deallocate(cs)
322 endif
323end subroutine rgc_tracer_end
324
325end module rgc_tracer