ISOMIP_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!> Routines used to set up and use a set of (one for now)
6!! dynamically passive tracers in the ISOMIP configuration.
7!!
8!! For now, just one passive tracer is injected in
9!! the sponge layer.
10module isomip_tracer
11
12! Original sample tracer package by Robert Hallberg, 2002
13! Adapted to the ISOMIP test case by Gustavo Marques, May 2016
14
15use mom_coms, only : max_across_pes
16use mom_coupler_types, only : set_coupler_type_data, atmos_ocn_coupler_flux
18use mom_error_handler, only : mom_error, fatal, warning
19use mom_file_parser, only : get_param, log_param, log_version, param_file_type
20use mom_forcing_type, only : forcing
22use mom_grid, only : ocean_grid_type
23use 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
27use mom_time_manager, only : time_type
28use mom_tracer_registry, only : register_tracer, tracer_registry_type
31use mom_variables, only : surface
33
34implicit none ; private
35
36#include <MOM_memory.h>
37
38!< Publicly available functions
41
42integer, parameter :: ntr = 1 !< ntr is the number of tracers in this module.
43
44!> ISOMIP tracer package control structure
45type, public :: isomip_tracer_cs ; private
46 logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
47 character(len = 200) :: tracer_ic_file !< The full path to the IC file, or " " to initialize internally.
48 type(time_type), pointer :: time !< A pointer to the ocean model's clock.
49 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the MOM tracer registry
50 real, pointer :: tr(:,:,:,:) => null() !< The array of tracers used in this package, in [conc] (g m-3)?
51 real :: land_val(ntr) = -1.0 !< The value of tr used where land is masked out [conc].
52 logical :: use_sponge !< If true, sponges may be applied somewhere in the domain.
53
54 integer, dimension(NTR) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux
55 !< if it is used and the surface tracer concentrations are to be
56 !< provided to the coupler.
57
58 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
59 !! timing of diagnostic output.
60
61 type(vardesc) :: tr_desc(ntr) !< Descriptions and metadata for the tracers in this package
62end type isomip_tracer_cs
63
64contains
65
66
67!> This subroutine is used to register tracer fields
68function register_isomip_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
69 type(hor_index_type), intent(in) :: hi !<A horizontal index type structure.
70 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
71 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
72 !! to parse for model parameter values.
73 type(isomip_tracer_cs), pointer :: cs !<A pointer that is set to point to the control
74 !! structure for this module (in/out).
75 type(tracer_registry_type), pointer :: tr_reg !<A pointer to the tracer registry.
76 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control struct
77
78 character(len=80) :: name, longname
79 ! This include declares and sets the variable "version".
80# include "version_variable.h"
81 character(len=40) :: mdl = "ISOMIP_tracer" ! This module's name.
82 character(len=200) :: inputdir
83 character(len=48) :: flux_units ! The units for tracer fluxes, usually
84 ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
85 real, pointer :: tr_ptr(:,:,:) => null() ! The tracer concentration [conc]
86 logical :: register_isomip_tracer
87 integer :: isd, ied, jsd, jed, nz, m
88 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
89
90 if (associated(cs)) then
91 call mom_error(fatal, "ISOMIP_register_tracer called with an "// &
92 "associated control structure.")
93 endif
94 allocate(cs)
95
96 ! Read all relevant parameters and write them to the model log.
97 call log_version(param_file, mdl, version, "")
98 call get_param(param_file, mdl, "ISOMIP_TRACER_IC_FILE", cs%tracer_IC_file, &
99 "The name of a file from which to read the initial "//&
100 "conditions for the ISOMIP tracers, or blank to initialize "//&
101 "them internally.", default=" ")
102 if (len_trim(cs%tracer_IC_file) >= 1) then
103 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
104 inputdir = slasher(inputdir)
105 cs%tracer_IC_file = trim(inputdir)//trim(cs%tracer_IC_file)
106 call log_param(param_file, mdl, "INPUTDIR/ISOMIP_TRACER_IC_FILE", &
107 cs%tracer_IC_file)
108 endif
109 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
110 "If true, sponges may be applied anywhere in the domain. "//&
111 "The exact location and properties of those sponges are "//&
112 "specified from MOM_initialization.F90.", default=.false.)
113
114 allocate(cs%tr(isd:ied,jsd:jed,nz,ntr), source=0.0)
115
116 do m=1,ntr
117 write(name,'("tr_D",I0)') m
118 write(longname,'("Concentration of ISOMIP Tracer ",I2.2)') m
119 cs%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mdl)
120 if (gv%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1"
121 else ; flux_units = "kg s-1" ; endif
122
123 ! This is needed to force the compiler not to do a copy in the registration
124 ! calls. Curses on the designers and implementers of Fortran90.
125 tr_ptr => cs%tr(:,:,:,m)
126 ! Register the tracer for horizontal advection, diffusion, and restarts.
127 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
128 name=name, longname=longname, units="kg kg-1", &
129 registry_diags=.true., flux_units=flux_units, &
130 restart_cs=restart_cs)
131
132 ! Set coupled_tracers to be true (hard-coded above) to provide the surface
133 ! values to the coupler (if any). This is meta-code and its arguments will
134 ! currently (deliberately) give fatal errors if it is used.
135 if (cs%coupled_tracers) &
136 cs%ind_tr(m) = atmos_ocn_coupler_flux(trim(name)//'_flux', &
137 flux_type=' ', implementation=' ', caller="register_ISOMIP_tracer")
138 enddo
139
140 cs%tr_Reg => tr_reg
142end function register_isomip_tracer
143
144!> Initializes the NTR tracer fields in tr(:,:,:,:)
145! and it sets up the tracer output.
146subroutine initialize_isomip_tracer(restart, day, G, GV, h, diag, OBC, CS, &
147 ALE_sponge_CSp)
148
149 type(ocean_grid_type), intent(in) :: g !< Grid structure.
150 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
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 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
155 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
156 !! diagnostic output.
157 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
158 !! whether, where, and what open boundary conditions
159 !! are used. This is not being used for now.
160 type(isomip_tracer_cs), pointer :: cs !< The control structure returned by a previous call
161 !! to ISOMIP_register_tracer.
162 type(ale_sponge_cs), pointer :: ale_sponge_csp !< A pointer to the control structure for
163 !! the sponges, if they are in use. Otherwise this
164 !! may be unassociated.
165
166 character(len=16) :: name ! A variable's name in a NetCDF file.
167 real :: h_neglect ! A thickness that is so small it is usually lost
168 ! in roundoff and can be neglected [H ~> m or kg m-2].
169 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
170 integer :: isdb, iedb, jsdb, jedb
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, "ISOMIP_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_ISOMIP_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 endif
198 endif ! restart
199
200! the following does not work in layer mode yet
201!! if ( CS%use_sponge ) then
202 ! If sponges are used, this example damps tracers in sponges in the
203 ! northern half of the domain to 1 and tracers in the southern half
204 ! to 0. For any tracers that are not damped in the sponge, the call
205 ! to set_up_sponge_field can simply be omitted.
206! if (.not.associated(ALE_sponge_CSp)) &
207! call MOM_error(FATAL, "ISOMIP_initialize_tracer: "// &
208! "The pointer to ALEsponge_CSp must be associated if SPONGE is defined.")
209
210! allocate(temp(G%isd:G%ied,G%jsd:G%jed,nz))
211
212! do j=js,je ; do i=is,ie
213! if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then
214! temp(i,j,:) = 1.0
215! else
216! temp(i,j,:) = 0.0
217! endif
218! enddo ; enddo
219
220 ! do m=1,NTR
221! do m=1,1
222 ! This is needed to force the compiler not to do a copy in the sponge
223 ! calls. Curses on the designers and implementers of Fortran90.
224! tr_ptr => CS%tr(:,:,:,m)
225! call set_up_ALE_sponge_field(temp, G, tr_ptr, ALE_sponge_CSp)
226! enddo
227! deallocate(temp)
228! endif
229
230end subroutine initialize_isomip_tracer
231
232!> This subroutine applies diapycnal diffusion, including the surface boundary
233!! conditions and any other column tracer physics or chemistry to the tracers from this file.
234subroutine isomip_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
235 evap_CFL_limit, minimum_forcing_depth)
236 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
237 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
238 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
239 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
240 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
241 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
242 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
243 intent(in) :: ea !< an array to which the amount of fluid entrained
244 !! from the layer above during this call will be
245 !! added [H ~> m or kg m-2].
246 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
247 intent(in) :: eb !< an array to which the amount of fluid entrained
248 !! from the layer below during this call will be
249 !! added [H ~> m or kg m-2].
250 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
251 !! and tracer forcing fields. Unused fields have NULL ptrs.
252 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
253 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
254 type(isomip_tracer_cs), pointer :: cs !< The control structure returned by a previous
255 !! call to ISOMIP_register_tracer.
256 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
257 !! be fluxed out of the top layer in a timestep [nondim]
258 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
259 !! fluxes can be applied [H ~> m or kg m-2]
260
261! The arguments to this subroutine are redundant in that
262! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
263
264 ! Local variables
265 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
266 real :: melt(szi_(g),szj_(g)) ! melt water (positive for melting, negative for freezing) [R Z T-1 ~> kg m-2 s-1]
267 real :: mmax ! The global maximum melting rate [R Z T-1 ~> kg m-2 s-1]
268 integer :: i, j, k, is, ie, js, je, nz, m
269 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
270
271 if (.not.associated(cs)) return
272
273 melt(:,:) = fluxes%iceshelf_melt(:,:)
274
275 ! max. melt
276 mmax = maxval(melt(is:ie,js:je))
277 call max_across_pes(mmax)
278 ! write(mesg,*) 'max melt = ', mmax
279 ! call MOM_mesg(mesg, 5)
280 ! dye melt water (m=1), dye = 1 if melt=max(melt)
281 do m=1,ntr
282 do j=js,je ; do i=is,ie
283 if (melt(i,j) > 0.0) then ! melting
284 cs%tr(i,j,1:2,m) = melt(i,j)/mmax ! inject dye in the ML
285 else ! freezing
286 cs%tr(i,j,1:2,m) = 0.0
287 endif
288 enddo ; enddo
289 enddo
290
291 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
292 do m=1,ntr
293 do k=1,nz ;do j=js,je ; do i=is,ie
294 h_work(i,j,k) = h_old(i,j,k)
295 enddo ; enddo ; enddo
296 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,m), dt, fluxes, h_work, &
297 evap_cfl_limit, minimum_forcing_depth)
298 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
299 enddo
300 else
301 do m=1,ntr
302 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
303 enddo
304 endif
305
306end subroutine isomip_tracer_column_physics
307
308!> This subroutine extracts the surface fields from this tracer package that
309!! are to be shared with the atmosphere in coupled configurations.
310!! This particular tracer package does not report anything back to the coupler.
311subroutine isomip_tracer_surface_state(sfc_state, h, G, GV, CS)
312 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
313 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
314 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
315 !! describe the surface state of the ocean.
316 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
317 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
318 type(isomip_tracer_cs), pointer :: cs !< The control structure returned by a previous
319 !! call to ISOMIP_register_tracer.
320
321 ! This particular tracer package does not report anything back to the coupler.
322 ! The code that is here is just a rough guide for packages that would.
323
324 integer :: m, is, ie, js, je, isd, ied, jsd, jed
325 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
326 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
327
328 if (.not.associated(cs)) return
329
330 if (cs%coupled_tracers) then
331 do m=1,ntr
332 ! This call loads the surface values into the appropriate array in the
333 ! coupler-type structure.
334 call set_coupler_type_data(cs%tr(:,:,1,m), cs%ind_tr(m), sfc_state%tr_fields, &
335 idim=(/isd, is, ie, ied/), jdim=(/jsd, js, je, jed/), turns=g%HI%turns)
336 enddo
337 endif
338
339end subroutine isomip_tracer_surface_state
340
341!> Deallocate any memory used by the ISOMIP tracer package
342subroutine isomip_tracer_end(CS)
343 type(isomip_tracer_cs), pointer :: cs !< The control structure returned by a previous
344 !! call to ISOMIP_register_tracer.
345
346 if (associated(cs)) then
347 if (associated(cs%tr)) deallocate(cs%tr)
348 deallocate(cs)
349 endif
350end subroutine isomip_tracer_end
351
352end module isomip_tracer