MOM6
MOM_ALE.F90
Go to the documentation of this file.
1 !> This module contains the main regridding routines.
2 !!
3 !! Regridding comprises two steps:
4 !! 1. Interpolation and creation of a new grid based on target interface
5 !! densities (or any other criterion).
6 !! 2. Remapping of quantities between old grid and new grid.
7 !!
8 !! Original module written by Laurent White, 2008.06.09
9 module mom_ale
10 
11 ! This file is part of MOM6. See LICENSE.md for the license.
12 
15 use mom_diag_mediator, only : time_type, diag_update_remap_grids
17 use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
18 use mom_eos, only : calculate_density
19 use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
20 use mom_error_handler, only : mom_error, fatal, warning
24 use mom_io, only : vardesc, var_desc, fieldtype, single_file
25 use mom_io, only : create_file, write_field, close_file
33 use mom_regridding, only : regriddingcoordinatemodedoc, default_coordinate_mode
48 use mom_variables, only : ocean_grid_type, thermo_var_ptrs
50 
51 !use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE
58 
59 
60 implicit none ; private
61 #include <MOM_memory.h>
62 
63 
64 !> ALE control structure
65 type, public :: ale_cs ; private
66  logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z"
67  !! method. If False, uses the new method that
68  !! remaps between grids described by h.
69 
70  real :: regrid_time_scale !< The time-scale used in blending between the current (old) grid
71  !! and the target (new) grid [T ~> s]
72 
73  type(regridding_cs) :: regridcs !< Regridding parameters and work arrays
74  type(remapping_cs) :: remapcs !< Remapping parameters and work arrays
75 
76  integer :: nk !< Used only for queries, not directly by this module
77 
78  logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state.
79 
80  logical :: answers_2018 !< If true, use the order of arithmetic and expressions for remapping
81  !! that recover the answers from the end of 2018. Otherwise, use more
82  !! robust and accurate forms of mathematically equivalent expressions.
83 
84  logical :: show_call_tree !< For debugging
85 
86  ! for diagnostics
87  type(diag_ctrl), pointer :: diag !< structure to regulate output
88  integer, dimension(:), allocatable :: id_tracer_remap_tendency !< diagnostic id
89  integer, dimension(:), allocatable :: id_htracer_remap_tendency !< diagnostic id
90  integer, dimension(:), allocatable :: id_htracer_remap_tendency_2d !< diagnostic id
91  logical, dimension(:), allocatable :: do_tendency_diag !< flag for doing diagnostics
92  integer :: id_dzregrid = -1 !< diagnostic id
93 
94  ! diagnostic for fields prior to applying ALE remapping
95  integer :: id_u_preale = -1 !< diagnostic id for zonal velocity before ALE.
96  integer :: id_v_preale = -1 !< diagnostic id for meridional velocity before ALE.
97  integer :: id_h_preale = -1 !< diagnostic id for layer thicknesses before ALE.
98  integer :: id_t_preale = -1 !< diagnostic id for temperatures before ALE.
99  integer :: id_s_preale = -1 !< diagnostic id for salinities before ALE.
100  integer :: id_e_preale = -1 !< diagnostic id for interface heights before ALE.
101  integer :: id_vert_remap_h = -1 !< diagnostic id for layer thicknesses used for remapping
102  integer :: id_vert_remap_h_tendency = -1 !< diagnostic id for layer thickness tendency due to ALE
103 
104 end type
105 
106 ! Publicly available functions
107 public ale_init
108 public ale_end
109 public ale_main
110 public ale_main_offline
111 public ale_offline_inputs
113 public ale_build_grid
115 public ale_remap_scalar
119 public ale_initregridding
120 public ale_getcoordinate
127 public ale_register_diags
128 
129 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
130 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
131 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
132 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
133 
134 contains
135 
136 !> This routine is typically called (from initialize_MOM in file MOM.F90)
137 !! before the main time integration loop to initialize the regridding stuff.
138 !! We read the MOM_input file to register the values of different
139 !! regridding/remapping parameters.
140 subroutine ale_init( param_file, GV, US, max_depth, CS)
141  type(param_file_type), intent(in) :: param_file !< Parameter file
142  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
143  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
144  real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
145  type(ale_cs), pointer :: cs !< Module control structure
146 
147  ! Local variables
148  real, dimension(:), allocatable :: dz
149  character(len=40) :: mdl = "MOM_ALE" ! This module's name.
150  character(len=80) :: string ! Temporary strings
151  real :: filter_shallow_depth, filter_deep_depth
152  logical :: default_2018_answers
153  logical :: check_reconstruction
154  logical :: check_remapping
155  logical :: force_bounds_in_subcell
156  logical :: local_logical
157  logical :: remap_boundary_extrap
158 
159  if (associated(cs)) then
160  call mom_error(warning, "ALE_init called with an associated "// &
161  "control structure.")
162  return
163  endif
164  allocate(cs)
165 
166  cs%show_call_tree = calltree_showquery()
167  if (cs%show_call_tree) call calltree_enter("ALE_init(), MOM_ALE.F90")
168 
169  call get_param(param_file, mdl, "REMAP_UV_USING_OLD_ALG", &
170  cs%remap_uv_using_old_alg, &
171  "If true, uses the old remapping-via-a-delta-z method for "//&
172  "remapping u and v. If false, uses the new method that remaps "//&
173  "between grids described by an old and new thickness.", &
174  default=.true.)
175 
176  ! Initialize and configure regridding
177  call ale_initregridding(gv, us, max_depth, param_file, mdl, cs%regridCS)
178 
179  ! Initialize and configure remapping
180  call get_param(param_file, mdl, "REMAPPING_SCHEME", string, &
181  "This sets the reconstruction scheme used "//&
182  "for vertical remapping for all variables. "//&
183  "It can be one of the following schemes: "//&
184  trim(remappingschemesdoc), default=remappingdefaultscheme)
185  call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
186  "If true, cell-by-cell reconstructions are checked for "//&
187  "consistency and if non-monotonicity or an inconsistency is "//&
188  "detected then a FATAL error is issued.", default=.false.)
189  call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", check_remapping, &
190  "If true, the results of remapping are checked for "//&
191  "conservation and new extrema and if an inconsistency is "//&
192  "detected then a FATAL error is issued.", default=.false.)
193  call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
194  "If true, the values on the intermediate grid used for remapping "//&
195  "are forced to be bounded, which might not be the case due to "//&
196  "round off.", default=.false.)
197  call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
198  "If true, values at the interfaces of boundary cells are "//&
199  "extrapolated instead of piecewise constant", default=.false.)
200  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
201  "This sets the default value for the various _2018_ANSWERS parameters.", &
202  default=.true.)
203  call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", cs%answers_2018, &
204  "If true, use the order of arithmetic and expressions that recover the "//&
205  "answers from the end of 2018. Otherwise, use updated and more robust "//&
206  "forms of the same expressions.", default=default_2018_answers)
207  call initialize_remapping( cs%remapCS, string, &
208  boundary_extrapolation=remap_boundary_extrap, &
209  check_reconstruction=check_reconstruction, &
210  check_remapping=check_remapping, &
211  force_bounds_in_subcell=force_bounds_in_subcell, &
212  answers_2018=cs%answers_2018)
213 
214  call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", cs%remap_after_initialization, &
215  "If true, applies regridding and remapping immediately after "//&
216  "initialization so that the state is ALE consistent. This is a "//&
217  "legacy step and should not be needed if the initialization is "//&
218  "consistent with the coordinate mode.", default=.true.)
219 
220  call get_param(param_file, mdl, "REGRID_TIME_SCALE", cs%regrid_time_scale, &
221  "The time-scale used in blending between the current (old) grid "//&
222  "and the target (new) grid. A short time-scale favors the target "//&
223  "grid (0. or anything less than DT_THERM) has no memory of the old "//&
224  "grid. A very long time-scale makes the model more Lagrangian.", &
225  units="s", default=0., scale=us%s_to_T)
226  call get_param(param_file, mdl, "REGRID_FILTER_SHALLOW_DEPTH", filter_shallow_depth, &
227  "The depth above which no time-filtering is applied. Above this depth "//&
228  "final grid exactly matches the target (new) grid.", &
229  units="m", default=0., scale=gv%m_to_H)
230  call get_param(param_file, mdl, "REGRID_FILTER_DEEP_DEPTH", filter_deep_depth, &
231  "The depth below which full time-filtering is applied with time-scale "//&
232  "REGRID_TIME_SCALE. Between depths REGRID_FILTER_SHALLOW_DEPTH and "//&
233  "REGRID_FILTER_SHALLOW_DEPTH the filter weights adopt a cubic profile.", &
234  units="m", default=0., scale=gv%m_to_H)
235  call set_regrid_params(cs%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &
236  depth_of_time_filter_deep=filter_deep_depth)
237  call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, &
238  "If true, the regridding ntegrates upwards from the bottom for "//&
239  "interface positions, much as the main model does. If false "//&
240  "regridding integrates downward, consistant with the remapping "//&
241  "code.", default=.true., do_not_log=.true.)
242  call set_regrid_params(cs%regridCS, integrate_downward_for_e=.not.local_logical)
243 
244  ! Keep a record of values for subsequent queries
245  cs%nk = gv%ke
246 
247  if (cs%show_call_tree) call calltree_leave("ALE_init()")
248 end subroutine ale_init
249 
250 !> Initialize diagnostics for the ALE module.
251 subroutine ale_register_diags(Time, G, GV, US, diag, CS)
252  type(time_type),target, intent(in) :: time !< Time structure
253  type(ocean_grid_type), intent(in) :: g !< Grid structure
254  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
255  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
256  type(diag_ctrl), target, intent(in) :: diag !< Diagnostics control structure
257  type(ale_cs), pointer :: cs !< Module control structure
258 
259  cs%diag => diag
260 
261  ! These diagnostics of the state variables before ALE are useful for
262  ! debugging the ALE code.
263  cs%id_u_preale = register_diag_field('ocean_model', 'u_preale', diag%axesCuL, time, &
264  'Zonal velocity before remapping', 'm s-1', conversion=us%L_T_to_m_s)
265  cs%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, time, &
266  'Meridional velocity before remapping', 'm s-1', conversion=us%L_T_to_m_s)
267  cs%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, time, &
268  'Layer Thickness before remapping', get_thickness_units(gv), &
269  conversion=gv%H_to_MKS, v_extensive=.true.)
270  cs%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, time, &
271  'Temperature before remapping', 'degC')
272  cs%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, time, &
273  'Salinity before remapping', 'PSU')
274  cs%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, time, &
275  'Interface Heights before remapping', 'm', conversion=us%Z_to_m)
276 
277  cs%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,time, &
278  'Change in interface height due to ALE regridding', 'm', &
279  conversion=gv%H_to_m)
280  cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', &
281  diag%axestl, time, 'layer thicknesses after ALE regridding and remapping', 'm', &
282  conversion=gv%H_to_m, v_extensive=.true.)
283  cs%id_vert_remap_h_tendency = register_diag_field('ocean_model','vert_remap_h_tendency',diag%axestl,time, &
284  'Layer thicknesses tendency due to ALE regridding and remapping', 'm', &
285  conversion=gv%H_to_m*us%s_to_T, v_extensive = .true.)
286 
287 end subroutine ale_register_diags
288 
289 !> Crudely adjust (initial) grid for integrity.
290 !! This routine is typically called (from initialize_MOM in file MOM.F90)
291 !! before the main time integration loop to initialize the regridding stuff.
292 !! We read the MOM_input file to register the values of different
293 !! regridding/remapping parameters.
294 subroutine adjustgridforintegrity( CS, G, GV, h )
295  type(ale_cs), pointer :: cs !< Regridding parameters and options
296  type(ocean_grid_type), intent(in) :: g !< Ocean grid informations
297  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
298  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid thickness that
299  !! are to be adjusted [H ~> m or kg-2]
300  call inflate_vanished_layers_old( cs%regridCS, g, gv, h(:,:,:) )
301 
302 end subroutine adjustgridforintegrity
303 
304 
305 !> End of regridding (memory deallocation).
306 !! This routine is typically called (from MOM_end in file MOM.F90)
307 !! after the main time integration loop to deallocate the regridding stuff.
308 subroutine ale_end(CS)
309  type(ale_cs), pointer :: cs !< module control structure
310 
311  ! Deallocate memory used for the regridding
312  call end_remapping( cs%remapCS )
313  call end_regridding( cs%regridCS )
314 
315  deallocate(cs)
316 
317 end subroutine ale_end
318 
319 !> Takes care of (1) building a new grid and (2) remapping all variables between
320 !! the old grid and the new grid. The creation of the new grid can be based
321 !! on z coordinates, target interface densities, sigma coordinates or any
322 !! arbitrary coordinate system.
323 subroutine ale_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
324  type(ocean_grid_type), intent(in) :: g !< Ocean grid informations
325  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
326  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
327  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the
328  !! last time step [H ~> m or kg m-2]
329  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocity field [L T-1 ~> m s-1]
330  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< Meridional velocity field [L T-1 ~> m s-1]
331  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
332  type(tracer_registry_type), pointer :: reg !< Tracer registry structure
333  type(ale_cs), pointer :: cs !< Regridding parameters and options
334  type(ocean_obc_type), pointer :: obc !< Open boundary structure
335  real, optional, intent(in) :: dt !< Time step between calls to ALE_main [T ~> s]
336  real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
337  ! Local variables
338  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
339  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale
340  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg-2]
341  integer :: nk, i, j, k, isc, iec, jsc, jec
342  logical :: ice_shelf
343 
344  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
345 
346  ice_shelf = .false.
347  if (present(frac_shelf_h)) then
348  if (associated(frac_shelf_h)) ice_shelf = .true.
349  endif
350 
351  if (cs%show_call_tree) call calltree_enter("ALE_main(), MOM_ALE.F90")
352 
353  ! These diagnostics of the state before ALE is applied are mostly used for debugging.
354  if (cs%id_u_preale > 0) call post_data(cs%id_u_preale, u, cs%diag)
355  if (cs%id_v_preale > 0) call post_data(cs%id_v_preale, v, cs%diag)
356  if (cs%id_h_preale > 0) call post_data(cs%id_h_preale, h, cs%diag)
357  if (cs%id_T_preale > 0) call post_data(cs%id_T_preale, tv%T, cs%diag)
358  if (cs%id_S_preale > 0) call post_data(cs%id_S_preale, tv%S, cs%diag)
359  if (cs%id_e_preale > 0) then
360  call find_eta(h, tv, g, gv, us, eta_preale)
361  call post_data(cs%id_e_preale, eta_preale, cs%diag)
362  endif
363 
364  if (present(dt)) then
365  call ale_update_regrid_weights( dt, cs )
366  endif
367  dzregrid(:,:,:) = 0.0
368 
369  ! Build new grid. The new grid is stored in h_new. The old grid is h.
370  ! Both are needed for the subsequent remapping of variables.
371  if (ice_shelf) then
372  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, frac_shelf_h)
373  else
374  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid)
375  endif
376 
377  call check_grid( g, gv, h, 0. )
378 
379  if (cs%show_call_tree) call calltree_waypoint("new grid generated (ALE_main)")
380 
381  ! The presence of dt is used for expediency to distinguish whether ALE_main is being called during init
382  ! or in the main loop. Tendency diagnostics in remap_all_state_vars also rely on this logic.
383  if (present(dt)) then
384  call diag_update_remap_grids(cs%diag)
385  endif
386  ! Remap all variables from old grid h onto new grid h_new
387  call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, obc, -dzregrid, &
388  u, v, cs%show_call_tree, dt )
389 
390  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_main)")
391 
392  ! Override old grid with new one. The new grid 'h_new' is built in
393  ! one of the 'build_...' routines above.
394  !$OMP parallel do default(shared)
395  do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1
396  h(i,j,k) = h_new(i,j,k)
397  enddo ; enddo ; enddo
398 
399  if (cs%show_call_tree) call calltree_leave("ALE_main()")
400 
401  if (cs%id_dzRegrid>0 .and. present(dt)) call post_data(cs%id_dzRegrid, dzregrid, cs%diag)
402 
403 
404 end subroutine ale_main
405 
406 !> Takes care of (1) building a new grid and (2) remapping all variables between
407 !! the old grid and the new grid. The creation of the new grid can be based
408 !! on z coordinates, target interface densities, sigma coordinates or any
409 !! arbitrary coordinate system.
410 subroutine ale_main_offline( G, GV, h, tv, Reg, CS, OBC, dt)
411  type(ocean_grid_type), intent(in) :: g !< Ocean grid informations
412  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
413  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the
414  !! last time step [H ~> m or kg-2]
415  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
416  type(tracer_registry_type), pointer :: reg !< Tracer registry structure
417  type(ale_cs), pointer :: cs !< Regridding parameters and options
418  type(ocean_obc_type), pointer :: obc !< Open boundary structure
419  real, optional, intent(in) :: dt !< Time step between calls to ALE_main [T ~> s]
420  ! Local variables
421  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
422  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg-2]
423  integer :: nk, i, j, k, isc, iec, jsc, jec
424 
425  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
426 
427  if (cs%show_call_tree) call calltree_enter("ALE_main_offline(), MOM_ALE.F90")
428 
429  if (present(dt)) then
430  call ale_update_regrid_weights( dt, cs )
431  endif
432  dzregrid(:,:,:) = 0.0
433 
434  ! Build new grid. The new grid is stored in h_new. The old grid is h.
435  ! Both are needed for the subsequent remapping of variables.
436  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid )
437 
438  call check_grid( g, gv, h, 0. )
439 
440  if (cs%show_call_tree) call calltree_waypoint("new grid generated (ALE_main)")
441 
442  ! Remap all variables from old grid h onto new grid h_new
443 
444  call remap_all_state_vars(cs%remapCS, cs, g, gv, h, h_new, reg, obc, &
445  debug=cs%show_call_tree, dt=dt )
446 
447  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_main)")
448 
449  ! Override old grid with new one. The new grid 'h_new' is built in
450  ! one of the 'build_...' routines above.
451  !$OMP parallel do default(shared)
452  do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1
453  h(i,j,k) = h_new(i,j,k)
454  enddo ; enddo ; enddo
455 
456  if (cs%show_call_tree) call calltree_leave("ALE_main()")
457  if (cs%id_dzRegrid>0 .and. present(dt)) call post_data(cs%id_dzRegrid, dzregrid, cs%diag)
458 
459 end subroutine ale_main_offline
460 
461 !> Regrid/remap stored fields used for offline tracer integrations. These input fields are assumed to have
462 !! the same layer thicknesses at the end of the last offline interval (which should be a Zstar grid). This
463 !! routine builds a grid on the runtime specified vertical coordinate
464 subroutine ale_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
465  type(ale_cs), pointer :: cs !< Regridding parameters and options
466  type(ocean_grid_type), intent(in ) :: g !< Ocean grid informations
467  type(verticalgrid_type), intent(in ) :: gv !< Ocean vertical grid structure
468  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses
469  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
470  type(tracer_registry_type), pointer :: reg !< Tracer registry structure
471  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Zonal mass fluxes
472  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Meridional mass fluxes
473  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: kd !< Input diffusivites
474  logical, intent(in ) :: debug !< If true, then turn checksums
475  type(ocean_obc_type), pointer :: obc !< Open boundary structure
476  ! Local variables
477  integer :: nk, i, j, k, isc, iec, jsc, jec
478  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! Layer thicknesses after regridding
479  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
480  real, dimension(SZK_(GV)) :: h_src
481  real, dimension(SZK_(GV)) :: h_dest, uh_dest
482  real, dimension(SZK_(GV)) :: temp_vec
483 
484  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
485  dzregrid(:,:,:) = 0.0
486  h_new(:,:,:) = 0.0
487 
488  if (debug) call mom_tracer_chkinv("Before ALE_offline_inputs", g, h, reg%Tr, reg%ntr)
489 
490  ! Build new grid from the Zstar state onto the requested vertical coordinate. The new grid is stored
491  ! in h_new. The old grid is h. Both are needed for the subsequent remapping of variables. Convective
492  ! adjustment right now is not used because it is unclear what to do with vanished layers
493  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, conv_adjust = .false. )
494  call check_grid( g, gv, h_new, 0. )
495  if (cs%show_call_tree) call calltree_waypoint("new grid generated (ALE_offline_inputs)")
496 
497  ! Remap all variables from old grid h onto new grid h_new
498  call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, obc, debug=cs%show_call_tree )
499  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_inputs)")
500 
501  ! Reintegrate mass transports from Zstar to the offline vertical coordinate
502  do j=jsc,jec ; do i=g%iscB,g%iecB
503  if (g%mask2dCu(i,j)>0.) then
504  h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
505  h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i+1,j,:))
506  call reintegrate_column(nk, h_src, uhtr(i,j,:), nk, h_dest, 0., temp_vec)
507  uhtr(i,j,:) = temp_vec
508  endif
509  enddo ; enddo
510  do j=g%jscB,g%jecB ; do i=isc,iec
511  if (g%mask2dCv(i,j)>0.) then
512  h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
513  h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i,j+1,:))
514  call reintegrate_column(nk, h_src, vhtr(i,j,:), nk, h_dest, 0., temp_vec)
515  vhtr(i,j,:) = temp_vec
516  endif
517  enddo ; enddo
518 
519  do j = jsc,jec ; do i=isc,iec
520  if (g%mask2dT(i,j)>0.) then
521  if (check_column_integrals(nk, h_src, nk, h_dest)) then
522  call mom_error(fatal, "ALE_offline_inputs: Kd interpolation columns do not match")
523  endif
524  call interpolate_column(nk, h(i,j,:), kd(i,j,:), nk, h_new(i,j,:), 0., kd(i,j,:))
525  endif
526  enddo ; enddo
527 
528  call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%T, h_new, tv%T)
529  call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%S, h_new, tv%S)
530 
531  if (debug) call mom_tracer_chkinv("After ALE_offline_inputs", g, h_new, reg%Tr, reg%ntr)
532 
533  ! Copy over the new layer thicknesses
534  do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1
535  h(i,j,k) = h_new(i,j,k)
536  enddo ; enddo ; enddo
537 
538  if (cs%show_call_tree) call calltree_leave("ALE_offline_inputs()")
539 end subroutine ale_offline_inputs
540 
541 
542 !> Remaps all tracers from h onto h_target. This is intended to be called when tracers
543 !! are done offline. In the case where transports don't quite conserve, we still want to
544 !! make sure that layer thicknesses offline do not drift too far away from the online model
545 subroutine ale_offline_tracer_final( G, GV, h, tv, h_target, Reg, CS, OBC)
546  type(ocean_grid_type), intent(in) :: g !< Ocean grid informations
547  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
548  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the
549  !! last time step [H ~> m or kg-2]
550  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
551  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_target !< Current 3D grid obtained after
552  !! last time step [H ~> m or kg-2]
553  type(tracer_registry_type), pointer :: reg !< Tracer registry structure
554  type(ale_cs), pointer :: cs !< Regridding parameters and options
555  type(ocean_obc_type), pointer :: obc !< Open boundary structure
556  ! Local variables
557 
558  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid !< The change in grid interface positions
559  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new !< Regridded target thicknesses
560  integer :: nk, i, j, k, isc, iec, jsc, jec
561 
562  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
563 
564  if (cs%show_call_tree) call calltree_enter("ALE_offline_tracer_final(), MOM_ALE.F90")
565  ! Need to make sure that h_target is consistent with the current offline ALE confiuration
566  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h_target, tv, h_new, dzregrid )
567  call check_grid( g, gv, h_target, 0. )
568 
569 
570  if (cs%show_call_tree) call calltree_waypoint("Source and target grids checked (ALE_offline_tracer_final)")
571 
572  ! Remap all variables from old grid h onto new grid h_new
573 
574  call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, obc, debug=cs%show_call_tree )
575 
576  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_offline_tracer_final)")
577 
578  ! Override old grid with new one. The new grid 'h_new' is built in
579  ! one of the 'build_...' routines above.
580  !$OMP parallel do default(shared)
581  do k = 1,nk
582  do j = jsc-1,jec+1 ; do i = isc-1,iec+1
583  h(i,j,k) = h_new(i,j,k)
584  enddo ; enddo
585  enddo
586  if (cs%show_call_tree) call calltree_leave("ALE_offline_tracer_final()")
587 end subroutine ale_offline_tracer_final
588 
589 !> Check grid for negative thicknesses
590 subroutine check_grid( G, GV, h, threshold )
591  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
592  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
593  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Current 3D grid obtained after the
594  !! last time step [H ~> m or kg m-2]
595  real, intent(in) :: threshold !< Value below which to flag issues,
596  !! [H ~> m or kg m-2]
597  ! Local variables
598  integer :: i, j
599 
600  do j = g%jsc,g%jec ; do i = g%isc,g%iec
601  if (g%mask2dT(i,j)>0.) then
602  if (minval(h(i,j,:)) < threshold) then
603  write(0,*) 'check_grid: i,j=',i,j,'h(i,j,:)=',h(i,j,:)
604  if (threshold <= 0.) then
605  call mom_error(fatal,"MOM_ALE, check_grid: negative thickness encountered.")
606  else
607  call mom_error(fatal,"MOM_ALE, check_grid: too tiny thickness encountered.")
608  endif
609  endif
610  endif
611  enddo ; enddo
612 
613 
614 end subroutine check_grid
615 
616 !> Generates new grid
617 subroutine ale_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h )
618  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
619  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
620  type(regridding_cs), intent(in) :: regridcs !< Regridding parameters and options
621  type(remapping_cs), intent(in) :: remapcs !< Remapping parameters and options
622  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variable structure
623  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the
624  !! last time step [H ~> m or kg-2]
625  logical, optional, intent(in) :: debug !< If true, show the call tree
626  real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
627  ! Local variables
628  integer :: nk, i, j, k
629  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
630  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! The new grid thicknesses
631  logical :: show_call_tree, use_ice_shelf
632 
633  show_call_tree = .false.
634  if (present(debug)) show_call_tree = debug
635  if (show_call_tree) call calltree_enter("ALE_build_grid(), MOM_ALE.F90")
636  use_ice_shelf = .false.
637  if (present(frac_shelf_h)) then
638  if (associated(frac_shelf_h)) use_ice_shelf = .true.
639  endif
640 
641  ! Build new grid. The new grid is stored in h_new. The old grid is h.
642  ! Both are needed for the subsequent remapping of variables.
643  if (use_ice_shelf) then
644  call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid, frac_shelf_h )
645  else
646  call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid )
647  endif
648 
649  ! Override old grid with new one. The new grid 'h_new' is built in
650  ! one of the 'build_...' routines above.
651 !$OMP parallel do default(none) shared(G,h,h_new)
652  do j = g%jsc,g%jec ; do i = g%isc,g%iec
653  if (g%mask2dT(i,j)>0.) h(i,j,:) = h_new(i,j,:)
654  enddo ; enddo
655 
656  if (show_call_tree) call calltree_leave("ALE_build_grid()")
657 end subroutine ale_build_grid
658 
659 !> For a state-based coordinate, accelerate the process of regridding by
660 !! repeatedly applying the grid calculation algorithm
661 subroutine ale_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzRegrid, initial)
662  type(ale_cs), pointer :: cs !< ALE control structure
663  type(ocean_grid_type), intent(inout) :: g !< Ocean grid
664  type(verticalgrid_type), intent(in) :: gv !< Vertical grid
665  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
666  intent(inout) :: h !< Original thicknesses [H ~> m or kg-2]
667  type(thermo_var_ptrs), intent(inout) :: tv !< Thermo vars (T/S/EOS)
668  integer, intent(in) :: n !< Number of times to regrid
669  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
670  intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
671  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
672  intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
673  type(ocean_obc_type), pointer :: obc !< Open boundary structure
674  type(tracer_registry_type), &
675  optional, pointer :: reg !< Tracer registry to remap onto new grid
676  real, optional, intent(in) :: dt !< Model timestep to provide a timescale for regridding [T ~> s]
677  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
678  optional, intent(inout) :: dzregrid !< Final change in interface positions
679  logical, optional, intent(in) :: initial !< Whether we're being called from an initialization
680  !! routine (and expect diagnostics to work)
681 
682  ! Local variables
683  integer :: i, j, k, nz
684  type(thermo_var_ptrs) :: tv_local ! local/intermediate temp/salt
685  type(group_pass_type) :: pass_t_s_h ! group pass if the coordinate has a stencil
686  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc, h_orig ! A working copy of layer thicknesses
687  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: t, s ! local temporary state
688  ! we have to keep track of the total dzInterface if for some reason
689  ! we're using the old remapping algorithm for u/v
690  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzinterface, dzinttotal
691 
692  nz = gv%ke
693 
694  ! initial total interface displacement due to successive regridding
695  dzinttotal(:,:,:) = 0.
696 
697  call create_group_pass(pass_t_s_h, t, g%domain)
698  call create_group_pass(pass_t_s_h, s, g%domain)
699  call create_group_pass(pass_t_s_h, h_loc, g%domain)
700 
701  ! copy original temp/salt and set our local tv_pointers to them
702  tv_local = tv
703  t(:,:,:) = tv%T(:,:,:)
704  s(:,:,:) = tv%S(:,:,:)
705  tv_local%T => t
706  tv_local%S => s
707 
708  ! get local copy of thickness and save original state for remapping
709  h_loc(:,:,:) = h(:,:,:)
710  h_orig(:,:,:) = h(:,:,:)
711 
712  ! Apply timescale to regridding (for e.g. filtered_grid_motion)
713  if (present(dt)) &
714  call ale_update_regrid_weights(dt, cs)
715 
716  do k = 1, n
717  call do_group_pass(pass_t_s_h, g%domain)
718 
719  ! generate new grid
720  call regridding_main(cs%remapCS, cs%regridCS, g, gv, h_loc, tv_local, h, dzinterface)
721  dzinttotal(:,:,:) = dzinttotal(:,:,:) + dzinterface(:,:,:)
722 
723  ! remap from original grid onto new grid
724  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
725  call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:))
726  call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:))
727  enddo ; enddo
728 
729  ! starting grid for next iteration
730  h_loc(:,:,:) = h(:,:,:)
731  enddo
732 
733  ! remap all state variables (including those that weren't needed for regridding)
734  call remap_all_state_vars(cs%remapCS, cs, g, gv, h_orig, h, reg, obc, dzinttotal, u, v)
735 
736  ! save total dzregrid for diags if needed?
737  if (present(dzregrid)) dzregrid(:,:,:) = dzinttotal(:,:,:)
738 end subroutine ale_regrid_accelerated
739 
740 !> This routine takes care of remapping all variable between the old and the
741 !! new grids. When velocity components need to be remapped, thicknesses at
742 !! velocity points are taken to be arithmetic averages of tracer thicknesses.
743 !! This routine is called during initialization of the model at time=0, to
744 !! remap initiali conditions to the model grid. It is also called during a
745 !! time step to update the state.
746 subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, OBC, &
747  dxInterface, u, v, debug, dt)
748  type(remapping_cs), intent(in) :: CS_remapping !< Remapping control structure
749  type(ale_cs), intent(in) :: CS_ALE !< ALE control structure
750  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
751  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
752  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid
753  !! [H ~> m or kg-2]
754  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid
755  !! [H ~> m or kg-2]
756  type(tracer_registry_type), pointer :: Reg !< Tracer registry structure
757  type(ocean_obc_type), pointer :: OBC !< Open boundary structure
758  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
759  optional, intent(in) :: dxInterface !< Change in interface position
760  !! [H ~> m or kg-2]
761  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
762  optional, intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
763  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
764  optional, intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
765  logical, optional, intent(in) :: debug !< If true, show the call tree
766  real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s]
767  ! Local variables
768  integer :: i, j, k, m
769  integer :: nz, ntr
770  real, dimension(GV%ke+1) :: dx
771  real, dimension(GV%ke) :: h1, u_column
772  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_conc
773  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_cont
774  real, dimension(SZI_(G), SZJ_(G)) :: work_2d
775  real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
776  real :: ppt2mks
777  real, dimension(GV%ke) :: h2
778  real :: h_neglect, h_neglect_edge
779  logical :: show_call_tree
780  type(tracer_type), pointer :: Tr => null()
781 
782  show_call_tree = .false.
783  if (present(debug)) show_call_tree = debug
784  if (show_call_tree) call calltree_enter("remap_all_state_vars(), MOM_ALE.F90")
785 
786  ! If remap_uv_using_old_alg is .true. and u or v is requested, then we must have dxInterface. Otherwise,
787  ! u and v can be remapped without dxInterface
788  if ( .not. present(dxinterface) .and. (cs_ale%remap_uv_using_old_alg .and. (present(u) .or. present(v))) ) then
789  call mom_error(fatal, "remap_all_state_vars: dxInterface must be present if using old algorithm "// &
790  "and u/v are to be remapped")
791  endif
792 
793  !### Try replacing both of these with GV%H_subroundoff
794  if (gv%Boussinesq) then
795  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
796  else
797  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
798  endif
799 
800  nz = gv%ke
801  ppt2mks = 0.001
802 
803  ntr = 0 ; if (associated(reg)) ntr = reg%ntr
804 
805  if (present(dt)) then
806  idt = 1.0/dt
807  work_conc(:,:,:) = 0.0
808  work_cont(:,:,:) = 0.0
809  endif
810 
811  ! Remap tracer
812  if (ntr>0) then
813  if (show_call_tree) call calltree_waypoint("remapping tracers (remap_all_state_vars)")
814  !$OMP parallel do default(shared) private(h1,h2,u_column,Tr)
815  do m=1,ntr ! For each tracer
816  tr => reg%Tr(m)
817  do j = g%jsc,g%jec ; do i = g%isc,g%iec ; if (g%mask2dT(i,j)>0.) then
818  ! Build the start and final grids
819  h1(:) = h_old(i,j,:)
820  h2(:) = h_new(i,j,:)
821  call remapping_core_h(cs_remapping, nz, h1, tr%t(i,j,:), nz, h2, &
822  u_column, h_neglect, h_neglect_edge)
823 
824  ! Intermediate steps for tendency of tracer concentration and tracer content.
825  if (present(dt)) then
826  if (tr%id_remap_conc > 0) then
827  do k=1,gv%ke
828  work_conc(i,j,k) = (u_column(k) - tr%t(i,j,k)) * idt
829  enddo
830  endif
831  if (tr%id_remap_cont > 0 .or. tr%id_remap_cont_2d > 0) then
832  do k=1,gv%ke
833  work_cont(i,j,k) = (u_column(k)*h2(k) - tr%t(i,j,k)*h1(k)) * idt
834  enddo
835  endif
836  endif
837  ! update tracer concentration
838  tr%t(i,j,:) = u_column(:)
839  endif ; enddo ; enddo
840 
841  ! tendency diagnostics.
842  if (present(dt)) then
843  if (tr%id_remap_conc > 0) then
844  call post_data(tr%id_remap_conc, work_conc, cs_ale%diag)
845  endif
846  if (tr%id_remap_cont > 0) then
847  call post_data(tr%id_remap_cont, work_cont, cs_ale%diag)
848  endif
849  if (tr%id_remap_cont_2d > 0) then
850  do j = g%jsc,g%jec ; do i = g%isc,g%iec
851  work_2d(i,j) = 0.0
852  do k = 1,gv%ke
853  work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k)
854  enddo
855  enddo ; enddo
856  call post_data(tr%id_remap_cont_2d, work_2d, cs_ale%diag)
857  endif
858  endif
859  enddo ! m=1,ntr
860 
861  endif ! endif for ntr > 0
862 
863  if (show_call_tree) call calltree_waypoint("tracers remapped (remap_all_state_vars)")
864 
865  ! Remap u velocity component
866  if ( present(u) ) then
867  !$OMP parallel do default(shared) private(h1,h2,dx,u_column)
868  do j = g%jsc,g%jec ; do i = g%iscB,g%iecB ; if (g%mask2dCu(i,j)>0.) then
869  ! Build the start and final grids
870  h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i+1,j,:) )
871  if (cs_ale%remap_uv_using_old_alg) then
872  dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i+1,j,:) )
873  do k = 1, nz
874  h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
875  enddo
876  else
877  h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i+1,j,:) )
878  endif
879  if (associated(obc)) then
880  if (obc%segnum_u(i,j) .ne. 0) then
881  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
882  h1(:) = h_old(i,j,:)
883  h2(:) = h_new(i,j,:)
884  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
885  h1(:) = h_old(i+1,j,:)
886  h2(:) = h_new(i+1,j,:)
887  endif
888  endif
889  endif
890  call remapping_core_h(cs_remapping, nz, h1, u(i,j,:), nz, h2, &
891  u_column, h_neglect, h_neglect_edge)
892  u(i,j,:) = u_column(:)
893  endif ; enddo ; enddo
894  endif
895 
896  if (show_call_tree) call calltree_waypoint("u remapped (remap_all_state_vars)")
897 
898  ! Remap v velocity component
899  if ( present(v) ) then
900  !$OMP parallel do default(shared) private(h1,h2,dx,u_column)
901  do j = g%jscB,g%jecB ; do i = g%isc,g%iec ; if (g%mask2dCv(i,j)>0.) then
902  ! Build the start and final grids
903  h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i,j+1,:) )
904  if (cs_ale%remap_uv_using_old_alg) then
905  dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i,j+1,:) )
906  do k = 1, nz
907  h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
908  enddo
909  else
910  h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i,j+1,:) )
911  endif
912  if (associated(obc)) then
913  if (obc%segnum_v(i,j) .ne. 0) then
914  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
915  h1(:) = h_old(i,j,:)
916  h2(:) = h_new(i,j,:)
917  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
918  h1(:) = h_old(i,j+1,:)
919  h2(:) = h_new(i,j+1,:)
920  endif
921  endif
922  endif
923  call remapping_core_h(cs_remapping, nz, h1, v(i,j,:), nz, h2, &
924  u_column, h_neglect, h_neglect_edge)
925  v(i,j,:) = u_column(:)
926  endif ; enddo ; enddo
927  endif
928 
929  if (cs_ale%id_vert_remap_h > 0) call post_data(cs_ale%id_vert_remap_h, h_old, cs_ale%diag)
930  if ((cs_ale%id_vert_remap_h_tendency > 0) .and. present(dt)) then
931  do k = 1, nz ; do j = g%jsc,g%jec ; do i = g%isc,g%iec
932  work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*idt
933  enddo ; enddo ; enddo
934  call post_data(cs_ale%id_vert_remap_h_tendency, work_cont, cs_ale%diag)
935  endif
936  if (show_call_tree) call calltree_waypoint("v remapped (remap_all_state_vars)")
937  if (show_call_tree) call calltree_leave("remap_all_state_vars()")
938 
939 end subroutine remap_all_state_vars
940 
941 
942 !> Remaps a single scalar between grids described by thicknesses h_src and h_dst.
943 !! h_dst must be dimensioned as a model array with GV%ke layers while h_src can
944 !! have an arbitrary number of layers specified by nk_src.
945 subroutine ale_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap )
946  type(remapping_cs), intent(in) :: cs !< Remapping control structure
947  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
948  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
949  integer, intent(in) :: nk_src !< Number of levels on source grid
950  real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: h_src !< Level thickness of source grid
951  !! [H ~> m or kg-2]
952  real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: s_src !< Scalar on source grid
953  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid
954  !! [H ~> m or kg-2]
955  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst !< Scalar on destination grid
956  logical, optional, intent(in) :: all_cells !< If false, only reconstruct for
957  !! non-vanished cells. Use all vanished
958  !! layers otherwise (default).
959  logical, optional, intent(in) :: old_remap !< If true, use the old "remapping_core_w"
960  !! method, otherwise use "remapping_core_h".
961  ! Local variables
962  integer :: i, j, k, n_points
963  real :: dx(gv%ke+1)
964  real :: h_neglect, h_neglect_edge
965  logical :: ignore_vanished_layers, use_remapping_core_w
966 
967  ignore_vanished_layers = .false.
968  if (present(all_cells)) ignore_vanished_layers = .not. all_cells
969  use_remapping_core_w = .false.
970  if (present(old_remap)) use_remapping_core_w = old_remap
971  n_points = nk_src
972 
973  !### Try replacing both of these with GV%H_subroundoff
974  if (gv%Boussinesq) then
975  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
976  else
977  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
978  endif
979 
980  !$OMP parallel do default(shared) firstprivate(n_points,dx)
981  do j = g%jsc,g%jec ; do i = g%isc,g%iec
982  if (g%mask2dT(i,j) > 0.) then
983  if (ignore_vanished_layers) then
984  n_points = 0
985  do k = 1, nk_src
986  if (h_src(i,j,k)>0.) n_points = n_points + 1
987  enddo
988  s_dst(i,j,:) = 0.
989  endif
990  if (use_remapping_core_w) then
991  call dzfromh1h2( n_points, h_src(i,j,1:n_points), gv%ke, h_dst(i,j,:), dx )
992  call remapping_core_w(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
993  gv%ke, dx, s_dst(i,j,:), h_neglect, h_neglect_edge)
994  else
995  call remapping_core_h(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
996  gv%ke, h_dst(i,j,:), s_dst(i,j,:), h_neglect, h_neglect_edge)
997  endif
998  else
999  s_dst(i,j,:) = 0.
1000  endif
1001  enddo ; enddo
1002 
1003 end subroutine ale_remap_scalar
1004 
1005 
1006 !> Use plm reconstruction for pressure gradient (determine edge values)
1007 !! By using a PLM (limited piecewise linear method) reconstruction, this
1008 !! routine determines the edge values for the salinity and temperature
1009 !! within each layer. These edge values are returned and are used to compute
1010 !! the pressure gradient (by computing the densities).
1011 subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
1012  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
1013  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1014  type(ale_cs), intent(inout) :: cs !< module control structure
1015  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1016  intent(inout) :: s_t !< Salinity at the top edge of each layer
1017  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1018  intent(inout) :: s_b !< Salinity at the bottom edge of each layer
1019  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1020  intent(inout) :: t_t !< Temperature at the top edge of each layer
1021  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1022  intent(inout) :: t_b !< Temperature at the bottom edge of each layer
1023  type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure
1024  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1025  intent(in) :: h !< layer thickness [H ~> m or kg m-2]
1026  logical, intent(in) :: bdry_extrap !< If true, use high-order boundary
1027  !! extrapolation within boundary cells
1028 
1029  ! Local variables
1030  integer :: i, j, k
1031  real :: htmp(gv%ke)
1032  real :: tmp(gv%ke)
1033  real, dimension(CS%nk,2) :: ppol_e !Edge value of polynomial
1034  real, dimension(CS%nk,2) :: ppol_coefs !Coefficients of polynomial
1035  real :: h_neglect
1036 
1037  !### Replace this with GV%H_subroundoff
1038  if (gv%Boussinesq) then
1039  h_neglect = gv%m_to_H*1.0e-30
1040  else
1041  h_neglect = gv%kg_m2_to_H*1.0e-30
1042  endif
1043 
1044  ! Determine reconstruction within each column
1045  !$OMP parallel do default(shared) private(hTmp,ppol_E,ppol_coefs,tmp)
1046  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1047  ! Build current grid
1048  htmp(:) = h(i,j,:)
1049  tmp(:) = tv%S(i,j,:)
1050  ! Reconstruct salinity profile
1051  ppol_e(:,:) = 0.0
1052  ppol_coefs(:,:) = 0.0
1053  call plm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1054  if (bdry_extrap) &
1055  call plm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1056 
1057  do k = 1,gv%ke
1058  s_t(i,j,k) = ppol_e(k,1)
1059  s_b(i,j,k) = ppol_e(k,2)
1060  enddo
1061 
1062  ! Reconstruct temperature profile
1063  ppol_e(:,:) = 0.0
1064  ppol_coefs(:,:) = 0.0
1065  tmp(:) = tv%T(i,j,:)
1066  call plm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1067  if (bdry_extrap) &
1068  call plm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1069 
1070  do k = 1,gv%ke
1071  t_t(i,j,k) = ppol_e(k,1)
1072  t_b(i,j,k) = ppol_e(k,2)
1073  enddo
1074 
1075  enddo ; enddo
1076 
1077 end subroutine pressure_gradient_plm
1078 
1079 
1080 !> Use ppm reconstruction for pressure gradient (determine edge values)
1081 !> By using a PPM (limited piecewise linear method) reconstruction, this
1082 !> routine determines the edge values for the salinity and temperature
1083 !> within each layer. These edge values are returned and are used to compute
1084 !> the pressure gradient (by computing the densities).
1085 subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
1086  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
1087  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1088  type(ale_cs), intent(inout) :: cs !< module control structure
1089  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1090  intent(inout) :: s_t !< Salinity at the top edge of each layer
1091  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1092  intent(inout) :: s_b !< Salinity at the bottom edge of each layer
1093  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1094  intent(inout) :: t_t !< Temperature at the top edge of each layer
1095  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1096  intent(inout) :: t_b !< Temperature at the bottom edge of each layer
1097  type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure
1098  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1099  intent(in) :: h !< layer thicknesses [H ~> m or kg m-2]
1100  logical, intent(in) :: bdry_extrap !< If true, use high-order boundary
1101  !! extrapolation within boundary cells
1102 
1103  ! Local variables
1104  integer :: i, j, k
1105  real :: htmp(gv%ke) ! A 1-d copy of h [H ~> m or kg m-2]
1106  real :: tmp(gv%ke) ! A 1-d copy of a column of temperature [degC] or salinity [ppt]
1107  real, dimension(CS%nk,2) :: &
1108  ppol_e ! Edge value of polynomial in [degC] or [ppt]
1109  real, dimension(CS%nk,3) :: &
1110  ppol_coefs ! Coefficients of polynomial, all in [degC] or [ppt]
1111  real :: h_neglect, h_neglect_edge ! Tiny thicknesses [H ~> m or kg m-2]
1112 
1113  !### Try replacing both of these with GV%H_subroundoff
1114  if (gv%Boussinesq) then
1115  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1116  else
1117  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1118  endif
1119 
1120  ! Determine reconstruction within each column
1121  !$OMP parallel do default(shared) private(hTmp,tmp,ppol_E,ppol_coefs)
1122  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1123 
1124  ! Build current grid
1125  htmp(:) = h(i,j,:)
1126  tmp(:) = tv%S(i,j,:)
1127 
1128  ! Reconstruct salinity profile
1129  ppol_e(:,:) = 0.0
1130  ppol_coefs(:,:) = 0.0
1131  !### Try to replace the following value of h_neglect with GV%H_subroundoff.
1132  call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=h_neglect_edge, &
1133  answers_2018=cs%answers_2018 )
1134  call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1135  if (bdry_extrap) &
1136  call ppm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1137 
1138  do k = 1,gv%ke
1139  s_t(i,j,k) = ppol_e(k,1)
1140  s_b(i,j,k) = ppol_e(k,2)
1141  enddo
1142 
1143  ! Reconstruct temperature profile
1144  ppol_e(:,:) = 0.0
1145  ppol_coefs(:,:) = 0.0
1146  tmp(:) = tv%T(i,j,:)
1147  !### Try to replace the following value of h_neglect with GV%H_subroundoff.
1148  call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=1.0e-10*gv%m_to_H, &
1149  answers_2018=cs%answers_2018 )
1150  call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1151  if (bdry_extrap) &
1152  call ppm_boundary_extrapolation(gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1153 
1154  do k = 1,gv%ke
1155  t_t(i,j,k) = ppol_e(k,1)
1156  t_b(i,j,k) = ppol_e(k,2)
1157  enddo
1158 
1159  enddo ; enddo
1160 
1161 end subroutine pressure_gradient_ppm
1162 
1163 
1164 !> Initializes regridding for the main ALE algorithm
1165 subroutine ale_initregridding(GV, US, max_depth, param_file, mdl, regridCS)
1166  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1167  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1168  real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
1169  type(param_file_type), intent(in) :: param_file !< parameter file
1170  character(len=*), intent(in) :: mdl !< Name of calling module
1171  type(regridding_cs), intent(out) :: regridcs !< Regridding parameters and work arrays
1172  ! Local variables
1173  character(len=30) :: coord_mode
1174 
1175  call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", coord_mode, &
1176  "Coordinate mode for vertical regridding. "//&
1177  "Choose among the following possibilities: "//&
1178  trim(regriddingcoordinatemodedoc), &
1179  default=default_coordinate_mode, fail_if_missing=.true.)
1180 
1181  call initialize_regridding(regridcs, gv, us, max_depth, param_file, mdl, coord_mode, '', '')
1182 
1183 end subroutine ale_initregridding
1184 
1185 !> Query the target coordinate interfaces positions
1186 function ale_getcoordinate( CS )
1187  type(ale_cs), pointer :: cs !< module control structure
1188 
1189  real, dimension(CS%nk+1) :: ale_getcoordinate
1190  ale_getcoordinate(:) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )
1191 
1192 end function ale_getcoordinate
1193 
1194 
1195 !> Query the target coordinate units
1196 function ale_getcoordinateunits( CS )
1197  type(ale_cs), pointer :: cs !< module control structure
1198 
1199  character(len=20) :: ale_getcoordinateunits
1200 
1201  ale_getcoordinateunits = getcoordinateunits( cs%regridCS )
1202 
1203 end function ale_getcoordinateunits
1204 
1205 
1206 !> Returns true if initial conditions should be regridded and remapped
1207 logical function ale_remap_init_conds( CS )
1208  type(ale_cs), pointer :: cs !< module control structure
1209 
1210  ale_remap_init_conds = .false.
1211  if (associated(cs)) ale_remap_init_conds = cs%remap_after_initialization
1212 end function ale_remap_init_conds
1213 
1214 !> Updates the weights for time filtering the new grid generated in regridding
1215 subroutine ale_update_regrid_weights( dt, CS )
1216  real, intent(in) :: dt !< Time-step used between ALE calls [T ~> s]
1217  type(ale_cs), pointer :: cs !< ALE control structure
1218  ! Local variables
1219  real :: w ! An implicit weighting estimate.
1220 
1221  if (associated(cs)) then
1222  w = 0.0
1223  if (cs%regrid_time_scale > 0.0) then
1224  w = cs%regrid_time_scale / (cs%regrid_time_scale + dt)
1225  endif
1226  call set_regrid_params(cs%regridCS, old_grid_weight=w)
1227  endif
1228 
1229 end subroutine ale_update_regrid_weights
1230 
1231 !> Update the vertical grid type with ALE information.
1232 !! This subroutine sets information in the verticalGrid_type to be
1233 !! consistent with the use of ALE mode.
1234 subroutine ale_updateverticalgridtype(CS, GV)
1235  type(ale_cs), pointer :: cs !< ALE control structure
1236  type(verticalgrid_type), pointer :: gv !< vertical grid information
1237 
1238  integer :: nk
1239 
1240  nk = gv%ke
1241  gv%sInterface(1:nk+1) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )
1242  gv%sLayer(1:nk) = 0.5*( gv%sInterface(1:nk) + gv%sInterface(2:nk+1) )
1243  gv%zAxisUnits = getcoordinateunits( cs%regridCS )
1244  gv%zAxisLongName = getcoordinateshortname( cs%regridCS )
1245  gv%direction = -1 ! Because of ferret in z* mode. Need method to set
1246  ! as function of coordinae mode.
1247 
1248 end subroutine ale_updateverticalgridtype
1249 
1250 
1251 !> Write the vertical coordinate information into a file.
1252 !! This subroutine writes out a file containing any available data related
1253 !! to the vertical grid used by the MOM ocean model when in ALE mode.
1254 subroutine ale_writecoordinatefile( CS, GV, directory )
1255  type(ale_cs), pointer :: cs !< module control structure
1256  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1257  character(len=*), intent(in) :: directory !< directory for writing grid info
1258 
1259  character(len=240) :: filepath
1260  type(vardesc) :: vars(2)
1261  type(fieldtype) :: fields(2)
1262  integer :: unit
1263  real :: ds(gv%ke), dsi(gv%ke+1)
1264 
1265  filepath = trim(directory) // trim("Vertical_coordinate")
1266  ds(:) = getcoordinateresolution( cs%regridCS, undo_scaling=.true. )
1267  dsi(1) = 0.5*ds(1)
1268  dsi(2:gv%ke) = 0.5*( ds(1:gv%ke-1) + ds(2:gv%ke) )
1269  dsi(gv%ke+1) = 0.5*ds(gv%ke)
1270 
1271  vars(1) = var_desc('ds', getcoordinateunits( cs%regridCS ), &
1272  'Layer Coordinate Thickness','1','L','1')
1273  vars(2) = var_desc('ds_interface', getcoordinateunits( cs%regridCS ), &
1274  'Layer Center Coordinate Separation','1','i','1')
1275 
1276  call create_file(unit, trim(filepath), vars, 2, fields, single_file, gv=gv)
1277  call write_field(unit, fields(1), ds)
1278  call write_field(unit, fields(2), dsi)
1279  call close_file(unit)
1280 
1281 end subroutine ale_writecoordinatefile
1282 
1283 
1284 !> Set h to coordinate values for fixed coordinate systems
1285 subroutine ale_initthicknesstocoord( CS, G, GV, h )
1286  type(ale_cs), intent(inout) :: cs !< module control structure
1287  type(ocean_grid_type), intent(in) :: g !< module grid structure
1288  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1289  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< layer thickness [H ~> m or kg m-2]
1290 
1291  ! Local variables
1292  integer :: i, j, k
1293 
1294  do j = g%jsd,g%jed ; do i = g%isd,g%ied
1295  h(i,j,:) = gv%Z_to_H * getstaticthickness( cs%regridCS, 0., g%bathyT(i,j) )
1296  enddo ; enddo
1297 
1298 end subroutine ale_initthicknesstocoord
1299 
1300 end module mom_ale
mom_ale::ale_offline_inputs
subroutine, public ale_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
Regrid/remap stored fields used for offline tracer integrations. These input fields are assumed to ha...
Definition: MOM_ALE.F90:465
ppm_functions::ppm_boundary_extrapolation
subroutine, public ppm_boundary_extrapolation(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Reconstruction by parabolas within boundary cells.
Definition: PPM_functions.F90:134
mom_ale::check_grid
subroutine check_grid(G, GV, h, threshold)
Check grid for negative thicknesses.
Definition: MOM_ALE.F90:591
plm_functions::plm_boundary_extrapolation
subroutine, public plm_boundary_extrapolation(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Reconstruction by linear polynomials within boundary cells.
Definition: PLM_functions.F90:206
ppm_functions
Provides functions used with the Piecewise-Parabolic-Method in the vertical ALE algorithm.
Definition: PPM_functions.F90:2
mom_io::var_desc
type(vardesc) function, public var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_name, cmor_units, cmor_longname, conversion, caller)
Returns a vardesc type whose elements have been filled with the provided fields. The argument name is...
Definition: MOM_io.F90:600
mom_regridding::set_regrid_params
subroutine, public set_regrid_params(CS, boundary_extrapolation, min_thickness, old_grid_weight, interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_to_interior, fix_haloclines, halocline_filt_len, halocline_strat_tol, integrate_downward_for_e, adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
Can be used to set any of the parameters for MOM_regridding.
Definition: MOM_regridding.F90:2218
mom_open_boundary::obc_direction_n
integer, parameter, public obc_direction_n
Indicates the boundary is an effective northern boundary.
Definition: MOM_open_boundary.F90:63
mom_open_boundary::obc_direction_s
integer, parameter, public obc_direction_s
Indicates the boundary is an effective southern boundary.
Definition: MOM_open_boundary.F90:64
mom_diag_vkernels
Provides kernels for single-column interpolation, re-integration (re-mapping of integrated quantities...
Definition: MOM_diag_vkernels.F90:3
mom_regridding::initialize_regridding
subroutine, public initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_mode, param_prefix, param_suffix)
Initialization and configures a regridding control structure based on customizable run-time parameter...
Definition: MOM_regridding.F90:177
mom_regridding::regridding_cs
Regridding control structure.
Definition: MOM_regridding.F90:45
mom_diag_vkernels::interpolate_column
subroutine, public interpolate_column(nsrc, h_src, u_src, ndest, h_dest, missing_value, u_dest)
Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest.
Definition: MOM_diag_vkernels.F90:17
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_ale::ale_remap_scalar
subroutine, public ale_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap)
Remaps a single scalar between grids described by thicknesses h_src and h_dst. h_dst must be dimensio...
Definition: MOM_ALE.F90:946
mom_regridding::getstaticthickness
real function, dimension(cs%nk), public getstaticthickness(CS, SSH, depth)
Return coordinate-derived thicknesses for fixed coordinate systems.
Definition: MOM_regridding.F90:2343
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
regrid_consts::state_dependent
Returns true if the coordinate is dependent on the state density, returns false otherwise.
Definition: regrid_consts.F90:44
mom_tracer_registry::tracer_type
The tracer type.
Definition: MOM_tracer_registry.F90:38
mom_regridding::regriddinginterpschemedoc
character(len= *), parameter, public regriddinginterpschemedoc
Documentation for regridding interpolation schemes.
Definition: MOM_regridding.F90:152
mom_ale::adjustgridforintegrity
subroutine, public adjustgridforintegrity(CS, G, GV, h)
Crudely adjust (initial) grid for integrity. This routine is typically called (from initialize_MOM in...
Definition: MOM_ALE.F90:295
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
p3m_functions::p3m_interpolation
subroutine, public p3m_interpolation(N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
Set up a piecewise cubic interpolation from cell averages and estimated edge slopes and values.
Definition: P3M_functions.F90:29
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_open_boundary::obc_direction_e
integer, parameter, public obc_direction_e
Indicates the boundary is an effective eastern boundary.
Definition: MOM_open_boundary.F90:65
mom_open_boundary::obc_direction_w
integer, parameter, public obc_direction_w
Indicates the boundary is an effective western boundary.
Definition: MOM_open_boundary.F90:66
p3m_functions
Cubic interpolation functions.
Definition: P3M_functions.F90:2
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_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
regrid_consts::coordinateunits
Returns a string with the coordinate units associated with the coordinate mode.
Definition: regrid_consts.F90:38
mom_verticalgrid::get_thickness_units
character(len=48) function, public get_thickness_units(GV)
Returns the model's thickness units, usually m or kg/m^2.
Definition: MOM_verticalGrid.F90:190
mom_ale::ale_writecoordinatefile
subroutine, public ale_writecoordinatefile(CS, GV, directory)
Write the vertical coordinate information into a file. This subroutine writes out a file containing a...
Definition: MOM_ALE.F90:1255
mom_regridding::regridding_main
subroutine, public regridding_main(remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
Dispatching regridding routine for orchestrating regridding & remapping.
Definition: MOM_regridding.F90:791
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_regridding::regriddingdefaultboundaryextrapolation
logical, parameter, public regriddingdefaultboundaryextrapolation
Default mode for boundary extrapolation.
Definition: MOM_regridding.F90:167
mom_remapping::remapping_cs
Container for remapping parameters.
Definition: MOM_remapping.F90:22
mom_tracer_registry
This module contains the tracer_registry_type and the subroutines that handle registration of tracers...
Definition: MOM_tracer_registry.F90:5
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::ale_regrid_accelerated
subroutine, public ale_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzRegrid, initial)
For a state-based coordinate, accelerate the process of regridding by repeatedly applying the grid ca...
Definition: MOM_ALE.F90:662
mom_ale::ale_initregridding
subroutine, public ale_initregridding(GV, US, max_depth, param_file, mdl, regridCS)
Initializes regridding for the main ALE algorithm.
Definition: MOM_ALE.F90:1166
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_regridding::getcoordinateunits
character(len=20) function, public getcoordinateunits(CS)
Query the target coordinate units.
Definition: MOM_regridding.F90:2158
regrid_edge_values
Edge value estimation for high-order resconstruction.
Definition: regrid_edge_values.F90:2
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_string_functions::uppercase
character(len=len(input_string)) function, public uppercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:42
mom_debugging::check_column_integrals
logical function, public check_column_integrals(nk_1, field_1, nk_2, field_2, missing_value)
Returns false if the column integrals of two given quantities are within roundoff of each other.
Definition: MOM_debugging.F90:823
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
p1m_functions
Linear interpolation functions.
Definition: P1M_functions.F90:2
regrid_consts
Contains constants for interpreting input parameters that control regridding.
Definition: regrid_consts.F90:2
mom_remapping
Provides column-wise vertical remapping functions.
Definition: MOM_remapping.F90:2
mom_regridding::getcoordinateinterfaces
real function, dimension(cs%nk+1), public getcoordinateinterfaces(CS, undo_scaling)
Query the target coordinate interface positions.
Definition: MOM_regridding.F90:2113
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_regridding::getcoordinateresolution
real function, dimension(cs%nk), public getcoordinateresolution(CS, undo_scaling)
Query the fixed resolution data.
Definition: MOM_regridding.F90:2095
mom_ale::ale_initthicknesstocoord
subroutine, public ale_initthicknesstocoord(CS, G, GV, h)
Set h to coordinate values for fixed coordinate systems.
Definition: MOM_ALE.F90:1286
mom_tracer_registry::mom_tracer_chkinv
subroutine, public mom_tracer_chkinv(mesg, G, h, Tr, ntr)
Calculates and prints the global inventory of all tracers in the registry.
Definition: MOM_tracer_registry.F90:821
mom_regridding::regriddingdefaultminthickness
real, parameter, public regriddingdefaultminthickness
Default minimum thickness for some coordinate generation modes.
Definition: MOM_regridding.F90:169
mom_remapping::remappingdefaultscheme
character(len=3), public remappingdefaultscheme
Default remapping method.
Definition: MOM_remapping.F90:69
mom_regridding::inflate_vanished_layers_old
subroutine, public inflate_vanished_layers_old(CS, G, GV, h)
Definition: MOM_regridding.F90:1803
mom_ale::ale_build_grid
subroutine, public ale_build_grid(G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h)
Generates new grid.
Definition: MOM_ALE.F90:618
mom_ale::ale_updateverticalgridtype
subroutine, public ale_updateverticalgridtype(CS, GV)
Update the vertical grid type with ALE information. This subroutine sets information in the verticalG...
Definition: MOM_ALE.F90:1235
mom_diag_mediator::diag_update_remap_grids
subroutine, public diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S)
Build/update vertical grids for diagnostic remapping.
Definition: MOM_diag_mediator.F90:3187
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_ale::ale_offline_tracer_final
subroutine, public ale_offline_tracer_final(G, GV, h, tv, h_target, Reg, CS, OBC)
Remaps all tracers from h onto h_target. This is intended to be called when tracers are done offline....
Definition: MOM_ALE.F90:546
mom_remapping::remappingschemesdoc
character(len=256), public remappingschemesdoc
Documentation for external callers.
Definition: MOM_remapping.F90:62
p1m_functions::p1m_boundary_extrapolation
subroutine, public p1m_boundary_extrapolation(N, h, u, ppoly_E, ppoly_coef)
Interpolation by linear polynomials within boundary cells.
Definition: P1M_functions.F90:70
mom_domains::do_group_pass
subroutine, public do_group_pass(group, MOM_dom, clock)
do_group_pass carries out a group halo update.
Definition: MOM_domains.F90:1090
mom_ale::ale_getcoordinateunits
character(len=20) function, public ale_getcoordinateunits(CS)
Query the target coordinate units.
Definition: MOM_ALE.F90:1197
regrid_edge_values::edge_values_implicit_h4
subroutine, public edge_values_implicit_h4(N, h, u, edge_val, h_neglect, answers_2018)
Compute ih4 edge values (implicit fourth order accurate) in the same units as h.
Definition: regrid_edge_values.F90:492
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_ale
This module contains the main regridding routines.
Definition: MOM_ALE.F90:9
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_ale::ale_init
subroutine, public ale_init(param_file, GV, US, max_depth, CS)
This routine is typically called (from initialize_MOM in file MOM.F90) before the main time integrati...
Definition: MOM_ALE.F90:141
mom_ale::pressure_gradient_plm
subroutine, public pressure_gradient_plm(CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap)
Use plm reconstruction for pressure gradient (determine edge values) By using a PLM (limited piecewis...
Definition: MOM_ALE.F90:1012
mom_open_boundary
Controls where open boundary conditions are applied.
Definition: MOM_open_boundary.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_ale::ale_cs
ALE control structure.
Definition: MOM_ALE.F90:65
mom_error_handler::calltree_leave
subroutine, public calltree_leave(mesg)
Writes a message about leaving a subroutine if call tree reporting is active.
Definition: MOM_error_handler.F90:151
mom_io::create_file
subroutine, public create_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV, checksums)
Routine creates a new NetCDF file. It also sets up structures that describe this file and variables t...
Definition: MOM_io.F90:93
mom_regridding::end_regridding
subroutine, public end_regridding(CS)
Deallocation of regridding memory.
Definition: MOM_regridding.F90:772
plm_functions::plm_reconstruction
subroutine, public plm_reconstruction(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Reconstruction by linear polynomials within each cell.
Definition: PLM_functions.F90:19
mom_regridding::regriddingcoordinatemodedoc
character(len= *), parameter, public regriddingcoordinatemodedoc
Documentation for coordinate options.
Definition: MOM_regridding.F90:141
mom_remapping::end_remapping
subroutine, public end_remapping(CS)
Destrcutor for remapping control structure.
Definition: MOM_remapping.F90:1603
mom_ale::ale_getcoordinate
real function, dimension(cs%nk+1), public ale_getcoordinate(CS)
Query the target coordinate interfaces positions.
Definition: MOM_ALE.F90:1187
mom_string_functions::extractword
character(len=120) function, public extractword(string, n)
Returns the string corresponding to the nth word in the argument or "" if the string is not long enou...
Definition: MOM_string_functions.F90:196
mom_tracer_registry::tracer_registry_type
Type to carry basic tracer information.
Definition: MOM_tracer_registry.F90:138
mom_interface_heights::find_eta
Calculates the heights of sruface or all interfaces from layer thicknesses.
Definition: MOM_interface_heights.F90:21
mom_open_boundary::ocean_obc_type
Open-boundary data.
Definition: MOM_open_boundary.F90:195
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
regrid_consts::coordinatemode
integer function coordinatemode(string)
Parse a string parameter specifying the coordinate mode and return the appropriate enumerated integer...
Definition: regrid_consts.F90:54
mom_ale::remap_all_state_vars
subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, OBC, dxInterface, u, v, debug, dt)
This routine takes care of remapping all variable between the old and the new grids....
Definition: MOM_ALE.F90:748
mom_regridding::set_target_densities
subroutine, public set_target_densities(CS, rho_int)
Set target densities based on vector of interface values.
Definition: MOM_regridding.F90:2015
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_io::vardesc
Type for describing a variable, typically a tracer.
Definition: MOM_io.F90:53
mom_regridding::set_target_densities_from_gv
subroutine, public set_target_densities_from_gv(GV, US, CS)
Set target densities based on the old Rlay variable.
Definition: MOM_regridding.F90:1997
mom_error_handler::calltree_showquery
logical function, public calltree_showquery()
Returns True, if the verbosity>=6 indicating to show the call tree.
Definition: MOM_error_handler.F90:123
mom_ale::ale_main_offline
subroutine, public ale_main_offline(G, GV, h, tv, Reg, CS, OBC, dt)
Takes care of (1) building a new grid and (2) remapping all variables between the old grid and the ne...
Definition: MOM_ALE.F90:411
mom_regridding::regriddingdefaultinterpscheme
character(len= *), parameter, public regriddingdefaultinterpscheme
Default interpolation scheme.
Definition: MOM_regridding.F90:165
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_error_handler::calltree_waypoint
subroutine, public calltree_waypoint(mesg, n)
Writes a message about reaching a milestone if call tree reporting is active.
Definition: MOM_error_handler.F90:161
p1m_functions::p1m_interpolation
subroutine, public p1m_interpolation(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Linearly interpolate between edge values.
Definition: P1M_functions.F90:28
mom_ale::ale_end
subroutine, public ale_end(CS)
End of regridding (memory deallocation). This routine is typically called (from MOM_end in file MOM....
Definition: MOM_ALE.F90:309
mom_ale::ale_update_regrid_weights
subroutine, public ale_update_regrid_weights(dt, CS)
Updates the weights for time filtering the new grid generated in regridding.
Definition: MOM_ALE.F90:1216
mom_regridding::uniformresolution
real function, dimension(nk), public uniformresolution(nk, coordMode, maxDepth, rhoLight, rhoHeavy)
Return a uniform resolution vector in the units of the coordinate.
Definition: MOM_regridding.F90:1908
mom_regridding
Generates vertical grids as part of the ALE algorithm.
Definition: MOM_regridding.F90:2
mom_regridding::getcoordinateshortname
character(len=20) function, public getcoordinateshortname(CS)
Query the short name of the coordinate.
Definition: MOM_regridding.F90:2182
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_remapping::remapping_core_w
subroutine, public remapping_core_w(CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_edge)
Remaps column of values u0 on grid h0 to implied grid h1 where the interfaces of h1 differ from those...
Definition: MOM_remapping.F90:266
mom_ale::pressure_gradient_ppm
subroutine, public pressure_gradient_ppm(CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap)
Use ppm reconstruction for pressure gradient (determine edge values) By using a PPM (limited piecewis...
Definition: MOM_ALE.F90:1086
mom_ale::ale_remap_init_conds
logical function, public ale_remap_init_conds(CS)
Returns true if initial conditions should be regridded and remapped.
Definition: MOM_ALE.F90:1208
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
p3m_functions::p3m_boundary_extrapolation
subroutine, public p3m_boundary_extrapolation(N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect, h_neglect_edge)
Calculate the edge values and slopes at boundary cells as part of building a piecewise cubic sub-grid...
Definition: P3M_functions.F90:190
mom_domains::create_group_pass
Set up a group of halo updates.
Definition: MOM_domains.F90:79
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_ale::ale_register_diags
subroutine, public ale_register_diags(Time, G, GV, US, diag, CS)
Initialize diagnostics for the ALE module.
Definition: MOM_ALE.F90:252
mom_interface_heights
Functions for calculating interface heights, including free surface height.
Definition: MOM_interface_heights.F90:2
mom_error_handler::calltree_enter
subroutine, public calltree_enter(mesg, n)
Writes a message about entering a subroutine if call tree reporting is active.
Definition: MOM_error_handler.F90:130
plm_functions
Piecewise linear reconstruction functions.
Definition: PLM_functions.F90:2
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_remapping::dzfromh1h2
subroutine, public dzfromh1h2(n1, h1, n2, h2, dx)
Calculates the change in interface positions based on h1 and h2.
Definition: MOM_remapping.F90:1522
mom_ale::ale_main
subroutine, public ale_main(G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
Takes care of (1) building a new grid and (2) remapping all variables between the old grid and the ne...
Definition: MOM_ALE.F90:324
ppm_functions::ppm_reconstruction
subroutine, public ppm_reconstruction(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Builds quadratic polynomials coefficients from cell mean and edge values.
Definition: PPM_functions.F90:29
mom_string_functions::extract_integer
integer function, public extract_integer(string, separators, n, missing_value)
Returns the integer corresponding to the nth word in the argument.
Definition: MOM_string_functions.F90:244
mom_diag_vkernels::reintegrate_column
subroutine, public reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, uh_dest)
Conservatively calculate integrated data, uh_dest, on grid h_dest, from layer-integrated data,...
Definition: MOM_diag_vkernels.F90:92
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60