MOM6
MOM_ALE_sponge.F90
Go to the documentation of this file.
1 !> This module contains the routines used to apply sponge layers when using
2 !! the ALE mode.
3 !!
4 !! Applying sponges requires the following:
5 !! 1. initialize_ALE_sponge
6 !! 2. set_up_ALE_sponge_field (tracers) and set_up_ALE_sponge_vel_field (vel)
7 !! 3. apply_ALE_sponge
8 !! 4. init_ALE_sponge_diags (not being used for now)
9 !! 5. ALE_sponge_end (not being used for now)
10 
12 
13 
14 ! This file is part of MOM6. See LICENSE.md for the license.
15 use mom_coms, only : sum_across_pes
17 use mom_diag_mediator, only : diag_ctrl
18 use mom_error_handler, only : mom_error, fatal, note, warning, is_root_pe
20 use mom_grid, only : ocean_grid_type
23 use mom_time_manager, only : time_type, init_external_field, get_external_field_size, time_interp_external_init
27 
28 implicit none ; private
29 
30 #include <MOM_memory.h>
31 
32 !> Store the reference profile at h points for a variable
34  module procedure set_up_ale_sponge_field_fixed
35  module procedure set_up_ale_sponge_field_varying
36 end interface
37 
38 !> This subroutine stores the reference profile at u and v points for a vector
40  module procedure set_up_ale_sponge_vel_field_fixed
42 end interface
43 
44 !> Ddetermine the number of points which are within sponges in this computational domain.
45 !!
46 !! Only points that have positive values of Iresttime and which mask2dT indicates are ocean
47 !! points are included in the sponges. It also stores the target interface heights.
49  module procedure initialize_ale_sponge_fixed
50  module procedure initialize_ale_sponge_varying
51 end interface
52 
53 ! Publicly available functions
57 
58 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
59 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
60 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
61 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
62 
63 !> A structure for creating arrays of pointers to 3D arrays with extra gridding information
64 type :: p3d
65  integer :: id !< id for FMS external time interpolator
66  integer :: nz_data !< The number of vertical levels in the input field.
67  integer :: num_tlevs !< The number of time records contained in the file
68  real, dimension(:,:,:), pointer :: mask_in => null() !< pointer to the data mask.
69  real, dimension(:,:,:), pointer :: p => null() !< pointer to the data.
70  real, dimension(:,:,:), pointer :: h => null() !< pointer to the data grid.
71 end type p3d
72 
73 !> A structure for creating arrays of pointers to 2D arrays with extra gridding information
74 type :: p2d
75  integer :: id !< id for FMS external time interpolator
76  integer :: nz_data !< The number of vertical levels in the input field
77  integer :: num_tlevs !< The number of time records contained in the file
78  real, dimension(:,:), pointer :: mask_in => null()!< pointer to the data mask.
79  real, dimension(:,:), pointer :: p => null() !< pointer the data.
80  real, dimension(:,:), pointer :: h => null() !< pointer the data grid.
81 end type p2d
82 
83 !> ALE sponge control structure
84 type, public :: ale_sponge_cs ; private
85  integer :: nz !< The total number of layers.
86  integer :: nz_data !< The total number of arbritary layers (used by older code).
87  integer :: isc !< The starting i-index of the computational domain at h.
88  integer :: iec !< The ending i-index of the computational domain at h.
89  integer :: jsc !< The starting j-index of the computational domain at h.
90  integer :: jec !< The ending j-index of the computational domain at h.
91  integer :: iscb !< The starting I-index of the computational domain at u/v.
92  integer :: iecb !< The ending I-index of the computational domain at u/v.
93  integer :: jscb !< The starting J-index of the computational domain at u/v.
94  integer :: jecb !< The ending J-index of the computational domain at h.
95  integer :: isd !< The starting i-index of the data domain at h.
96  integer :: ied !< The ending i-index of the data domain at h.
97  integer :: jsd !< The starting j-index of the data domain at h.
98  integer :: jed !< The ending j-index of the data domain at h.
99  integer :: num_col !< The number of sponge tracer points within the computational domain.
100  integer :: num_col_u !< The number of sponge u-points within the computational domain.
101  integer :: num_col_v !< The number of sponge v-points within the computational domain.
102  integer :: fldno = 0 !< The number of fields which have already been
103  !! registered by calls to set_up_sponge_field
104  logical :: sponge_uv !< Control whether u and v are included in sponge
105  integer, pointer :: col_i(:) => null() !< Array of the i-indicies of each tracer columns being damped.
106  integer, pointer :: col_j(:) => null() !< Array of the j-indicies of each tracer columns being damped.
107  integer, pointer :: col_i_u(:) => null() !< Array of the i-indicies of each u-columns being damped.
108  integer, pointer :: col_j_u(:) => null() !< Array of the j-indicies of each u-columns being damped.
109  integer, pointer :: col_i_v(:) => null() !< Array of the i-indicies of each v-columns being damped.
110  integer, pointer :: col_j_v(:) => null() !< Array of the j-indicies of each v-columns being damped.
111 
112  real, pointer :: iresttime_col(:) => null() !< The inverse restoring time of each tracer column [T-1 ~> s-1].
113  real, pointer :: iresttime_col_u(:) => null() !< The inverse restoring time of each u-column [T-1 ~> s-1].
114  real, pointer :: iresttime_col_v(:) => null() !< The inverse restoring time of each v-column [T-1 ~> s-1].
115 
116  type(p3d) :: var(max_fields_) !< Pointers to the fields that are being damped.
117  type(p2d) :: ref_val(max_fields_) !< The values to which the fields are damped.
118  type(p2d) :: ref_val_u !< The values to which the u-velocities are damped.
119  type(p2d) :: ref_val_v !< The values to which the v-velocities are damped.
120  type(p3d) :: var_u !< Pointer to the u velocities. that are being damped.
121  type(p3d) :: var_v !< Pointer to the v velocities. that are being damped.
122  type(p2d) :: ref_h !< Grid on which reference data is provided (older code).
123  type(p2d) :: ref_hu !< u-point grid on which reference data is provided (older code).
124  type(p2d) :: ref_hv !< v-point grid on which reference data is provided (older code).
125 
126  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
127  !! timing of diagnostic output.
128 
129  type(remapping_cs) :: remap_cs !< Remapping parameters and work arrays
130 
131  logical :: time_varying_sponges !< True if using newer sponge code
132  logical :: spongedataongrid !< True if the sponge data are on the model horizontal grid
133 end type ale_sponge_cs
134 
135 contains
136 
137 !> This subroutine determines the number of points which are within sponges in this computational
138 !! domain. Only points that have positive values of Iresttime and which mask2dT indicates are ocean
139 !! points are included in the sponges. It also stores the target interface heights.
140 subroutine initialize_ale_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_data)
141 
142  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
143  integer, intent(in) :: nz_data !< The total number of sponge input layers.
144  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Iresttime !< The inverse of the restoring time [s-1].
145  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
146  !! to parse for model parameter values.
147  type(ale_sponge_cs), pointer :: CS !< A pointer that is set to point to the control
148  !! structure for this module (in/out).
149  real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_h !< The thicknesses of the sponge
150  !! input layers [H ~> m or kg m-2].
151 
152 
153 ! This include declares and sets the variable "version".
154 #include "version_variable.h"
155  character(len=40) :: mdl = "MOM_sponge" ! This module's name.
156  logical :: use_sponge
157  real, allocatable, dimension(:,:,:) :: data_hu !< thickness at u points [H ~> m or kg m-2]
158  real, allocatable, dimension(:,:,:) :: data_hv !< thickness at v points [H ~> m or kg m-2]
159  real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [s-1]
160  real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [s-1]
161  logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
162  integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v
163  character(len=10) :: remapScheme
164  if (associated(cs)) then
165  call mom_error(warning, "initialize_sponge called with an associated "// &
166  "control structure.")
167  return
168  endif
169 
170 ! Set default, read and log parameters
171  call log_version(param_file, mdl, version, "")
172  call get_param(param_file, mdl, "SPONGE", use_sponge, &
173  "If true, sponges may be applied anywhere in the domain. "//&
174  "The exact location and properties of those sponges are "//&
175  "specified from MOM_initialization.F90.", default=.false.)
176 
177  if (.not.use_sponge) return
178 
179  allocate(cs)
180 
181  call get_param(param_file, mdl, "SPONGE_UV", cs%sponge_uv, &
182  "Apply sponges in u and v, in addition to tracers.", &
183  default=.false.)
184 
185  call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
186  "This sets the reconstruction scheme used "//&
187  " for vertical remapping for all variables.", &
188  default="PLM", do_not_log=.true.)
189 
190  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
191  "When defined, a proper high-order reconstruction "//&
192  "scheme is used within boundary cells rather "//&
193  "than PCM. E.g., if PPM is used for remapping, a "//&
194  "PPM reconstruction will also be used within boundary cells.", &
195  default=.false., do_not_log=.true.)
196 
197  cs%time_varying_sponges = .false.
198  cs%nz = g%ke
199  cs%isc = g%isc ; cs%iec = g%iec ; cs%jsc = g%jsc ; cs%jec = g%jec
200  cs%isd = g%isd ; cs%ied = g%ied ; cs%jsd = g%jsd ; cs%jed = g%jed
201  cs%iscB = g%iscB ; cs%iecB = g%iecB; cs%jscB = g%jscB ; cs%jecB = g%jecB
202 
203  ! number of columns to be restored
204  cs%num_col = 0 ; cs%fldno = 0
205  do j=g%jsc,g%jec ; do i=g%isc,g%iec
206  if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) &
207  cs%num_col = cs%num_col + 1
208  enddo ; enddo
209 
210  if (cs%num_col > 0) then
211  allocate(cs%Iresttime_col(cs%num_col)) ; cs%Iresttime_col = 0.0
212  allocate(cs%col_i(cs%num_col)) ; cs%col_i = 0
213  allocate(cs%col_j(cs%num_col)) ; cs%col_j = 0
214  ! pass indices, restoring time to the CS structure
215  col = 1
216  do j=g%jsc,g%jec ; do i=g%isc,g%iec
217  if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) then
218  cs%col_i(col) = i ; cs%col_j(col) = j
219  cs%Iresttime_col(col) = g%US%T_to_s*iresttime(i,j)
220  col = col +1
221  endif
222  enddo ; enddo
223  ! same for total number of arbritary layers and correspondent data
224  cs%nz_data = nz_data
225  allocate(cs%Ref_h%p(cs%nz_data,cs%num_col))
226  do col=1,cs%num_col ; do k=1,cs%nz_data
227  cs%Ref_h%p(k,col) = data_h(cs%col_i(col),cs%col_j(col),k)
228  enddo ; enddo
229  endif
230 
231  total_sponge_cols = cs%num_col
232  call sum_across_pes(total_sponge_cols)
233 
234 ! Call the constructor for remapping control structure
235  call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation)
236 
237  call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, &
238  "The total number of columns where sponges are applied at h points.")
239 
240  if (cs%sponge_uv) then
241 
242  allocate(data_hu(g%isdB:g%iedB,g%jsd:g%jed,nz_data)); data_hu(:,:,:)=0.0
243  allocate(data_hv(g%isd:g%ied,g%jsdB:g%jedB,nz_data)); data_hv(:,:,:)=0.0
244  allocate(iresttime_u(g%isdB:g%iedB,g%jsd:g%jed)); iresttime_u(:,:)=0.0
245  allocate(iresttime_v(g%isd:g%ied,g%jsdB:g%jedB)); iresttime_v(:,:)=0.0
246 
247  ! u points
248  cs%num_col_u = 0 ; !CS%fldno_u = 0
249  do j=cs%jsc,cs%jec; do i=cs%iscB,cs%iecB
250  data_hu(i,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:))
251  iresttime_u(i,j) = 0.5 * (iresttime(i,j) + iresttime(i+1,j))
252  if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) &
253  cs%num_col_u = cs%num_col_u + 1
254  enddo ; enddo
255 
256  if (cs%num_col_u > 0) then
257 
258  allocate(cs%Iresttime_col_u(cs%num_col_u)) ; cs%Iresttime_col_u = 0.0
259  allocate(cs%col_i_u(cs%num_col_u)) ; cs%col_i_u = 0
260  allocate(cs%col_j_u(cs%num_col_u)) ; cs%col_j_u = 0
261 
262  ! pass indices, restoring time to the CS structure
263  col = 1
264  do j=cs%jsc,cs%jec ; do i=cs%iscB,cs%iecB
265  if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) then
266  cs%col_i_u(col) = i ; cs%col_j_u(col) = j
267  cs%Iresttime_col_u(col) = g%US%T_to_s*iresttime_u(i,j)
268  col = col +1
269  endif
270  enddo ; enddo
271 
272  ! same for total number of arbritary layers and correspondent data
273 
274  allocate(cs%Ref_hu%p(cs%nz_data,cs%num_col_u))
275  do col=1,cs%num_col_u ; do k=1,cs%nz_data
276  cs%Ref_hu%p(k,col) = data_hu(cs%col_i_u(col),cs%col_j_u(col),k)
277  enddo ; enddo
278  endif
279  total_sponge_cols_u = cs%num_col_u
280  call sum_across_pes(total_sponge_cols_u)
281  call log_param(param_file, mdl, "!Total sponge columns at u points", total_sponge_cols_u, &
282  "The total number of columns where sponges are applied at u points.")
283 
284  ! v points
285  cs%num_col_v = 0 ; !CS%fldno_v = 0
286  do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec
287  data_hv(i,j,:) = 0.5 * (data_h(i,j,:) + data_h(i,j+1,:))
288  iresttime_v(i,j) = 0.5 * (iresttime(i,j) + iresttime(i,j+1))
289  if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) &
290  cs%num_col_v = cs%num_col_v + 1
291  enddo ; enddo
292 
293  if (cs%num_col_v > 0) then
294 
295  allocate(cs%Iresttime_col_v(cs%num_col_v)) ; cs%Iresttime_col_v = 0.0
296  allocate(cs%col_i_v(cs%num_col_v)) ; cs%col_i_v = 0
297  allocate(cs%col_j_v(cs%num_col_v)) ; cs%col_j_v = 0
298 
299  ! pass indices, restoring time to the CS structure
300  col = 1
301  do j=cs%jscB,cs%jecB ; do i=cs%isc,cs%iec
302  if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) then
303  cs%col_i_v(col) = i ; cs%col_j_v(col) = j
304  cs%Iresttime_col_v(col) = g%US%T_to_s*iresttime_v(i,j)
305  col = col +1
306  endif
307  enddo ; enddo
308 
309  ! same for total number of arbritary layers and correspondent data
310  allocate(cs%Ref_hv%p(cs%nz_data,cs%num_col_v))
311  do col=1,cs%num_col_v ; do k=1,cs%nz_data
312  cs%Ref_hv%p(k,col) = data_hv(cs%col_i_v(col),cs%col_j_v(col),k)
313  enddo ; enddo
314  endif
315  total_sponge_cols_v = cs%num_col_v
316  call sum_across_pes(total_sponge_cols_v)
317  call log_param(param_file, mdl, "!Total sponge columns at v points", total_sponge_cols_v, &
318  "The total number of columns where sponges are applied at v points.")
319  endif
320 
321 end subroutine initialize_ale_sponge_fixed
322 
323 !> Return the number of layers in the data with a fixed ALE sponge, or 0 if there are
324 !! no sponge columns on this PE.
325 function get_ale_sponge_nz_data(CS)
326  type(ale_sponge_cs), pointer :: cs !< A pointer that is set to point to the control
327  !! structure for the ALE_sponge module.
328  integer :: get_ale_sponge_nz_data !< The number of layers in the fixed sponge data.
329 
330  if (associated(cs)) then
331  get_ale_sponge_nz_data = cs%nz_data
332  else
334  endif
335 end function get_ale_sponge_nz_data
336 
337 !> Return the thicknesses used for the data with a fixed ALE sponge
338 subroutine get_ale_sponge_thicknesses(G, data_h, sponge_mask, CS)
339  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure (in).
340  real, allocatable, dimension(:,:,:), &
341  intent(inout) :: data_h !< The thicknesses of the sponge input layers [H ~> m or kg m-2].
342  logical, dimension(SZI_(G),SZJ_(G)), &
343  intent(out) :: sponge_mask !< A logical mask that is true where
344  !! sponges are being applied.
345  type(ale_sponge_cs), pointer :: cs !< A pointer that is set to point to the control
346  !! structure for the ALE_sponge module.
347  integer :: c, i, j, k
348 
349  if (allocated(data_h)) call mom_error(fatal, &
350  "get_ALE_sponge_thicknesses called with an allocated data_h.")
351 
352  if (.not.associated(cs)) then
353  ! There are no sponge points on this PE.
354  allocate(data_h(g%isd:g%ied,g%jsd:g%jed,1)) ; data_h(:,:,:) = -1.0
355  sponge_mask(:,:) = .false.
356  return
357  endif
358 
359  allocate(data_h(g%isd:g%ied,g%jsd:g%jed,cs%nz_data)) ; data_h(:,:,:) = -1.0
360  sponge_mask(:,:) = .false.
361 
362  do c=1,cs%num_col
363  i = cs%col_i(c) ; j = cs%col_j(c)
364  sponge_mask(i,j) = .true.
365  do k=1,cs%nz_data
366  data_h(i,j,k) = cs%Ref_h%p(k,c)
367  enddo
368  enddo
369 
370 end subroutine get_ale_sponge_thicknesses
371 
372 !> This subroutine determines the number of points which are to be restoref in the computational
373 !! domain. Only points that have positive values of Iresttime and which mask2dT indicates are ocean
374 !! points are included in the sponges.
375 subroutine initialize_ale_sponge_varying(Iresttime, G, param_file, CS)
376 
377  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
378  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Iresttime !< The inverse of the restoring time [s-1].
379  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to parse
380  !! for model parameter values.
381  type(ale_sponge_cs), pointer :: CS !< A pointer that is set to point to the control
382  !! structure for this module (in/out).
383 
384 ! This include declares and sets the variable "version".
385 #include "version_variable.h"
386  character(len=40) :: mdl = "MOM_sponge" ! This module's name.
387  logical :: use_sponge
388  real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [s-1]
389  real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [s-1]
390  logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
391  logical :: spongeDataOngrid = .false.
392  integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v
393  character(len=10) :: remapScheme
394 
395  if (associated(cs)) then
396  call mom_error(warning, "initialize_sponge called with an associated "// &
397  "control structure.")
398  return
399  endif
400 ! Set default, read and log parameters
401  call log_version(param_file, mdl, version, "")
402  call get_param(param_file, mdl, "SPONGE", use_sponge, &
403  "If true, sponges may be applied anywhere in the domain. "//&
404  "The exact location and properties of those sponges are "//&
405  "specified from MOM_initialization.F90.", default=.false.)
406  if (.not.use_sponge) return
407  allocate(cs)
408  call get_param(param_file, mdl, "SPONGE_UV", cs%sponge_uv, &
409  "Apply sponges in u and v, in addition to tracers.", &
410  default=.false.)
411  call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
412  "This sets the reconstruction scheme used "//&
413  " for vertical remapping for all variables.", &
414  default="PLM", do_not_log=.true.)
415  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
416  "When defined, a proper high-order reconstruction "//&
417  "scheme is used within boundary cells rather "//&
418  "than PCM. E.g., if PPM is used for remapping, a "//&
419  "PPM reconstruction will also be used within boundary cells.", &
420  default=.false., do_not_log=.true.)
421  call get_param(param_file, mdl, "SPONGE_DATA_ONGRID", cs%spongeDataOngrid, &
422  "When defined, the incoming sponge data are "//&
423  "assumed to be on the model grid " , &
424  default=.false.)
425  cs%time_varying_sponges = .true.
426  cs%nz = g%ke
427  cs%isc = g%isc ; cs%iec = g%iec ; cs%jsc = g%jsc ; cs%jec = g%jec
428  cs%isd = g%isd ; cs%ied = g%ied ; cs%jsd = g%jsd ; cs%jed = g%jed
429  cs%iscB = g%iscB ; cs%iecB = g%iecB; cs%jscB = g%jscB ; cs%jecB = g%jecB
430 
431  ! number of columns to be restored
432  cs%num_col = 0 ; cs%fldno = 0
433  do j=g%jsc,g%jec ; do i=g%isc,g%iec
434  if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) &
435  cs%num_col = cs%num_col + 1
436  enddo ; enddo
437  if (cs%num_col > 0) then
438  allocate(cs%Iresttime_col(cs%num_col)) ; cs%Iresttime_col = 0.0
439  allocate(cs%col_i(cs%num_col)) ; cs%col_i = 0
440  allocate(cs%col_j(cs%num_col)) ; cs%col_j = 0
441  ! pass indices, restoring time to the CS structure
442  col = 1
443  do j=g%jsc,g%jec ; do i=g%isc,g%iec
444  if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) then
445  cs%col_i(col) = i ; cs%col_j(col) = j
446  cs%Iresttime_col(col) = g%US%T_to_s*iresttime(i,j)
447  col = col +1
448  endif
449  enddo ; enddo
450  endif
451  total_sponge_cols = cs%num_col
452  call sum_across_pes(total_sponge_cols)
453 
454 ! Call the constructor for remapping control structure
455  call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation)
456  call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, &
457  "The total number of columns where sponges are applied at h points.")
458  if (cs%sponge_uv) then
459  allocate(iresttime_u(g%isdB:g%iedB,g%jsd:g%jed)); iresttime_u(:,:)=0.0
460  allocate(iresttime_v(g%isd:g%ied,g%jsdB:g%jedB)); iresttime_v(:,:)=0.0
461  ! u points
462  cs%num_col_u = 0 ; !CS%fldno_u = 0
463  do j=cs%jsc,cs%jec; do i=cs%iscB,cs%iecB
464  iresttime_u(i,j) = 0.5 * (iresttime(i,j) + iresttime(i+1,j))
465  if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) &
466  cs%num_col_u = cs%num_col_u + 1
467  enddo ; enddo
468  if (cs%num_col_u > 0) then
469  allocate(cs%Iresttime_col_u(cs%num_col_u)) ; cs%Iresttime_col_u = 0.0
470  allocate(cs%col_i_u(cs%num_col_u)) ; cs%col_i_u = 0
471  allocate(cs%col_j_u(cs%num_col_u)) ; cs%col_j_u = 0
472  ! pass indices, restoring time to the CS structure
473  col = 1
474  do j=cs%jsc,cs%jec ; do i=cs%iscB,cs%iecB
475  if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) then
476  cs%col_i_u(col) = i ; cs%col_j_u(col) = j
477  cs%Iresttime_col_u(col) = g%US%T_to_s*iresttime_u(i,j)
478  col = col +1
479  endif
480  enddo ; enddo
481  ! same for total number of arbritary layers and correspondent data
482  endif
483  total_sponge_cols_u = cs%num_col_u
484  call sum_across_pes(total_sponge_cols_u)
485  call log_param(param_file, mdl, "!Total sponge columns at u points", total_sponge_cols_u, &
486  "The total number of columns where sponges are applied at u points.")
487  ! v points
488  cs%num_col_v = 0 ; !CS%fldno_v = 0
489  do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec
490  iresttime_v(i,j) = 0.5 * (iresttime(i,j) + iresttime(i,j+1))
491  if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) &
492  cs%num_col_v = cs%num_col_v + 1
493  enddo ; enddo
494  if (cs%num_col_v > 0) then
495  allocate(cs%Iresttime_col_v(cs%num_col_v)) ; cs%Iresttime_col_v = 0.0
496  allocate(cs%col_i_v(cs%num_col_v)) ; cs%col_i_v = 0
497  allocate(cs%col_j_v(cs%num_col_v)) ; cs%col_j_v = 0
498  ! pass indices, restoring time to the CS structure
499  col = 1
500  do j=cs%jscB,cs%jecB ; do i=cs%isc,cs%iec
501  if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) then
502  cs%col_i_v(col) = i ; cs%col_j_v(col) = j
503  cs%Iresttime_col_v(col) = g%US%T_to_s*iresttime_v(i,j)
504  col = col +1
505  endif
506  enddo ; enddo
507  endif
508  total_sponge_cols_v = cs%num_col_v
509  call sum_across_pes(total_sponge_cols_v)
510  call log_param(param_file, mdl, "!Total sponge columns at v points", total_sponge_cols_v, &
511  "The total number of columns where sponges are applied at v points.")
512  endif
513 
514 end subroutine initialize_ale_sponge_varying
515 
516 !> Initialize diagnostics for the ALE_sponge module.
517 ! GMM: this routine is not being used for now.
518 subroutine init_ale_sponge_diags(Time, G, diag, CS)
519  type(time_type), target, intent(in) :: time !< The current model time
520  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
521  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
522  !! output.
523  type(ale_sponge_cs), pointer :: cs !< ALE sponge control structure
524 
525  if (.not.associated(cs)) return
526 
527  cs%diag => diag
528 
529 end subroutine init_ale_sponge_diags
530 
531 !> This subroutine stores the reference profile at h points for the variable
532 !! whose address is given by f_ptr.
533 subroutine set_up_ale_sponge_field_fixed(sp_val, G, f_ptr, CS)
534  type(ocean_grid_type), intent(in) :: G !< Grid structure
535  type(ale_sponge_cs), pointer :: CS !< ALE sponge control structure (in/out).
536  real, dimension(SZI_(G),SZJ_(G),CS%nz_data), &
537  intent(in) :: sp_val !< Field to be used in the sponge, it has arbitrary number of layers.
538  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
539  target, intent(in) :: f_ptr !< Pointer to the field to be damped
540 
541  integer :: j, k, col
542  character(len=256) :: mesg ! String for error messages
543 
544  if (.not.associated(cs)) return
545 
546  cs%fldno = cs%fldno + 1
547 
548  if (cs%fldno > max_fields_) then
549  write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease &
550  &the number of fields to be damped in the call to &
551  &initialize_sponge." )') cs%fldno
552  call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
553  endif
554 
555  ! stores the reference profile
556  allocate(cs%Ref_val(cs%fldno)%p(cs%nz_data,cs%num_col))
557  cs%Ref_val(cs%fldno)%p(:,:) = 0.0
558  do col=1,cs%num_col
559  do k=1,cs%nz_data
560  cs%Ref_val(cs%fldno)%p(k,col) = sp_val(cs%col_i(col),cs%col_j(col),k)
561  enddo
562  enddo
563 
564  cs%var(cs%fldno)%p => f_ptr
565 
566 end subroutine set_up_ale_sponge_field_fixed
567 
568 !> This subroutine stores the reference profile at h points for the variable
569 !! whose address is given by filename and fieldname.
570 subroutine set_up_ale_sponge_field_varying(filename, fieldname, Time, G, GV, US, f_ptr, CS)
571  character(len=*), intent(in) :: filename !< The name of the file with the
572  !! time varying field data
573  character(len=*), intent(in) :: fieldname !< The name of the field in the file
574  !! with the time varying field data
575  type(time_type), intent(in) :: Time !< The current model time
576  type(ocean_grid_type), intent(in) :: G !< Grid structure (in).
577  type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
578  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
579  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
580  target, intent(in) :: f_ptr !< Pointer to the field to be damped (in).
581  type(ale_sponge_cs), pointer :: CS !< Sponge control structure (in/out).
582 
583  ! Local variables
584  real, allocatable, dimension(:,:,:) :: sp_val !< Field to be used in the sponge
585  real, allocatable, dimension(:,:,:) :: mask_z !< Field mask for the sponge data
586  real, allocatable, dimension(:), target :: z_in, z_edges_in ! Heights [Z ~> m].
587  real :: missing_value
588  integer :: j, k, col
589  integer :: isd,ied,jsd,jed
590  integer :: nPoints
591  integer, dimension(4) :: fld_sz
592  integer :: nz_data !< the number of vertical levels in this input field
593  character(len=256) :: mesg ! String for error messages
594  ! Local variables for ALE remapping
595  real, dimension(:), allocatable :: tmpT1d
596  real :: zTopOfCell, zBottomOfCell ! Heights [Z ~> m].
597  type(remapping_cs) :: remapCS ! Remapping parameters and work arrays
598 
599  if (.not.associated(cs)) return
600  ! initialize time interpolator module
601  call time_interp_external_init()
602  isd = g%isd; ied = g%ied; jsd = g%jsd; jed = g%jed
603  cs%fldno = cs%fldno + 1
604  if (cs%fldno > max_fields_) then
605  write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease &
606  &the number of fields to be damped in the call to &
607  &initialize_sponge." )') cs%fldno
608  call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
609  endif
610  ! get a unique time interp id for this field. If sponge data is ongrid, then setup
611  ! to only read on the computational domain
612  if (cs%spongeDataOngrid) then
613  cs%Ref_val(cs%fldno)%id = init_external_field(filename, fieldname,domain=g%Domain%mpp_domain)
614  else
615  cs%Ref_val(cs%fldno)%id = init_external_field(filename, fieldname)
616  endif
617  fld_sz(1:4)=-1
618  fld_sz = get_external_field_size(cs%Ref_val(cs%fldno)%id)
619  nz_data = fld_sz(3)
620  cs%Ref_val(cs%fldno)%nz_data = nz_data !< individual sponge fields may reside on a different vertical grid
621  cs%Ref_val(cs%fldno)%num_tlevs = fld_sz(4)
622  ! initializes the target profile array for this field
623  ! for all columns which will be masked
624  allocate(cs%Ref_val(cs%fldno)%p(nz_data,cs%num_col))
625  cs%Ref_val(cs%fldno)%p(:,:) = 0.0
626  allocate( cs%Ref_val(cs%fldno)%h(nz_data,cs%num_col) )
627  cs%Ref_val(cs%fldno)%h(:,:) = 0.0
628  cs%var(cs%fldno)%p => f_ptr
629 
630 
631 end subroutine set_up_ale_sponge_field_varying
632 
633 !> This subroutine stores the reference profile at u and v points for the variable
634 !! whose address is given by u_ptr and v_ptr.
635 subroutine set_up_ale_sponge_vel_field_fixed(u_val, v_val, G, u_ptr, v_ptr, CS)
636  type(ocean_grid_type), intent(in) :: G !< Grid structure (in).
637  type(ale_sponge_cs), pointer :: CS !< Sponge structure (in/out).
638  real, dimension(SZIB_(G),SZJ_(G),CS%nz_data), &
639  intent(in) :: u_val !< u field to be used in the sponge, it has arbritary number of layers.
640  real, dimension(SZI_(G),SZJB_(G),CS%nz_data), &
641  intent(in) :: v_val !< v field to be used in the sponge, it has arbritary number of layers.
642  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(in) :: u_ptr !< u pointer to the field to be damped
643  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(in) :: v_ptr !< v pointer to the field to be damped
644 
645  integer :: j, k, col
646  character(len=256) :: mesg ! String for error messages
647 
648  if (.not.associated(cs)) return
649 
650  ! stores the reference profile
651  allocate(cs%Ref_val_u%p(cs%nz_data,cs%num_col_u))
652  cs%Ref_val_u%p(:,:) = 0.0
653  do col=1,cs%num_col_u
654  do k=1,cs%nz_data
655  cs%Ref_val_u%p(k,col) = u_val(cs%col_i_u(col),cs%col_j_u(col),k)
656  enddo
657  enddo
658  cs%var_u%p => u_ptr
659  allocate(cs%Ref_val_v%p(cs%nz_data,cs%num_col_v))
660  cs%Ref_val_v%p(:,:) = 0.0
661  do col=1,cs%num_col_v
662  do k=1,cs%nz_data
663  cs%Ref_val_v%p(k,col) = v_val(cs%col_i_v(col),cs%col_j_v(col),k)
664  enddo
665  enddo
666  cs%var_v%p => v_ptr
667 
669 
670 !> This subroutine stores the reference profile at uand v points for the variable
671 !! whose address is given by u_ptr and v_ptr.
672 subroutine set_up_ale_sponge_vel_field_varying(filename_u, fieldname_u, filename_v, fieldname_v, &
673  Time, G, US, CS, u_ptr, v_ptr)
674  character(len=*), intent(in) :: filename_u !< File name for u field
675  character(len=*), intent(in) :: fieldname_u !< Name of u variable in file
676  character(len=*), intent(in) :: filename_v !< File name for v field
677  character(len=*), intent(in) :: fieldname_v !< Name of v variable in file
678  type(time_type), intent(in) :: Time !< Model time
679  type(ocean_grid_type), intent(inout) :: G !< Ocean grid (in)
680  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
681  type(ale_sponge_cs), pointer :: CS !< Sponge structure (in/out).
682  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(in) :: u_ptr !< u pointer to the field to be damped (in).
683  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(in) :: v_ptr !< v pointer to the field to be damped (in).
684  ! Local variables
685  real, allocatable, dimension(:,:,:) :: u_val !< U field to be used in the sponge.
686  real, allocatable, dimension(:,:,:) :: mask_u !< U field mask for the sponge data.
687  real, allocatable, dimension(:,:,:) :: v_val !< V field to be used in the sponge.
688  real, allocatable, dimension(:,:,:) :: mask_v !< V field mask for the sponge data.
689 
690  real, allocatable, dimension(:), target :: z_in, z_edges_in
691  real :: missing_value
692 
693  integer :: j, k, col
694  integer :: isd, ied, jsd, jed
695  integer :: isdB, iedB, jsdB, jedB
696  integer, dimension(4) :: fld_sz
697  character(len=256) :: mesg ! String for error messages
698 
699  if (.not.associated(cs)) return
700 
701  isd = g%isd; ied = g%ied; jsd = g%jsd; jed = g%jed
702  isdb = g%isdB; iedb = g%iedB; jsdb = g%jsdB; jedb = g%jedB
703  ! get a unique id for this field which will allow us to return an array
704  ! containing time-interpolated values from an external file corresponding
705  ! to the current model date.
706  cs%Ref_val_u%id = init_external_field(filename_u, fieldname_u)
707  fld_sz(1:4)=-1
708  fld_sz = get_external_field_size(cs%Ref_val_u%id)
709  cs%Ref_val_u%nz_data = fld_sz(3)
710  cs%Ref_val_u%num_tlevs = fld_sz(4)
711  cs%Ref_val_v%id = init_external_field(filename_v, fieldname_v)
712  fld_sz(1:4)=-1
713  fld_sz = get_external_field_size(cs%Ref_val_v%id)
714  cs%Ref_val_v%nz_data = fld_sz(3)
715  cs%Ref_val_v%num_tlevs = fld_sz(4)
716  allocate( u_val(isdb:iedb,jsd:jed, fld_sz(3)) )
717  allocate( mask_u(isdb:iedb,jsd:jed, fld_sz(3)) )
718  allocate( v_val(isd:ied,jsdb:jedb, fld_sz(3)) )
719  allocate( mask_v(isd:ied,jsdb:jedb, fld_sz(3)) )
720  ! Interpolate external file data to the model grid
721  ! I am hard-wiring this call to assume that the input grid is zonally re-entrant
722  ! In the future, this should be generalized using an interface to return the
723  ! modulo attribute of the zonal axis (mjh).
724  call horiz_interp_and_extrap_tracer(cs%Ref_val_u%id,time, 1.0,g,u_val,mask_u,z_in,z_edges_in,&
725  missing_value,.true.,.false.,.false., m_to_z=us%m_to_Z)
726  !!! TODO: add a velocity interface! (mjh)
727  ! Interpolate external file data to the model grid
728  ! I am hard-wiring this call to assume that the input grid is zonally re-entrant
729  ! In the future, this should be generalized using an interface to return the
730  ! modulo attribute of the zonal axis (mjh).
731  call horiz_interp_and_extrap_tracer(cs%Ref_val_v%id,time, 1.0,g,v_val,mask_v,z_in,z_edges_in, &
732  missing_value,.true.,.false.,.false., m_to_z=us%m_to_Z)
733  ! stores the reference profile
734  allocate(cs%Ref_val_u%p(fld_sz(3),cs%num_col_u))
735  cs%Ref_val_u%p(:,:) = 0.0
736  do col=1,cs%num_col_u
737  do k=1,fld_sz(3)
738  cs%Ref_val_u%p(k,col) = u_val(cs%col_i_u(col),cs%col_j_u(col),k)
739  enddo
740  enddo
741  cs%var_u%p => u_ptr
742  allocate(cs%Ref_val_v%p(fld_sz(3),cs%num_col_v))
743  cs%Ref_val_v%p(:,:) = 0.0
744  do col=1,cs%num_col_v
745  do k=1,fld_sz(3)
746  cs%Ref_val_v%p(k,col) = v_val(cs%col_i_v(col),cs%col_j_v(col),k)
747  enddo
748  enddo
749  cs%var_v%p => v_ptr
750 
752 
753 !> This subroutine applies damping to the layers thicknesses, temp, salt and a variety of tracers
754 !! for every column where there is damping.
755 subroutine apply_ale_sponge(h, dt, G, GV, US, CS, Time)
756  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure (in).
757  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
758  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
759  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
760  intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in)
761  real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
762  type(ale_sponge_cs), pointer :: cs !< A pointer to the control structure for this module
763  !! that is set by a previous call to initialize_sponge (in).
764  type(time_type), optional, intent(in) :: time !< The current model date
765 
766  real :: damp ! The timestep times the local damping coefficient [nondim].
767  real :: i1pdamp ! I1pdamp is 1/(1 + damp). [nondim].
768  real :: m_to_z ! A unit conversion factor from m to Z.
769  real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid
770  real, dimension(SZK_(G)) :: tmp_val1 ! data values remapped to model grid
771  real :: hu(szib_(g), szj_(g), szk_(g)) ! A temporary array for h at u pts
772  real :: hv(szi_(g), szjb_(g), szk_(g)) ! A temporary array for h at v pts
773  real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields
774  real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts
775  real, dimension(:), allocatable :: hsrc ! Source thicknesses [Z ~> m].
776  ! Local variables for ALE remapping
777  real, dimension(:), allocatable :: tmpt1d
778  integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz, nz_data
779  integer :: col, total_sponge_cols
780  real, allocatable, dimension(:), target :: z_in, z_edges_in
781  real :: missing_value
782  real :: h_neglect, h_neglect_edge
783  real :: ztopofcell, zbottomofcell ! Heights [Z ~> m].
784  integer :: npoints
785 
786  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
787  if (.not.associated(cs)) return
788 
789  if (gv%Boussinesq) then
790  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
791  else
792  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
793  endif
794 
795  if (cs%time_varying_sponges) then
796  if (.not. present(time)) &
797  call mom_error(fatal,"apply_ALE_sponge: No time information provided")
798  do m=1,cs%fldno
799  nz_data = cs%Ref_val(m)%nz_data
800  allocate(sp_val(g%isd:g%ied,g%jsd:g%jed,1:nz_data))
801  allocate(mask_z(g%isd:g%ied,g%jsd:g%jed,1:nz_data))
802  sp_val(:,:,:)=0.0
803  mask_z(:,:,:)=0.0
804  call horiz_interp_and_extrap_tracer(cs%Ref_val(m)%id,time, 1.0,g,sp_val,mask_z,z_in,z_edges_in, &
805  missing_value,.true., .false.,.false.,spongeongrid=cs%SpongeDataOngrid, m_to_z=us%m_to_Z)
806  allocate( hsrc(nz_data) )
807  allocate( tmpt1d(nz_data) )
808  do c=1,cs%num_col
809  i = cs%col_i(c) ; j = cs%col_j(c)
810  cs%Ref_val(m)%p(1:nz_data,c) = sp_val(i,j,1:nz_data)
811  ! Build the source grid
812  ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0; hsrc(:) = 0.0; tmpt1d(:) = -99.9
813  do k=1,nz_data
814  if (mask_z(cs%col_i(c),cs%col_j(c),k) == 1.0) then
815  zbottomofcell = -min( z_edges_in(k+1), g%bathyT(cs%col_i(c),cs%col_j(c)) )
816  tmpt1d(k) = sp_val(cs%col_i(c),cs%col_j(c),k)
817  elseif (k>1) then
818  zbottomofcell = -g%bathyT(cs%col_i(c),cs%col_j(c))
819  tmpt1d(k) = tmpt1d(k-1)
820  else ! This next block should only ever be reached over land
821  tmpt1d(k) = -99.9
822  endif
823  hsrc(k) = ztopofcell - zbottomofcell
824  if (hsrc(k)>0.) npoints = npoints + 1
825  ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
826  enddo
827  ! In case data is deeper than model
828  hsrc(nz_data) = hsrc(nz_data) + ( ztopofcell + g%bathyT(cs%col_i(c),cs%col_j(c)) )
829  cs%Ref_val(cs%fldno)%h(1:nz_data,c) = gv%Z_to_H*hsrc(1:nz_data)
830  cs%Ref_val(cs%fldno)%p(1:nz_data,c) = tmpt1d(1:nz_data)
831  do k=2,nz_data
832  ! if (mask_z(i,j,k)==0.) &
833  if (cs%Ref_val(m)%h(k,c) <= 0.001*gv%m_to_H) &
834  ! some confusion here about why the masks are not correct returning from horiz_interp
835  ! reverting to using a minimum thickness criteria
836  cs%Ref_val(m)%p(k,c) = cs%Ref_val(m)%p(k-1,c)
837  enddo
838  enddo
839  deallocate(sp_val, mask_z, hsrc, tmpt1d)
840  enddo
841  else
842  nz_data = cs%nz_data
843  endif
844 
845  allocate(tmp_val2(nz_data))
846  do m=1,cs%fldno
847  do c=1,cs%num_col
848 ! c is an index for the next 3 lines but a multiplier for the rest of the loop
849 ! Therefore we use c as per C code and increment the index where necessary.
850  i = cs%col_i(c) ; j = cs%col_j(c)
851  damp = dt * cs%Iresttime_col(c)
852  i1pdamp = 1.0 / (1.0 + damp)
853  tmp_val2(1:nz_data) = cs%Ref_val(m)%p(1:nz_data,c)
854  if (cs%time_varying_sponges) then
855  call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_val(m)%h(1:nz_data,c), tmp_val2, &
856  cs%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
857  else
858  call remapping_core_h(cs%remap_cs,nz_data, cs%Ref_h%p(1:nz_data,c), tmp_val2, &
859  cs%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
860  endif
861  !Backward Euler method
862  cs%var(m)%p(i,j,1:cs%nz) = i1pdamp * (cs%var(m)%p(i,j,1:cs%nz) + tmp_val1 * damp)
863  enddo
864  enddo
865 
866  ! for debugging
867  !c=CS%num_col
868  !do m=1,CS%fldno
869  ! write(*,*) 'APPLY SPONGE,m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1',&
870  ! m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1
871  !enddo
872 
873  if (cs%sponge_uv) then
874  ! u points
875  do j=cs%jsc,cs%jec; do i=cs%iscB,cs%iecB; do k=1,nz
876  hu(i,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
877  enddo ; enddo ; enddo
878  if (cs%time_varying_sponges) then
879  if (.not. present(time)) &
880  call mom_error(fatal,"apply_ALE_sponge: No time information provided")
881 
882  nz_data = cs%Ref_val_u%nz_data
883  allocate(sp_val(g%isdB:g%iedB,g%jsd:g%jed,1:nz_data))
884  allocate(mask_z(g%isdB:g%iedB,g%jsd:g%jed,1:nz_data))
885  ! Interpolate from the external horizontal grid and in time
886  call horiz_interp_and_extrap_tracer(cs%Ref_val_u%id,time, 1.0,g,sp_val,mask_z,z_in,z_edges_in, &
887  missing_value, .true., .false., .false., m_to_z=us%m_to_Z)
888 ! call pass_var(sp_val,G%Domain)
889 ! call pass_var(mask_z,G%Domain)
890  do c=1,cs%num_col
891  ! c is an index for the next 3 lines but a multiplier for the rest of the loop
892  ! Therefore we use c as per C code and increment the index where necessary.
893  i = cs%col_i(c) ; j = cs%col_j(c)
894  cs%Ref_val_u%p(1:nz_data,c) = sp_val(i,j,1:nz_data)
895  enddo
896 
897  deallocate (sp_val, mask_z)
898 
899  nz_data = cs%Ref_val_v%nz_data
900  allocate(sp_val(g%isd:g%ied,g%jsdB:g%jedB,1:nz_data))
901  allocate(mask_z(g%isd:g%ied,g%jsdB:g%jedB,1:nz_data))
902  ! Interpolate from the external horizontal grid and in time
903  call horiz_interp_and_extrap_tracer(cs%Ref_val_v%id,time, 1.0,g,sp_val,mask_z,z_in,z_edges_in, &
904  missing_value, .true., .false., .false., m_to_z=us%m_to_Z)
905 
906 ! call pass_var(sp_val,G%Domain)
907 ! call pass_var(mask_z,G%Domain)
908 
909  do c=1,cs%num_col
910  ! c is an index for the next 3 lines but a multiplier for the rest of the loop
911  ! Therefore we use c as per C code and increment the index where necessary.
912  i = cs%col_i(c) ; j = cs%col_j(c)
913  cs%Ref_val_v%p(1:nz_data,c) = sp_val(i,j,1:nz_data)
914  enddo
915 
916  deallocate (sp_val, mask_z)
917 
918  else
919  nz_data = cs%nz_data
920  endif
921 
922  do c=1,cs%num_col_u
923  i = cs%col_i_u(c) ; j = cs%col_j_u(c)
924  damp = dt * cs%Iresttime_col_u(c)
925  i1pdamp = 1.0 / (1.0 + damp)
926  if (cs%time_varying_sponges) nz_data = cs%Ref_val(m)%nz_data
927  tmp_val2(1:nz_data) = cs%Ref_val_u%p(1:nz_data,c)
928  if (cs%time_varying_sponges) then
929  call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_val_u%h(:,c), tmp_val2, &
930  cs%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
931  else
932  call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_hu%p(:,c), tmp_val2, &
933  cs%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
934  endif
935  !Backward Euler method
936  cs%var_u%p(i,j,:) = i1pdamp * (cs%var_u%p(i,j,:) + tmp_val1 * damp)
937  enddo
938 
939  ! v points
940  do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec; do k=1,nz
941  hv(i,j,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
942  enddo ; enddo ; enddo
943 
944  do c=1,cs%num_col_v
945  i = cs%col_i_v(c) ; j = cs%col_j_v(c)
946  damp = dt * cs%Iresttime_col_v(c)
947  i1pdamp = 1.0 / (1.0 + damp)
948  tmp_val2(1:nz_data) = cs%Ref_val_v%p(1:nz_data,c)
949  if (cs%time_varying_sponges) then
950  call remapping_core_h(cs%remap_cs, cs%nz_data, cs%Ref_val_v%h(:,c), tmp_val2, &
951  cs%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
952  else
953  call remapping_core_h(cs%remap_cs, cs%nz_data, cs%Ref_hv%p(:,c), tmp_val2, &
954  cs%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
955  endif
956  !Backward Euler method
957  cs%var_v%p(i,j,:) = i1pdamp * (cs%var_v%p(i,j,:) + tmp_val1 * damp)
958  enddo
959 
960  endif
961 
962  deallocate(tmp_val2)
963 
964 end subroutine apply_ale_sponge
965 
966 ! GMM: I could not find where sponge_end is being called, but I am keeping
967 ! ALE_sponge_end here so we can add that if needed.
968 !> This subroutine deallocates any memory associated with the ALE_sponge module.
969 subroutine ale_sponge_end(CS)
970  type(ale_sponge_cs), pointer :: cs !< A pointer to the control structure that is
971  !! set by a previous call to initialize_sponge.
972 
973  integer :: m
974 
975  if (.not.associated(cs)) return
976 
977  if (associated(cs%col_i)) deallocate(cs%col_i)
978  if (associated(cs%col_i_u)) deallocate(cs%col_i_u)
979  if (associated(cs%col_i_v)) deallocate(cs%col_i_v)
980  if (associated(cs%col_j)) deallocate(cs%col_j)
981  if (associated(cs%col_j_u)) deallocate(cs%col_j_u)
982  if (associated(cs%col_j_v)) deallocate(cs%col_j_v)
983 
984  if (associated(cs%Iresttime_col)) deallocate(cs%Iresttime_col)
985  if (associated(cs%Iresttime_col_u)) deallocate(cs%Iresttime_col_u)
986  if (associated(cs%Iresttime_col_v)) deallocate(cs%Iresttime_col_v)
987 
988  do m=1,cs%fldno
989  if (associated(cs%Ref_val(cs%fldno)%p)) deallocate(cs%Ref_val(cs%fldno)%p)
990  enddo
991 
992  deallocate(cs)
993 
994 end subroutine ale_sponge_end
995 
996 end module mom_ale_sponge
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
mom_spatial_means::global_i_mean
subroutine, public global_i_mean(array, i_mean, G, mask, scale, tmp_scale)
Determine the global mean of a field along rows of constant i, returning it in a 1-d array using the ...
Definition: MOM_spatial_means.F90:186
mom_spatial_means
Functions and routines to take area, volume, mass-weighted, layerwise, zonal or meridional means.
Definition: MOM_spatial_means.F90:2
mom_remapping::remapping_core_h
subroutine, public remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned.
Definition: MOM_remapping.F90:188
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Call this subroutine to determine whether the averaging is currently enabled. .true....
Definition: MOM_diag_mediator.F90:1850
mom_horizontal_regridding::horiz_interp_and_extrap_tracer
Extrapolate and interpolate data.
Definition: MOM_horizontal_regridding.F90:52
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_ale_sponge::ale_sponge_end
subroutine, public ale_sponge_end(CS)
This subroutine deallocates any memory associated with the ALE_sponge module.
Definition: MOM_ALE_sponge.F90:970
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_ale_sponge::initialize_ale_sponge
Ddetermine the number of points which are within sponges in this computational domain.
Definition: MOM_ALE_sponge.F90:48
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_ale_sponge::set_up_ale_sponge_vel_field_varying
subroutine set_up_ale_sponge_vel_field_varying(filename_u, fieldname_u, filename_v, fieldname_v, Time, G, US, CS, u_ptr, v_ptr)
This subroutine stores the reference profile at uand v points for the variable whose address is given...
Definition: MOM_ALE_sponge.F90:674
mom_ale_sponge::ale_sponge_cs
ALE sponge control structure.
Definition: MOM_ALE_sponge.F90:84
mom_remapping::remapping_cs
Container for remapping parameters.
Definition: MOM_remapping.F90:22
mom_ale_sponge::get_ale_sponge_nz_data
integer function, public get_ale_sponge_nz_data(CS)
Return the number of layers in the data with a fixed ALE sponge, or 0 if there are no sponge columns ...
Definition: MOM_ALE_sponge.F90:326
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_ale_sponge::set_up_ale_sponge_vel_field_fixed
subroutine set_up_ale_sponge_vel_field_fixed(u_val, v_val, G, u_ptr, v_ptr, CS)
This subroutine stores the reference profile at u and v points for the variable whose address is give...
Definition: MOM_ALE_sponge.F90:636
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_ale_sponge
This module contains the routines used to apply sponge layers when using the ALE mode.
Definition: MOM_ALE_sponge.F90:11
mom_horizontal_regridding
Horizontal interpolation.
Definition: MOM_horizontal_regridding.F90:2
mom_ale_sponge::p2d
A structure for creating arrays of pointers to 2D arrays with extra gridding information.
Definition: MOM_ALE_sponge.F90:74
mom_ale_sponge::apply_ale_sponge
subroutine, public apply_ale_sponge(h, dt, G, GV, US, CS, Time)
This subroutine applies damping to the layers thicknesses, temp, salt and a variety of tracers for ev...
Definition: MOM_ALE_sponge.F90:756
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_coms
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
Definition: MOM_coms.F90:3
mom_remapping
Provides column-wise vertical remapping functions.
Definition: MOM_remapping.F90:2
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_ale_sponge::init_ale_sponge_diags
subroutine, public init_ale_sponge_diags(Time, G, diag, CS)
Initialize diagnostics for the ALE_sponge module.
Definition: MOM_ALE_sponge.F90:519
mom_ale_sponge::set_up_ale_sponge_field
Store the reference profile at h points for a variable.
Definition: MOM_ALE_sponge.F90:33
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_ale_sponge::set_up_ale_sponge_field_varying
subroutine set_up_ale_sponge_field_varying(filename, fieldname, Time, G, GV, US, f_ptr, CS)
This subroutine stores the reference profile at h points for the variable whose address is given by f...
Definition: MOM_ALE_sponge.F90:571
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_ale_sponge::initialize_ale_sponge_fixed
subroutine initialize_ale_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_data)
This subroutine determines the number of points which are within sponges in this computational domain...
Definition: MOM_ALE_sponge.F90:141
mom_ale_sponge::p3d
A structure for creating arrays of pointers to 3D arrays with extra gridding information.
Definition: MOM_ALE_sponge.F90:64
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_ale_sponge::get_ale_sponge_thicknesses
subroutine, public get_ale_sponge_thicknesses(G, data_h, sponge_mask, CS)
Return the thicknesses used for the data with a fixed ALE sponge.
Definition: MOM_ALE_sponge.F90:339
mom_remapping::initialize_remapping
subroutine, public initialize_remapping(CS, remapping_scheme, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
Constructor for remapping control structure.
Definition: MOM_remapping.F90:1547
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_ale_sponge::initialize_ale_sponge_varying
subroutine initialize_ale_sponge_varying(Iresttime, G, param_file, CS)
This subroutine determines the number of points which are to be restoref in the computational domain....
Definition: MOM_ALE_sponge.F90:376
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_ale_sponge::set_up_ale_sponge_vel_field
This subroutine stores the reference profile at u and v points for a vector.
Definition: MOM_ALE_sponge.F90:39
mom_ale_sponge::set_up_ale_sponge_field_fixed
subroutine set_up_ale_sponge_field_fixed(sp_val, G, f_ptr, CS)
This subroutine stores the reference profile at h points for the variable whose address is given by f...
Definition: MOM_ALE_sponge.F90:534