MOM6
MOM_regridding.F90
Go to the documentation of this file.
1 !> Generates vertical grids as part of the ALE algorithm
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_error_handler, only : mom_error, fatal, warning
8 use mom_io, only : file_exists, field_exists, field_size, mom_read_data
9 use mom_io, only : slasher
11 use mom_variables, only : ocean_grid_type, thermo_var_ptrs
15 
16 use mom_remapping, only : remapping_cs
24 
32 
33 use netcdf ! Used by check_grid_def()
34 
35 implicit none ; private
36 
37 #include <MOM_memory.h>
38 
39 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
40 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
41 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
42 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
43 
44 !> Regridding control structure
45 type, public :: regridding_cs ; private
46 
47  !> This array is set by function setCoordinateResolution()
48  !! It contains the "resolution" or delta coordinate of the target
49  !! coordinate. It has the units of the target coordinate, e.g.
50  !! [Z ~> m] for z*, non-dimensional for sigma, etc.
51  real, dimension(:), allocatable :: coordinateresolution
52 
53  !> This is a scaling factor that restores coordinateResolution to values in
54  !! the natural units for output.
55  real :: coord_scale = 1.0
56 
57  !> This array is set by function set_target_densities()
58  !! This array is the nominal coordinate of interfaces and is the
59  !! running sum of coordinateResolution, in [R ~> kg m-3]. i.e.
60  !! target_density(k+1) = coordinateResolution(k) + coordinateResolution(k)
61  !! It is only used in "rho", "SLight" or "Hycom" mode.
62  real, dimension(:), allocatable :: target_density
63 
64  !> A flag to indicate that the target_density arrays has been filled with data.
65  logical :: target_density_set = .false.
66 
67  !> This array is set by function set_regrid_max_depths()
68  !! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
69  real, dimension(:), allocatable :: max_interface_depths
70 
71  !> This array is set by function set_regrid_max_thickness()
72  !! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
73  real, dimension(:), allocatable :: max_layer_thickness
74 
75  integer :: nk !< Number of layers/levels in generated grid
76 
77  !> Indicates which grid to use in the vertical (z*, sigma, target interface
78  !! densities)
79  integer :: regridding_scheme
80 
81  !> Interpolation control structure
82  type(interp_cs_type) :: interp_cs
83 
84  !> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2].
85  real :: min_thickness
86 
87  !> Reference pressure for potential density calculations (Pa)
88  real :: ref_pressure = 2.e7
89 
90  !> Weight given to old coordinate when blending between new and old grids [nondim]
91  !! Used only below depth_of_time_filter_shallow, with a cubic variation
92  !! from zero to full effect between depth_of_time_filter_shallow and
93  !! depth_of_time_filter_deep.
94  real :: old_grid_weight = 0.
95 
96  !> Depth above which no time-filtering of grid is applied [H ~> m or kg m-2]
97  real :: depth_of_time_filter_shallow = 0.
98 
99  !> Depth below which time-filtering of grid is applied at full effect [H ~> m or kg m-2]
100  real :: depth_of_time_filter_deep = 0.
101 
102  !> Fraction (between 0 and 1) of compressibility to add to potential density
103  !! profiles when interpolating for target grid positions. [nondim]
104  real :: compressibility_fraction = 0.
105 
106  !> If true, each interface is given a maximum depth based on a rescaling of
107  !! the indexing of coordinateResolution.
108  logical :: set_maximum_depths = .false.
109 
110  !> A scaling factor (> 1) of the rate at which the coordinateResolution list
111  !! is traversed to set the minimum depth of interfaces.
112  real :: max_depth_index_scale = 2.0
113 
114  !> If true, integrate for interface positions from the top downward.
115  !! If false, integrate from the bottom upward, as does the rest of the model.
116  logical :: integrate_downward_for_e = .true.
117 
118  type(zlike_cs), pointer :: zlike_cs => null() !< Control structure for z-like coordinate generator
119  type(sigma_cs), pointer :: sigma_cs => null() !< Control structure for sigma coordinate generator
120  type(rho_cs), pointer :: rho_cs => null() !< Control structure for rho coordinate generator
121  type(hycom_cs), pointer :: hycom_cs => null() !< Control structure for hybrid coordinate generator
122  type(slight_cs), pointer :: slight_cs => null() !< Control structure for Slight-coordinate generator
123  type(adapt_cs), pointer :: adapt_cs => null() !< Control structure for adaptive coordinate generator
124 
125 end type
126 
127 ! The following routines are visible to the outside world
132 public build_rho_column
137 public default_coordinate_mode
139 
140 !> Documentation for coordinate options
141 character(len=*), parameter, public :: regriddingcoordinatemodedoc = &
142  " LAYER - Isopycnal or stacked shallow water layers\n"//&
143  " ZSTAR, Z* - stretched geopotential z*\n"//&
144  " SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"//&
145  " SIGMA - terrain following coordinates\n"//&
146  " RHO - continuous isopycnal\n"//&
147  " HYCOM1 - HyCOM-like hybrid coordinate\n"//&
148  " SLIGHT - stretched coordinates above continuous isopycnal\n"//&
149  " ADAPTIVE - optimize for smooth neutral density surfaces"
150 
151 !> Documentation for regridding interpolation schemes
152 character(len=*), parameter, public :: regriddinginterpschemedoc = &
153  " P1M_H2 (2nd-order accurate)\n"//&
154  " P1M_H4 (2nd-order accurate)\n"//&
155  " P1M_IH4 (2nd-order accurate)\n"//&
156  " PLM (2nd-order accurate)\n"//&
157  " PPM_H4 (3rd-order accurate)\n"//&
158  " PPM_IH4 (3rd-order accurate)\n"//&
159  " P3M_IH4IH3 (4th-order accurate)\n"//&
160  " P3M_IH6IH5 (4th-order accurate)\n"//&
161  " PQM_IH4IH3 (4th-order accurate)\n"//&
162  " PQM_IH6IH5 (5th-order accurate)"
163 
164 !> Default interpolation scheme
165 character(len=*), parameter, public :: regriddingdefaultinterpscheme = "P1M_H2"
166 !> Default mode for boundary extrapolation
167 logical, parameter, public :: regriddingdefaultboundaryextrapolation = .false.
168 !> Default minimum thickness for some coordinate generation modes
169 real, parameter, public :: regriddingdefaultminthickness = 1.e-3
170 
171 #undef __DO_SAFETY_CHECKS__
172 
173 contains
174 
175 !> Initialization and configures a regridding control structure based on customizable run-time parameters
176 subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_mode, param_prefix, param_suffix)
177  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
178  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
179  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
180  real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
181  type(param_file_type), intent(in) :: param_file !< Parameter file
182  character(len=*), intent(in) :: mdl !< Name of calling module.
183  character(len=*), intent(in) :: coord_mode !< Coordinate mode
184  character(len=*), intent(in) :: param_prefix !< String to prefix to parameter names.
185  !! If empty, causes main model parameters to be used.
186  character(len=*), intent(in) :: param_suffix !< String to append to parameter names.
187 
188  ! Local variables
189  integer :: ke ! Number of levels
190  character(len=80) :: string, string2, varname ! Temporary strings
191  character(len=40) :: coord_units, param_name, coord_res_param ! Temporary strings
192  character(len=200) :: inputdir, filename
193  character(len=320) :: message ! Temporary strings
194  character(len=12) :: expected_units ! Temporary strings
195  logical :: tmplogical, fix_haloclines, set_max, do_sum, main_parameters
196  logical :: coord_is_state_dependent, ierr
197  real :: filt_len, strat_tol, index_scale, tmpreal
198  real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
199  real :: dz_fixed_sfc, rho_avg_depth, nlay_sfc_int
200  real :: adapttimeratio, adaptzoom, adaptzoomcoeff, adaptbuoycoeff, adaptalpha
201  integer :: nz_fixed_sfc, k, nzf(4)
202  real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m]
203  ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
204  real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2]
205  real, dimension(:), allocatable :: z_max ! Maximum interface depths [H ~> m or kg m-2] or other
206  ! units depending on the coordinate
207  real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths
208  ! [H ~> m or kg m-2] or other units
209  real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3]
210  ! Thicknesses [m] that give level centers corresponding to table 2 of WOA09
211  real, dimension(40) :: woa09_dz = (/ 5., 10., 10., 15., 22.5, 25., 25., 25., &
212  37.5, 50., 50., 75., 100., 100., 100., 100., &
213  100., 100., 100., 100., 100., 100., 100., 175., &
214  250., 375., 500., 500., 500., 500., 500., 500., &
215  500., 500., 500., 500., 500., 500., 500., 500. /)
216 
217  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
218  inputdir = slasher(inputdir)
219 
220  main_parameters=.false.
221  if (len_trim(param_prefix)==0) main_parameters=.true.
222  if (main_parameters .and. len_trim(param_suffix)>0) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
223  'Suffix provided without prefix for parameter names!')
224 
225  cs%nk = 0
226  cs%regridding_scheme = coordinatemode(coord_mode)
227  coord_is_state_dependent = state_dependent(coord_mode)
228  maximum_depth = us%Z_to_m*max_depth
229 
230  if (main_parameters) then
231  ! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS)
232  call get_param(param_file, mdl, "REGRIDDING_COORDINATE_UNITS", coord_units, &
233  "Units of the regridding coordinate.", default=coordinateunits(coord_mode))
234  else
235  coord_units=coordinateunits(coord_mode)
236  endif
237 
238  if (coord_is_state_dependent) then
239  if (main_parameters) then
240  param_name = "INTERPOLATION_SCHEME"
242  else
243  param_name = trim(param_prefix)//"_INTERP_SCHEME_"//trim(param_suffix)
244  string2 = 'PPM_H4' ! Default for diagnostics
245  endif
246  call get_param(param_file, mdl, "INTERPOLATION_SCHEME", string, &
247  "This sets the interpolation scheme to use to "//&
248  "determine the new grid. These parameters are "//&
249  "only relevant when REGRIDDING_COORDINATE_MODE is "//&
250  "set to a function of state. Otherwise, it is not "//&
251  "used. It can be one of the following schemes: "//&
252  trim(regriddinginterpschemedoc), default=trim(string2))
253  call set_regrid_params(cs, interp_scheme=string)
254  endif
255 
256  if (main_parameters .and. coord_is_state_dependent) then
257  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", tmplogical, &
258  "When defined, a proper high-order reconstruction "//&
259  "scheme is used within boundary cells rather "//&
260  "than PCM. E.g., if PPM is used for remapping, a "//&
261  "PPM reconstruction will also be used within "//&
262  "boundary cells.", default=regriddingdefaultboundaryextrapolation)
263  call set_regrid_params(cs, boundary_extrapolation=tmplogical)
264  else
265  call set_regrid_params(cs, boundary_extrapolation=.false.)
266  endif
267 
268  ! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG)
269  if (main_parameters) then
270  param_name = "ALE_COORDINATE_CONFIG"
271  coord_res_param = "ALE_RESOLUTION"
272  string2 = 'UNIFORM'
273  else
274  param_name = trim(param_prefix)//"_DEF_"//trim(param_suffix)
275  coord_res_param = trim(param_prefix)//"_RES_"//trim(param_suffix)
276  string2 = 'UNIFORM'
277  if (maximum_depth>3000.) string2='WOA09' ! For convenience
278  endif
279  call get_param(param_file, mdl, param_name, string, &
280  "Determines how to specify the coordinate "//&
281  "resolution. Valid options are:\n"//&
282  " PARAM - use the vector-parameter "//trim(coord_res_param)//"\n"//&
283  " UNIFORM[:N] - uniformly distributed\n"//&
284  " FILE:string - read from a file. The string specifies\n"//&
285  " the filename and variable name, separated\n"//&
286  " by a comma or space, e.g. FILE:lev.nc,dz\n"//&
287  " or FILE:lev.nc,interfaces=zw\n"//&
288  " WOA09[:N] - the WOA09 vertical grid (approximately)\n"//&
289  " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
290  " HYBRID:string - read from a file. The string specifies\n"//&
291  " the filename and two variable names, separated\n"//&
292  " by a comma or space, for sigma-2 and dz. e.g.\n"//&
293  " HYBRID:vgrid.nc,sigma2,dz",&
294  default=trim(string2))
295  message = "The distribution of vertical resolution for the target\n"//&
296  "grid used for Eulerian-like coordinates. For example,\n"//&
297  "in z-coordinate mode, the parameter is a list of level\n"//&
298  "thicknesses (in m). In sigma-coordinate mode, the list\n"//&
299  "is of non-dimensional fractions of the water column."
300  if (index(trim(string),'UNIFORM')==1) then
301  if (len_trim(string)==7) then
302  ke = gv%ke ! Use model nk by default
303  tmpreal = maximum_depth
304  elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then
305  ! Format is "UNIFORM:N" or "UNIFORM:N,dz"
306  ke = extract_integer(string(9:len_trim(string)),'',1)
307  tmpreal = extract_real(string(9:len_trim(string)),',',2,missing_value=maximum_depth)
308  else
309  call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
310  'Unable to interpret "'//trim(string)//'".')
311  endif
312  allocate(dz(ke))
313  dz(:) = uniformresolution(ke, coord_mode, tmpreal, &
314  us%R_to_kg_m3*(gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(min(2,ke)))), &
315  us%R_to_kg_m3*(gv%Rlay(ke) + 0.5*(gv%Rlay(ke)-gv%Rlay(max(ke-1,1)))) )
316  if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
317  trim(message), units=trim(coord_units))
318  elseif (trim(string)=='PARAM') then
319  ! Read coordinate resolution (main model = ALE_RESOLUTION)
320  ke = gv%ke ! Use model nk by default
321  allocate(dz(ke))
322  call get_param(param_file, mdl, coord_res_param, dz, &
323  trim(message), units=trim(coord_units), fail_if_missing=.true.)
324  elseif (index(trim(string),'FILE:')==1) then
325  ! FILE:filename,var_name is assumed to be reading level thickness variables
326  ! FILE:filename,interfaces=var_name reads positions
327  if (string(6:6)=='.' .or. string(6:6)=='/') then
328  ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
329  filename = trim( extractword(trim(string(6:80)), 1) )
330  else
331  ! Otherwise assume we should look for the file in INPUTDIR
332  filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
333  endif
334  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
335  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
336 
337  varname = trim( extractword(trim(string(6:)), 2) )
338  if (len_trim(varname)==0) then
339  if (field_exists(filename,'dz')) then; varname = 'dz'
340  elseif (field_exists(filename,'dsigma')) then; varname = 'dsigma'
341  elseif (field_exists(filename,'ztest')) then; varname = 'ztest'
342  else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
343  "Coordinate variable not specified and none could be guessed.")
344  endif
345  endif
346  ! This check fails when the variable is a dimension variable! -AJA
347  !if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
348  ! "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
349  if (cs%regridding_scheme == regridding_sigma) then
350  expected_units = 'nondim'
351  elseif (cs%regridding_scheme == regridding_rho) then
352  expected_units = 'kg m-3'
353  else
354  expected_units = 'meters'
355  endif
356  if (index(trim(varname),'interfaces=')==1) then
357  varname=trim(varname(12:))
358  call check_grid_def(filename, varname, expected_units, message, ierr)
359  if (ierr) call mom_error(fatal,trim(mdl)//", initialize_regridding: "//&
360  "Unsupported format in grid definition '"//trim(filename)//"'. Error message "//trim(message))
361  call field_size(trim(filename), trim(varname), nzf)
362  ke = nzf(1)-1
363  if (cs%regridding_scheme == regridding_rho) then
364  allocate(rho_target(ke+1))
365  call mom_read_data(trim(filename), trim(varname), rho_target)
366  else
367  allocate(dz(ke))
368  allocate(z_max(ke+1))
369  call mom_read_data(trim(filename), trim(varname), z_max)
370  dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
371  deallocate(z_max)
372  endif
373  else
374  ! Assume reading resolution
375  call field_size(trim(filename), trim(varname), nzf)
376  ke = nzf(1)
377  allocate(dz(ke))
378  call mom_read_data(trim(filename), trim(varname), dz)
379  endif
380  if (main_parameters .and. ke/=gv%ke) then
381  call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
382  'Mismatch in number of model levels and "'//trim(string)//'".')
383  endif
384  if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
385  trim(message), units=coordinateunits(coord_mode))
386  elseif (index(trim(string),'FNC1:')==1) then
387  ke = gv%ke; allocate(dz(ke))
388  call dz_function1( trim(string(6:)), dz )
389  if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
390  trim(message), units=coordinateunits(coord_mode))
391  elseif (index(trim(string),'RFNC1:')==1) then
392  ! Function used for set target interface densities
393  ke = rho_function1( trim(string(7:)), rho_target )
394  elseif (index(trim(string),'HYBRID:')==1) then
395  ke = gv%ke; allocate(dz(ke))
396  ! The following assumes the FILE: syntax of above but without "FILE:" in the string
397  allocate(rho_target(ke+1))
398  filename = trim( extractword(trim(string(8:)), 1) )
399  if (filename(1:1)/='.' .and. filename(1:1)/='/') filename = trim(inputdir) // trim( filename )
400  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
401  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
402  varname = trim( extractword(trim(string(8:)), 2) )
403  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
404  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
405  call mom_read_data(trim(filename), trim(varname), rho_target)
406  varname = trim( extractword(trim(string(8:)), 3) )
407  if (varname(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz
408  call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
409  else ! Read dz from file
410  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
411  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
412  call mom_read_data(trim(filename), trim(varname), dz)
413  endif
414  if (main_parameters) then
415  call log_param(param_file, mdl, "!"//coord_res_param, dz, &
416  trim(message), units=coordinateunits(coord_mode))
417  call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, &
418  'HYBRID target densities for interfaces', units=coordinateunits(coord_mode))
419  endif
420  elseif (index(trim(string),'WOA09')==1) then
421  if (len_trim(string)==5) then
422  tmpreal = 0. ; ke = 0
423  do while (tmpreal<maximum_depth)
424  ke = ke + 1
425  tmpreal = tmpreal + woa09_dz(ke)
426  enddo
427  elseif (index(trim(string),'WOA09:')==1) then
428  if (len_trim(string)==6) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
429  'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
430  ke = extract_integer(string(7:len_trim(string)),'',1)
431  endif
432  if (ke>40 .or. ke<1) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
433  'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
434  allocate(dz(ke))
435  dz(1:ke) = woa09_dz(1:ke)
436  else
437  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
438  "Unrecognized coordinate configuration"//trim(string))
439  endif
440 
441  if (main_parameters) then
442  ! This is a work around to apparently needed to work with the from_Z initialization... ???
443  if (coordinatemode(coord_mode) == regridding_zstar .or. &
444  coordinatemode(coord_mode) == regridding_hycom1 .or. &
445  coordinatemode(coord_mode) == regridding_slight .or. &
446  coordinatemode(coord_mode) == regridding_adaptive) then
447  ! Adjust target grid to be consistent with maximum_depth
448  tmpreal = sum( dz(:) )
449  if (tmpreal < maximum_depth) then
450  dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
451  elseif (tmpreal > maximum_depth) then
452  if ( dz(ke) + ( maximum_depth - tmpreal ) > 0. ) then
453  dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
454  else
455  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
456  "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
457  endif
458  endif
459  endif
460  endif
461 
462  cs%nk=ke
463 
464  ! Target resolution (for fixed coordinates)
465  allocate( cs%coordinateResolution(cs%nk) ); cs%coordinateResolution(:) = -1.e30
466  if (state_dependent(cs%regridding_scheme)) then
467  ! Target values
468  allocate( cs%target_density(cs%nk+1) ); cs%target_density(:) = -1.e30*us%kg_m3_to_R
469  endif
470 
471  if (allocated(dz)) then
472  if (coordinatemode(coord_mode) == regridding_sigma) then
473  call setcoordinateresolution(dz, cs, scale=1.0)
474  elseif (coordinatemode(coord_mode) == regridding_rho) then
475  call setcoordinateresolution(dz, cs, scale=us%kg_m3_to_R)
476  cs%coord_scale = us%R_to_kg_m3
477  elseif (coordinatemode(coord_mode) == regridding_adaptive) then
478  call setcoordinateresolution(dz, cs, scale=gv%m_to_H)
479  cs%coord_scale = gv%H_to_m
480  else
481  call setcoordinateresolution(dz, cs, scale=us%m_to_Z)
482  cs%coord_scale = us%Z_to_m
483  endif
484  endif
485 
486  if (allocated(rho_target)) then
487  call set_target_densities(cs, us%kg_m3_to_R*rho_target)
488  deallocate(rho_target)
489 
490  ! \todo This line looks like it would overwrite the target densities set just above?
491  elseif (coordinatemode(coord_mode) == regridding_rho) then
492  call set_target_densities_from_gv(gv, us, cs)
493  call log_param(param_file, mdl, "!TARGET_DENSITIES", us%R_to_kg_m3*cs%target_density(:), &
494  'RHO target densities for interfaces', units=coordinateunits(coord_mode))
495  endif
496 
497  ! initialise coordinate-specific control structure
498  call initcoord(cs, gv, us, coord_mode)
499 
500  if (main_parameters .and. coord_is_state_dependent) then
501  call get_param(param_file, mdl, "REGRID_COMPRESSIBILITY_FRACTION", tmpreal, &
502  "When interpolating potential density profiles we can add "//&
503  "some artificial compressibility solely to make homogeneous "//&
504  "regions appear stratified.", default=0.)
505  call set_regrid_params(cs, compress_fraction=tmpreal)
506  endif
507 
508  if (main_parameters) then
509  call get_param(param_file, mdl, "MIN_THICKNESS", tmpreal, &
510  "When regridding, this is the minimum layer "//&
511  "thickness allowed.", units="m", scale=gv%m_to_H, &
513  call set_regrid_params(cs, min_thickness=tmpreal)
514  else
515  call set_regrid_params(cs, min_thickness=0.)
516  endif
517 
518  if (coordinatemode(coord_mode) == regridding_slight) then
519  ! Set SLight-specific regridding parameters.
520  call get_param(param_file, mdl, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, &
521  "The nominal thickness of fixed thickness near-surface "//&
522  "layers with the SLight coordinate.", units="m", default=1.0, scale=gv%m_to_H)
523  call get_param(param_file, mdl, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, &
524  "The number of fixed-depth surface layers with the SLight "//&
525  "coordinate.", units="nondimensional", default=2)
526  call get_param(param_file, mdl, "SLIGHT_SURFACE_AVG_DEPTH", rho_avg_depth, &
527  "The thickness of the surface region over which to average "//&
528  "when calculating the density to use to define the interior "//&
529  "with the SLight coordinate.", units="m", default=1.0, scale=gv%m_to_H)
530  call get_param(param_file, mdl, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, &
531  "The number of layers to offset the surface density when "//&
532  "defining where the interior ocean starts with SLight.", &
533  units="nondimensional", default=2.0)
534  call get_param(param_file, mdl, "SLIGHT_FIX_HALOCLINES", fix_haloclines, &
535  "If true, identify regions above the reference pressure "//&
536  "where the reference pressure systematically underestimates "//&
537  "the stratification and use this in the definition of the "//&
538  "interior with the SLight coordinate.", default=.false.)
539 
540  call set_regrid_params(cs, dz_min_surface=dz_fixed_sfc, &
541  nz_fixed_surface=nz_fixed_sfc, rho_ml_avg_depth=rho_avg_depth, &
542  nlay_ml_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines)
543  if (fix_haloclines) then
544  ! Set additional parameters related to SLIGHT_FIX_HALOCLINES.
545  call get_param(param_file, mdl, "HALOCLINE_FILTER_LENGTH", filt_len, &
546  "A length scale over which to smooth the temperature and "//&
547  "salinity before identifying erroneously unstable haloclines.", &
548  units="m", default=2.0)
549  call get_param(param_file, mdl, "HALOCLINE_STRAT_TOL", strat_tol, &
550  "A tolerance for the ratio of the stratification of the "//&
551  "apparent coordinate stratification to the actual value "//&
552  "that is used to identify erroneously unstable haloclines. "//&
553  "This ratio is 1 when they are equal, and sensible values "//&
554  "are between 0 and 0.5.", units="nondimensional", default=0.2)
555  call set_regrid_params(cs, halocline_filt_len=filt_len, &
556  halocline_strat_tol=strat_tol)
557  endif
558 
559  endif
560 
561  if (coordinatemode(coord_mode) == regridding_adaptive) then
562  call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adapttimeratio, &
563  "Ratio of ALE timestep to grid timescale.", units="s", default=1e-1) !### Should the units be "nondim"?
564  call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptzoom, &
565  "Depth of near-surface zooming region.", units="m", default=200.0, scale=gv%m_to_H)
566  call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptzoomcoeff, &
567  "Coefficient of near-surface zooming diffusivity.", &
568  units="nondim", default=0.2)
569  call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptbuoycoeff, &
570  "Coefficient of buoyancy diffusivity.", &
571  units="nondim", default=0.8)
572  call get_param(param_file, mdl, "ADAPT_ALPHA", adaptalpha, &
573  "Scaling on optimization tendency.", &
574  units="nondim", default=1.0)
575  call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmplogical, &
576  "If true, make a HyCOM-like mixed layer by preventing interfaces "//&
577  "from being shallower than the depths specified by the regridding coordinate.", &
578  default=.false.)
579 
580  call set_regrid_params(cs, adapttimeratio=adapttimeratio, adaptzoom=adaptzoom, &
581  adaptzoomcoeff=adaptzoomcoeff, adaptbuoycoeff=adaptbuoycoeff, adaptalpha=adaptalpha, &
582  adaptdomin=tmplogical)
583  endif
584 
585  if (main_parameters .and. coord_is_state_dependent) then
586  call get_param(param_file, mdl, "MAXIMUM_INT_DEPTH_CONFIG", string, &
587  "Determines how to specify the maximum interface depths.\n"//&
588  "Valid options are:\n"//&
589  " NONE - there are no maximum interface depths\n"//&
590  " PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
591  " FILE:string - read from a file. The string specifies\n"//&
592  " the filename and variable name, separated\n"//&
593  " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
594  " FNC1:string - FNC1:dz_min,H_total,power,precision",&
595  default='NONE')
596  message = "The list of maximum depths for each interface."
597  allocate(z_max(ke+1))
598  allocate(dz_max(ke))
599  if ( trim(string) == "NONE") then
600  ! Do nothing.
601  elseif ( trim(string) == "PARAM") then
602  call get_param(param_file, mdl, "MAXIMUM_INTERFACE_DEPTHS", z_max, &
603  trim(message), units="m", scale=gv%m_to_H, fail_if_missing=.true.)
604  call set_regrid_max_depths(cs, z_max)
605  elseif (index(trim(string),'FILE:')==1) then
606  if (string(6:6)=='.' .or. string(6:6)=='/') then
607  ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
608  filename = trim( extractword(trim(string(6:80)), 1) )
609  else
610  ! Otherwise assume we should look for the file in INPUTDIR
611  filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
612  endif
613  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
614  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
615 
616  do_sum = .false.
617  varname = trim( extractword(trim(string(6:)), 2) )
618  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
619  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
620  if (len_trim(varname)==0) then
621  if (field_exists(filename,'z_max')) then; varname = 'z_max'
622  elseif (field_exists(filename,'dz')) then; varname = 'dz' ; do_sum = .true.
623  elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max' ; do_sum = .true.
624  else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
625  "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
626  endif
627  endif
628  if (do_sum) then
629  call mom_read_data(trim(filename), trim(varname), dz_max)
630  z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
631  else
632  call mom_read_data(trim(filename), trim(varname), z_max)
633  endif
634  call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
635  trim(message), units=coordinateunits(coord_mode))
636  call set_regrid_max_depths(cs, z_max, gv%m_to_H)
637  elseif (index(trim(string),'FNC1:')==1) then
638  call dz_function1( trim(string(6:)), dz_max )
639  if ((coordinatemode(coord_mode) == regridding_slight) .and. &
640  (dz_fixed_sfc > 0.0)) then
641  do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo
642  endif
643  z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
644  call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
645  trim(message), units=coordinateunits(coord_mode))
646  call set_regrid_max_depths(cs, z_max, gv%m_to_H)
647  else
648  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
649  "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
650  endif
651  deallocate(z_max)
652  deallocate(dz_max)
653 
654  ! Optionally specify maximum thicknesses for each layer, enforced by moving
655  ! the interface below a layer downward.
656  call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", string, &
657  "Determines how to specify the maximum layer thicknesses.\n"//&
658  "Valid options are:\n"//&
659  " NONE - there are no maximum layer thicknesses\n"//&
660  " PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
661  " FILE:string - read from a file. The string specifies\n"//&
662  " the filename and variable name, separated\n"//&
663  " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
664  " FNC1:string - FNC1:dz_min,H_total,power,precision",&
665  default='NONE')
666  message = "The list of maximum thickness for each layer."
667  allocate(h_max(ke))
668  if ( trim(string) == "NONE") then
669  ! Do nothing.
670  elseif ( trim(string) == "PARAM") then
671  call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, &
672  trim(message), units="m", fail_if_missing=.true., scale=gv%m_to_H)
673  call set_regrid_max_thickness(cs, h_max)
674  elseif (index(trim(string),'FILE:')==1) then
675  if (string(6:6)=='.' .or. string(6:6)=='/') then
676  ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
677  filename = trim( extractword(trim(string(6:80)), 1) )
678  else
679  ! Otherwise assume we should look for the file in INPUTDIR
680  filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
681  endif
682  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
683  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
684 
685  varname = trim( extractword(trim(string(6:)), 2) )
686  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
687  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
688  if (len_trim(varname)==0) then
689  if (field_exists(filename,'h_max')) then; varname = 'h_max'
690  elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max'
691  else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
692  "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
693  endif
694  endif
695  call mom_read_data(trim(filename), trim(varname), h_max)
696  call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
697  trim(message), units=coordinateunits(coord_mode))
698  call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
699  elseif (index(trim(string),'FNC1:')==1) then
700  call dz_function1( trim(string(6:)), h_max )
701  call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
702  trim(message), units=coordinateunits(coord_mode))
703  call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
704  else
705  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
706  "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
707  endif
708  deallocate(h_max)
709  endif
710 
711  if (allocated(dz)) deallocate(dz)
712 end subroutine initialize_regridding
713 
714 !> Do some basic checks on the vertical grid definition file, variable
715 subroutine check_grid_def(filename, varname, expected_units, msg, ierr)
716  character(len=*), intent(in) :: filename !< File name
717  character(len=*), intent(in) :: varname !< Variable name
718  character(len=*), intent(in) :: expected_units !< Expected units of variable
719  character(len=*), intent(inout) :: msg !< Message to use for errors
720  logical, intent(out) :: ierr !< True if an error occurs
721  ! Local variables
722  character (len=200) :: units, long_name
723  integer :: ncid, status, intid, vid
724  integer :: i
725 
726  ierr = .false.
727  status = nf90_open(trim(filename), nf90_nowrite, ncid)
728  if (status /= nf90_noerr) then
729  ierr = .true.
730  msg = 'File not found: '//trim(filename)
731  return
732  endif
733 
734  status = nf90_inq_varid(ncid, trim(varname), vid)
735  if (status /= nf90_noerr) then
736  ierr = .true.
737  msg = 'Var not found: '//trim(varname)
738  return
739  endif
740 
741  status = nf90_get_att(ncid, vid, "units", units)
742  if (status /= nf90_noerr) then
743  ierr = .true.
744  msg = 'Attribute not found: units'
745  return
746  endif
747  ! NF90_GET_ATT can return attributes with null characters, which TRIM will not truncate.
748  ! This loop replaces any null characters with a space so that the following check between
749  ! the read units and the expected units will pass
750  do i=1,len_trim(units)
751  if (units(i:i) == char(0)) units(i:i) = " "
752  enddo
753 
754  if (trim(units) /= trim(expected_units)) then
755  if (trim(expected_units) == "meters") then
756  if (trim(units) /= "m") then
757  ierr = .true.
758  endif
759  else
760  ierr = .true.
761  endif
762  endif
763 
764  if (ierr) then
765  msg = 'Units incorrect: '//trim(units)//' /= '//trim(expected_units)
766  endif
767 
768 end subroutine check_grid_def
769 
770 !> Deallocation of regridding memory
771 subroutine end_regridding(CS)
772  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
773 
774  if (associated(cs%zlike_CS)) call end_coord_zlike(cs%zlike_CS)
775  if (associated(cs%sigma_CS)) call end_coord_sigma(cs%sigma_CS)
776  if (associated(cs%rho_CS)) call end_coord_rho(cs%rho_CS)
777  if (associated(cs%hycom_CS)) call end_coord_hycom(cs%hycom_CS)
778  if (associated(cs%slight_CS)) call end_coord_slight(cs%slight_CS)
779  if (associated(cs%adapt_CS)) call end_coord_adapt(cs%adapt_CS)
780 
781  deallocate( cs%coordinateResolution )
782  if (allocated(cs%target_density)) deallocate( cs%target_density )
783  if (allocated(cs%max_interface_depths) ) deallocate( cs%max_interface_depths )
784  if (allocated(cs%max_layer_thickness) ) deallocate( cs%max_layer_thickness )
785 
786 end subroutine end_regridding
787 
788 !------------------------------------------------------------------------------
789 !> Dispatching regridding routine for orchestrating regridding & remapping
790 subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
791 !------------------------------------------------------------------------------
792 ! This routine takes care of (1) building a new grid and (2) remapping between
793 ! the old grid and the new grid. The creation of the new grid can be based
794 ! on z coordinates, target interface densities, sigma coordinates or any
795 ! arbitrary coordinate system.
796 ! The MOM6 interface positions are always calculated from the bottom up by
797 ! accumulating the layer thicknesses starting at z=-G%bathyT. z increases
798 ! upwards (decreasing k-index).
799 ! The new grid is defined by the change in position of those interfaces in z
800 ! dzInterface = zNew - zOld.
801 ! Thus, if the regridding inflates the top layer, hNew(1) > hOld(1), then the
802 ! second interface moves downward, zNew(2) < zOld(2), and dzInterface(2) < 0.
803 ! hNew(k) = hOld(k) - dzInterface(k+1) + dzInterface(k)
804 ! IMPORTANT NOTE:
805 ! This is the converse of the sign convention used in the remapping code!
806 !------------------------------------------------------------------------------
807 
808  ! Arguments
809  type(remapping_cs), intent(in) :: remapcs !< Remapping parameters and options
810  type(regridding_cs), intent(in) :: cs !< Regridding control structure
811  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
812  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
813  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after
814  !! the last time step
815  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...)
816  real, dimension(SZI_(G),SZJ_(G), CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate
817  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzinterface !< The change in position of each interface
818  real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
819  logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment
820  ! Local variables
821  real :: trickgnucompiler
822  logical :: use_ice_shelf
823  logical :: do_convective_adjustment
824 
825  do_convective_adjustment = .true.
826  if (present(conv_adjust)) do_convective_adjustment = conv_adjust
827 
828  use_ice_shelf = .false.
829  if (present(frac_shelf_h)) then
830  if (associated(frac_shelf_h)) use_ice_shelf = .true.
831  endif
832 
833  select case ( cs%regridding_scheme )
834 
835  case ( regridding_zstar )
836  if (use_ice_shelf) then
837  call build_zstar_grid( cs, g, gv, h, dzinterface, frac_shelf_h )
838  else
839  call build_zstar_grid( cs, g, gv, h, dzinterface )
840  endif
841  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
842 
843  case ( regridding_sigma_shelf_zstar)
844  call build_zstar_grid( cs, g, gv, h, dzinterface )
845  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
846 
847  case ( regridding_sigma )
848  call build_sigma_grid( cs, g, gv, h, dzinterface )
849  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
850 
851  case ( regridding_rho )
852  if (do_convective_adjustment) call convective_adjustment(g, gv, h, tv)
853  call build_rho_grid( g, gv, h, tv, dzinterface, remapcs, cs )
854  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
855 
856  case ( regridding_arbitrary )
857  call build_grid_arbitrary( g, gv, h, dzinterface, trickgnucompiler, cs )
858  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
859 
860  case ( regridding_hycom1 )
861  call build_grid_hycom1( g, gv, h, tv, h_new, dzinterface, cs )
862 
863  case ( regridding_slight )
864  call build_grid_slight( g, gv, h, tv, dzinterface, cs )
865  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
866 
867  case ( regridding_adaptive )
868  call build_grid_adaptive(g, gv, h, tv, dzinterface, remapcs, cs)
869  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
870 
871  case default
872  call mom_error(fatal,'MOM_regridding, regridding_main: '//&
873  'Unknown regridding scheme selected!')
874 
875  end select ! type of grid
876 
877 #ifdef __DO_SAFETY_CHECKS__
878  call check_remapping_grid(g, gv, h, dzinterface,'in regridding_main')
879 #endif
880 
881 end subroutine regridding_main
882 
883 !> Calculates h_new from h + delta_k dzInterface
884 subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
885  type(regridding_cs), intent(in) :: CS !< Regridding control structure
886  type(ocean_grid_type), intent(in) :: G !< Grid structure
887  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
888  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Old layer thicknesses (arbitrary units)
889  real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: dzInterface !< Change in interface positions (same as h)
890  real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses (same as h)
891  ! Local variables
892  integer :: i, j, k, nki
893 
894  nki = min(cs%nk, gv%ke)
895 
896  !$OMP parallel do default(shared)
897  do j = g%jsc-1,g%jec+1
898  do i = g%isc-1,g%iec+1
899  if (g%mask2dT(i,j)>0.) then
900  do k=1,nki
901  h_new(i,j,k) = max( 0., h(i,j,k) + ( dzinterface(i,j,k) - dzinterface(i,j,k+1) ) )
902  enddo
903  if (cs%nk > gv%ke) then
904  do k=nki+1, cs%nk
905  h_new(i,j,k) = max( 0., dzinterface(i,j,k) - dzinterface(i,j,k+1) )
906  enddo
907  endif
908  else
909  h_new(i,j,1:nki) = h(i,j,1:nki)
910  if (cs%nk > gv%ke) h_new(i,j,nki+1:cs%nk) = 0.
911  ! On land points, why are we keeping the original h rather than setting to zero? -AJA
912  endif
913  enddo
914  enddo
915 
916 end subroutine calc_h_new_by_dz
917 
918 !> Check that the total thickness of two grids match
919 subroutine check_remapping_grid( G, GV, h, dzInterface, msg )
920  type(ocean_grid_type), intent(in) :: g !< Grid structure
921  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
922  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
923  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzinterface !< Change in interface positions
924  !! [H ~> m or kg m-2]
925  character(len=*), intent(in) :: msg !< Message to append to errors
926  ! Local variables
927  integer :: i, j
928 
929  !$OMP parallel do default(shared)
930  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
931  if (g%mask2dT(i,j)>0.) call check_grid_column( gv%ke, gv%Z_to_H*g%bathyT(i,j), h(i,j,:), dzinterface(i,j,:), msg )
932  enddo ; enddo
933 
934 end subroutine check_remapping_grid
935 
936 !> Check that the total thickness of new and old grids are consistent
937 subroutine check_grid_column( nk, depth, h, dzInterface, msg )
938  integer, intent(in) :: nk !< Number of cells
939  real, intent(in) :: depth !< Depth of bottom [Z ~> m] or arbitrary units
940  real, dimension(nk), intent(in) :: h !< Cell thicknesses [Z ~> m] or arbitrary units
941  real, dimension(nk+1), intent(in) :: dzinterface !< Change in interface positions (same units as h)
942  character(len=*), intent(in) :: msg !< Message to append to errors
943  ! Local variables
944  integer :: k
945  real :: eps, total_h_old, total_h_new, h_new, z_old, z_new
946 
947  eps =1. ; eps = epsilon(eps)
948 
949  ! Total thickness of grid h
950  total_h_old = 0.
951  do k = 1,nk
952  total_h_old = total_h_old + h(k)
953  enddo
954 
955  ! Integrate upwards for the interfaces consistent with the rest of MOM6
956  z_old = - depth
957  if (depth == 0.) z_old = - total_h_old
958  total_h_new = 0.
959  do k = nk,1,-1
960  z_old = z_old + h(k) ! Old interface position above layer k
961  z_new = z_old + dzinterface(k) ! New interface position based on dzInterface
962  h_new = h(k) + ( dzinterface(k) - dzinterface(k+1) ) ! New thickness
963  if (h_new<0.) then
964  write(0,*) 'k,h,hnew=',k,h(k),h_new
965  write(0,*) 'dzI(k+1),dzI(k)=',dzinterface(k+1),dzinterface(k)
966  call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
967  'Negative layer thickness implied by re-gridding, '//trim(msg))
968  endif
969  total_h_new = total_h_new + h_new
970 
971  enddo
972 
973  ! Conservation by implied h_new
974  if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps) then
975  write(0,*) 'nk=',nk
976  do k = 1,nk
977  write(0,*) 'k,h,hnew=',k,h(k),h(k)+(dzinterface(k)-dzinterface(k+1))
978  enddo
979  write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
980  write(0,*) 'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
981  call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
982  'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
983  endif
984 
985  ! Check that the top and bottom are intentionally moving
986  if (dzinterface(1) /= 0.) call mom_error( fatal, &
987  'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
988  if (dzinterface(nk+1) /= 0.) call mom_error( fatal, &
989  'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
990 
991 end subroutine check_grid_column
992 
993 !> Returns the change in interface position motion after filtering and
994 !! assuming the top and bottom interfaces do not move. The filtering is
995 !! a function of depth, and is applied as the integrated average filtering
996 !! over the trajectory of the interface. By design, this code can not give
997 !! tangled interfaces provided that z_old and z_new are not already tangled.
998 subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
999  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1000  integer, intent(in) :: nk !< Number of cells in source grid
1001  real, dimension(nk+1), intent(in) :: z_old !< Old grid position [m]
1002  real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position [m]
1003  real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [m]
1004  ! Local variables
1005  real :: sgn ! The sign convention for downward.
1006  real :: dz_tgt, zr1, z_old_k
1007  real :: Aq, Bq, dz0, z0, F0
1008  real :: zs, zd, dzwt, Idzwt
1009  real :: wtd, Iwtd
1010  real :: Int_zs, Int_zd, dInt_zs_zd
1011 ! For debugging:
1012  real, dimension(nk+1) :: z_act
1013 ! real, dimension(nk+1) :: ddz_g_s, ddz_g_d
1014  logical :: debug = .false.
1015  integer :: k
1016 
1017  if ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) < 0.0) then
1018  call mom_error(fatal, "filtered_grid_motion: z_old and z_new use different sign conventions.")
1019  elseif ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) == 0.0) then
1020  ! This is a massless column, so do nothing and return.
1021  do k=1,cs%nk+1 ; dz_g(k) = 0.0 ; enddo ; return
1022  elseif ((z_old(nk+1) - z_old(1)) + (z_new(cs%nk+1) - z_new(1)) > 0.0) then
1023  sgn = 1.0
1024  else
1025  sgn = -1.0
1026  endif
1027 
1028  if (debug) then
1029  do k=2,cs%nk+1
1030  if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
1031  call mom_error(fatal, "filtered_grid_motion: z_new is tangled.")
1032  enddo
1033  do k=2,nk+1
1034  if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
1035  call mom_error(fatal, "filtered_grid_motion: z_old is tangled.")
1036  enddo
1037  ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0
1038  endif
1039 
1040  zs = cs%depth_of_time_filter_shallow
1041  zd = cs%depth_of_time_filter_deep
1042  wtd = 1.0 - cs%old_grid_weight
1043  iwtd = 1.0 / wtd
1044 
1045  dzwt = (zd - zs)
1046  idzwt = 0.0 ; if (abs(zd - zs) > 0.0) idzwt = 1.0 / (zd - zs)
1047  dint_zs_zd = 0.5*(1.0 + iwtd) * (zd - zs)
1048  aq = 0.5*(iwtd - 1.0)
1049 
1050  dz_g(1) = 0.0
1051  z_old_k = z_old(1)
1052  do k = 2,cs%nk+1
1053  if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1054  ! zr1 is positive and increases with depth, and dz_tgt is positive downward.
1055  dz_tgt = sgn*(z_new(k) - z_old_k)
1056  zr1 = sgn*(z_old_k - z_old(1))
1057 
1058  ! First, handle the two simple and common cases that do not pass through
1059  ! the adjustment rate transition zone.
1060  if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then
1061  dz_g(k) = sgn * wtd * dz_tgt
1062  elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then
1063  dz_g(k) = sgn * dz_tgt
1064  else
1065  ! Find the new value by inverting the equation
1066  ! integral(0 to dz_new) Iwt(z) dz = dz_tgt
1067  ! This is trivial where Iwt is a constant, and agrees with the two limits above.
1068 
1069  ! Take test values at the transition points to figure out which segment
1070  ! the new value will be found in.
1071  if (zr1 >= zd) then
1072  int_zd = iwtd*(zd - zr1)
1073  int_zs = int_zd - dint_zs_zd
1074  elseif (zr1 <= zs) then
1075  int_zs = (zs - zr1)
1076  int_zd = dint_zs_zd + (zs - zr1)
1077  else
1078 ! Int_zd = (zd - zr1) * (Iwtd + 0.5*(1.0 - Iwtd) * (zd - zr1) / (zd - zs))
1079  int_zd = (zd - zr1) * (iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * idzwt
1080  int_zs = (zs - zr1) * (0.5*iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * idzwt
1081  ! It has been verified that Int_zs = Int_zd - dInt_zs_zd to within roundoff.
1082  endif
1083 
1084  if (dz_tgt >= int_zd) then ! The new location is in the deep, slow region.
1085  dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - int_zd))
1086  elseif (dz_tgt <= int_zs) then ! The new location is in the shallow region.
1087  dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - int_zs))
1088  else ! We need to solve a quadratic equation for z_new.
1089  ! For accuracy, do the integral from the starting depth or the nearest
1090  ! edge of the transition region. The results with each choice are
1091  ! mathematically equivalent, but differ in roundoff, and this choice
1092  ! should minimize the likelihood of inadvertently overlapping interfaces.
1093  if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; f0 = dz_tgt - int_zs
1094  elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; f0 = dz_tgt - int_zd
1095  else ; dz0 = 0.0 ; z0 = zr1 ; f0 = dz_tgt ; endif
1096 
1097  bq = (dzwt + 2.0*aq*(z0-zs))
1098  ! Solve the quadratic: Aq*(zn-z0)**2 + Bq*(zn-z0) - F0*dzwt = 0
1099  ! Note that b>=0, and the two terms in the standard form cancel for the right root.
1100  dz_g(k) = sgn * (dz0 + 2.0*f0*dzwt / (bq + sqrt(bq**2 + 4.0*aq*f0*dzwt) ))
1101 
1102 ! if (debug) then
1103 ! dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0*Aq*(z0-zs))
1104 ! ddz_g_s(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
1105 ! dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0*Aq*(z0-zs))
1106 ! ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
1107 !
1108 ! if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
1109 ! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).")
1110 ! if (abs(ddz_g_d(k) - ddz_g_s(k)) > 1e-12*(abs(dz_g(k)+ddz_g_d(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
1111 ! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.")
1112 ! endif
1113  endif
1114 
1115  endif
1116  enddo
1117  !dz_g(CS%nk+1) = 0.0
1118 
1119  if (debug) then
1120  z_old_k = z_old(1)
1121  do k=1,cs%nk+1
1122  if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1123  z_act(k) = z_old_k + dz_g(k)
1124  enddo
1125  do k=2,cs%nk+1
1126  if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) &
1127  call mom_error(fatal, "filtered_grid_motion: z_output is tangled.")
1128  enddo
1129  endif
1130 
1131 end subroutine filtered_grid_motion
1132 
1133 !> Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004).
1134 !! z* is defined as
1135 !! z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H .
1136 subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)
1138  ! Arguments
1139  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1140  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1141  type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
1142  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1143  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1144  !! [H ~> m or kg m-2].
1145  real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage.
1146  ! Local variables
1147  integer :: i, j, k
1148  integer :: nz
1149  real :: nominalDepth, totalThickness, dh
1150  real, dimension(SZK_(GV)+1) :: zOld, zNew
1151  real :: minThickness
1152  logical :: ice_shelf
1153 
1154  nz = gv%ke
1155  minthickness = cs%min_thickness
1156  ice_shelf = .false.
1157  if (present(frac_shelf_h)) then
1158  if (associated(frac_shelf_h)) ice_shelf = .true.
1159  endif
1160 
1161 !$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, &
1162 !$OMP ice_shelf,minThickness) &
1163 !$OMP private(nominalDepth,totalThickness, &
1164 !$OMP zNew,dh,zOld)
1165  do j = g%jsc-1,g%jec+1
1166  do i = g%isc-1,g%iec+1
1167 
1168  if (g%mask2dT(i,j)==0.) then
1169  dzinterface(i,j,:) = 0.
1170  cycle
1171  endif
1172 
1173  ! Local depth (G%bathyT is positive)
1174  nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1175 
1176  ! Determine water column thickness
1177  totalthickness = 0.0
1178  do k = 1,nz
1179  totalthickness = totalthickness + h(i,j,k)
1180  enddo
1181 
1182  zold(nz+1) = - nominaldepth
1183  do k = nz,1,-1
1184  zold(k) = zold(k+1) + h(i,j,k)
1185  enddo
1186 
1187  if (ice_shelf) then
1188  if (frac_shelf_h(i,j) > 0.) then ! under ice shelf
1189  call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, znew, &
1190  z_rigid_top = totalthickness-nominaldepth, &
1191  eta_orig=zold(1), zscale=gv%Z_to_H)
1192  else
1193  call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1194  znew, zscale=gv%Z_to_H)
1195  endif
1196  else
1197  call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1198  znew, zscale=gv%Z_to_H)
1199  endif
1200 
1201  ! Calculate the final change in grid position after blending new and old grids
1202  call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1203 
1204 #ifdef __DO_SAFETY_CHECKS__
1205  dh=max(nominaldepth,totalthickness)
1206  if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
1207  write(0,*) 'min_thickness=',minthickness
1208  write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1209  write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz
1210  do k=1,nz+1
1211  write(0,*) k,zold(k),znew(k)
1212  enddo
1213  do k=1,nz
1214  write(0,*) k,h(i,j,k),znew(k)-znew(k+1),cs%coordinateResolution(k)
1215  enddo
1216  call mom_error( fatal, &
1217  'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
1218  endif
1219 #endif
1220 
1221  call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
1222 
1223  enddo
1224  enddo
1225 
1226 end subroutine build_zstar_grid
1227 
1228 !------------------------------------------------------------------------------
1229 ! Build sigma grid
1230 !> This routine builds a grid based on terrain-following coordinates.
1231 subroutine build_sigma_grid( CS, G, GV, h, dzInterface )
1232 !------------------------------------------------------------------------------
1233 ! This routine builds a grid based on terrain-following coordinates.
1234 ! The module parameter coordinateResolution(:) determines the resolution in
1235 ! sigma coordinate, dSigma(:). sigma-coordinates are defined by
1236 ! sigma = (eta-z)/(H+eta) s.t. sigma=0 at z=eta and sigma=1 at z=-H .
1237 !------------------------------------------------------------------------------
1238 
1239  ! Arguments
1240  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1241  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1242  type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
1243  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1244  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1245  !! [H ~> m or kg m-2]
1246 
1247  ! Local variables
1248  integer :: i, j, k
1249  integer :: nz
1250  real :: nominalDepth, totalThickness, dh
1251  real, dimension(SZK_(GV)+1) :: zOld, zNew
1252 
1253  nz = gv%ke
1254 
1255  do i = g%isc-1,g%iec+1
1256  do j = g%jsc-1,g%jec+1
1257 
1258  if (g%mask2dT(i,j)==0.) then
1259  dzinterface(i,j,:) = 0.
1260  cycle
1261  endif
1262 
1263  ! The rest of the model defines grids integrating up from the bottom
1264  nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1265 
1266  ! Determine water column height
1267  totalthickness = 0.0
1268  do k = 1,nz
1269  totalthickness = totalthickness + h(i,j,k)
1270  enddo
1271 
1272  call build_sigma_column(cs%sigma_CS, nominaldepth, totalthickness, znew)
1273 
1274  ! Calculate the final change in grid position after blending new and old grids
1275  zold(nz+1) = -nominaldepth
1276  do k = nz,1,-1
1277  zold(k) = zold(k+1) + h(i, j, k)
1278  enddo
1279 
1280  call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1281 
1282 #ifdef __DO_SAFETY_CHECKS__
1283  dh=max(nominaldepth,totalthickness)
1284  if (abs(znew(1)-zold(1))>(cs%nk-1)*0.5*epsilon(dh)*dh) then
1285  write(0,*) 'min_thickness=',cs%min_thickness
1286  write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1287  write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz,cs%nk
1288  do k=1,nz+1
1289  write(0,*) k,zold(k),znew(k)
1290  enddo
1291  do k=1,cs%nk
1292  write(0,*) k,h(i,j,k),znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k),cs%coordinateResolution(k)
1293  enddo
1294  call mom_error( fatal, &
1295  'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
1296  endif
1297  dzinterface(i,j,1) = 0.
1298  dzinterface(i,j,cs%nk+1) = 0.
1299 #endif
1300 
1301  enddo
1302  enddo
1303 
1304 end subroutine build_sigma_grid
1305 
1306 !------------------------------------------------------------------------------
1307 ! Build grid based on target interface densities
1308 !------------------------------------------------------------------------------
1309 !> This routine builds a new grid based on a given set of target interface densities.
1310 subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
1311 !------------------------------------------------------------------------------
1312 ! This routine builds a new grid based on a given set of target interface
1313 ! densities (these target densities are computed by taking the mean value
1314 ! of given layer densities). The algorithn operates as follows within each
1315 ! column:
1316 ! 1. Given T & S within each layer, the layer densities are computed.
1317 ! 2. Based on these layer densities, a global density profile is reconstructed
1318 ! (this profile is monotonically increasing and may be discontinuous)
1319 ! 3. The new grid interfaces are determined based on the target interface
1320 ! densities.
1321 ! 4. T & S are remapped onto the new grid.
1322 ! 5. Return to step 1 until convergence or until the maximum number of
1323 ! iterations is reached, whichever comes first.
1324 !------------------------------------------------------------------------------
1325 
1326  ! Arguments
1327  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1328  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1329  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1330  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1331  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
1332  !! [H ~> m or kg m-2]
1333  type(remapping_cs), intent(in) :: remapCS !< The remapping control structure
1334  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1335 
1336  ! Local variables
1337  integer :: nz
1338  integer :: i, j, k
1339  real :: nominalDepth, totalThickness
1340  real, dimension(SZK_(GV)+1) :: zOld, zNew
1341  real :: h_neglect, h_neglect_edge
1342 #ifdef __DO_SAFETY_CHECKS__
1343  real :: dh
1344 #endif
1345 
1346  !### Try replacing both of these with GV%H_subroundoff
1347  if (gv%Boussinesq) then
1348  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1349  else
1350  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1351  endif
1352 
1353  nz = gv%ke
1354 
1355  if (.not.cs%target_density_set) call mom_error(fatal, "build_rho_grid: "//&
1356  "Target densities must be set before build_rho_grid is called.")
1357 
1358  ! Build grid based on target interface densities
1359  do j = g%jsc-1,g%jec+1
1360  do i = g%isc-1,g%iec+1
1361 
1362  if (g%mask2dT(i,j)==0.) then
1363  dzinterface(i,j,:) = 0.
1364  cycle
1365  endif
1366 
1367 
1368  ! Local depth (G%bathyT is positive)
1369  nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1370 
1371  call build_rho_column(cs%rho_CS, nz, nominaldepth, h(i, j, :), &
1372  tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, znew, &
1373  h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1374 
1375  if (cs%integrate_downward_for_e) then
1376  zold(1) = 0.
1377  do k = 1,nz
1378  zold(k+1) = zold(k) - h(i,j,k)
1379  enddo
1380  else
1381  ! The rest of the model defines grids integrating up from the bottom
1382  zold(nz+1) = - nominaldepth
1383  do k = nz,1,-1
1384  zold(k) = zold(k+1) + h(i,j,k)
1385  enddo
1386  endif
1387 
1388  ! Calculate the final change in grid position after blending new and old grids
1389  call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1390 
1391 #ifdef __DO_SAFETY_CHECKS__
1392  do k = 2,nz
1393  if (znew(k) > zold(1)) then
1394  write(0,*) 'zOld=',zold
1395  write(0,*) 'zNew=',znew
1396  call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
1397  'interior interface above surface!' )
1398  endif
1399  if (znew(k) > znew(k-1)) then
1400  write(0,*) 'zOld=',zold
1401  write(0,*) 'zNew=',znew
1402  call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
1403  'interior interfaces cross!' )
1404  endif
1405  enddo
1406 
1407  totalthickness = 0.0
1408  do k = 1,nz
1409  totalthickness = totalthickness + h(i,j,k)
1410  enddo
1411 
1412  dh=max(nominaldepth,totalthickness)
1413  if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
1414  write(0,*) 'min_thickness=',cs%min_thickness
1415  write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1416  write(0,*) 'zNew(1)-zOld(1) = ',znew(1)-zold(1),epsilon(dh),nz
1417  do k=1,nz+1
1418  write(0,*) k,zold(k),znew(k)
1419  enddo
1420  do k=1,nz
1421  write(0,*) k,h(i,j,k),znew(k)-znew(k+1)
1422  enddo
1423  call mom_error( fatal, &
1424  'MOM_regridding, build_rho_grid: top surface has moved!!!' )
1425  endif
1426 #endif
1427 
1428  enddo ! end loop on i
1429  enddo ! end loop on j
1430 
1431 end subroutine build_rho_grid
1432 
1433 !> Builds a simple HyCOM-like grid with the deepest location of potential
1434 !! density interpolated from the column profile and a clipping of depth for
1435 !! each interface to a fixed z* or p* grid. This should probably be (optionally?)
1436 !! changed to find the nearest location of the target density.
1437 !! \remark { Based on Bleck, 2002: An oceanice general circulation model framed in
1438 !! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88.
1439 !! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 }
1440 subroutine build_grid_hycom1( G, GV, h, tv, h_new, dzInterface, CS )
1441  type(ocean_grid_type), intent(in) :: G !< Grid structure
1442  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1443  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2]
1444  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1445  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1446  real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
1447  real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position
1448 
1449  ! Local variables
1450  real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
1451  real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2]
1452  real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
1453  real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [Pa]
1454  integer :: i, j, k, nki
1455  real :: depth
1456  real :: h_neglect, h_neglect_edge
1457 
1458  !### Try replacing both of these with GV%H_subroundoff
1459  if (gv%Boussinesq) then
1460  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1461  else
1462  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1463  endif
1464 
1465  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_HyCOM1 : "//&
1466  "Target densities must be set before build_grid_HyCOM1 is called.")
1467 
1468  nki = min(gv%ke, cs%nk)
1469 
1470  ! Build grid based on target interface densities
1471  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1472  if (g%mask2dT(i,j)>0.) then
1473 
1474  depth = g%bathyT(i,j) * gv%Z_to_H
1475 
1476  z_col(1) = 0. ! Work downward rather than bottom up
1477  do k = 1, gv%ke
1478  z_col(k+1) = z_col(k) + h(i,j,k)
1479  p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1480  ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1481  enddo
1482 
1483  call build_hycom1_column(cs%hycom_CS, tv%eqn_of_state, gv%ke, depth, &
1484  h(i, j, :), tv%T(i, j, :), tv%S(i, j, :), p_col, &
1485  z_col, z_col_new, zscale=gv%Z_to_H, &
1486  h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1487 
1488  ! Calculate the final change in grid position after blending new and old grids
1489  call filtered_grid_motion( cs, gv%ke, z_col, z_col_new, dz_col )
1490 
1491  ! This adjusts things robust to round-off errors
1492  dz_col(:) = -dz_col(:)
1493  call adjust_interface_motion( cs, gv%ke, h(i,j,:), dz_col(:) )
1494 
1495  dzinterface(i,j,1:nki+1) = dz_col(1:nki+1)
1496  if (nki<cs%nk) dzinterface(i,j,nki+2:cs%nk+1) = 0.
1497 
1498  else ! on land
1499  dzinterface(i,j,:) = 0.
1500  endif ! mask2dT
1501  enddo ; enddo ! i,j
1502 
1503  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1504 
1505 end subroutine build_grid_hycom1
1506 
1507 !> This subroutine builds an adaptive grid that follows density surfaces where
1508 !! possible, subject to constraints on the smoothness of interface heights.
1509 subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS)
1510  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1511  type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1512  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1513  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
1514  !! thermodynamic variables
1515  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
1516  !! [H ~> m or kg m-2]
1517  type(remapping_cs), intent(in) :: remapCS !< The remapping control structure
1518  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1519 
1520  ! local variables
1521  integer :: i, j, k, nz ! indices and dimension lengths
1522  ! temperature, salinity and pressure on interfaces
1523  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt
1524  ! current interface positions and after tendency term is applied
1525  ! positive downward
1526  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt
1527  real, dimension(SZK_(GV)+1) :: zNext
1528 
1529  nz = gv%ke
1530 
1531  ! position surface at z = 0.
1532  zint(:,:,1) = 0.
1533 
1534  ! work on interior interfaces
1535  do k = 2, nz ; do j = g%jsc-2,g%jec+2 ; do i = g%isc-2,g%iec+2
1536  tint(i,j,k) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
1537  sint(i,j,k) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
1538  zint(i,j,k) = zint(i,j,k-1) + h(i,j,k-1) ! zInt in [H]
1539  enddo ; enddo ; enddo
1540 
1541  ! top and bottom temp/salt interfaces are just the layer
1542  ! average values
1543  tint(:,:,1) = tv%T(:,:,1) ; tint(:,:,nz+1) = tv%T(:,:,nz)
1544  sint(:,:,1) = tv%S(:,:,1) ; sint(:,:,nz+1) = tv%S(:,:,nz)
1545 
1546  ! set the bottom interface depth
1547  zint(:,:,nz+1) = zint(:,:,nz) + h(:,:,nz)
1548 
1549  ! calculate horizontal density derivatives (alpha/beta)
1550  ! between cells in a 5-point stencil, columnwise
1551  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1552  if (g%mask2dT(i,j) < 0.5) then
1553  dzinterface(i,j,:) = 0. ! land point, don't move interfaces, and skip
1554  cycle
1555  endif
1556 
1557  call build_adapt_column(cs%adapt_CS, g, gv, tv, i, j, zint, tint, sint, h, znext)
1558 
1559  call filtered_grid_motion(cs, nz, zint(i,j,:), znext, dzinterface(i,j,:))
1560  ! convert from depth to z
1561  do k = 1, nz+1 ; dzinterface(i,j,k) = -dzinterface(i,j,k) ; enddo
1562  call adjust_interface_motion(cs, nz, h(i,j,:), dzinterface(i,j,:))
1563  enddo ; enddo
1564 end subroutine build_grid_adaptive
1565 
1566 !> Builds a grid that tracks density interfaces for water that is denser than
1567 !! the surface density plus an increment of some number of layers, and uses all
1568 !! lighter layers uniformly above this location. Note that this amounts to
1569 !! interpolating to find the depth of an arbitrary (non-integer) interface index
1570 !! which should make the results vary smoothly in space to the extent that the
1571 !! surface density and interior stratification vary smoothly in space. Over
1572 !! shallow topography, this will tend to give a uniform sigma-like coordinate.
1573 !! For sufficiently shallow water, a minimum grid spacing is used to avoid
1574 !! certain instabilities.
1575 subroutine build_grid_slight(G, GV, h, tv, dzInterface, CS)
1576  type(ocean_grid_type), intent(in) :: G !< Grid structure
1577  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1578  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2]
1579  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1580  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position
1581  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1582 
1583  real, dimension(SZK_(GV)+1) :: z_col ! Interface positions relative to the surface [H ~> m or kg m-2]
1584  real, dimension(SZK_(GV)+1) :: z_col_new ! Interface positions relative to the surface [H ~> m or kg m-2]
1585  real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
1586  real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [Pa]
1587  real :: depth
1588  integer :: i, j, k, nz
1589  real :: h_neglect, h_neglect_edge
1590 
1591  !### Try replacing both of these with GV%H_subroundoff
1592  if (gv%Boussinesq) then
1593  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1594  else
1595  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1596  endif
1597 
1598  nz = gv%ke
1599 
1600  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_SLight : "//&
1601  "Target densities must be set before build_grid_SLight is called.")
1602 
1603  ! Build grid based on target interface densities
1604  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1605  if (g%mask2dT(i,j)>0.) then
1606 
1607  depth = g%bathyT(i,j) * gv%Z_to_H
1608  z_col(1) = 0. ! Work downward rather than bottom up
1609  do k=1,nz
1610  z_col(k+1) = z_col(k) + h(i,j,k)
1611  p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1612  ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1613  enddo
1614 
1615  call build_slight_column(cs%slight_CS, tv%eqn_of_state, gv%H_to_Pa, &
1616  gv%H_subroundoff, nz, depth, h(i, j, :), &
1617  tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, &
1618  h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1619 
1620  ! Calculate the final change in grid position after blending new and old grids
1621  call filtered_grid_motion( cs, nz, z_col, z_col_new, dz_col )
1622  do k=1,nz+1 ; dzinterface(i,j,k) = -dz_col(k) ; enddo
1623 #ifdef __DO_SAFETY_CHECKS__
1624  if (dzinterface(i,j,1) /= 0.) stop 'build_grid_SLight: Surface moved?!'
1625  if (dzinterface(i,j,nz+1) /= 0.) stop 'build_grid_SLight: Bottom moved?!'
1626 #endif
1627 
1628  ! This adjusts things robust to round-off errors
1629  call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
1630 
1631  else ! on land
1632  dzinterface(i,j,:) = 0.
1633  endif ! mask2dT
1634  enddo ; enddo ! i,j
1635 
1636 end subroutine build_grid_slight
1637 
1638 !> Adjust dz_Interface to ensure non-negative future thicknesses
1639 subroutine adjust_interface_motion( CS, nk, h_old, dz_int )
1640  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1641  integer, intent(in) :: nk !< Number of layers in h_old
1642  real, dimension(nk), intent(in) :: h_old !< Minium allowed thickness of h [H ~> m or kg m-2]
1643  real, dimension(CS%nk+1), intent(inout) :: dz_int !< Minium allowed thickness of h [H ~> m or kg m-2]
1644  ! Local variables
1645  integer :: k
1646  real :: h_new, eps, h_total, h_err
1647 
1648  eps = 1. ; eps = epsilon(eps)
1649 
1650  h_total = 0. ; h_err = 0.
1651  do k = 1, min(cs%nk,nk)
1652  h_total = h_total + h_old(k)
1653  h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1654  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1655  if (h_new < -3.0*h_err) then
1656  write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
1657  'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
1658  'h_new=',h_new,'h_err=',h_err
1659  call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
1660  'implied h<0 is larger than roundoff!')
1661  endif
1662  enddo
1663  if (cs%nk>nk) then
1664  do k = nk+1, cs%nk
1665  h_err = h_err + max( abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1666  h_new = ( dz_int(k) - dz_int(k+1) )
1667  if (h_new < -3.0*h_err) then
1668  write(0,*) 'h<0 at k=',k,'h_old was empty',&
1669  'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
1670  'h_new=',h_new,'h_err=',h_err
1671  call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
1672  'implied h<0 is larger than roundoff!')
1673  endif
1674  enddo
1675  endif
1676  do k = min(cs%nk,nk),2,-1
1677  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1678  if (h_new<cs%min_thickness) &
1679  dz_int(k) = ( dz_int(k+1) - h_old(k) ) + cs%min_thickness ! Implies next h_new = min_thickness
1680  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1681  if (h_new<0.) &
1682  dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) ) ! Backup in case min_thickness==0
1683  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1684  if (h_new<0.) then
1685  write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
1686  'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
1687  'h_new=',h_new
1688  stop 'Still did not work!'
1689  call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
1690  'Repeated adjustment for roundoff h<0 failed!')
1691  endif
1692  enddo
1693  !if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'
1694 
1695 end subroutine adjust_interface_motion
1696 
1697 !------------------------------------------------------------------------------
1698 ! Build arbitrary grid
1699 !------------------------------------------------------------------------------
1700 subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS )
1701 !------------------------------------------------------------------------------
1702 ! This routine builds a grid based on arbitrary rules
1703 !------------------------------------------------------------------------------
1704 
1705  ! Arguments
1706  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1707  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1708  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(in) :: h !< Original layer thicknesses [H ~> m or kg m-2]
1709  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface
1710  !! depth [H ~> m or kg m-2]
1711  real, intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
1712  type(regridding_cs), intent(in) :: CS !< Regridding control structure
1713 
1714  ! Local variables
1715  integer :: i, j, k
1716  integer :: nz
1717  real :: z_inter(SZK_(GV)+1)
1718  real :: total_height
1719  real :: delta_h
1720  real :: max_depth
1721  real :: eta ! local elevation
1722  real :: local_depth
1723  real :: x1, y1, x2, y2
1724  real :: x, t
1725 
1726  nz = gv%ke
1727  max_depth = g%max_depth*gv%Z_to_H
1728 
1729  do j = g%jsc-1,g%jec+1
1730  do i = g%isc-1,g%iec+1
1731 
1732  ! Local depth
1733  local_depth = g%bathyT(i,j)*gv%Z_to_H
1734 
1735  ! Determine water column height
1736  total_height = 0.0
1737  do k = 1,nz
1738  total_height = total_height + h(i,j,k)
1739  enddo
1740 
1741  eta = total_height - local_depth
1742 
1743  ! Compute new thicknesses based on stretched water column
1744  delta_h = (max_depth + eta) / nz
1745 
1746  ! Define interfaces
1747  z_inter(1) = eta
1748  do k = 1,nz
1749  z_inter(k+1) = z_inter(k) - delta_h
1750  enddo
1751 
1752  ! Refine grid in the middle
1753  do k = 1,nz+1
1754  x1 = 0.35; y1 = 0.45; x2 = 0.65; y2 = 0.55
1755 
1756  x = - ( z_inter(k) - eta ) / max_depth
1757 
1758  if ( x <= x1 ) then
1759  t = y1*x/x1
1760  elseif ( (x > x1 ) .and. ( x < x2 )) then
1761  t = y1 + (y2-y1) * (x-x1) / (x2-x1)
1762  else
1763  t = y2 + (1.0-y2) * (x-x2) / (1.0-x2)
1764  endif
1765 
1766  z_inter(k) = -t * max_depth + eta
1767 
1768  enddo
1769 
1770  ! Modify interface heights to account for topography
1771  z_inter(nz+1) = - local_depth
1772 
1773  ! Modify interface heights to avoid layers of zero thicknesses
1774  do k = nz,1,-1
1775  if ( z_inter(k) < (z_inter(k+1) + cs%min_thickness) ) then
1776  z_inter(k) = z_inter(k+1) + cs%min_thickness
1777  endif
1778  enddo
1779 
1780  ! Chnage in interface position
1781  x = 0. ! Left boundary at x=0
1782  dzinterface(i,j,1) = 0.
1783  do k = 2,nz
1784  x = x + h(i,j,k)
1785  dzinterface(i,j,k) = z_inter(k) - x
1786  enddo
1787  dzinterface(i,j,nz+1) = 0.
1788 
1789  enddo
1790  enddo
1791 
1792 stop 'OOOOOOPS' ! For some reason the gnu compiler will not let me delete this
1793  ! routine????
1794 
1795 end subroutine build_grid_arbitrary
1796 
1797 
1798 
1799 !------------------------------------------------------------------------------
1800 ! Check grid integrity
1801 !------------------------------------------------------------------------------
1802 subroutine inflate_vanished_layers_old( CS, G, GV, h )
1803 !------------------------------------------------------------------------------
1804 ! This routine is called when initializing the regridding options. The
1805 ! objective is to make sure all layers are at least as thick as the minimum
1806 ! thickness allowed for regridding purposes (this parameter is set in the
1807 ! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they
1808 ! are inflated up to the minmum thickness.
1809 !------------------------------------------------------------------------------
1810 
1811  ! Arguments
1812  type(regridding_cs), intent(in) :: cs !< Regridding control structure
1813  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1814  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1815  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
1816 
1817  ! Local variables
1818  integer :: i, j, k
1819  real :: htmp(gv%ke)
1820 
1821  do i = g%isc-1,g%iec+1
1822  do j = g%jsc-1,g%jec+1
1823 
1824  ! Build grid for current column
1825  do k = 1,gv%ke
1826  htmp(k) = h(i,j,k)
1827  enddo
1828 
1829  call old_inflate_layers_1d( cs%min_thickness, gv%ke, htmp )
1830 
1831  ! Save modified grid
1832  do k = 1,gv%ke
1833  h(i,j,k) = htmp(k)
1834  enddo
1835 
1836  enddo
1837  enddo
1838 
1839 end subroutine inflate_vanished_layers_old
1840 
1841 !------------------------------------------------------------------------------
1842 !> Achieve convective adjustment by swapping layers
1843 subroutine convective_adjustment(G, GV, h, tv)
1844  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1845  type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1846  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1847  intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
1848  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1849 !------------------------------------------------------------------------------
1850 ! Check each water column to see whether it is stratified. If not, sort the
1851 ! layers by successive swappings of water masses (bubble sort algorithm)
1852 !------------------------------------------------------------------------------
1853 
1854  ! Local variables
1855  integer :: i, j, k
1856  real :: T0, T1 ! temperatures
1857  real :: S0, S1 ! salinities
1858  real :: r0, r1 ! densities
1859  real :: h0, h1
1860  logical :: stratified
1861  real, dimension(GV%ke) :: p_col, densities
1862 
1863  p_col(:) = 0.
1864 
1865  ! Loop on columns
1866  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1867 
1868  ! Compute densities within current water column
1869  call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, &
1870  densities, 1, gv%ke, tv%eqn_of_state )
1871 
1872  ! Repeat restratification until complete
1873  do
1874  stratified = .true.
1875  do k = 1,gv%ke-1
1876  ! Gather information of current and next cells
1877  t0 = tv%T(i,j,k) ; t1 = tv%T(i,j,k+1)
1878  s0 = tv%S(i,j,k) ; s1 = tv%S(i,j,k+1)
1879  r0 = densities(k) ; r1 = densities(k+1)
1880  h0 = h(i,j,k) ; h1 = h(i,j,k+1)
1881  ! If the density of the current cell is larger than the density
1882  ! below it, we swap the cells and recalculate the densitiies
1883  ! within the swapped cells
1884  if ( r0 > r1 ) then
1885  tv%T(i,j,k) = t1 ; tv%T(i,j,k+1) = t0
1886  tv%S(i,j,k) = s1 ; tv%S(i,j,k+1) = s0
1887  h(i,j,k) = h1 ; h(i,j,k+1) = h0
1888  ! Recompute densities at levels k and k+1
1889  call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), &
1890  densities(k), tv%eqn_of_state )
1891  call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
1892  densities(k+1), tv%eqn_of_state )
1893  stratified = .false.
1894  endif
1895  enddo ! k
1896 
1897  if ( stratified ) exit
1898  enddo
1899 
1900  enddo ; enddo ! i & j
1901 
1902 end subroutine convective_adjustment
1903 
1904 
1905 !------------------------------------------------------------------------------
1906 !> Return a uniform resolution vector in the units of the coordinate
1907 function uniformresolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy)
1908 !------------------------------------------------------------------------------
1909 ! Calculate a vector of uniform resolution in the units of the coordinate
1910 !------------------------------------------------------------------------------
1911  ! Arguments
1912  integer, intent(in) :: nk !< Number of cells in source grid
1913  character(len=*), intent(in) :: coordmode !< A string indicating the coordinate mode.
1914  !! See the documenttion for regrid_consts
1915  !! for the recognized values.
1916  real, intent(in) :: maxdepth !< The range of the grid values in some modes
1917  real, intent(in) :: rholight !< The minimum value of the grid in RHO mode
1918  real, intent(in) :: rhoheavy !< The maximum value of the grid in RHO mode
1919 
1920  real :: uniformresolution(nk) !< The returned uniform resolution grid.
1921 
1922  ! Local variables
1923  integer :: scheme
1924 
1925  scheme = coordinatemode(coordmode)
1926  select case ( scheme )
1927 
1928  case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_sigma_shelf_zstar, &
1930  uniformresolution(:) = maxdepth / real(nk)
1931 
1932  case ( regridding_rho )
1933  uniformresolution(:) = (rhoheavy - rholight) / real(nk)
1934 
1935  case ( regridding_sigma )
1936  uniformresolution(:) = 1. / real(nk)
1937 
1938  case default
1939  call mom_error(fatal, "MOM_regridding, uniformResolution: "//&
1940  "Unrecognized choice for coordinate mode ("//trim(coordmode)//").")
1941 
1942  end select ! type of grid
1943 
1944 end function uniformresolution
1945 
1946 !> Initialize the coordinate resolutions by calling the appropriate initialization
1947 !! routine for the specified coordinate mode.
1948 subroutine initcoord(CS, GV, US, coord_mode)
1949  type(regridding_cs), intent(inout) :: CS !< Regridding control structure
1950  character(len=*), intent(in) :: coord_mode !< A string indicating the coordinate mode.
1951  !! See the documenttion for regrid_consts
1952  !! for the recognized values.
1953  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1954  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1955 
1956  select case (coordinatemode(coord_mode))
1957  case (regridding_zstar)
1958  call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
1959  case (regridding_sigma_shelf_zstar)
1960  call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
1961  case (regridding_sigma)
1962  call init_coord_sigma(cs%sigma_CS, cs%nk, cs%coordinateResolution)
1963  case (regridding_rho)
1964  call init_coord_rho(cs%rho_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS, &
1965  rho_scale=us%kg_m3_to_R)
1966  case (regridding_hycom1)
1967  call init_coord_hycom(cs%hycom_CS, cs%nk, cs%coordinateResolution, cs%target_density, &
1968  cs%interp_CS, rho_scale=us%kg_m3_to_R)
1969  case (regridding_slight)
1970  call init_coord_slight(cs%slight_CS, cs%nk, cs%ref_pressure, cs%target_density, &
1971  cs%interp_CS, gv%m_to_H, rho_scale=us%kg_m3_to_R)
1972  case (regridding_adaptive)
1973  call init_coord_adapt(cs%adapt_CS, cs%nk, cs%coordinateResolution, gv%m_to_H)
1974  end select
1975 end subroutine initcoord
1976 
1977 !------------------------------------------------------------------------------
1978 !> Set the fixed resolution data
1979 subroutine setcoordinateresolution( dz, CS, scale )
1980  real, dimension(:), intent(in) :: dz !< A vector of vertical grid spacings
1981  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
1982  real, optional, intent(in) :: scale !< A scaling factor converting dz to coordRes
1983 
1984  if (size(dz)/=cs%nk) call mom_error( fatal, &
1985  'setCoordinateResolution: inconsistent number of levels' )
1986 
1987  if (present(scale)) then
1988  cs%coordinateResolution(:) = scale*dz(:)
1989  else
1990  cs%coordinateResolution(:) = dz(:)
1991  endif
1992 
1993 end subroutine setcoordinateresolution
1994 
1995 !> Set target densities based on the old Rlay variable
1996 subroutine set_target_densities_from_gv( GV, US, CS )
1997  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1998  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1999  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2000  ! Local variables
2001  integer :: k, nz
2002 
2003  nz = cs%nk
2004  cs%target_density(1) = (gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)))
2005  cs%target_density(nz+1) = (gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)))
2006  do k = 2,nz
2007  cs%target_density(k) = cs%target_density(k-1) + cs%coordinateResolution(k)
2008  enddo
2009  cs%target_density_set = .true.
2010 
2011 end subroutine set_target_densities_from_gv
2012 
2013 !> Set target densities based on vector of interface values
2014 subroutine set_target_densities( CS, rho_int )
2015  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2016  real, dimension(CS%nk+1), intent(in) :: rho_int !< Interface densities [R ~> kg m-3]
2017 
2018  if (size(cs%target_density)/=size(rho_int)) then
2019  call mom_error(fatal, "set_target_densities inconsistent args!")
2020  endif
2021 
2022  cs%target_density(:) = rho_int(:)
2023  cs%target_density_set = .true.
2024 
2025 end subroutine set_target_densities
2026 
2027 !> Set maximum interface depths based on a vector of input values.
2028 subroutine set_regrid_max_depths( CS, max_depths, units_to_H )
2029  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2030  real, dimension(CS%nk+1), intent(in) :: max_depths !< Maximum interface depths, in arbitrary units
2031  real, optional, intent(in) :: units_to_h !< A conversion factor for max_depths into H units
2032  ! Local variables
2033  real :: val_to_h
2034  integer :: k
2035 
2036  if (.not.allocated(cs%max_interface_depths)) allocate(cs%max_interface_depths(1:cs%nk+1))
2037 
2038  val_to_h = 1.0 ; if (present(units_to_h)) val_to_h = units_to_h
2039  if (max_depths(cs%nk+1) < max_depths(1)) val_to_h = -1.0*val_to_h
2040 
2041  ! Check for sign reversals in the depths.
2042  if (max_depths(cs%nk+1) < max_depths(1)) then
2043  do k=1,cs%nk ; if (max_depths(k+1) > max_depths(k)) &
2044  call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths!")
2045  enddo
2046  else
2047  do k=1,cs%nk ; if (max_depths(k+1) < max_depths(k)) &
2048  call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths.")
2049  enddo
2050  endif
2051 
2052  do k=1,cs%nk+1
2053  cs%max_interface_depths(k) = val_to_h * max_depths(k)
2054  enddo
2055 
2056  ! set max depths for coordinate
2057  select case (cs%regridding_scheme)
2058  case (regridding_hycom1)
2059  call set_hycom_params(cs%hycom_CS, max_interface_depths=cs%max_interface_depths)
2060  case (regridding_slight)
2061  call set_slight_params(cs%slight_CS, max_interface_depths=cs%max_interface_depths)
2062  end select
2063 end subroutine set_regrid_max_depths
2064 
2065 !> Set maximum layer thicknesses based on a vector of input values.
2066 subroutine set_regrid_max_thickness( CS, max_h, units_to_H )
2067  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2068  real, dimension(CS%nk+1), intent(in) :: max_h !< Maximum interface depths, in arbitrary units
2069  real, optional, intent(in) :: units_to_h !< A conversion factor for max_h into H units
2070  ! Local variables
2071  real :: val_to_h
2072  integer :: k
2073 
2074  if (.not.allocated(cs%max_layer_thickness)) allocate(cs%max_layer_thickness(1:cs%nk))
2075 
2076  val_to_h = 1.0 ; if (present( units_to_h)) val_to_h = units_to_h
2077 
2078  do k=1,cs%nk
2079  cs%max_layer_thickness(k) = val_to_h * max_h(k)
2080  enddo
2081 
2082  ! set max thickness for coordinate
2083  select case (cs%regridding_scheme)
2084  case (regridding_hycom1)
2085  call set_hycom_params(cs%hycom_CS, max_layer_thickness=cs%max_layer_thickness)
2086  case (regridding_slight)
2087  call set_slight_params(cs%slight_CS, max_layer_thickness=cs%max_layer_thickness)
2088  end select
2089 end subroutine set_regrid_max_thickness
2090 
2091 
2092 !------------------------------------------------------------------------------
2093 !> Query the fixed resolution data
2094 function getcoordinateresolution( CS, undo_scaling )
2095  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2096  logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2097  !! rescaling of the resolution data.
2098  real, dimension(CS%nk) :: getcoordinateresolution
2099 
2100  logical :: unscale
2101  unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2102 
2103  if (unscale) then
2104  getcoordinateresolution(:) = cs%coord_scale * cs%coordinateResolution(:)
2105  else
2106  getcoordinateresolution(:) = cs%coordinateResolution(:)
2107  endif
2108 
2109 end function getcoordinateresolution
2110 
2111 !> Query the target coordinate interface positions
2112 function getcoordinateinterfaces( CS, undo_scaling )
2113  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2114  logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2115  !! rescaling of the resolution data.
2116  real, dimension(CS%nk+1) :: getcoordinateinterfaces !< Interface positions in target coordinate
2117 
2118  integer :: k
2119  logical :: unscale
2120  unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2121 
2122  ! When using a coordinate with target densities, we need to get the actual
2123  ! densities, rather than computing the interfaces based on resolution
2124  if (cs%regridding_scheme == regridding_rho) then
2125  if (.not. cs%target_density_set) &
2126  call mom_error(fatal, 'MOM_regridding, getCoordinateInterfaces: '//&
2127  'target densities not set!')
2128 
2129  if (unscale) then
2130  getcoordinateinterfaces(:) = cs%coord_scale * cs%target_density(:)
2131  else
2132  getcoordinateinterfaces(:) = cs%target_density(:)
2133  endif
2134  else
2135  if (unscale) then
2136  getcoordinateinterfaces(1) = 0.
2137  do k = 1, cs%nk
2139  cs%coord_scale * cs%coordinateResolution(k)
2140  enddo
2141  else
2142  getcoordinateinterfaces(1) = 0.
2143  do k = 1, cs%nk
2145  cs%coordinateResolution(k)
2146  enddo
2147  endif
2148  ! The following line has an "abs()" to allow ferret users to reference
2149  ! data by index. It is a temporary work around... :( -AJA
2151  endif
2152 
2153 end function getcoordinateinterfaces
2154 
2155 !------------------------------------------------------------------------------
2156 !> Query the target coordinate units
2157 function getcoordinateunits( CS )
2158  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2159  character(len=20) :: getcoordinateunits
2160 
2161  select case ( cs%regridding_scheme )
2162  case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
2163  getcoordinateunits = 'meter'
2164  case ( regridding_sigma_shelf_zstar )
2165  getcoordinateunits = 'meter/fraction'
2166  case ( regridding_sigma )
2167  getcoordinateunits = 'fraction'
2168  case ( regridding_rho )
2169  getcoordinateunits = 'kg/m3'
2170  case ( regridding_arbitrary )
2171  getcoordinateunits = 'unknown'
2172  case default
2173  call mom_error(fatal,'MOM_regridding, getCoordinateUnits: '//&
2174  'Unknown regridding scheme selected!')
2175  end select ! type of grid
2176 
2177 end function getcoordinateunits
2178 
2179 !------------------------------------------------------------------------------
2180 !> Query the short name of the coordinate
2181 function getcoordinateshortname( CS )
2182  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2183  character(len=20) :: getcoordinateshortname
2184 
2185  select case ( cs%regridding_scheme )
2186  case ( regridding_zstar )
2187  !getCoordinateShortName = 'z*'
2188  ! The following line is a temporary work around... :( -AJA
2189  getcoordinateshortname = 'pseudo-depth, -z*'
2190  case ( regridding_sigma_shelf_zstar )
2191  getcoordinateshortname = 'pseudo-depth, -z*/sigma'
2192  case ( regridding_sigma )
2193  getcoordinateshortname = 'sigma'
2194  case ( regridding_rho )
2195  getcoordinateshortname = 'rho'
2196  case ( regridding_arbitrary )
2197  getcoordinateshortname = 'coordinate'
2198  case ( regridding_hycom1 )
2199  getcoordinateshortname = 'z-rho'
2200  case ( regridding_slight )
2201  getcoordinateshortname = 's-rho'
2202  case ( regridding_adaptive )
2203  getcoordinateshortname = 'adaptive'
2204  case default
2205  call mom_error(fatal,'MOM_regridding, getCoordinateShortName: '//&
2206  'Unknown regridding scheme selected!')
2207  end select ! type of grid
2208 
2209 end function getcoordinateshortname
2210 
2211 !> Can be used to set any of the parameters for MOM_regridding.
2212 subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, &
2213  interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, &
2214  compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, &
2215  nlay_ML_to_interior, fix_haloclines, halocline_filt_len, &
2216  halocline_strat_tol, integrate_downward_for_e, &
2217  adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
2218  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2219  logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
2220  real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the
2221  !! new grid [H ~> m or kg m-2]
2222  real, optional, intent(in) :: old_grid_weight !< Weight given to old coordinate when time-filtering grid
2223  character(len=*), optional, intent(in) :: interp_scheme !< Interpolation method for state-dependent coordinates
2224  real, optional, intent(in) :: depth_of_time_filter_shallow !< Depth to start cubic [H ~> m or kg m-2]
2225  real, optional, intent(in) :: depth_of_time_filter_deep !< Depth to end cubic [H ~> m or kg m-2]
2226  real, optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density
2227  real, optional, intent(in) :: dz_min_surface !< The fixed resolution in the topmost
2228  !! SLight_nkml_min layers [H ~> m or kg m-2]
2229  integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickness layers at the top of the model
2230  real, optional, intent(in) :: rho_ml_avg_depth !< Averaging depth over which to determine mixed layer potential
2231  !! density [H ~> m or kg m-2]
2232  real, optional, intent(in) :: nlay_ml_to_interior !< Number of layers to offset the mixed layer density to find
2233  !! resolved stratification [nondim]
2234  logical, optional, intent(in) :: fix_haloclines !< Detect regions with much weaker stratification in the coordinate
2235  real, optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for
2236  !! spuriously unstable water mass profiles [m]
2237  real, optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic
2238  !! halocline region.
2239  logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward
2240  !! from the top.
2241  real, optional, intent(in) :: adapttimeratio !< Ratio of the ALE timestep to the grid timescale [nondim].
2242  real, optional, intent(in) :: adaptzoom !< Depth of near-surface zooming region [H ~> m or kg m-2].
2243  real, optional, intent(in) :: adaptzoomcoeff !< Coefficient of near-surface zooming diffusivity [nondim].
2244  real, optional, intent(in) :: adaptbuoycoeff !< Coefficient of buoyancy diffusivity [nondim].
2245  real, optional, intent(in) :: adaptalpha !< Scaling factor on optimization tendency [nondim].
2246  logical, optional, intent(in) :: adaptdomin !< If true, make a HyCOM-like mixed layer by
2247  !! preventing interfaces from being shallower than
2248  !! the depths specified by the regridding coordinate.
2249 
2250  if (present(interp_scheme)) call set_interp_scheme(cs%interp_CS, interp_scheme)
2251  if (present(boundary_extrapolation)) call set_interp_extrap(cs%interp_CS, boundary_extrapolation)
2252 
2253  if (present(old_grid_weight)) then
2254  if (old_grid_weight<0. .or. old_grid_weight>1.) &
2255  call mom_error(fatal,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
2256  cs%old_grid_weight = old_grid_weight
2257  endif
2258  if (present(depth_of_time_filter_shallow)) cs%depth_of_time_filter_shallow = depth_of_time_filter_shallow
2259  if (present(depth_of_time_filter_deep)) cs%depth_of_time_filter_deep = depth_of_time_filter_deep
2260  if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then
2261  if (cs%depth_of_time_filter_deep<cs%depth_of_time_filter_shallow) call mom_error(fatal,'MOM_regridding, '//&
2262  'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
2263  endif
2264 
2265  if (present(min_thickness)) cs%min_thickness = min_thickness
2266  if (present(compress_fraction)) cs%compressibility_fraction = compress_fraction
2267  if (present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
2268 
2269  select case (cs%regridding_scheme)
2270  case (regridding_zstar)
2271  if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2272  case (regridding_sigma_shelf_zstar)
2273  if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2274  case (regridding_sigma)
2275  if (present(min_thickness)) call set_sigma_params(cs%sigma_CS, min_thickness=min_thickness)
2276  case (regridding_rho)
2277  if (present(min_thickness)) call set_rho_params(cs%rho_CS, min_thickness=min_thickness)
2278  if (present(integrate_downward_for_e)) &
2279  call set_rho_params(cs%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
2280  if (associated(cs%rho_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
2281  call set_rho_params(cs%rho_CS, interp_cs=cs%interp_CS)
2282  case (regridding_hycom1)
2283  if (associated(cs%hycom_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
2284  call set_hycom_params(cs%hycom_CS, interp_cs=cs%interp_CS)
2285  case (regridding_slight)
2286  if (present(min_thickness)) call set_slight_params(cs%slight_CS, min_thickness=min_thickness)
2287  if (present(dz_min_surface)) call set_slight_params(cs%slight_CS, dz_ml_min=dz_min_surface)
2288  if (present(nz_fixed_surface)) call set_slight_params(cs%slight_CS, nz_fixed_surface=nz_fixed_surface)
2289  if (present(rho_ml_avg_depth)) call set_slight_params(cs%slight_CS, rho_ml_avg_depth=rho_ml_avg_depth)
2290  if (present(nlay_ml_to_interior)) call set_slight_params(cs%slight_CS, nlay_ml_offset=nlay_ml_to_interior)
2291  if (present(fix_haloclines)) call set_slight_params(cs%slight_CS, fix_haloclines=fix_haloclines)
2292  if (present(halocline_filt_len)) call set_slight_params(cs%slight_CS, halocline_filter_length=halocline_filt_len)
2293  if (present(halocline_strat_tol)) call set_slight_params(cs%slight_CS, halocline_strat_tol=halocline_strat_tol)
2294  if (present(compress_fraction)) call set_slight_params(cs%slight_CS, compressibility_fraction=compress_fraction)
2295  if (associated(cs%slight_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
2296  call set_slight_params(cs%slight_CS, interp_cs=cs%interp_CS)
2297  case (regridding_adaptive)
2298  if (present(adapttimeratio)) call set_adapt_params(cs%adapt_CS, adapttimeratio=adapttimeratio)
2299  if (present(adaptzoom)) call set_adapt_params(cs%adapt_CS, adaptzoom=adaptzoom)
2300  if (present(adaptzoomcoeff)) call set_adapt_params(cs%adapt_CS, adaptzoomcoeff=adaptzoomcoeff)
2301  if (present(adaptbuoycoeff)) call set_adapt_params(cs%adapt_CS, adaptbuoycoeff=adaptbuoycoeff)
2302  if (present(adaptalpha)) call set_adapt_params(cs%adapt_CS, adaptalpha=adaptalpha)
2303  if (present(adaptdomin)) call set_adapt_params(cs%adapt_CS, adaptdomin=adaptdomin)
2304  end select
2305 
2306 end subroutine set_regrid_params
2307 
2308 !> Returns the number of levels/layers in the regridding control structure
2309 integer function get_regrid_size(CS)
2310  type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2311 
2312  get_regrid_size = cs%nk
2313 
2314 end function get_regrid_size
2315 
2316 !> This returns a copy of the zlike_CS stored in the regridding control structure.
2317 function get_zlike_cs(CS)
2318  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2319  type(zlike_cs) :: get_zlike_cs
2320 
2321  get_zlike_cs = cs%zlike_CS
2322 end function get_zlike_cs
2323 
2324 !> This returns a copy of the sigma_CS stored in the regridding control structure.
2325 function get_sigma_cs(CS)
2326  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2327  type(sigma_cs) :: get_sigma_cs
2328 
2329  get_sigma_cs = cs%sigma_CS
2330 end function get_sigma_cs
2331 
2332 !> This returns a copy of the rho_CS stored in the regridding control structure.
2333 function get_rho_cs(CS)
2334  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2335  type(rho_cs) :: get_rho_cs
2336 
2337  get_rho_cs = cs%rho_CS
2338 end function get_rho_cs
2339 
2340 !------------------------------------------------------------------------------
2341 !> Return coordinate-derived thicknesses for fixed coordinate systems
2342 function getstaticthickness( CS, SSH, depth )
2343  type(regridding_cs), intent(in) :: cs !< Regridding control structure
2344  real, intent(in) :: ssh !< The sea surface height, in the same units as depth
2345  real, intent(in) :: depth !< The maximum depth of the grid, often [Z ~> m]
2346  real, dimension(CS%nk) :: getstaticthickness !< The returned thicknesses in the units of depth
2347  ! Local
2348  integer :: k
2349  real :: z, dz
2350 
2351  select case ( cs%regridding_scheme )
2352  case ( regridding_zstar, regridding_sigma_shelf_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
2353  if (depth>0.) then
2354  z = ssh
2355  do k = 1, cs%nk
2356  dz = cs%coordinateResolution(k) * ( 1. + ssh/depth ) ! Nominal dz*
2357  dz = max(dz, 0.) ! Avoid negative incase ssh=-depth
2358  dz = min(dz, depth - z) ! Clip if below topography
2359  z = z + dz ! Bottom of layer
2360  getstaticthickness(k) = dz
2361  enddo
2362  else
2363  getstaticthickness(:) = 0. ! On land ...
2364  endif
2365  case ( regridding_sigma )
2366  getstaticthickness(:) = cs%coordinateResolution(:) * ( depth + ssh )
2367  case ( regridding_rho )
2368  getstaticthickness(:) = 0. ! Not applicable
2369  case ( regridding_arbitrary )
2370  getstaticthickness(:) = 0. ! Not applicable
2371  case default
2372  call mom_error(fatal,'MOM_regridding, getStaticThickness: '//&
2373  'Unknown regridding scheme selected!')
2374  end select ! type of grid
2375 
2376 end function getstaticthickness
2377 
2378 !> Parses a string and generates a dz(:) profile that goes like k**power.
2379 subroutine dz_function1( string, dz )
2380  character(len=*), intent(in) :: string !< String with list of parameters in form
2381  !! dz_min, H_total, power, precision
2382  real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses
2383  ! Local variables
2384  integer :: nk, k
2385  real :: dz_min, power, prec, H_total
2386 
2387  nk = size(dz) ! Number of cells
2388  prec = -1024.
2389  read( string, *) dz_min, h_total, power, prec
2390  if (prec == -1024.) call mom_error(fatal,"dz_function1: "// &
2391  "Problem reading FNC1: string ="//trim(string))
2392  ! Create profile of ( dz - dz_min )
2393  do k = 1, nk
2394  dz(k) = (real(k-1)/real(nk-1))**power
2395  enddo
2396  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2397  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2398  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2399  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2400  dz(nk) = dz(nk) + ( h_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer
2401  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2402  dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min
2403 
2404 end subroutine dz_function1
2405 
2406 !> Parses a string and generates a rho_target(:) profile with refined resolution downward
2407 !! and returns the number of levels
2408 integer function rho_function1( string, rho_target )
2409  character(len=*), intent(in) :: string !< String with list of parameters in form
2410  !! dz_min, H_total, power, precision
2411  real, dimension(:), allocatable, intent(inout) :: rho_target !< Profile of interface densities [kg m-3]
2412  ! Local variables
2413  integer :: nki, k, nk
2414  real :: ddx, dx, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2415 
2416  read( string, *) nk, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2417  allocate(rho_target(nk+1))
2418  nki = nk + 1 - 4 ! Number of interfaces minus 4 specified values
2419  rho_target(1) = rho_1
2420  rho_target(2) = rho_2
2421  dx = 0.
2422  do k = 0, nki
2423  ddx = max( drho_min, real(nki-k)/real(nki*nki) )
2424  dx = dx + ddx
2425  rho_target(3+k) = rho_3 + (2. * drho) * dx
2426  enddo
2427  rho_target(nki+4) = rho_4
2428 
2429  rho_function1 = nk
2430 
2431 end function rho_function1
2432 
2433 !> \namespace mom_regridding
2434 !!
2435 !! A vertical grid is defined solely by the cell thicknesses, \f$h\f$.
2436 !! Most calculations in this module start with the coordinate at the bottom
2437 !! of the column set to -depth, and use a increasing value of coordinate with
2438 !! decreasing k. This is consistent with the rest of MOM6 that uses position,
2439 !! \f$z\f$ which is a negative quantity for most of the ocean.
2440 !!
2441 !! A change in grid is define through a change in position of the interfaces:
2442 !! \f[
2443 !! z^n_{k+1/2} = z^{n-1}_{k+1/2} + \Delta z_{k+1/2}
2444 !! \f]
2445 !! with the positive upward coordinate convention
2446 !! \f[
2447 !! z_{k-1/2} = z_{k+1/2} + h_k
2448 !! \f]
2449 !! so that
2450 !! \f[
2451 !! h^n_k = h^{n-1}_k + ( \Delta z_{k-1/2} - \Delta z_{k+1/2} )
2452 !! \f]
2453 !!
2454 !! Original date of creation: 2008.06.09 by L. White
2455 
2456 end module mom_regridding
mom_regridding::check_grid_column
subroutine, public check_grid_column(nk, depth, h, dzInterface, msg)
Check that the total thickness of new and old grids are consistent.
Definition: MOM_regridding.F90:938
coord_slight::slight_cs
Control structure containing required parameters for the SLight coordinate.
Definition: coord_slight.F90:15
mom_regridding::setcoordinateresolution
subroutine, public setcoordinateresolution(dz, CS, scale)
Set the fixed resolution data.
Definition: MOM_regridding.F90:1980
mom_regridding::get_sigma_cs
type(sigma_cs) function, public get_sigma_cs(CS)
This returns a copy of the sigma_CS stored in the regridding control structure.
Definition: MOM_regridding.F90:2326
coord_rho::set_rho_params
subroutine, public set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS)
This subroutine can be used to set the parameters for the coord_rho module.
Definition: coord_rho.F90:81
coord_rho
Regrid columns for the continuous isopycnal (rho) coordinate.
Definition: coord_rho.F90:2
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_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
coord_slight::end_coord_slight
subroutine, public end_coord_slight(CS)
This subroutine deallocates memory in the control structure for the coord_slight module.
Definition: coord_slight.F90:110
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_interp::set_interp_scheme
subroutine, public set_interp_scheme(CS, interp_scheme)
Store the interpolation_scheme value in the interp_CS based on the input string.
Definition: regrid_interp.F90:518
regrid_consts::state_dependent
Returns true if the coordinate is dependent on the state density, returns false otherwise.
Definition: regrid_consts.F90:44
coord_hycom::build_hycom1_column
subroutine, public build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, z_col, z_col_new, zScale, h_neglect, h_neglect_edge)
Build a HyCOM coordinate column.
Definition: coord_hycom.F90:105
mom_regridding::regriddinginterpschemedoc
character(len= *), parameter, public regriddinginterpschemedoc
Documentation for regridding interpolation schemes.
Definition: MOM_regridding.F90:152
coord_sigma::init_coord_sigma
subroutine, public init_coord_sigma(CS, nk, coordinateResolution)
Initialise a sigma_CS with pointers to parameters.
Definition: coord_sigma.F90:29
mom_regridding::build_grid_adaptive
subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS)
This subroutine builds an adaptive grid that follows density surfaces where possible,...
Definition: MOM_regridding.F90:1510
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
regrid_consts::regridding_arbitrary
integer, parameter regridding_arbitrary
Arbitrary coordinates identifier.
Definition: regrid_consts.F90:17
mom_regridding::rho_function1
integer function rho_function1(string, rho_target)
Parses a string and generates a rho_target(:) profile with refined resolution downward and returns th...
Definition: MOM_regridding.F90:2409
coord_zlike::zlike_cs
Control structure containing required parameters for a z-like coordinate.
Definition: coord_zlike.F90:11
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
coord_adapt::set_adapt_params
subroutine, public set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptDrho0, adaptDoMin)
This subtroutine can be used to set the parameters for coord_adapt module.
Definition: coord_adapt.F90:93
regrid_consts::coordinateunits
Returns a string with the coordinate units associated with the coordinate mode.
Definition: regrid_consts.F90:38
coord_zlike::build_zstar_column
subroutine, public build_zstar_column(CS, depth, total_thickness, zInterface, z_rigid_top, eta_orig, zScale)
Builds a z* coordinate with a minimum thickness.
Definition: coord_zlike.F90:65
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_string_functions::extract_real
real function, public extract_real(string, separators, n, missing_value)
Returns the real corresponding to the nth word in the argument.
Definition: MOM_string_functions.F90:268
coord_rho::rho_cs
Control structure containing required parameters for the rho coordinate.
Definition: coord_rho.F90:14
coord_adapt::init_coord_adapt
subroutine, public init_coord_adapt(CS, nk, coordinateResolution, m_to_H)
Initialise an adapt_CS with parameters.
Definition: coord_adapt.F90:53
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_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_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
mom_regridding::build_grid_arbitrary
subroutine build_grid_arbitrary(G, GV, h, dzInterface, h_new, CS)
Definition: MOM_regridding.F90:1701
mom_regridding::convective_adjustment
subroutine convective_adjustment(G, GV, h, tv)
Achieve convective adjustment by swapping layers.
Definition: MOM_regridding.F90:1844
regrid_consts::regridding_layer
integer, parameter regridding_layer
Layer mode identifier.
Definition: regrid_consts.F90:13
mom_regridding::initcoord
subroutine initcoord(CS, GV, US, coord_mode)
Initialize the coordinate resolutions by calling the appropriate initialization routine for the speci...
Definition: MOM_regridding.F90:1949
mom_regridding::dz_function1
subroutine dz_function1(string, dz)
Parses a string and generates a dz(:) profile that goes like k**power.
Definition: MOM_regridding.F90:2380
regrid_interp::set_interp_extrap
subroutine, public set_interp_extrap(CS, extrap)
Store the boundary_extrapolation value in the interp_CS.
Definition: regrid_interp.F90:528
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_regridding::check_remapping_grid
subroutine, public check_remapping_grid(G, GV, h, dzInterface, msg)
Check that the total thickness of two grids match.
Definition: MOM_regridding.F90:920
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
coord_rho::init_coord_rho
subroutine, public init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS, rho_scale)
Initialise a rho_CS with pointers to parameters.
Definition: coord_rho.F90:50
mom_regridding::set_regrid_max_thickness
subroutine, public set_regrid_max_thickness(CS, max_h, units_to_H)
Set maximum layer thicknesses based on a vector of input values.
Definition: MOM_regridding.F90:2067
regrid_consts::regridding_sigma_shelf_zstar
integer, parameter regridding_sigma_shelf_zstar
Identifiered for z* coordinates at the bottom, sigma-near the top.
Definition: regrid_consts.F90:21
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
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
coord_hycom
Regrid columns for the HyCOM coordinate.
Definition: coord_hycom.F90:2
mom_regridding::filtered_grid_motion
subroutine filtered_grid_motion(CS, nk, z_old, z_new, dz_g)
Returns the change in interface position motion after filtering and assuming the top and bottom inter...
Definition: MOM_regridding.F90:999
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
regrid_consts::regridding_hycom1
integer, parameter regridding_hycom1
Simple HyCOM coordinates without BBL.
Definition: regrid_consts.F90:18
regrid_interp
Vertical interpolation for regridding.
Definition: regrid_interp.F90:2
mom_regridding::getcoordinateresolution
real function, dimension(cs%nk), public getcoordinateresolution(CS, undo_scaling)
Query the fixed resolution data.
Definition: MOM_regridding.F90:2095
mom_regridding::build_grid_slight
subroutine build_grid_slight(G, GV, h, tv, dzInterface, CS)
Builds a grid that tracks density interfaces for water that is denser than the surface density plus a...
Definition: MOM_regridding.F90:1576
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
regrid_consts::regridding_slight
integer, parameter regridding_slight
Identifier for stretched coordinates in the lightest water, isopycnal below.
Definition: regrid_consts.F90:19
mom_regridding::regriddingdefaultminthickness
real, parameter, public regriddingdefaultminthickness
Default minimum thickness for some coordinate generation modes.
Definition: MOM_regridding.F90:169
mom_regridding::build_zstar_grid
subroutine build_zstar_grid(CS, G, GV, h, dzInterface, frac_shelf_h)
Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004). z* is defined as z* = (z-e...
Definition: MOM_regridding.F90:1137
coord_zlike::end_coord_zlike
subroutine, public end_coord_zlike(CS)
Deallocates the zlike control structure.
Definition: coord_zlike.F90:44
mom_regridding::get_regrid_size
integer function, public get_regrid_size(CS)
Returns the number of levels/layers in the regridding control structure.
Definition: MOM_regridding.F90:2310
mom_regridding::inflate_vanished_layers_old
subroutine, public inflate_vanished_layers_old(CS, G, GV, h)
Definition: MOM_regridding.F90:1803
coord_adapt::end_coord_adapt
subroutine, public end_coord_adapt(CS)
Clean up the coordinate control structure.
Definition: coord_adapt.F90:82
regrid_consts::regridding_sigma
integer, parameter regridding_sigma
Sigma coordinates identifier.
Definition: regrid_consts.F90:16
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_regridding::build_rho_grid
subroutine build_rho_grid(G, GV, h, tv, dzInterface, remapCS, CS)
This routine builds a new grid based on a given set of target interface densities.
Definition: MOM_regridding.F90:1311
coord_zlike
Regrid columns for a z-like coordinate (z-star, z-level)
Definition: coord_zlike.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
coord_slight::set_slight_params
subroutine, public set_slight_params(CS, max_interface_depths, max_layer_thickness, min_thickness, compressibility_fraction, dz_ml_min, nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_offset, fix_haloclines, halocline_filter_length, halocline_strat_tol, interp_CS)
This subroutine can be used to set the parameters for the coord_slight module.
Definition: coord_slight.F90:123
coord_sigma::build_sigma_column
subroutine, public build_sigma_column(CS, depth, totalThickness, zInterface)
Build a sigma coordinate column.
Definition: coord_sigma.F90:64
mom_regridding::check_grid_def
subroutine check_grid_def(filename, varname, expected_units, msg, ierr)
Do some basic checks on the vertical grid definition file, variable.
Definition: MOM_regridding.F90:716
coord_slight::build_slight_column
subroutine, public build_slight_column(CS, eqn_of_state, H_to_Pa, H_subroundoff, nz, depth, h_col, T_col, S_col, p_col, z_col, z_col_new, h_neglect, h_neglect_edge)
Build a SLight coordinate column.
Definition: coord_slight.F90:188
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
coord_adapt::adapt_cs
Control structure for adaptive coordinates (coord_adapt).
Definition: coord_adapt.F90:16
coord_hycom::set_hycom_params
subroutine, public set_hycom_params(CS, max_interface_depths, max_layer_thickness, interp_CS)
This subroutine can be used to set the parameters for the coord_hycom module.
Definition: coord_hycom.F90:78
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
regrid_consts::regridding_zstar
integer, parameter regridding_zstar
z* coordinates identifier
Definition: regrid_consts.F90:14
mom_regridding::set_regrid_max_depths
subroutine, public set_regrid_max_depths(CS, max_depths, units_to_H)
Set maximum interface depths based on a vector of input values.
Definition: MOM_regridding.F90:2029
mom_regridding::end_regridding
subroutine, public end_regridding(CS)
Deallocation of regridding memory.
Definition: MOM_regridding.F90:772
mom_regridding::regriddingcoordinatemodedoc
character(len= *), parameter, public regriddingcoordinatemodedoc
Documentation for coordinate options.
Definition: MOM_regridding.F90:141
coord_sigma::sigma_cs
Control structure containing required parameters for the sigma coordinate.
Definition: coord_sigma.F90:11
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_regridding::build_grid_hycom1
subroutine build_grid_hycom1(G, GV, h, tv, h_new, dzInterface, CS)
Builds a simple HyCOM-like grid with the deepest location of potential density interpolated from the ...
Definition: MOM_regridding.F90:1441
mom_regridding::build_sigma_grid
subroutine build_sigma_grid(CS, G, GV, h, dzInterface)
This routine builds a grid based on terrain-following coordinates.
Definition: MOM_regridding.F90:1232
coord_hycom::end_coord_hycom
subroutine, public end_coord_hycom(CS)
This subroutine deallocates memory in the control structure for the coord_hycom module.
Definition: coord_hycom.F90:65
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
coord_zlike::set_zlike_params
subroutine, public set_zlike_params(CS, min_thickness)
Set parameters in the zlike structure.
Definition: coord_zlike.F90:54
mom_regridding::calc_h_new_by_dz
subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
Calculates h_new from h + delta_k dzInterface.
Definition: MOM_regridding.F90:885
coord_slight::init_coord_slight
subroutine, public init_coord_slight(CS, nk, ref_pressure, target_density, interp_CS, m_to_H, rho_scale)
Initialise a slight_CS with pointers to parameters.
Definition: coord_slight.F90:76
mom_regridding::get_rho_cs
type(rho_cs) function, public get_rho_cs(CS)
This returns a copy of the rho_CS stored in the regridding control structure.
Definition: MOM_regridding.F90:2334
coord_adapt::build_adapt_column
subroutine, public build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
Definition: coord_adapt.F90:118
regrid_interp::interp_cs_type
Control structure for regrid_interp module.
Definition: regrid_interp.F90:23
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
regrid_consts::default_coordinate_mode
character(len= *), parameter default_coordinate_mode
Default coordinate mode.
Definition: regrid_consts.F90:35
mom_io::file_exists
Indicate whether a file exists, perhaps with domain decomposition.
Definition: MOM_io.F90:68
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
coord_rho::build_rho_column
subroutine, public build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, h_neglect, h_neglect_edge)
Build a rho coordinate column.
Definition: coord_rho.F90:102
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
coord_rho::end_coord_rho
subroutine, public end_coord_rho(CS)
This subroutine deallocates memory in the control structure for the coord_rho module.
Definition: coord_rho.F90:71
coord_sigma::set_sigma_params
subroutine, public set_sigma_params(CS, min_thickness)
This subroutine can be used to set the parameters for the coord_sigma module.
Definition: coord_sigma.F90:53
regrid_consts::regridding_rho
integer, parameter regridding_rho
Density coordinates identifier.
Definition: regrid_consts.F90:15
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_regridding::adjust_interface_motion
subroutine adjust_interface_motion(CS, nk, h_old, dz_int)
Adjust dz_Interface to ensure non-negative future thicknesses.
Definition: MOM_regridding.F90:1640
coord_zlike::init_coord_zlike
subroutine, public init_coord_zlike(CS, nk, coordinateResolution)
Initialise a zlike_CS with pointers to parameters.
Definition: coord_zlike.F90:30
mom_regridding::get_zlike_cs
type(zlike_cs) function, public get_zlike_cs(CS)
This returns a copy of the zlike_CS stored in the regridding control structure.
Definition: MOM_regridding.F90:2318
coord_sigma::end_coord_sigma
subroutine, public end_coord_sigma(CS)
This subroutine deallocates memory in the control structure for the coord_sigma module.
Definition: coord_sigma.F90:43
coord_sigma
Regrid columns for the sigma coordinate.
Definition: coord_sigma.F90:2
regrid_consts::regridding_adaptive
integer, parameter regridding_adaptive
Adaptive coordinate mode identifier.
Definition: regrid_consts.F90:23
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
coord_hycom::hycom_cs
Control structure containing required parameters for the HyCOM coordinate.
Definition: coord_hycom.F90:13
coord_slight
Regrid columns for the SLight coordinate.
Definition: coord_slight.F90:2
coord_adapt
Regrid columns for the adaptive coordinate.
Definition: coord_adapt.F90:2
coord_hycom::init_coord_hycom
subroutine, public init_coord_hycom(CS, nk, coordinateResolution, target_density, interp_CS, rho_scale)
Initialise a hycom_CS with pointers to parameters.
Definition: coord_hycom.F90:43
coord_rho::old_inflate_layers_1d
subroutine, public old_inflate_layers_1d(min_thickness, nk, h)
Inflate vanished layers to finite (nonzero) width.
Definition: coord_rho.F90:357
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_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60