MOM_ALE_sponge.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 apply sponge layers when using
6!! the ALE mode.
7!!
8!! Applying sponges requires the following:
9!! 1. initialize_ALE_sponge
10!! 2. set_up_ALE_sponge_field (tracers) and set_up_ALE_sponge_vel_field (vel)
11!! 3. apply_ALE_sponge
12!! 4. init_ALE_sponge_diags (not being used for now)
13!! 5. ALE_sponge_end (not being used for now)
14
15module mom_ale_sponge
16
17
18use mom_array_transform, only: rotate_array
19use mom_coms, only : sum_across_pes
22use mom_domains, only : pass_var, to_all, omit_corners
23use mom_error_handler, only : mom_error, fatal, note, warning, is_root_pe
24use mom_file_parser, only : get_param, log_param, log_version, param_file_type
25use mom_grid, only : ocean_grid_type
26use mom_horizontal_regridding, only : horiz_interp_and_extrap_tracer
27use mom_interpolate, only : init_external_field, time_interp_external_init
29use mom_interpolate, only : external_field
30use mom_io, only : axis_info
33use mom_time_manager, only : time_type
37
38implicit none ; private
39
40#include <MOM_memory.h>
41
42!> Store the reference profile at h points for a variable
43interface set_up_ale_sponge_field
44 module procedure set_up_ale_sponge_field_fixed
45 module procedure set_up_ale_sponge_field_varying
46end interface
47
48!> This subroutine stores the reference profile at u and v points for a vector
49interface set_up_ale_sponge_vel_field
50 module procedure set_up_ale_sponge_vel_field_fixed
51 module procedure set_up_ale_sponge_vel_field_varying
52end interface
53
54!> Determine the number of points which are within sponges in this computational domain.
55!!
56!! Only points that have positive values of Iresttime and which mask2dT indicates are ocean
57!! points are included in the sponges. It also stores the target interface heights.
58interface initialize_ale_sponge
59 module procedure initialize_ale_sponge_fixed
60 module procedure initialize_ale_sponge_varying
61end interface
62
63! Publicly available functions
64public set_up_ale_sponge_field, set_up_ale_sponge_vel_field
65public get_ale_sponge_thicknesses, get_ale_sponge_nz_data
66public initialize_ale_sponge, apply_ale_sponge, ale_sponge_end, init_ale_sponge_diags
67public rotate_ale_sponge, update_ale_sponge_field
68
69! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
70! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
71! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
72! vary with the Boussinesq approximation, the Boussinesq variant is given first.
73
74!> A structure for creating arrays of pointers to 3D arrays with extra gridding information
75type :: p3d ; private
76 !integer :: id !< id for FMS external time interpolator
77 integer :: nz_data !< The number of vertical levels in the input field.
78 integer :: num_tlevs !< The number of time records contained in the file
79 real, dimension(:,:,:), pointer :: p => null() !< pointer to the data [various]
80 real, dimension(:,:,:), pointer :: dz => null() !< pointer to the data grid spacing [Z ~> m]
81end type p3d
82
83!> A structure for creating arrays of pointers to 2D arrays with extra gridding information
84type :: p2d ; private
85 type(external_field) :: field !< Time interpolator field handle
86 integer :: nz_data !< The number of vertical levels in the input field
87 integer :: num_tlevs !< The number of time records contained in the file
88 real :: scale = 1.0 !< A multiplicative factor by which to rescale input data [various]
89 real, dimension(:,:), pointer :: p => null() !< pointer to the data [various]
90 real, dimension(:,:), pointer :: dz => null() !< pointer to the data grid spacing [Z ~> m]
91 character(len=:), allocatable :: name !< The name of the input field
92 character(len=:), allocatable :: long_name !< The long name of the input field
93 character(len=:), allocatable :: unit !< The unit of the input field
94 type(axis_info), allocatable :: axes_data(:) !< Axis types for the input field
95end type p2d
96
97!> ALE sponge control structure
98type, public :: ale_sponge_cs ; private
99 integer :: nz !< The total number of layers.
100 integer :: nz_data !< The total number of arbitrary layers (used by older code).
101 integer :: num_col !< The number of sponge tracer points within the computational domain.
102 integer :: num_col_u !< The number of sponge u-points within the computational domain.
103 integer :: num_col_v !< The number of sponge v-points within the computational domain.
104 integer :: fldno = 0 !< The number of fields which have already been
105 !! registered by calls to set_up_sponge_field
106 logical :: sponge_uv !< Control whether u and v are included in sponge
107 integer, allocatable :: col_i(:) !< Array of the i-indices of each tracer column being damped
108 integer, allocatable :: col_j(:) !< Array of the j-indices of each tracer column being damped
109 integer, allocatable :: col_i_u(:) !< Array of the i-indices of each u-column being damped
110 integer, allocatable :: col_j_u(:) !< Array of the j-indices of each u-column being damped
111 integer, allocatable :: col_i_v(:) !< Array of the i-indices of each v-column being damped
112 integer, allocatable :: col_j_v(:) !< Array of the j-indices of each v-column being damped
113
114 real, allocatable :: iresttime_col(:) !< The inverse restoring time of each tracer column [T-1 ~> s-1]
115 real, allocatable :: iresttime_col_u(:) !< The inverse restoring time of each u-column [T-1 ~> s-1]
116 real, allocatable :: iresttime_col_v(:) !< The inverse restoring time of each v-column [T-1 ~> s-1]
117
118 type(p3d) :: var(max_fields_) !< Pointers to the fields that are being damped.
119 type(p2d) :: ref_val(max_fields_) !< The values to which the fields are damped.
120 type(p2d) :: ref_val_u !< The values to which the u-velocities are damped.
121 type(p2d) :: ref_val_v !< The values to which the v-velocities are damped.
122 type(p3d) :: var_u !< Pointer to the u velocities that are being damped.
123 type(p3d) :: var_v !< Pointer to the v velocities that are being damped.
124 type(p2d) :: ref_dz !< Grid on which reference data is provided (older code).
125 type(p2d) :: ref_dzu !< u-point grid on which reference data is provided (older code).
126 type(p2d) :: ref_dzv !< v-point grid on which reference data is provided (older code).
127
128 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
129 !! timing of diagnostic output.
130
131 type(remapping_cs) :: remap_cs !< Remapping parameters and work arrays
132 integer :: remap_answer_date !< The vintage of the order of arithmetic and expressions to use
133 !! for remapping. Values below 20190101 recover the remapping
134 !! answers from 2018, while higher values use more robust
135 !! forms of the same remapping expressions.
136 integer :: hor_regrid_answer_date !< The vintage of the order of arithmetic and expressions to use
137 !! for horizontal regridding. Values below 20190101 recover the
138 !! answers from 2018, while higher values use expressions that have
139 !! been rearranged for rotational invariance.
140
141 logical :: time_varying_sponges !< True if using newer sponge code
142 logical :: spongedataongrid !< True if the sponge data are on the model horizontal grid
143 real :: varying_input_dz_mask !< An input file thickness below which the target values with time-varying
144 !! sponges are replaced by the value above [Z ~> m].
145 !! It is not clear why this needs to be greater than 0.
146
147 !>@{ Diagnostic IDs
148 integer, dimension(MAX_FIELDS_) :: id_sp_tendency = reshape([-1], [max_fields_], [-1]) !< Diagnostic ids for tracer
149 !! tendencies due to sponges.
150 !! Init all to -1.
151 integer :: id_sp_u_tendency !< Diagnostic id for zonal momentum tendency due to
152 !! Rayleigh damping
153 integer :: id_sp_v_tendency !< Diagnostic id for meridional momentum tendency due to
154 !! Rayleigh damping
155end type ale_sponge_cs
156
157contains
158
159!> This subroutine determines the number of points which are within sponges in this computational
160!! domain. Only points that have positive values of Iresttime and which mask2dT indicates are ocean
161!! points are included in the sponges. It also stores the target interface heights.
162subroutine initialize_ale_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, nz_data, &
163 Iresttime_u_in, Iresttime_v_in, data_h_is_Z)
164
165 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
166 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
167 integer, intent(in) :: nz_data !< The total number of sponge input layers.
168 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: Iresttime !< The inverse of the restoring time [T-1 ~> s-1].
169 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
170 !! to parse for model parameter values.
171 type(ale_sponge_cs), pointer :: CS !< A pointer that is set to point to the control
172 !! structure for this module (in/out).
173 real, dimension(SZI_(G),SZJ_(G),nz_data), intent(inout) :: data_h !< The thicknesses of the sponge
174 !! input layers, in [H ~> m or kg m-2] or [Z ~> m]
175 !! depending on data_h_is_Z.
176 real, dimension(SZIB_(G),SZJ_(G)), optional, intent(in) :: Iresttime_u_in !< The inverse of the restoring
177 !! time at U-points [T-1 ~> s-1].
178 real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: Iresttime_v_in !< The inverse of the restoring
179 !! time at v-points [T-1 ~> s-1].
180 logical, optional, intent(in) :: data_h_is_Z !< If present and true data_h is already in
181 !! depth units. Omitting this is the same as setting
182 !! it to false.
183
184 ! Local variables
185 character(len=40) :: mdl = "MOM_sponge" ! This module's name.
186 real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1]
187 real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1]
188 real, dimension(SZI_(G),SZJ_(G),nz_data) :: data_dz !< The vertical extent of the sponge
189 !! input layers [Z ~> m].
190 real :: data_h_to_Z_scale ! A scaling factor to convert data_h into the right units, often [Z H-1 ~> 1 or m3 kg-1]
191 ! This include declares and sets the variable "version".
192# include "version_variable.h"
193 character(len=64) :: remapScheme
194 logical :: use_sponge
195 logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
196 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
197 logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
198 integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v
199
200 if (associated(cs)) then
201 call mom_error(warning, "initialize_ALE_sponge_fixed called with an associated "// &
202 "control structure.")
203 return
204 endif
205
206! Set default, read and log parameters
207 call log_version(param_file, mdl, version, "")
208 call get_param(param_file, mdl, "SPONGE", use_sponge, &
209 "If true, sponges may be applied anywhere in the domain. "//&
210 "The exact location and properties of those sponges are "//&
211 "specified from MOM_initialization.F90.", default=.false.)
212
213 if (.not.use_sponge) return
214
215 allocate(cs)
216
217 call get_param(param_file, mdl, "SPONGE_UV", cs%sponge_uv, &
218 "Apply sponges in u and v, in addition to tracers.", &
219 default=.false.)
220
221 call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
222 default="PLM", do_not_log=.true.)
223 call get_param(param_file, mdl, "SPONGE_REMAPPING_SCHEME", remapscheme, &
224 "This sets the reconstruction scheme used "//&
225 "for vertical remapping for all SPONGE variables.", default=remapscheme)
226
227 !This default should be from REMAP_BOUNDARY_EXTRAP
228 call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
229 default=.false., do_not_log=.true.)
230 call get_param(param_file, mdl, "SPONGE_BOUNDARY_EXTRAP", bndextrapolation, &
231 "If true, values at the interfaces of SPONGE boundary cells are "//&
232 "extrapolated instead of piecewise constant", default=bndextrapolation)
233 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
234 "This sets the default value for the various _ANSWER_DATE parameters.", &
235 default=99991231)
236 call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", cs%remap_answer_date, &
237 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
238 "Values below 20190101 result in the use of older, less accurate expressions "//&
239 "that were in use at the end of 2018. Higher values result in the use of more "//&
240 "robust and accurate forms of mathematically equivalent expressions.", &
241 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
242 if (.not.gv%Boussinesq) cs%remap_answer_date = max(cs%remap_answer_date, 20230701)
243 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
244 do_not_log=.true., default=.true.)
245 call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
246 "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//&
247 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
248 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
249
250 call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", cs%hor_regrid_answer_date, &
251 "The vintage of the order of arithmetic for horizontal regridding. "//&
252 "Dates before 20190101 give the same answers as the code did in late 2018, "//&
253 "while later versions add parentheses for rotational symmetry. "//&
254 "Dates after 20230101 use reproducing sums for global averages.", &
255 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
256 if (.not.gv%Boussinesq) cs%hor_regrid_answer_date = max(cs%hor_regrid_answer_date, 20230701)
257
258 cs%time_varying_sponges = .false.
259 cs%nz = gv%ke
260
261 data_h_to_z_scale = gv%H_to_Z ; if (present(data_h_is_z)) data_h_to_z_scale = 1.0
262
263 do k=1,nz_data ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
264 data_dz(i,j,k) = data_h_to_z_scale * data_h(i,j,k)
265 enddo ; enddo ; enddo
266 ! number of columns to be restored
267 cs%num_col = 0 ; cs%fldno = 0
268 do j=g%jsc,g%jec ; do i=g%isc,g%iec
269 if ((iresttime(i,j) > 0.0) .and. (g%mask2dT(i,j) > 0.0)) &
270 cs%num_col = cs%num_col + 1
271 enddo ; enddo
272
273 if (cs%num_col > 0) then
274 allocate(cs%Iresttime_col(cs%num_col), source=0.0)
275 allocate(cs%col_i(cs%num_col), source=0)
276 allocate(cs%col_j(cs%num_col), source=0)
277 ! pass indices, restoring time to the CS structure
278 col = 1
279 do j=g%jsc,g%jec ; do i=g%isc,g%iec
280 if ((iresttime(i,j) > 0.0) .and. (g%mask2dT(i,j) > 0.0)) then
281 cs%col_i(col) = i ; cs%col_j(col) = j
282 cs%Iresttime_col(col) = iresttime(i,j)
283 col = col + 1
284 endif
285 enddo ; enddo
286 ! same for total number of arbitrary layers and correspondent data
287 cs%nz_data = nz_data
288 allocate(cs%Ref_dz%p(cs%nz_data,cs%num_col), source=0.0)
289 do col=1,cs%num_col ; do k=1,cs%nz_data
290 cs%Ref_dz%p(k,col) = data_dz(cs%col_i(col),cs%col_j(col),k)
291 enddo ; enddo
292 endif
293
294 total_sponge_cols = cs%num_col
295 call sum_across_pes(total_sponge_cols)
296
297! Call the constructor for remapping control structure
298 call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation, &
299 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
300 answer_date=cs%remap_answer_date)
301
302 call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, &
303 "The total number of columns where sponges are applied at h points.", like_default=.true.)
304
305 if (cs%sponge_uv) then
306 allocate(iresttime_u(g%isdB:g%iedB,g%jsd:g%jed), source=0.0)
307 allocate(iresttime_v(g%isd:g%ied,g%jsdB:g%jedB), source=0.0)
308
309 call pass_var(iresttime, g%Domain, to_all+omit_corners, halo=1)
310 call pass_var(data_dz, g%Domain, to_all+omit_corners, halo=1)
311
312 ! u points
313 cs%num_col_u = 0
314 if (present(iresttime_u_in)) then
315 iresttime_u(:,:) = iresttime_u_in(:,:)
316 else
317 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
318 iresttime_u(i,j) = 0.5 * (iresttime(i,j) + iresttime(i+1,j))
319 enddo ; enddo
320 endif
321 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
322 if ((iresttime_u(i,j) > 0.0) .and. (g%mask2dCu(i,j) > 0.0)) &
323 cs%num_col_u = cs%num_col_u + 1
324 enddo ; enddo
325
326 if (cs%num_col_u > 0) then
327
328 allocate(cs%Iresttime_col_u(cs%num_col_u), source=0.0)
329 allocate(cs%col_i_u(cs%num_col_u), source=0)
330 allocate(cs%col_j_u(cs%num_col_u), source=0)
331
332 ! Store the column indices and restoring rates in the CS structure
333 col = 1
334 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
335 if ((iresttime_u(i,j) > 0.0) .and. (g%mask2dCu(i,j) > 0.0)) then
336 cs%col_i_u(col) = i ; cs%col_j_u(col) = j
337 cs%Iresttime_col_u(col) = iresttime_u(i,j)
338 col = col + 1
339 endif
340 enddo ; enddo
341
342 ! same for total number of arbitrary layers and correspondent data
343 allocate(cs%Ref_dzu%p(cs%nz_data,cs%num_col_u), source=0.0)
344 do col=1,cs%num_col_u
345 i = cs%col_i_u(col) ; j = cs%col_j_u(col)
346 do k=1,cs%nz_data
347 cs%Ref_dzu%p(k,col) = 0.5 * (data_dz(i,j,k) + data_dz(i+1,j,k))
348 enddo
349 enddo
350 endif
351 total_sponge_cols_u = cs%num_col_u
352 call sum_across_pes(total_sponge_cols_u)
353 call log_param(param_file, mdl, "!Total sponge columns at u points", total_sponge_cols_u, &
354 "The total number of columns where sponges are applied at u points.", like_default=.true.)
355
356 ! v points
357 cs%num_col_v = 0
358 if (present(iresttime_v_in)) then
359 iresttime_v(:,:) = iresttime_v_in(:,:)
360 else
361 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
362 iresttime_v(i,j) = 0.5 * (iresttime(i,j) + iresttime(i,j+1))
363 enddo ; enddo
364 endif
365 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
366 if ((iresttime_v(i,j) > 0.0) .and. (g%mask2dCv(i,j) > 0.0)) &
367 cs%num_col_v = cs%num_col_v + 1
368 enddo ; enddo
369
370 if (cs%num_col_v > 0) then
371
372 allocate(cs%Iresttime_col_v(cs%num_col_v), source=0.0)
373 allocate(cs%col_i_v(cs%num_col_v), source=0)
374 allocate(cs%col_j_v(cs%num_col_v), source=0)
375
376 ! pass indices, restoring time to the CS structure
377 col = 1
378 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
379 if ((iresttime_v(i,j) > 0.0) .and. (g%mask2dCv(i,j) > 0.0)) then
380 cs%col_i_v(col) = i ; cs%col_j_v(col) = j
381 cs%Iresttime_col_v(col) = iresttime_v(i,j)
382 col = col + 1
383 endif
384 enddo ; enddo
385
386 ! same for total number of arbitrary layers and correspondent data
387 allocate(cs%Ref_dzv%p(cs%nz_data,cs%num_col_v), source=0.0)
388 do col=1,cs%num_col_v
389 i = cs%col_i_v(col) ; j = cs%col_j_v(col)
390 do k=1,cs%nz_data
391 cs%Ref_dzv%p(k,col) = 0.5 * (data_dz(i,j,k) + data_dz(i,j+1,k))
392 enddo
393 enddo
394 endif
395 total_sponge_cols_v = cs%num_col_v
396 call sum_across_pes(total_sponge_cols_v)
397 call log_param(param_file, mdl, "!Total sponge columns at v points", total_sponge_cols_v, &
398 "The total number of columns where sponges are applied at v points.", like_default=.true.)
399 endif
400
401end subroutine initialize_ale_sponge_fixed
402
403!> Return the number of layers in the data with a fixed ALE sponge, or 0 if there are
404!! no sponge columns on this PE.
405function get_ale_sponge_nz_data(CS)
406 type(ale_sponge_cs), intent(in) :: cs !< ALE sponge control struct
407 integer :: get_ale_sponge_nz_data !< The number of layers in the fixed sponge data.
408
409 get_ale_sponge_nz_data = cs%nz_data
410end function get_ale_sponge_nz_data
411
412!> Return the thicknesses used for the data with a fixed ALE sponge
413subroutine get_ale_sponge_thicknesses(G, GV, data_h, sponge_mask, CS, data_h_in_Z)
414 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure (in).
415 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
416 real, allocatable, dimension(:,:,:), &
417 intent(inout) :: data_h !< The thicknesses of the sponge input layers expressed
418 !! as vertical extents [Z ~> m] or in thickness units
419 !! [H ~> m or kg m-2], depending on the value of data_h_in_Z.
420 logical, dimension(SZI_(G),SZJ_(G)), &
421 intent(out) :: sponge_mask !< A logical mask that is true where
422 !! sponges are being applied.
423 type(ale_sponge_cs), pointer :: cs !< A pointer that is set to point to the control
424 !! structure for the ALE_sponge module.
425 logical, optional, intent(in) :: data_h_in_z !< If present and true data_h is returned in
426 !! depth units. Omitting this is the same as setting
427 !! it to false.
428 real :: z_to_data_h_units ! A scaling factor to return data_h in the right units, often [H Z-1 ~> 1 or kg m-3]
429 integer :: c, i, j, k
430
431 if (allocated(data_h)) call mom_error(fatal, &
432 "get_ALE_sponge_thicknesses called with an allocated data_h.")
433
434 if (.not.associated(cs)) then
435 ! There are no sponge points on this PE.
436 allocate(data_h(g%isd:g%ied,g%jsd:g%jed,1), source=-1.0)
437 sponge_mask(:,:) = .false.
438 return
439 endif
440
441 allocate(data_h(g%isd:g%ied,g%jsd:g%jed,cs%nz_data), source=-1.0)
442 sponge_mask(:,:) = .false.
443
444 z_to_data_h_units = gv%Z_to_H ; if (present(data_h_in_z)) z_to_data_h_units = 1.0
445
446 do c=1,cs%num_col
447 i = cs%col_i(c) ; j = cs%col_j(c)
448 sponge_mask(i,j) = .true.
449 do k=1,cs%nz_data
450 data_h(i,j,k) = z_to_data_h_units*cs%Ref_dz%p(k,c)
451 enddo
452 enddo
453
454end subroutine get_ale_sponge_thicknesses
455
456!> This subroutine determines the number of points which are to be restored in the computational
457!! domain. Only points that have positive values of Iresttime and which mask2dT indicates are ocean
458!! points are included in the sponges.
459subroutine initialize_ale_sponge_varying(Iresttime, G, GV, US, param_file, CS, Iresttime_u_in, Iresttime_v_in)
460
461 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
462 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
463 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
464 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: Iresttime !< The inverse of the restoring time [T-1 ~> s-1].
465 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to parse
466 !! for model parameter values.
467 type(ale_sponge_cs), pointer :: CS !< A pointer that is set to point to the control
468 !! structure for this module (in/out).
469 real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: Iresttime_u_in !< The inverse of the restoring time
470 !! for u [T-1 ~> s-1].
471 real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: Iresttime_v_in !< The inverse of the restoring time
472 !! for v [T-1 ~> s-1].
473
474 ! Local variables
475 character(len=40) :: mdl = "MOM_sponge" ! This module's name.
476 real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1]
477 real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1]
478 real :: dz_neglect, dz_neglect_edge ! Negligible layer extents [Z ~> m]
479 ! This include declares and sets the variable "version".
480# include "version_variable.h"
481 character(len=64) :: remapScheme
482 logical :: use_sponge
483 logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
484 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
485 logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
486 integer :: i, j, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v
487
488 if (associated(cs)) then
489 call mom_error(warning, "initialize_ALE_sponge_varying called with an associated "// &
490 "control structure.")
491 return
492 endif
493! Set default, read and log parameters
494 call log_version(param_file, mdl, version, "")
495 call get_param(param_file, mdl, "SPONGE", use_sponge, &
496 "If true, sponges may be applied anywhere in the domain. "//&
497 "The exact location and properties of those sponges are "//&
498 "specified from MOM_initialization.F90.", default=.false.)
499 if (.not.use_sponge) return
500 allocate(cs)
501 call get_param(param_file, mdl, "SPONGE_UV", cs%sponge_uv, &
502 "Apply sponges in u and v, in addition to tracers.", &
503 default=.false.)
504 call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
505 default="PLM", do_not_log=.true.)
506 call get_param(param_file, mdl, "SPONGE_REMAPPING_SCHEME", remapscheme, &
507 "This sets the reconstruction scheme used "//&
508 "for vertical remapping for all SPONGE variables.", default=remapscheme)
509 !This default should be from REMAP_BOUNDARY_EXTRAP
510 call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
511 default=.false., do_not_log=.true.)
512 call get_param(param_file, mdl, "SPONGE_BOUNDARY_EXTRAP", bndextrapolation, &
513 "If true, values at the interfaces of SPONGE boundary cells are "//&
514 "extrapolated instead of piecewise constant", default=bndextrapolation)
515 call get_param(param_file, mdl, "VARYING_SPONGE_MASK_THICKNESS", cs%varying_input_dz_mask, &
516 "An input file thickness below which the target values with "//&
517 "time-varying sponges are replaced by the value above.", &
518 units="m", default=0.001, scale=us%m_to_Z)
519
520 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
521 "This sets the default value for the various _ANSWER_DATE parameters.", &
522 default=99991231)
523 call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", cs%remap_answer_date, &
524 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
525 "Values below 20190101 result in the use of older, less accurate expressions "//&
526 "that were in use at the end of 2018. Higher values result in the use of more "//&
527 "robust and accurate forms of mathematically equivalent expressions.", &
528 default=default_answer_date)
529 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
530 do_not_log=.true., default=.true.)
531 call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
532 "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//&
533 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
534 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
535 call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", cs%hor_regrid_answer_date, &
536 "The vintage of the order of arithmetic for horizontal regridding. "//&
537 "Dates before 20190101 give the same answers as the code did in late 2018, "//&
538 "while later versions add parentheses for rotational symmetry. "//&
539 "Dates after 20230101 use reproducing sums for global averages.", &
540 default=default_answer_date)
541 call get_param(param_file, mdl, "SPONGE_DATA_ONGRID", cs%spongeDataOngrid, &
542 "When defined, the incoming sponge data are "//&
543 "assumed to be on the model grid " , &
544 default=.false.)
545
546 cs%time_varying_sponges = .true.
547 cs%nz = gv%ke
548
549 ! number of columns to be restored
550 cs%num_col = 0 ; cs%fldno = 0
551 do j=g%jsc,g%jec ; do i=g%isc,g%iec
552 if ((iresttime(i,j) > 0.0) .and. (g%mask2dT(i,j) > 0.0)) &
553 cs%num_col = cs%num_col + 1
554 enddo ; enddo
555 if (cs%num_col > 0) then
556 allocate(cs%Iresttime_col(cs%num_col), source=0.0)
557 allocate(cs%col_i(cs%num_col), source=0)
558 allocate(cs%col_j(cs%num_col), source=0)
559 ! pass indices, restoring time to the CS structure
560 col = 1
561 do j=g%jsc,g%jec ; do i=g%isc,g%iec
562 if ((iresttime(i,j) > 0.0) .and. (g%mask2dT(i,j) > 0.0)) then
563 cs%col_i(col) = i ; cs%col_j(col) = j
564 cs%Iresttime_col(col) = iresttime(i,j)
565 col = col + 1
566 endif
567 enddo ; enddo
568 endif
569 total_sponge_cols = cs%num_col
570 call sum_across_pes(total_sponge_cols)
571
572! Call the constructor for remapping control structure
573 if (cs%remap_answer_date >= 20190101) then
574 dz_neglect = gv%dZ_subroundoff ; dz_neglect_edge = gv%dZ_subroundoff
575 elseif (gv%Boussinesq) then
576 dz_neglect = us%m_to_Z*1.0e-30 ; dz_neglect_edge = us%m_to_Z*1.0e-10
577 elseif (gv%semi_Boussinesq) then
578 dz_neglect = gv%kg_m2_to_H*gv%H_to_Z*1.0e-30 ; dz_neglect_edge = gv%kg_m2_to_H*gv%H_to_Z*1.0e-10
579 else
580 dz_neglect = gv%dZ_subroundoff ; dz_neglect_edge = gv%dZ_subroundoff
581 endif
582 call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation, &
583 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
584 answer_date=cs%remap_answer_date, &
585 h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge)
586 call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, &
587 "The total number of columns where sponges are applied at h points.", like_default=.true.)
588 if (cs%sponge_uv) then
589 allocate(iresttime_u(g%isdB:g%iedB,g%jsd:g%jed), source=0.0)
590 allocate(iresttime_v(g%isd:g%ied,g%jsdB:g%jedB), source=0.0)
591
592 call pass_var(iresttime, g%Domain, to_all+omit_corners, halo=1)
593 ! u points
594 if (present(iresttime_u_in)) then
595 iresttime_u(:,:) = iresttime_u_in(:,:)
596 else
597 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
598 iresttime_u(i,j) = 0.5 * (iresttime(i,j) + iresttime(i+1,j))
599 enddo ; enddo
600 endif
601 cs%num_col_u = 0
602 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
603 if ((iresttime_u(i,j) > 0.0) .and. (g%mask2dCu(i,j) > 0.0)) &
604 cs%num_col_u = cs%num_col_u + 1
605 enddo ; enddo
606 if (cs%num_col_u > 0) then
607 allocate(cs%Iresttime_col_u(cs%num_col_u), source=0.0)
608 allocate(cs%col_i_u(cs%num_col_u), source=0)
609 allocate(cs%col_j_u(cs%num_col_u), source=0)
610 ! pass indices, restoring time to the CS structure
611 col = 1
612 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
613 if ((iresttime_u(i,j) > 0.0) .and. (g%mask2dCu(i,j) > 0.0)) then
614 cs%col_i_u(col) = i ; cs%col_j_u(col) = j
615 cs%Iresttime_col_u(col) = iresttime_u(i,j)
616 col = col + 1
617 endif
618 enddo ; enddo
619 ! same for total number of arbitrary layers and correspondent data
620 endif
621 total_sponge_cols_u = cs%num_col_u
622 call sum_across_pes(total_sponge_cols_u)
623 call log_param(param_file, mdl, "!Total sponge columns at u points", total_sponge_cols_u, &
624 "The total number of columns where sponges are applied at u points.", like_default=.true.)
625 ! v points
626 if (present(iresttime_v_in)) then
627 iresttime_v(:,:) = iresttime_v_in(:,:)
628 else
629 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
630 iresttime_v(i,j) = 0.5 * (iresttime(i,j) + iresttime(i,j+1))
631 enddo ; enddo
632 endif
633 cs%num_col_v = 0
634 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
635 if ((iresttime_v(i,j) > 0.0) .and. (g%mask2dCv(i,j) > 0.0)) &
636 cs%num_col_v = cs%num_col_v + 1
637 enddo ; enddo
638 if (cs%num_col_v > 0) then
639 allocate(cs%Iresttime_col_v(cs%num_col_v), source=0.0)
640 allocate(cs%col_i_v(cs%num_col_v), source=0)
641 allocate(cs%col_j_v(cs%num_col_v), source=0)
642 ! pass indices, restoring time to the CS structure
643 col = 1
644 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
645 if ((iresttime_v(i,j) > 0.0) .and. (g%mask2dCv(i,j) > 0.0)) then
646 cs%col_i_v(col) = i ; cs%col_j_v(col) = j
647 cs%Iresttime_col_v(col) = iresttime_v(i,j)
648 col = col + 1
649 endif
650 enddo ; enddo
651 endif
652 total_sponge_cols_v = cs%num_col_v
653 call sum_across_pes(total_sponge_cols_v)
654 call log_param(param_file, mdl, "!Total sponge columns at v points", total_sponge_cols_v, &
655 "The total number of columns where sponges are applied at v points.", like_default=.true.)
656 endif
657
658end subroutine initialize_ale_sponge_varying
659
660!> Initialize diagnostics for the ALE_sponge module.
661! GMM: this routine is not being used for now.
662subroutine init_ale_sponge_diags(Time, G, diag, CS, US)
663 type(time_type), target, intent(in) :: time !< The current model time
664 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
665 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
666 !! output.
667 type(ale_sponge_cs), intent(inout) :: cs !< ALE sponge control structure
668 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
669 ! Local Variables
670 character(len=:), allocatable :: tend_unit ! The units for a sponge tendency diagnostic.
671 real :: tend_conv ! The conversion factor use for the sponge tendency [A T-1 ~> a s-1]
672 integer :: m
673
674 cs%diag => diag
675
676 do m=1,cs%fldno
677 if ((trim(cs%Ref_val(m)%unit) == 'none') .or. (len_trim(cs%Ref_val(m)%unit) == 0)) then
678 tend_unit = "s-1"
679 else
680 tend_unit = trim(cs%Ref_val(m)%unit)//" s-1"
681 endif
682 tend_conv = us%s_to_T ; if (cs%Ref_val(m)%scale /= 0.0) tend_conv = us%s_to_T / cs%Ref_val(m)%scale
683 cs%id_sp_tendency(m) = register_diag_field('ocean_model', 'sp_tendency_'//cs%Ref_val(m)%name, &
684 diag%axesTL, time, long_name='Time tendency due to restoring '//cs%Ref_val(m)%long_name, &
685 units=tend_unit, conversion=tend_conv)
686 enddo
687
688 cs%id_sp_u_tendency = -1
689 cs%id_sp_u_tendency = register_diag_field('ocean_model', 'sp_tendency_u', diag%axesCuL, time, &
690 'Zonal acceleration due to sponges', 'm s-2', conversion=us%L_T2_to_m_s2)
691 cs%id_sp_v_tendency = -1
692 cs%id_sp_v_tendency = register_diag_field('ocean_model', 'sp_tendency_v', diag%axesCvL, time, &
693 'Meridional acceleration due to sponges', 'm s-2', conversion=us%L_T2_to_m_s2)
694
695end subroutine init_ale_sponge_diags
696
697!> This subroutine stores the reference profile at h points for the variable
698!! whose address is given by f_ptr.
699subroutine set_up_ale_sponge_field_fixed(sp_val, G, GV, f_ptr, CS, &
700 sp_name, sp_long_name, sp_unit, scale)
701 type(ocean_grid_type), intent(in) :: G !< Grid structure
702 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
703 type(ale_sponge_cs), pointer :: CS !< ALE sponge control structure (in/out).
704 real, dimension(SZI_(G),SZJ_(G),CS%nz_data), &
705 intent(in) :: sp_val !< Field to be used in the sponge, it can have an
706 !! arbitrary number of layers [various]
707 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
708 target, intent(in) :: f_ptr !< Pointer to the field to be damped [various]
709 character(len=*), intent(in) :: sp_name !< The name of the tracer field
710 character(len=*), optional, &
711 intent(in) :: sp_long_name !< The long name of the tracer field
712 !! if not given, use the sp_name
713 character(len=*), optional, &
714 intent(in) :: sp_unit !< The unit of the tracer field
715 !! if not given, use the none
716 real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any
717 !! contributions due to dimensional rescaling [various ~> 1].
718 !! The default is 1.
719
720 real :: scale_fac ! A factor by which to scale sp_val before storing it [various ~> 1]
721 integer :: k, col
722 character(len=256) :: mesg ! String for error messages
723 character(len=256) :: long_name ! The long name of the tracer field
724 character(len=256) :: unit ! The unit of the tracer field
725
726 if (.not.associated(cs)) return
727
728 scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
729 long_name = sp_name ; if (present(sp_long_name)) long_name = sp_long_name
730 unit = 'none' ; if (present(sp_unit)) unit = sp_unit
731
732 cs%fldno = cs%fldno + 1
733 if (cs%fldno > max_fields_) then
734 write(mesg,'("Increase MAX_FIELDS_ to at least ",I0," in MOM_memory.h or decrease &
735 &the number of fields to be damped in the call to &
736 &initialize_ALE_sponge." )') cs%fldno
737 call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
738 endif
739
740 ! stores the reference profile
741 cs%Ref_val(cs%fldno)%nz_data = cs%nz_data
742 cs%Ref_val(cs%fldno)%name = sp_name
743 cs%Ref_val(cs%fldno)%long_name = long_name
744 cs%Ref_val(cs%fldno)%unit = unit
745 cs%Ref_val(cs%fldno)%scale = scale_fac
746 allocate(cs%Ref_val(cs%fldno)%p(cs%nz_data,cs%num_col), source=0.0)
747 do col=1,cs%num_col
748 do k=1,cs%nz_data
749 cs%Ref_val(cs%fldno)%p(k,col) = scale_fac*sp_val(cs%col_i(col),cs%col_j(col),k)
750 enddo
751 enddo
752
753 cs%var(cs%fldno)%p => f_ptr
754
755end subroutine set_up_ale_sponge_field_fixed
756
757!> This subroutine stores the reference profile at h points for the variable
758!! whose address is given by filename and fieldname.
759subroutine set_up_ale_sponge_field_varying(filename, fieldname, Time, G, GV, US, f_ptr, CS, &
760 sp_name, sp_long_name, sp_unit, scale)
761 character(len=*), intent(in) :: filename !< The name of the file with the
762 !! time varying field data
763 character(len=*), intent(in) :: fieldname !< The name of the field in the file
764 !! with the time varying field data
765 type(time_type), intent(in) :: Time !< The current model time
766 type(ocean_grid_type), intent(in) :: G !< Grid structure (in).
767 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
768 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
769 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
770 target, intent(in) :: f_ptr !< Pointer to the field to be damped (in) [various].
771 type(ale_sponge_cs), pointer :: CS !< Sponge control structure (in/out).
772 character(len=*), intent(in) :: sp_name !< The name of the tracer field
773 character(len=*), optional, &
774 intent(in) :: sp_long_name !< The long name of the tracer field
775 !! if not given, use the sp_name
776 character(len=*), optional, &
777 intent(in) :: sp_unit !< The unit of the tracer field
778 !! if not given, use 'none'
779 real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any
780 !! contributions due to dimensional rescaling [various ~> 1].
781
782 ! Local variables
783 integer :: isd, ied, jsd, jed
784 integer, dimension(4) :: fld_sz
785 integer :: nz_data !< the number of vertical levels in this input field
786 character(len=256) :: mesg ! String for error messages
787 character(len=256) :: long_name ! The long name of the tracer field
788 character(len=256) :: unit ! The unit of the tracer field
789 long_name = sp_name ; if (present(sp_long_name)) long_name = sp_long_name
790 unit = 'none' ; if (present(sp_unit)) unit = sp_unit
791
792 ! Local variables for ALE remapping
793
794 if (.not.associated(cs)) return
795 ! initialize time interpolator module
796 call time_interp_external_init()
797 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
798 cs%fldno = cs%fldno + 1
799 if (cs%fldno > max_fields_) then
800 write(mesg, '("Increase MAX_FIELDS_ to at least ",I0," in MOM_memory.h or decrease "//&
801 &"the number of fields to be damped in the call to initialize_ALE_sponge." )') cs%fldno
802 call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
803 endif
804 ! get a unique time interp id for this field. If sponge data is on-grid, then setup
805 ! to only read on the computational domain
806 if (cs%spongeDataOngrid) then
807 cs%Ref_val(cs%fldno)%field = init_external_field(filename, fieldname, mom_domain=g%Domain)
808 else
809 cs%Ref_val(cs%fldno)%field = init_external_field(filename, fieldname)
810 endif
811 cs%Ref_val(cs%fldno)%name = sp_name
812 cs%Ref_val(cs%fldno)%long_name = long_name
813 cs%Ref_val(cs%fldno)%unit = unit
814 fld_sz(1:4) = -1
815 call get_external_field_info(cs%Ref_val(cs%fldno)%field, size=fld_sz, axes=cs%Ref_val(cs%fldno)%axes_data)
816 nz_data = fld_sz(3)
817 cs%Ref_val(cs%fldno)%nz_data = nz_data !< individual sponge fields may reside on a different vertical grid
818 cs%Ref_val(cs%fldno)%num_tlevs = fld_sz(4)
819 cs%Ref_val(cs%fldno)%scale = 1.0 ; if (present(scale)) cs%Ref_val(cs%fldno)%scale = scale
820 ! initializes the target profile array for this field
821 ! for all columns which will be masked
822 allocate(cs%Ref_val(cs%fldno)%p(nz_data,cs%num_col), source=0.0)
823 allocate(cs%Ref_val(cs%fldno)%dz(nz_data,cs%num_col), source=0.0)
824 cs%var(cs%fldno)%p => f_ptr
825
826end subroutine set_up_ale_sponge_field_varying
827
828!> This subroutine stores the reference profile at u and v points for the variable
829!! whose address is given by u_ptr and v_ptr.
830subroutine set_up_ale_sponge_vel_field_fixed(u_val, v_val, G, GV, u_ptr, v_ptr, CS, scale)
831 type(ocean_grid_type), intent(in) :: G !< Grid structure (in).
832 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
833 type(ale_sponge_cs), pointer :: CS !< Sponge structure (in/out).
834 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
835 intent(in) :: u_val !< u field to be used in the sponge [L T-1 ~> m s-1],
836 !! it is provided on its own vertical grid that may
837 !! have fewer layers than the model itself, but not more.
838 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
839 intent(in) :: v_val !< v field to be used in the sponge [L T-1 ~> m s-1],
840 !! it is provided on its own vertical grid that may
841 !! have fewer layers than the model itself, but not more.
842 real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1]
843 real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1]
844 real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any
845 !! contributions due to dimensional rescaling [various ~> 1].
846 !! The default is 1.
847
848 real :: scale_fac ! A dimensional rescaling factor [various ~> 1]
849 integer :: k, col
850
851 if (.not.associated(cs)) return
852
853 scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
854
855 ! stores the reference profile
856 allocate(cs%Ref_val_u%p(cs%nz_data,cs%num_col_u), source=0.0)
857 do col=1,cs%num_col_u
858 do k=1,cs%nz_data
859 cs%Ref_val_u%p(k,col) = scale_fac*u_val(cs%col_i_u(col),cs%col_j_u(col),k)
860 enddo
861 enddo
862 cs%var_u%p => u_ptr
863 allocate(cs%Ref_val_v%p(cs%nz_data,cs%num_col_v), source=0.0)
864 do col=1,cs%num_col_v
865 do k=1,cs%nz_data
866 cs%Ref_val_v%p(k,col) = scale_fac*v_val(cs%col_i_v(col),cs%col_j_v(col),k)
867 enddo
868 enddo
869 cs%var_v%p => v_ptr
870
871end subroutine set_up_ale_sponge_vel_field_fixed
872
873!> This subroutine stores the reference profile at u and v points for the variable
874!! whose address is given by u_ptr and v_ptr.
875subroutine set_up_ale_sponge_vel_field_varying(filename_u, fieldname_u, filename_v, fieldname_v, &
876 Time, G, GV, US, CS, u_ptr, v_ptr, scale)
877 character(len=*), intent(in) :: filename_u !< File name for u field
878 character(len=*), intent(in) :: fieldname_u !< Name of u variable in file
879 character(len=*), intent(in) :: filename_v !< File name for v field
880 character(len=*), intent(in) :: fieldname_v !< Name of v variable in file
881 type(time_type), intent(in) :: Time !< Model time
882 type(ocean_grid_type), intent(in) :: G !< Ocean grid (in)
883 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
884 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
885 type(ale_sponge_cs), pointer :: CS !< Sponge structure (in/out).
886 real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1]
887 real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1]
888 real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any
889 !! contributions due to dimensional rescaling, often in
890 !! [L s T-1 m-1 ~> 1]. For varying velocities the
891 !! default is the same as using US%m_s_to_L_T.
892
893 ! Local variables
894 logical :: override
895 integer :: isd, ied, jsd, jed
896 integer :: isdB, iedB, jsdB, jedB
897 integer, dimension(4) :: fld_sz
898 if (.not.associated(cs)) return
899
900 override =.true.
901
902 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
903 isdb = g%isdB ; iedb = g%iedB ; jsdb = g%jsdB ; jedb = g%jedB
904 ! get a unique id for this field which will allow us to return an array
905 ! containing time-interpolated values from an external file corresponding
906 ! to the current model date.
907 if (cs%spongeDataOngrid) then
908 cs%Ref_val_u%field = init_external_field(filename_u, fieldname_u, domain=g%Domain%mpp_domain)
909 else
910 cs%Ref_val_u%field = init_external_field(filename_u, fieldname_u)
911 endif
912 fld_sz(1:4) = -1
913 call get_external_field_info(cs%Ref_val_u%field, size=fld_sz, axes=cs%Ref_val_u%axes_data)
914 cs%Ref_val_u%nz_data = fld_sz(3)
915 cs%Ref_val_u%num_tlevs = fld_sz(4)
916 cs%Ref_val_u%scale = us%m_s_to_L_T ; if (present(scale)) cs%Ref_val_u%scale = scale
917
918 if (cs%spongeDataOngrid) then
919 cs%Ref_val_v%field = init_external_field(filename_v, fieldname_v, domain=g%Domain%mpp_domain)
920 else
921 cs%Ref_val_v%field = init_external_field(filename_v, fieldname_v)
922 endif
923 fld_sz(1:4) = -1
924 call get_external_field_info(cs%Ref_val_v%field, size=fld_sz, axes=cs%Ref_val_v%axes_data)
925 cs%Ref_val_v%nz_data = fld_sz(3)
926 cs%Ref_val_v%num_tlevs = fld_sz(4)
927 cs%Ref_val_v%scale = us%m_s_to_L_T ; if (present(scale)) cs%Ref_val_v%scale = scale
928
929 ! stores the reference profile
930 allocate(cs%Ref_val_u%p(fld_sz(3),cs%num_col_u), source=0.0)
931 allocate(cs%Ref_val_u%dz(fld_sz(3),cs%num_col_u), source=0.0)
932 cs%var_u%p => u_ptr
933 allocate(cs%Ref_val_v%p(fld_sz(3),cs%num_col_v), source=0.0)
934 allocate(cs%Ref_val_v%dz(fld_sz(3),cs%num_col_v), source=0.0)
935 cs%var_v%p => v_ptr
936
937end subroutine set_up_ale_sponge_vel_field_varying
938
939!> This subroutine applies damping to the layers thicknesses, temp, salt and a variety of tracers
940!! for every column where there is damping.
941subroutine apply_ale_sponge(h, tv, dt, G, GV, US, CS, Time)
942 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure (in).
943 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
944 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
945 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
946 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in)
947 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
948 !! thermodynamic variables
949 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
950 type(ale_sponge_cs), pointer :: cs !< A pointer to the control structure for this module
951 !! that is set by a previous call to initialize_ALE_sponge (in).
952 type(time_type), intent(in) :: time !< The current model date
953
954 ! Local variables
955 real :: damp ! The timestep times the local damping coefficient [nondim].
956 real :: i1pdamp ! I1pdamp is 1/(1 + damp). [nondim].
957 real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid [various]
958 real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid [various]
959 real, dimension(SZK_(GV)) :: dz_col ! A column of thicknesses at h, u or v points [Z ~> m]
960 real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields [various]
961 real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts [nondim]
962 real, allocatable, dimension(:,:,:) :: mask_u ! A temporary array for field mask at u pts [nondim]
963 real, allocatable, dimension(:,:,:) :: mask_v ! A temporary array for field mask at v pts [nondim]
964 real, allocatable, dimension(:,:,:) :: tmp !< A temporary array for thermodynamic sponge tendency
965 !! diagnostics [various] then in [various T-1 ~> various s-1]
966 real, allocatable, dimension(:,:,:) :: tmp_u !< A temporary array for u sponge acceleration diagnostics
967 !! first in [L T-1 ~> m s-1] then in [L T-2 ~> m s-2]
968 real, allocatable, dimension(:,:,:) :: tmp_v !< A temporary array for v sponge acceleration diagnostics
969 !! first in [L T-1 ~> m s-1] then in [L T-2 ~> m s-2]
970 real, dimension(:), allocatable :: dz_src ! Source thicknesses [Z ~> m].
971 real :: dz_model(szi_(g),szj_(g),szk_(gv)) ! Vertical distance across model layers [Z ~> m]
972
973 ! Local variables for ALE remapping
974 real, dimension(:), allocatable :: tmpt1d ! A temporary variable for ALE remapping [various]
975 integer :: c, m, i, j, k, is, ie, js, je, nz, nz_data
976 real, allocatable, dimension(:), target :: z_in ! The depths (positive downward) in the input file [Z ~> m]
977 real, allocatable, dimension(:), target :: z_edges_in ! The depths (positive downward) of the
978 ! edges in the input file [Z ~> m]
979 real :: missing_value ! The missing value in the input data field [various]
980 real :: idt ! The inverse of the timestep [T-1 ~> s-1]
981 real :: ztopofcell, zbottomofcell ! Interface heights (positive upward) in the input dataset [Z ~> m].
982 real :: sp_val_u ! Interpolation of sp_val to u-points, often a velocity in [L T-1 ~> m s-1]
983 real :: sp_val_v ! Interpolation of sp_val to v-points, often a velocity in [L T-1 ~> m s-1]
984 integer :: npoints
985
986 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
987 if (.not.associated(cs)) return
988
989 idt = 1.0/dt
990
991 if (cs%time_varying_sponges) then
992 do m=1,cs%fldno
993 nz_data = cs%Ref_val(m)%nz_data
994 call horiz_interp_and_extrap_tracer(cs%Ref_val(m)%field, time, g, sp_val, &
995 mask_z, z_in, z_edges_in, missing_value, &
996 scale=cs%Ref_val(m)%scale, spongeongrid=cs%SpongeDataOngrid, m_to_z=us%m_to_Z, &
997 answer_date=cs%hor_regrid_answer_date, axes=cs%Ref_val(m)%axes_data)
998 allocate( dz_src(nz_data) )
999 allocate( tmpt1d(nz_data) )
1000 do c=1,cs%num_col
1001 ! Set i and j to the structured indices of column c.
1002 i = cs%col_i(c) ; j = cs%col_j(c)
1003 cs%Ref_val(m)%p(1:nz_data,c) = sp_val(i,j,1:nz_data)
1004 ! Build the source grid
1005 ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0 ; dz_src(:) = 0.0 ; tmpt1d(:) = -99.9
1006 do k=1,nz_data
1007 if (mask_z(cs%col_i(c),cs%col_j(c),k) == 1.0) then
1008 zbottomofcell = -min( z_edges_in(k+1) - g%Z_ref, g%bathyT(cs%col_i(c),cs%col_j(c)) )
1009 tmpt1d(k) = sp_val(cs%col_i(c),cs%col_j(c),k)
1010 elseif (k>1) then
1011 zbottomofcell = -g%bathyT(cs%col_i(c),cs%col_j(c))
1012 tmpt1d(k) = tmpt1d(k-1)
1013 else ! This next block should only ever be reached over land
1014 tmpt1d(k) = -99.9
1015 endif
1016 dz_src(k) = ztopofcell - zbottomofcell
1017 if (dz_src(k) > 0.) npoints = npoints + 1
1018 ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
1019 enddo
1020 ! In case data is deeper than model
1021 dz_src(nz_data) = dz_src(nz_data) + ( ztopofcell + g%bathyT(cs%col_i(c),cs%col_j(c)) )
1022 cs%Ref_val(m)%dz(1:nz_data,c) = dz_src(1:nz_data)
1023 cs%Ref_val(m)%p(1:nz_data,c) = tmpt1d(1:nz_data)
1024 do k=2,nz_data
1025 if (cs%Ref_val(m)%dz(k,c) <= cs%varying_input_dz_mask) &
1026 ! some confusion here about why the masks are not correct returning from horiz_interp
1027 ! reverting to using a minimum thickness criteria
1028 cs%Ref_val(m)%p(k,c) = cs%Ref_val(m)%p(k-1,c)
1029 enddo
1030 enddo
1031 deallocate(sp_val, mask_z, dz_src, tmpt1d)
1032 enddo
1033 endif
1034
1035 tmp_val1(:) = 0.0 ; dz_col(:) = 0.0
1036 do m=1,cs%fldno
1037 nz_data = cs%Ref_val(m)%nz_data
1038 allocate(tmp_val2(cs%Ref_val(m)%nz_data))
1039 if (cs%id_sp_tendency(m) > 0) then
1040 allocate(tmp(g%isd:g%ied,g%jsd:g%jed,nz), source=0.0)
1041 endif
1042 do c=1,cs%num_col
1043 ! Set i and j to the structured indices of column c.
1044 i = cs%col_i(c) ; j = cs%col_j(c)
1045 damp = dt * cs%Iresttime_col(c)
1046 i1pdamp = 1.0 / (1.0 + damp)
1047 tmp_val2(1:nz_data) = cs%Ref_val(m)%p(1:nz_data,c)
1048 if ((.not.gv%Boussinesq) .and. allocated(tv%SpV_avg)) then
1049 do k=1,nz
1050 dz_col(k) = gv%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
1051 enddo
1052 else
1053 do k=1,nz
1054 dz_col(k) = gv%H_to_Z * h(i,j,k)
1055 enddo
1056 endif
1057 if (cs%time_varying_sponges) then
1058 call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_val(m)%dz(1:nz_data,c), tmp_val2, &
1059 cs%nz, dz_col, tmp_val1)
1060 else
1061 call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_dz%p(1:nz_data,c), tmp_val2, &
1062 cs%nz, dz_col, tmp_val1)
1063 endif
1064 !Backward Euler method
1065 if (cs%id_sp_tendency(m) > 0) tmp(i,j,1:nz) = cs%var(m)%p(i,j,1:nz)
1066 cs%var(m)%p(i,j,1:nz) = i1pdamp * (cs%var(m)%p(i,j,1:nz) + tmp_val1(1:nz) * damp)
1067 if (cs%id_sp_tendency(m) > 0) &
1068 tmp(i,j,1:cs%nz) = idt*(cs%var(m)%p(i,j,1:nz) - tmp(i,j,1:nz))
1069 enddo
1070
1071 if (cs%id_sp_tendency(m) > 0) then
1072 call post_data(cs%id_sp_tendency(m), tmp, cs%diag)
1073 deallocate(tmp)
1074 endif
1075 deallocate(tmp_val2)
1076 enddo
1077
1078 if (cs%sponge_uv) then
1079
1080 if (cs%time_varying_sponges) then
1081 nz_data = cs%Ref_val_u%nz_data
1082 ! Interpolate from the external horizontal grid and in time
1083 call horiz_interp_and_extrap_tracer(cs%Ref_val_u%field, time, g, sp_val, &
1084 mask_z, z_in, z_edges_in, missing_value, &
1085 scale=cs%Ref_val_u%scale, spongeongrid=cs%SpongeDataOngrid, m_to_z=us%m_to_Z, &
1086 answer_date=cs%hor_regrid_answer_date, axes=cs%Ref_val_u%axes_data)
1087
1088 ! Initialize mask_z halos to zero before pass_var, in case of no update
1089 mask_z(g%isc-1, g%jsc:g%jec, :) = 0.
1090 mask_z(g%iec+1, g%jsc:g%jec, :) = 0.
1091 call pass_var(sp_val, g%Domain, to_all+omit_corners, halo=1)
1092 call pass_var(mask_z, g%Domain, to_all+omit_corners, halo=1)
1093
1094 allocate(mask_u(g%isdB:g%iedB,g%jsd:g%jed,1:nz_data))
1095 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
1096 mask_u(i,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data))
1097 enddo ; enddo
1098
1099 allocate( dz_src(nz_data) )
1100 do c=1,cs%num_col_u
1101 ! Set i and j to the structured indices of column c.
1102 i = cs%col_i_u(c) ; j = cs%col_j_u(c)
1103 if (mask_u(i,j,1) == 1.0) then
1104 do k=1,nz_data
1105 sp_val_u = 0.5 * (sp_val(i,j,k) + sp_val(i+1,j,k))
1106 cs%Ref_val_u%p(k,c) = sp_val_u
1107 enddo
1108 else
1109 cs%Ref_val_u%p(1:nz_data,c) = 0.0
1110 endif
1111 ! Build the source grid
1112 ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0 ; dz_src(:) = 0.0
1113 do k=1,nz_data
1114 if (mask_u(i,j,k) == 1.0) then
1115 zbottomofcell = -min( z_edges_in(k+1) - g%Z_ref, g%bathyT(i,j) )
1116 elseif (k>1) then
1117 zbottomofcell = -g%bathyT(i,j)
1118 else ! This next block should only ever be reached over land
1119 endif
1120 dz_src(k) = ztopofcell - zbottomofcell
1121 if (dz_src(k) > 0.) npoints = npoints + 1
1122 ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
1123 enddo
1124 ! In case data is deeper than model
1125 dz_src(nz_data) = dz_src(nz_data) + ( ztopofcell + g%bathyT(i,j) )
1126 cs%Ref_val_u%dz(1:nz_data,c) = dz_src(1:nz_data)
1127 enddo
1128 deallocate(sp_val, mask_u, mask_z, dz_src)
1129 nz_data = cs%Ref_val_v%nz_data
1130 ! Interpolate from the external horizontal grid and in time
1131 call horiz_interp_and_extrap_tracer(cs%Ref_val_v%field, time, g, sp_val, &
1132 mask_z, z_in, z_edges_in, missing_value, &
1133 scale=cs%Ref_val_v%scale, spongeongrid=cs%SpongeDataOngrid, m_to_z=us%m_to_Z,&
1134 answer_date=cs%hor_regrid_answer_date, axes=cs%Ref_val_v%axes_data)
1135 ! Initialize mask_z halos to zero before pass_var, in case of no update
1136 mask_z(g%isc:g%iec, g%jsc-1, :) = 0.
1137 mask_z(g%isc:g%iec, g%jec+1, :) = 0.
1138 call pass_var(sp_val, g%Domain, to_all+omit_corners, halo=1)
1139 call pass_var(mask_z, g%Domain, to_all+omit_corners, halo=1)
1140
1141 allocate(mask_v(g%isd:g%ied,g%jsdB:g%jedB,1:nz_data))
1142 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
1143 mask_v(i,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data))
1144 enddo ; enddo
1145
1146 allocate( dz_src(nz_data) )
1147 do c=1,cs%num_col_v
1148 ! Set i and j to the structured indices of column c.
1149 i = cs%col_i_v(c) ; j = cs%col_j_v(c)
1150 if (mask_v(i,j,1) == 1.0) then
1151 do k=1,nz_data
1152 sp_val_v = 0.5 * (sp_val(i,j,k) + sp_val(i,j+1,k))
1153 cs%Ref_val_v%p(k,c) = sp_val_v
1154 enddo
1155 else
1156 cs%Ref_val_v%p(1:nz_data,c) = 0.0
1157 endif
1158 ! Build the source grid
1159 ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0 ; dz_src(:) = 0.0
1160 do k=1,nz_data
1161 if (mask_v(i,j,k) == 1.0) then
1162 zbottomofcell = -min( z_edges_in(k+1) - g%Z_ref, g%bathyT(i,j) )
1163 elseif (k>1) then
1164 zbottomofcell = -g%bathyT(i,j)
1165 else ! This next block should only ever be reached over land
1166 endif
1167 dz_src(k) = ztopofcell - zbottomofcell
1168 if (dz_src(k) > 0.) npoints = npoints + 1
1169 ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
1170 enddo
1171 ! In case data is deeper than model
1172 dz_src(nz_data) = dz_src(nz_data) + ( ztopofcell + g%bathyT(i,j) )
1173 cs%Ref_val_v%dz(1:nz_data,c) = dz_src(1:nz_data)
1174 enddo
1175 deallocate(sp_val, mask_v, mask_z, dz_src)
1176 endif
1177
1178 ! Because we can not be certain whether there are velocity points at the processor
1179 ! boundaries, and the thicknesses might not have been updated there, we need to
1180 ! calculate the tracer point layer vertical extents and then do a halo update.
1181 if ((.not.gv%Boussinesq) .and. allocated(tv%SpV_avg)) then
1182 do k=1,nz ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
1183 dz_model(i,j,k) = gv%H_to_RZ * (h(i,j,k)*tv%SpV_avg(i,j,k))
1184 enddo ; enddo ; enddo
1185 else
1186 do k=1,nz ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
1187 dz_model(i,j,k) = gv%H_to_Z * h(i,j,k)
1188 enddo ; enddo ; enddo
1189 endif
1190 call pass_var(dz_model, g%Domain, to_all+omit_corners, halo=1)
1191
1192 nz_data = cs%Ref_val_u%nz_data
1193 allocate(tmp_val2(nz_data))
1194 if (cs%id_sp_u_tendency > 0) then
1195 allocate(tmp_u(g%isdB:g%iedB,g%jsd:g%jed,nz), source=0.0)
1196 endif
1197 ! u points
1198 do c=1,cs%num_col_u
1199 i = cs%col_i_u(c) ; j = cs%col_j_u(c)
1200 damp = dt * cs%Iresttime_col_u(c)
1201 i1pdamp = 1.0 / (1.0 + damp)
1202 tmp_val2(1:nz_data) = cs%Ref_val_u%p(1:nz_data,c)
1203 do k=1,nz
1204 dz_col(k) = 0.5 * (dz_model(i,j,k) + dz_model(i+1,j,k))
1205 enddo
1206 if (cs%time_varying_sponges) then
1207 call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_val_u%dz(1:nz_data,c), tmp_val2, &
1208 cs%nz, dz_col, tmp_val1)
1209 else
1210 call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_dzu%p(1:nz_data,c), tmp_val2, &
1211 cs%nz, dz_col, tmp_val1)
1212 endif
1213 if (cs%id_sp_u_tendency > 0) tmp_u(i,j,1:nz) = cs%var_u%p(i,j,1:nz)
1214 !Backward Euler method
1215 cs%var_u%p(i,j,1:nz) = i1pdamp * (cs%var_u%p(i,j,1:nz) + tmp_val1 * damp)
1216 if (cs%id_sp_u_tendency > 0) tmp_u(i,j,1:nz) = idt*(cs%var_u%p(i,j,1:nz) - tmp_u(i,j,1:nz))
1217 enddo
1218 deallocate(tmp_val2)
1219 if (cs%id_sp_u_tendency > 0) then
1220 call post_data(cs%id_sp_u_tendency, tmp_u, cs%diag)
1221 deallocate(tmp_u)
1222 endif
1223 ! v points
1224 if (cs%id_sp_v_tendency > 0) then
1225 allocate(tmp_v(g%isd:g%ied,g%jsdB:g%jedB,nz), source=0.0)
1226 endif
1227 nz_data = cs%Ref_val_v%nz_data
1228 allocate(tmp_val2(nz_data))
1229
1230 do c=1,cs%num_col_v
1231 i = cs%col_i_v(c) ; j = cs%col_j_v(c)
1232 damp = dt * cs%Iresttime_col_v(c)
1233 i1pdamp = 1.0 / (1.0 + damp)
1234 if (cs%time_varying_sponges) nz_data = cs%Ref_val_v%nz_data
1235 tmp_val2(1:nz_data) = cs%Ref_val_v%p(1:nz_data,c)
1236 do k=1,nz
1237 dz_col(k) = 0.5 * (dz_model(i,j,k) + dz_model(i,j+1,k))
1238 enddo
1239 if (cs%time_varying_sponges) then
1240 call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_val_v%dz(1:nz_data,c), tmp_val2, &
1241 cs%nz, dz_col, tmp_val1)
1242 else
1243 call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_dzv%p(1:nz_data,c), tmp_val2, &
1244 cs%nz, dz_col, tmp_val1)
1245 endif
1246 if (cs%id_sp_v_tendency > 0) tmp_v(i,j,1:nz) = cs%var_v%p(i,j,1:nz)
1247 !Backward Euler method
1248 cs%var_v%p(i,j,1:nz) = i1pdamp * (cs%var_v%p(i,j,1:nz) + tmp_val1 * damp)
1249 if (cs%id_sp_v_tendency > 0) tmp_v(i,j,1:nz) = idt*(cs%var_v%p(i,j,1:nz) - tmp_v(i,j,1:nz))
1250 enddo
1251 if (cs%id_sp_v_tendency > 0) then
1252 call post_data(cs%id_sp_v_tendency, tmp_v, cs%diag)
1253 deallocate(tmp_v)
1254 endif
1255 deallocate(tmp_val2)
1256 endif
1257
1258end subroutine apply_ale_sponge
1259
1260!> Rotate the ALE sponge fields from the input to the model index map.
1261subroutine rotate_ale_sponge(sponge_in, G_in, sponge, G, GV, US, turns, param_file)
1262 type(ale_sponge_cs), intent(in) :: sponge_in !< The control structure for this module with the
1263 !! original grid rotation
1264 type(ocean_grid_type), intent(in) :: g_in !< The ocean's grid structure with the original rotation.
1265 type(ale_sponge_cs), pointer :: sponge !< A pointer to the control that will be set up with
1266 !! the new grid rotation
1267 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure with the new rotation.
1268 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1269 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1270 integer, intent(in) :: turns !< The number of 90-degree turns between grids
1271 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
1272 !! to parse for model parameter values.
1273
1274 ! First part: Index construction
1275 ! 1. Reconstruct Iresttime(:,:) from sponge_in
1276 ! 2. rotate Iresttime(:,:)
1277 ! 3. Call initialize_ALE_sponge using new grid and rotated Iresttime(:,:)
1278 ! All the index adjustment should follow from the Iresttime rotation
1279
1280 real, dimension(:,:), allocatable :: iresttime_in ! Restoring rate on the input sponges [T-1 ~> s-1]
1281 real, dimension(:,:), allocatable :: iresttime ! Restoring rate on the output sponges [T-1 ~> s-1]
1282 real, dimension(:,:,:), allocatable :: data_dz_in ! Grid for the input sponges [Z ~> m]
1283 real, dimension(:,:,:), allocatable :: data_dz ! Grid for the output sponges [Z ~> m]
1284 real, dimension(:,:,:), allocatable :: sp_val_in ! Target data for the input sponges [various]
1285 real, dimension(:,:,:), allocatable :: sp_val ! Target data for the output sponges [various]
1286 real, dimension(:,:,:), pointer :: sp_ptr => null() ! Target data for the input sponges [various]
1287 integer :: c, c_i, c_j
1288 integer :: k, nz_data
1289 integer :: n
1290 logical :: fixed_sponge
1291
1292 fixed_sponge = .not. sponge_in%time_varying_sponges
1293 ! NOTE: nz_data is only conditionally set when fixed_sponge is true.
1294
1295 allocate(iresttime_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed), source=0.0)
1296 allocate(iresttime(g%isd:g%ied, g%jsd:g%jed))
1297
1298 if (fixed_sponge) then
1299 nz_data = sponge_in%nz_data
1300 allocate(data_dz_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz_data), source=0.0)
1301 allocate(data_dz(g%isd:g%ied, g%jsd:g%jed, nz_data))
1302 endif
1303
1304 ! Re-populate the 2D Iresttime and data_dz arrays on the original grid
1305 do c=1,sponge_in%num_col
1306 c_i = sponge_in%col_i(c)
1307 c_j = sponge_in%col_j(c)
1308 iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c)
1309 if (fixed_sponge) then
1310 do k = 1, nz_data
1311 data_dz_in(c_i, c_j, k) = sponge_in%Ref_dz%p(k,c)
1312 enddo
1313 endif
1314 enddo
1315
1316 call rotate_array(iresttime_in, turns, iresttime)
1317 if (fixed_sponge) then
1318 call rotate_array(data_dz_in, turns, data_dz)
1319 call initialize_ale_sponge_fixed(iresttime, g, gv, param_file, sponge, &
1320 data_dz, nz_data, data_h_is_z=.true.)
1321 else
1322 call initialize_ale_sponge_varying(iresttime, g, gv, us, param_file, sponge)
1323 endif
1324
1325 deallocate(iresttime_in)
1326 deallocate(iresttime)
1327 if (fixed_sponge) then
1328 deallocate(data_dz_in)
1329 deallocate(data_dz)
1330 endif
1331
1332 ! Second part: Provide rotated fields for which relaxation is applied
1333
1334 if (fixed_sponge) then
1335 allocate(sp_val_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz_data))
1336 allocate(sp_val(g%isd:g%ied, g%jsd:g%jed, nz_data))
1337 ! For a fixed sponge, sponge%fldno is incremented from 0 in the calls to set_up_ALE_sponge_field.
1338 sponge%fldno = 0
1339 else
1340 sponge%fldno = sponge_in%fldno
1341 endif
1342
1343 do n=1,sponge_in%fldno
1344 ! Assume that tracers are pointers and are remapped in other functions(?)
1345 sp_ptr => sponge_in%var(n)%p
1346 if (fixed_sponge) then
1347 sp_val_in(:,:,:) = 0.0
1348 do c = 1, sponge_in%num_col
1349 c_i = sponge_in%col_i(c)
1350 c_j = sponge_in%col_j(c)
1351 do k = 1, nz_data
1352 sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c)
1353 enddo
1354 enddo
1355
1356 call rotate_array(sp_val_in, turns, sp_val)
1357
1358 ! NOTE: This points sp_val with the unrotated field. See note below.
1359 call set_up_ale_sponge_field(sp_val, g, gv, sp_ptr, sponge, &
1360 sponge_in%Ref_val(n)%name, sp_long_name=sponge_in%Ref_val(n)%long_name, &
1361 sp_unit=sponge_in%Ref_val(n)%unit)
1362
1363 deallocate(sp_val_in)
1364 else
1365 ! We don't want to repeat FMS init in set_up_ALE_sponge_field_varying()
1366 ! (time_interp_external_init, init_external_field, etc), so we manually
1367 ! do a portion of this function below.
1368 sponge%Ref_val(n)%field = sponge_in%Ref_val(n)%field
1369 sponge%Ref_val(n)%num_tlevs = sponge_in%Ref_val(n)%num_tlevs
1370
1371 nz_data = sponge_in%Ref_val(n)%nz_data
1372 sponge%Ref_val(n)%nz_data = nz_data
1373
1374 allocate(sponge%Ref_val(n)%p(nz_data, sponge_in%num_col), source=0.0)
1375 allocate(sponge%Ref_val(n)%dz(nz_data, sponge_in%num_col), source=0.0)
1376
1377 ! TODO: There is currently no way to associate a generic field pointer to
1378 ! its rotated equivalent without introducing a new data structure which
1379 ! explicitly tracks the pairing.
1380 !
1381 ! As a temporary fix, we store the pointer to the unrotated field in
1382 ! the rotated sponge, and use this reference to replace the pointer
1383 ! to the rotated field update_ALE_sponge field.
1384 !
1385 ! This makes a lot of unverifiable assumptions, and should not be
1386 ! considered the final solution.
1387 sponge%var(n)%p => sp_ptr
1388 endif
1389 enddo
1390
1391 ! TODO: var_u and var_v sponge damping is not yet supported.
1392 if (associated(sponge_in%var_u%p) .or. associated(sponge_in%var_v%p)) &
1393 call mom_error(fatal, "Rotation of ALE sponge velocities is not yet implemented.")
1394
1395 ! Transfer any existing diag_CS reference pointer
1396 sponge%diag => sponge_in%diag
1397
1398 ! NOTE: initialize_ALE_sponge_* resolves remap_cs
1399end subroutine rotate_ale_sponge
1400
1401
1402!> Scan the ALE sponge variables and replace a prescribed pointer to a new value.
1403! TODO: This function solely exists to replace field pointers in the sponge
1404! after rotation. This function is part of a temporary solution until
1405! something more robust is developed.
1406subroutine update_ale_sponge_field(sponge, p_old, G, GV, p_new)
1407 type(ale_sponge_cs), intent(inout) :: sponge !< ALE sponge control struct
1408 real, dimension(:,:,:), &
1409 target, intent(in) :: p_old !< The previous array of target values [various]
1410 type(ocean_grid_type), intent(in) :: g !< The updated ocean grid structure
1411 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1412 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1413 target, intent(in) :: p_new !< The new array of target values [various]
1414
1415 integer :: n
1416
1417 do n=1,sponge%fldno
1418 if (associated(sponge%var(n)%p, p_old)) sponge%var(n)%p => p_new
1419 enddo
1420end subroutine update_ale_sponge_field
1421
1422
1423! GMM: I could not find where sponge_end is being called, but I am keeping
1424! ALE_sponge_end here so we can add that if needed.
1425!> This subroutine deallocates any memory associated with the ALE_sponge module.
1426subroutine ale_sponge_end(CS)
1427 type(ale_sponge_cs), pointer :: cs !< A pointer to the control structure that is
1428 !! set by a previous call to initialize_ALE_sponge.
1429
1430 integer :: m
1431
1432 if (.not.associated(cs)) return
1433
1434 if (allocated(cs%col_i)) deallocate(cs%col_i)
1435 if (allocated(cs%col_i_u)) deallocate(cs%col_i_u)
1436 if (allocated(cs%col_i_v)) deallocate(cs%col_i_v)
1437 if (allocated(cs%col_j)) deallocate(cs%col_j)
1438 if (allocated(cs%col_j_u)) deallocate(cs%col_j_u)
1439 if (allocated(cs%col_j_v)) deallocate(cs%col_j_v)
1440
1441 if (allocated(cs%Iresttime_col)) deallocate(cs%Iresttime_col)
1442 if (allocated(cs%Iresttime_col_u)) deallocate(cs%Iresttime_col_u)
1443 if (allocated(cs%Iresttime_col_v)) deallocate(cs%Iresttime_col_v)
1444
1445 do m=1,cs%fldno
1446 if (associated(cs%Ref_val(m)%p)) deallocate(cs%Ref_val(m)%p)
1447 enddo
1448
1449 deallocate(cs)
1450
1451end subroutine ale_sponge_end
1452
1453end module mom_ale_sponge