MOM_regridding.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Generates vertical grids as part of the ALE algorithm
7
8use mom_error_handler, only : mom_error, fatal, warning, note, assert
9use mom_file_parser, only : param_file_type, get_param, log_param
10use mom_io, only : file_exists, field_exists, field_size, mom_read_data
11use mom_io, only : read_variable
12use mom_io, only : vardesc, var_desc, single_file
13use mom_io, only : mom_netcdf_file, mom_field
14use mom_io, only : create_mom_file, mom_write_field
15use mom_io, only : verify_variable_units, slasher
17use mom_variables, only : ocean_grid_type, thermo_var_ptrs
19use mom_eos, only : eos_type, calculate_density
20use mom_domains, only : max_across_pes, pass_var
22
24use regrid_consts, only : state_dependent, coordinateunits
25use regrid_consts, only : coordinatemode, default_coordinate_mode
26use regrid_consts, only : regridding_layer, regridding_zstar
27use regrid_consts, only : regridding_rho, regridding_sigma
28use regrid_consts, only : regridding_arbitrary, regridding_sigma_shelf_zstar
29use regrid_consts, only : regridding_hycom1, regridding_hybgen, regridding_adaptive
32
33use coord_zlike, only : zlike_cs
35use coord_sigma, only : sigma_cs
39use coord_hycom, only : hycom_cs
42use coord_adapt, only : adapt_cs
46
47implicit none ; private
48
49#include <MOM_memory.h>
50
51! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
52! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
53! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
54! vary with the Boussinesq approximation, the Boussinesq variant is given first.
55
56!> Regridding control structure
57type, public :: regridding_cs ; private
58
59 !> This array is set by function setCoordinateResolution()
60 !! It contains the "resolution" or delta coordinate of the target
61 !! coordinate. It has the units of the target coordinate, e.g.
62 !! [Z ~> m] for z*, [nondim] for sigma, etc.
63 real, dimension(:), allocatable :: coordinateresolution
64
65 !> This is a scaling factor that restores coordinateResolution to values in
66 !! the natural units for output, perhaps [nondim]
67 real :: coord_scale = 1.0
68
69 !> This array is set by function set_target_densities()
70 !! This array is the nominal coordinate of interfaces and is the
71 !! running sum of coordinateResolution, in [R ~> kg m-3]. i.e.
72 !! target_density(k+1) = coordinateResolution(k) + coordinateResolution(k)
73 !! It is only used in "rho" or "Hycom" mode.
74 real, dimension(:), allocatable :: target_density
75
76 !> A flag to indicate that the target_density arrays has been filled with data.
77 logical :: target_density_set = .false.
78
79 !> Nominal HYCOM1 3D near-surface resolution [Z ~> m]
80 real, allocatable, dimension(:,:,:) :: coordinateresolution_3d
81
82 !> Nominal HYCOM1 3D density of interfaces [R ~> kg m-3]
83 real, allocatable, dimension(:,:,:) :: target_density_3d
84
85 !> This array is set by function set_regrid_max_depths()
86 !! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
87 real, dimension(:), allocatable :: max_interface_depths
88
89 !> This array is set by function set_regrid_max_thickness()
90 !! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
91 real, dimension(:), allocatable :: max_layer_thickness
92
93 integer :: nk !< Number of layers/levels in generated grid
94
95 !> Indicates which grid to use in the vertical (z*, sigma, target interface
96 !! densities)
97 integer :: regridding_scheme
98
99 !> Interpolation control structure
100 type(interp_cs_type) :: interp_cs
101
102 !> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2].
103 real :: min_thickness
104
105 !> If true, call adjust_interface_motion() after initial grid generation
106 logical :: use_adjust_interface_motion
107
108 !> Reference pressure for potential density calculations [R L2 T-2 ~> Pa]
109 real :: ref_pressure = 2.e7
110
111 !> If true, always pass through the depth-based time filtering that uses CS%old_grid_weight
112 !! If false, allows bypassing of the call if CS%old_grid_weight==0
113 logical :: use_depth_based_time_filter
114
115 !> Weight given to old coordinate when blending between new and old grids [nondim]
116 !! Used only below depth_of_time_filter_shallow, with a cubic variation
117 !! from zero to full effect between depth_of_time_filter_shallow and
118 !! depth_of_time_filter_deep.
119 real :: old_grid_weight = 0.
120
121 !> Depth above which no time-filtering of grid is applied [H ~> m or kg m-2]
122 real :: depth_of_time_filter_shallow = 0.
123
124 !> Depth below which time-filtering of grid is applied at full effect [H ~> m or kg m-2]
125 real :: depth_of_time_filter_deep = 0.
126
127 !> Fraction (between 0 and 1) of compressibility to add to potential density
128 !! profiles when interpolating for target grid positions [nondim]
129 real :: compressibility_fraction = 0.
130
131 !> If true, each interface is given a maximum depth based on a rescaling of
132 !! the indexing of coordinateResolution.
133 logical :: set_maximum_depths = .false.
134
135 !> If true, integrate for interface positions from the top downward.
136 !! If false, integrate from the bottom upward, as does the rest of the model.
137 logical :: integrate_downward_for_e = .true.
138
139 !> The vintage of the order of arithmetic and expressions to use for remapping.
140 !! Values below 20190101 recover the remapping answers from 2018.
141 !! Higher values use more robust forms of the same remapping expressions.
142 integer :: remap_answer_date = 99991231
143
144 logical :: use_hybgen_unmix = .false. !< If true, use the hybgen unmixing code before remapping
145
146 type(zlike_cs), pointer :: zlike_cs => null() !< Control structure for z-like coordinate generator
147 type(sigma_cs), pointer :: sigma_cs => null() !< Control structure for sigma coordinate generator
148 type(rho_cs), pointer :: rho_cs => null() !< Control structure for rho coordinate generator
149 type(hycom_cs), pointer :: hycom_cs => null() !< Control structure for hybrid coordinate generator
150 type(adapt_cs), pointer :: adapt_cs => null() !< Control structure for adaptive coordinate generator
151 type(hybgen_regrid_cs), pointer :: hybgen_cs => null() !< Control structure for hybgen regridding
152
153end type
154
155! The following routines are visible to the outside world
165public default_coordinate_mode
168
169!> Documentation for coordinate options
170character(len=*), parameter, public :: regriddingcoordinatemodedoc = &
171 " LAYER - Isopycnal or stacked shallow water layers\n"//&
172 " ZSTAR, Z* - stretched geopotential z*\n"//&
173 " SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"//&
174 " SIGMA - terrain following coordinates\n"//&
175 " RHO - continuous isopycnal\n"//&
176 " HYCOM1 - HyCOM-like hybrid coordinate\n"//&
177 " HYBGEN - Hybrid coordinate from the Hycom hybgen code\n"//&
178 " ADAPTIVE - optimize for smooth neutral density surfaces"
179
180!> Documentation for regridding interpolation schemes
181character(len=*), parameter, public :: regriddinginterpschemedoc = &
182 " P1M_H2 (2nd-order accurate)\n"//&
183 " P1M_H4 (2nd-order accurate)\n"//&
184 " P1M_IH4 (2nd-order accurate)\n"//&
185 " PLM (2nd-order accurate)\n"//&
186 " PPM_CW (3rd-order accurate)\n"//&
187 " PPM_H4 (3rd-order accurate)\n"//&
188 " PPM_IH4 (3rd-order accurate)\n"//&
189 " P3M_IH4IH3 (4th-order accurate)\n"//&
190 " P3M_IH6IH5 (4th-order accurate)\n"//&
191 " PQM_IH4IH3 (4th-order accurate)\n"//&
192 " PQM_IH6IH5 (5th-order accurate)"
193
194!> Default interpolation scheme
195character(len=*), parameter, public :: regriddingdefaultinterpscheme = "P1M_H2"
196!> Default mode for boundary extrapolation
197logical, parameter, public :: regriddingdefaultboundaryextrapolation = .false.
198!> Default minimum thickness for some coordinate generation modes [m]
199real, parameter, public :: regriddingdefaultminthickness = 1.e-3
200
201!> Maximum length of parameters
202integer, parameter :: max_param_length = 120
203
204#undef __DO_SAFETY_CHECKS__
205
206contains
207
208!> Initialization and configures a regridding control structure based on customizable run-time parameters
209subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, &
210 coord_mode, param_prefix, param_suffix)
211 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
212 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
213 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
214 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
215 real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
216 type(param_file_type), intent(in) :: param_file !< Parameter file
217 character(len=*), intent(in) :: mdl !< Name of calling module.
218 character(len=*), intent(in) :: coord_mode !< Coordinate mode
219 character(len=*), intent(in) :: param_prefix !< String to prefix to parameter names.
220 !! If empty, causes main model parameters to be used.
221 character(len=*), intent(in) :: param_suffix !< String to append to parameter names.
222
223 ! Local variables
224 integer :: ke ! Number of levels
225 integer :: n_sigma ! Number of shallow dz's, for HYBRID_MAP or HYBRID_3D
226 integer :: np ! Number of profiles, for HYBRID_MAP
227 integer :: nceiling ! ceiling of map index, for HYBRID_MAP
228 integer :: nfloor ! floor of map index, for HYBRID_MAP
229 real :: nfrac ! fraction of map index, for HYBRID_MAP [nondim]
230 character(len=80) :: string, string2, varname ! Temporary strings
231 character(len=40) :: coord_units, coord_res_param ! Temporary strings
232 character(len=MAX_PARAM_LENGTH) :: param_name
233 character(len=200) :: inputdir, filename, longstring
234 character(len=320) :: message ! Temporary strings
235 character(len=12) :: expected_units, alt_units ! Temporary strings
236 logical :: tmplogical, do_sum, main_parameters
237 logical :: coord_is_state_dependent, ierr
238 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
239 integer :: remap_answer_date ! The vintage of the remapping expressions to use.
240 integer :: regrid_answer_date ! The vintage of the regridding expressions to use.
241 real :: tmpreal ! A temporary variable used in setting other variables [various]
242 real :: p_ref ! The coordinate variable reference pression [R L2 T-2 ~> Pa]
243 real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
244 real :: dz_extra ! The thickness of an added layer to append to the woa09_dz profile when
245 ! maximum_depth is large [m] (not in Z).
246 real :: nominaldepth ! Depth of ocean bottom in thickness units (positive downward) [H ~> m or kg m-2]
247 real :: depth_q ! A depth scale factor [nondim]
248 real :: depth_s ! The end of the shallow Z regime [m]
249 real :: depth_d ! The start of the deep Z regime [m]
250 real :: adapttimeratio, adaptzoomcoeff ! Temporary variables for input parameters [nondim]
251 real :: adaptbuoycoeff, adaptalpha ! Temporary variables for input parameters [nondim]
252 real :: adaptzoom ! The thickness of the near-surface zooming region with the adaptive coordinate [H ~> m or kg m-2]
253 real :: adaptdrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3]
254 integer :: i, j, k, nzf(4)
255 real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m]
256 ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
257 real, dimension(:,:), allocatable :: dz_2d ! 2D resolution (thickness) in units of coordinate, which may be [m]
258 ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
259 real, dimension(:,:,:), allocatable :: dz_3d ! 3D resolution (thickness) in units of coordinate, which may be [m]
260 ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
261 real, dimension(:), allocatable :: dz_shallow ! Shallow resolution (thickness), for HYBRID_MAP or HYBRID_3D [m]
262 real, dimension(:,:), allocatable :: rho_target_2d ! 2D target density used in HYBRID mode [kg m-3]
263 real, dimension(:,:,:), allocatable :: rho_target_3d ! 3D target density used in HYBRID mode [kg m-3]
264 real, dimension(:,:), allocatable :: index_map ! Region array of indexes for HYBRID_MAP [nondim]
265 real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2]
266 real, dimension(:), allocatable :: z_max ! Maximum interface depths [H ~> m or kg m-2] or other
267 ! units depending on the coordinate
268 real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths
269 ! [H ~> m or kg m-2] or other units
270 real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3]
271 ! Thicknesses [m] that give level centers approximately corresponding to table 2 of WOA09
272 ! These are approximate because the WOA09 depths are not smoothly spaced. Levels
273 ! 1, 4, 5, 9, 12, 24, and 36 are 2.5, 2.5, 1.25 12.5, 37.5 and 62.5 m deeper than WOA09
274 ! but all others are identical.
275 real, dimension(40) :: woa09_dz_approx = (/ 5., 10., 10., 15., 22.5, 25., 25., 25., &
276 37.5, 50., 50., 75., 100., 100., 100., 100., &
277 100., 100., 100., 100., 100., 100., 100., 175., &
278 250., 375., 500., 500., 500., 500., 500., 500., &
279 500., 500., 500., 500., 500., 500., 500., 500. /)
280 ! These are the actual spacings [m] between WOA09 depths which, if used for layer thickness, places
281 ! the interfaces at the WOA09 depths.
282 real, dimension(39) :: woa09_dzi = (/ 10., 10., 10., 20., 25., 25., 25., 25., &
283 50., 50., 50., 100., 100., 100., 100., 100., &
284 100., 100., 100., 100., 100., 100., 100., 250., &
285 250., 500., 500., 500., 500., 500., 500., 500., &
286 500., 500., 500., 500., 500., 500., 500. /)
287 ! These are the spacings [m] between WOA23 depths from table 3 of
288 ! https://www.ncei.noaa.gov/data/oceans/woa/WOA13/DOC/woa13documentation.pdf
289 real, dimension(136) :: woa23_dzi = (/ 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., &
290 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., &
291 25., 25., 25., 25., 25., 25., 25., 25., 25., 25., &
292 25., 25., 25., 25., 25., 25., 50., 50., 50., 50., &
293 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., &
294 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., &
295 50., 50., 50., 50., 50., 50., 100., 100., 100., 100., &
296 100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
297 100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
298 100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
299 100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
300 100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
301 100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
302 100., 100., 100., 100., 100., 100. /)
303
304 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
305 inputdir = slasher(inputdir)
306
307 main_parameters=.false.
308 if (len_trim(param_prefix)==0) main_parameters=.true.
309 if (main_parameters .and. len_trim(param_suffix)>0) call mom_error(fatal,trim(mdl)//&
310 ' initialize_regridding: Suffix provided without prefix for parameter names!')
311
312 cs%nk = 0
313 cs%regridding_scheme = coordinatemode(coord_mode)
314 coord_is_state_dependent = state_dependent(coord_mode)
315 maximum_depth = us%Z_to_m*max_depth
316
317 if (main_parameters) then
318 ! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS)
319 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_UNITS", coord_units, &
320 "Units of the regridding coordinate.", default=coordinateunits(coord_mode))
321 else
322 coord_units=coordinateunits(coord_mode)
323 endif
324
325 if (coord_is_state_dependent) then
326 if (main_parameters) then
327 param_name = "INTERPOLATION_SCHEME"
328 string2 = regriddingdefaultinterpscheme
329 else
330 param_name = create_coord_param(param_prefix, "INTERP_SCHEME", param_suffix)
331 string2 = 'PPM_H4' ! Default for diagnostics
332 endif
333 call get_param(param_file, mdl, param_name, string, &
334 "This sets the interpolation scheme to use to "//&
335 "determine the new grid. These parameters are "//&
336 "only relevant when REGRIDDING_COORDINATE_MODE is "//&
337 "set to a function of state. Otherwise, it is not "//&
338 "used. It can be one of the following schemes: \n"//&
339 trim(regriddinginterpschemedoc), default=trim(string2))
340 call set_regrid_params(cs, interp_scheme=string)
341
342 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
343 "This sets the default value for the various _ANSWER_DATE parameters.", &
344 default=99991231)
345 call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
346 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
347 "Values below 20190101 result in the use of older, less accurate expressions "//&
348 "that were in use at the end of 2018. Higher values result in the use of more "//&
349 "robust and accurate forms of mathematically equivalent expressions.", &
350 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
351 if (.not.gv%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
352 call set_regrid_params(cs, remap_answer_date=remap_answer_date)
353 call get_param(param_file, mdl, "REGRIDDING_ANSWER_DATE", regrid_answer_date, &
354 "The vintage of the expressions and order of arithmetic to use for regridding. "//&
355 "Values below 20190101 result in the use of older, less accurate expressions "//&
356 "that were in use at the end of 2018. Higher values result in the use of more "//&
357 "robust and accurate forms of mathematically equivalent expressions.", &
358 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
359 if (.not.gv%Boussinesq) regrid_answer_date = max(regrid_answer_date, 20230701)
360 call set_regrid_params(cs, regrid_answer_date=regrid_answer_date)
361 endif
362
363 if (main_parameters .and. coord_is_state_dependent) then
364 call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", tmplogical, &
365 "When defined, a proper high-order reconstruction "//&
366 "scheme is used within boundary cells rather "//&
367 "than PCM. E.g., if PPM is used for remapping, a "//&
368 "PPM reconstruction will also be used within "//&
369 "boundary cells.", default=regriddingdefaultboundaryextrapolation)
370 call set_regrid_params(cs, boundary_extrapolation=tmplogical)
371 else
372 call set_regrid_params(cs, boundary_extrapolation=.false.)
373 endif
374
375 ! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG)
376 if (main_parameters) then
377 param_name = "ALE_COORDINATE_CONFIG"
378 coord_res_param = "ALE_RESOLUTION"
379 string2 = 'UNIFORM'
380 else
381 param_name = create_coord_param(param_prefix, "DEF", param_suffix)
382 coord_res_param = create_coord_param(param_prefix, "RES", param_suffix)
383 string2 = 'UNIFORM'
384 if ((maximum_depth>3000.) .and. (maximum_depth<9250.)) string2='WOA09' ! For convenience
385 endif
386 call get_param(param_file, mdl, param_name, string, &
387 "Determines how to specify the coordinate "//&
388 "resolution. Valid options are:\n"//&
389 " PARAM - use the vector-parameter "//trim(coord_res_param)//"\n"//&
390 " UNIFORM[:N] - uniformly distributed\n"//&
391 " FILE:string - read from a file. The string specifies\n"//&
392 " the filename and variable name, separated\n"//&
393 " by a comma or space, e.g. FILE:lev.nc,dz\n"//&
394 " or FILE:lev.nc,interfaces=zw\n"//&
395 " WOA09[:N] - the WOA09 vertical grid (approximately)\n"//&
396 " WOA09INT[:N] - layers spanned by the WOA09 depths\n"//&
397 " WOA23INT[:N] - layers spanned by the WOA23 depths\n"//&
398 " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
399 " HYBRID:string - read from a file. The string specifies\n"//&
400 " the filename and two variable names, separated\n"//&
401 " by a comma or space, for sigma-2 and dz.\n"//&
402 " e.g. HYBRID:vgrid.nc,sigma2,dz\n"//&
403 " HYBRID_3D:string - read from a file. The string specifies\n"//&
404 " the filename and two 3D variable names, separated\n"//&
405 " by a comma or space, for sigma-2 and dz. The\n"//&
406 " latter can be FNC1:string which is used everywhere.\n"//&
407 " e.g. HYBRID_3D:vgrid.nc,sigma2,dz\n"//&
408 " HYBRID_MAP:string - read from a file. The string specifies\n"//&
409 " the filename and three variable names, separated\n"//&
410 " by a comma or space, for map, sigma-2 and dz.\n"//&
411 " Map is a spatial index array with, maxval(map)=N,\n"//&
412 " and the others are 2D arrays containing N profiles.\n"//&
413 " Map typically contains integer values, but it can\n"//&
414 " contain real values, I+w, which imply using\n"//&
415 " the weighted sum of profiles I and I+1.\n"//&
416 " Dz can be FNC1:string which is used everywhere.\n"//&
417 " e.g. HYBRID_MAP:vgrid.nc,map,sigma2,dz",&
418 default=trim(string2))
419 message = "The distribution of vertical resolution for the target\n"//&
420 "grid used for Eulerian-like coordinates. For example,\n"//&
421 "in z-coordinate mode, the parameter is a list of level\n"//&
422 "thicknesses (in m). In sigma-coordinate mode, the list\n"//&
423 "is of non-dimensional fractions of the water column."
424 if (index(trim(string),'UNIFORM')==1) then
425 if (len_trim(string)==7) then
426 ke = gv%ke ! Use model nk by default
427 tmpreal = maximum_depth
428 elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then
429 ! Format is "UNIFORM:N" or "UNIFORM:N,dz"
430 ke = extract_integer(string(9:len_trim(string)),'',1)
431 tmpreal = extract_real(string(9:len_trim(string)),',',2,missing_value=maximum_depth)
432 else
433 call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
434 'Unable to interpret "'//trim(string)//'".')
435 endif
436 allocate(dz(ke))
437 dz(:) = uniformresolution(ke, coord_mode, tmpreal, &
438 us%R_to_kg_m3*(gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(min(2,ke)))), &
439 us%R_to_kg_m3*(gv%Rlay(ke) + 0.5*(gv%Rlay(ke)-gv%Rlay(max(ke-1,1)))) )
440 if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
441 trim(message), units=trim(coord_units))
442 elseif (trim(string)=='PARAM') then
443 ! Read coordinate resolution (main model = ALE_RESOLUTION)
444 allocate(dz(1001))
445 dz(:) = -1. ! Setting to <0 allows detection of unset elements
446 call get_param(param_file, mdl, coord_res_param, dz, "Scan", units="", do_not_log=.true.)
447 if (dz(1001)>=0.) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
448 "PARAM specification is limited to 1000 values. Hack the code to use more!")
449 do ke=1,1000 ! Find number of defined levels
450 if (dz(ke+1)<0.) exit
451 enddo
452 deallocate(dz)
453 allocate(dz(ke)) ! Allocate with the correct number of levels, and re-read thicknesses
454 call get_param(param_file, mdl, coord_res_param, dz, &
455 trim(message), units=trim(coord_units), fail_if_missing=.true.)
456 elseif (index(trim(string),'FILE:')==1) then
457 ! FILE:filename,var_name is assumed to be reading level thickness variables
458 ! FILE:filename,interfaces=var_name reads positions
459 if (string(6:6)=='.' .or. string(6:6)=='/') then
460 ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
461 filename = trim( extractword(trim(string(6:80)), 1) )
462 else
463 ! Otherwise assume we should look for the file in INPUTDIR
464 filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
465 endif
466 if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
467 "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
468
469 varname = trim( extractword(trim(string(6:)), 2) )
470 if (len_trim(varname)==0) then
471 if (field_exists(filename,'dz')) then ; varname = 'dz'
472 elseif (field_exists(filename,'dsigma')) then ; varname = 'dsigma'
473 elseif (field_exists(filename,'ztest')) then ; varname = 'ztest'
474 else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
475 "Coordinate variable not specified and none could be guessed.")
476 endif
477 endif
478 ! This check fails when the variable is a dimension variable! -AJA
479 !if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
480 ! "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
481 if (cs%regridding_scheme == regridding_sigma) then
482 expected_units = 'nondim' ; alt_units = expected_units
483 elseif (cs%regridding_scheme == regridding_rho) then
484 expected_units = 'kg m-3' ; alt_units = expected_units
485 else
486 expected_units = 'meters' ; alt_units = 'm'
487 endif
488 if (index(trim(varname),'interfaces=')==1) then
489 varname=trim(varname(12:))
490 call verify_variable_units(filename, varname, expected_units, message, ierr, alt_units)
491 if (ierr) call mom_error(fatal, trim(mdl)//", initialize_regridding: "//&
492 "Unsupported format in grid definition '"//trim(filename)//&
493 "'. Error message "//trim(message))
494 call field_size(trim(filename), trim(varname), nzf)
495 ke = nzf(1)-1
496 if (ke < 1) call mom_error(fatal, trim(mdl)//" initialize_regridding via Var "//&
497 trim(varname)//"in FILE "//trim(filename)//&
498 " requires at least 2 target interface values.")
499 if (cs%regridding_scheme == regridding_rho) then
500 allocate(rho_target(ke+1))
501 call mom_read_data(trim(filename), trim(varname), rho_target)
502 else
503 allocate(dz(ke))
504 allocate(z_max(ke+1))
505 call mom_read_data(trim(filename), trim(varname), z_max)
506 dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
507 deallocate(z_max)
508 endif
509 else
510 ! Assume reading resolution
511 call field_size(trim(filename), trim(varname), nzf)
512 ke = nzf(1)
513 allocate(dz(ke))
514 call mom_read_data(trim(filename), trim(varname), dz)
515 endif
516 if (main_parameters .and. (ke/=gv%ke)) then
517 call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
518 'Mismatch in number of model levels and "'//trim(string)//'".')
519 endif
520 if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
521 trim(message), units=coordinateunits(coord_mode))
522 elseif (index(trim(string),'FNC1:')==1) then
523 ke = gv%ke ; allocate(dz(ke))
524 call dz_function1( trim(string(6:)), dz )
525 if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
526 trim(message), units=coordinateunits(coord_mode))
527 elseif (index(trim(string),'RFNC1:')==1) then
528 ! Function used for set target interface densities
529 ke = rho_function1( trim(string(7:)), rho_target )
530 elseif (index(trim(string),'HYBRID:')==1) then
531 ke = gv%ke
532 allocate(dz(ke))
533 allocate(rho_target(ke+1))
534 ! The following assumes the FILE: syntax of above but without "FILE:" in the string
535 varname = trim( extractword(trim(string(8:)), 3) )
536 if (varname == " ") call mom_error(fatal, &
537 trim(mdl)//", initialize_regridding: HYBRID "// &
538 "Too few arguments in ("//trim(string)//")")
539 filename = trim( extractword(trim(string(8:)), 1) )
540 if (filename(1:1)/='.' .and. filename(1:1)/='/') filename = trim(inputdir) // trim( filename )
541 if (.not. file_exists(filename)) call mom_error(fatal, &
542 trim(mdl)//", initialize_regridding: HYBRID "// &
543 "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
544 varname = trim( extractword(trim(string(8:)), 2) )
545 if (.not. field_exists(filename,varname)) call mom_error(fatal, &
546 trim(mdl)//", initialize_regridding: HYBRID "// &
547 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
548 call mom_read_data(trim(filename), trim(varname), rho_target)
549 varname = trim( extractword(trim(string(8:)), 3) )
550 if (varname(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz
551 call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
552 else ! Read dz from file
553 if (.not. field_exists(filename,varname)) call mom_error(fatal, &
554 trim(mdl)//", initialize_regridding: HYBRID "// &
555 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
556 call mom_read_data(trim(filename), trim(varname), dz)
557 endif
558 if (main_parameters) then
559 call log_param(param_file, mdl, "!"//coord_res_param, dz, &
560 trim(message), units=coordinateunits(coord_mode))
561 call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, &
562 'HYBRID target densities for interfaces', units="kg m-3")
563 endif
564 elseif (index(trim(string),'HYBRID_3D:')==1) then
565 ke = gv%ke
566 allocate(dz_3d(szi_(g),szj_(g),ke), source=0.0)
567 allocate(rho_target_3d(szi_(g),szj_(g),ke+1), source=0.0)
568 ! The following assumes the FILE: syntax of above but without "FILE:" in the string
569 varname = trim( extractword(trim(string(11:)), 3) )
570 if (varname == " ") call mom_error(fatal, &
571 trim(mdl)//", initialize_regridding: HYBRID_3D "// &
572 "Too few arguments in ("//trim(string)//")")
573 filename = trim( extractword(trim(string(11:)), 1) )
574 if (filename(1:1)/='.' .and. filename(1:1)/='/') filename = trim(inputdir) // trim( filename )
575 if (.not. file_exists(filename)) call mom_error(fatal, &
576 trim(mdl)//", initialize_regridding: HYBRID_3D "// &
577 "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
578 varname = trim( extractword(trim(string(11:)), 2) )
579 if (.not. field_exists(filename,varname)) call mom_error(fatal, &
580 trim(mdl)//", initialize_regridding: HYBRID_3D "// &
581 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
582 call mom_read_data(trim(filename), trim(varname), rho_target_3d, g%Domain)
583 call pass_var(rho_target_3d, g%Domain, halo=1)
584 varname = trim( extractword(trim(string(11:)), 3) )
585 if (varname(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz_3d
586 allocate(dz(ke))
587 call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
588 ! Adjust target grid to be consistent with maximum_depth
589 tmpreal = sum( dz(:) )
590 if (tmpreal < maximum_depth) then
591 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
592 endif
593 do i=g%isc-1,g%iec+1 ; do j=g%jsc-1,g%jec+1
594 if (g%mask2dT(i,j)>0.) then
595 do k=1,ke
596 dz_3d(i,j,k) = dz(k)
597 enddo
598 endif !mask2dT
599 enddo ; enddo
600 if (main_parameters) then
601 call log_param(param_file, mdl, "!"//coord_res_param, dz, &
602 trim(message), units=coordinateunits(coord_mode))
603 endif
604 else ! Read dz from file
605 if (.not. field_exists(filename,varname)) call mom_error(fatal, &
606 trim(mdl)//", initialize_regridding: HYBRID_3D "// &
607 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
608 call mom_read_data(trim(filename), trim(varname), dz_3d, g%Domain)
609 call pass_var(dz_3d, g%Domain, halo=1)
610 ! set nominal 1-d dz to UNIFORM
611 allocate(dz(ke))
612 dz(:) = uniformresolution(ke, coord_mode, maximum_depth, &
613 us%R_to_kg_m3*(gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(min(2,ke)))), &
614 us%R_to_kg_m3*(gv%Rlay(ke) + 0.5*(gv%Rlay(ke)-gv%Rlay(max(ke-1,1)))) )
615 endif !dz
616 elseif (index(trim(string),'HYBRID_MAP:')==1) then
617 ke = gv%ke
618 allocate(dz_3d(szi_(g),szj_(g),ke), source=0.0)
619 allocate(rho_target_3d(szi_(g),szj_(g),ke+1), source=0.0)
620 allocate(index_map(szi_(g),szj_(g)), source=1.0)
621 ! The following assumes the FILE: syntax of above but without "FILE:" in the string
622 varname = trim( extractword(trim(string(12:)), 4) )
623 if (varname == " ") call mom_error(fatal, &
624 trim(mdl)//", initialize_regridding: HYBRID_3D "// &
625 "Too few arguments in ("//trim(string)//")")
626 filename = trim( extractword(trim(string(12:)), 1) )
627 if (filename(1:1)/='.' .and. filename(1:1)/='/') filename = trim(inputdir) // trim( filename )
628 if (.not. file_exists(filename)) call mom_error(fatal, &
629 trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
630 "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
631 varname = trim( extractword(trim(string(12:)), 2) )
632 if (.not. field_exists(filename,varname)) call mom_error(fatal, &
633 trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
634 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
635 call mom_read_data(trim(filename), trim(varname), index_map, g%Domain)
636 call pass_var(index_map, g%Domain, halo=1)
637 !find maximum index
638 np = 1
639 do j=g%jsc, g%jec ; do i=g%isc, g%iec
640 np = max(np,ceiling(index_map(i,j)))
641 enddo ; enddo
642 call max_across_pes(np)
643 write(string2,"(i3)") np
644 call mom_error(note, &
645 trim(mdl)//", initialize_regridding: HYBRID_MAP NP="//trim(string2))
646 if (np<1) call mom_error(fatal, &
647 trim(mdl)//", initialize_regridding: HYBRID_MAP to small NP from "//trim(varname))
648 allocate(dz_2d(ke,np))
649 allocate(rho_target_2d(ke+1,np))
650 varname = trim( extractword(trim(string(12:)), 3) )
651 if (.not. field_exists(filename,varname)) call mom_error(fatal, &
652 trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
653 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
654 ! MOM_read_data can't handle this array
655 call read_variable(trim(filename), trim(varname), rho_target_2d)
656 if (main_parameters) then
657 call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target_2d(:,1), &
658 'HYBRID target densities for interfaces', units="kg m-3")
659 endif
660 do i=g%isc-1,g%iec+1 ; do j=g%jsc-1,g%jec+1
661 if (g%mask2dT(i,j)>0.) then
662 nfloor = floor(index_map(i,j))
663 nceiling = ceiling(index_map(i,j))
664 if (nfloor<1 .or. nceiling>np) then
665 write(0,'(a,2i5,a,g20.6)') 'HYBRID_MAP: i,j=',i,j,'index_map(i,j)=', index_map(i,j)
666 call mom_error(fatal, trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
667 "index_map out of range")
668 endif
669 if (nfloor == nceiling) then
670 do k=1,ke+1
671 rho_target_3d(i,j,k) = rho_target_2d(k,nfloor)
672 enddo
673 else
674 nfrac = index_map(i,j) - nfloor !between 0.0 and 1.0
675 do k=1,ke+1
676 rho_target_3d(i,j,k) = (1.0-nfrac)*rho_target_2d(k,nfloor) + &
677 nfrac *rho_target_2d(k,nceiling)
678 enddo
679 endif !integer:else
680 endif !mask2dT
681 enddo ; enddo
682 varname = trim( extractword(trim(string(12:)), 4) )
683 if (varname(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz_3d
684 allocate(dz(ke))
685 call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
686 ! Adjust target grid to be consistent with maximum_depth
687 tmpreal = sum( dz(:) )
688 if (tmpreal < maximum_depth) then
689 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
690 endif
691 do i=g%isc-1,g%iec+1 ; do j=g%jsc-1,g%jec+1
692 if (g%mask2dT(i,j)>0.) then
693 do k=1,ke
694 dz_3d(i,j,k) = dz(k)
695 enddo
696 endif !mask2dT
697 enddo ; enddo
698 if (main_parameters) then
699 call log_param(param_file, mdl, "!"//coord_res_param, dz, &
700 trim(message), units=coordinateunits(coord_mode))
701 endif
702 else ! Read dz from file
703 if (.not. field_exists(filename,varname)) call mom_error(fatal, &
704 trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
705 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
706 ! MOM_read_data can't handle this array
707 call read_variable(trim(filename), trim(varname), dz_2d)
708 if (main_parameters) then
709 call log_param(param_file, mdl, "!"//coord_res_param, dz_2d(:,1), &
710 trim(message), units=coordinateunits(coord_mode))
711 endif
712 do i=1,np
713 tmpreal = sum( dz_2d(:,i) )
714 if (tmpreal < maximum_depth) then
715 dz_2d(ke,i) = dz_2d(ke,i) + ( maximum_depth - tmpreal )
716 endif
717 enddo
718 allocate(dz(ke))
719 dz(:) = dz_2d(:,1)
720 if (main_parameters) then
721 call log_param(param_file, mdl, "!"//coord_res_param, dz, &
722 trim(message), units=coordinateunits(coord_mode))
723 endif
724 do i=g%isc-1,g%iec+1 ; do j=g%jsc-1,g%jec+1
725 if (g%mask2dT(i,j)>0.) then
726 nfloor = floor(index_map(i,j))
727 nceiling = ceiling(index_map(i,j))
728 if (nfloor == nceiling) then
729 do k=1,ke
730 dz_3d(i,j,k) = dz_2d(k,nfloor)
731 enddo
732 else
733 nfrac = index_map(i,j) - nfloor !between 0.0 and 1.0
734 do k=1,ke
735 dz_3d(i,j,k) = (1.0-nfrac)*dz_2d(k,nfloor) + &
736 nfrac *dz_2d(k,nceiling)
737 enddo
738 endif !integer:else
739 endif !mask2dT
740 enddo ; enddo
741 endif !dz
742 deallocate(index_map)
743 deallocate(rho_target_2d)
744 deallocate(dz_2d)
745 elseif (index(trim(string),'WOA09INT')==1) then
746 if (len_trim(string)==8) then ! string=='WOA09INT'
747 tmpreal = 0. ; ke = 0 ; dz_extra = 0.
748 do while (tmpreal<maximum_depth)
749 ke = ke + 1
750 if (ke > size(woa09_dzi)) then
751 dz_extra = maximum_depth - tmpreal
752 exit
753 endif
754 tmpreal = tmpreal + woa09_dzi(ke)
755 enddo
756 elseif (index(trim(string),'WOA09INT:')==1) then ! string starts with 'WOA09INT:'
757 if (len_trim(string)==9) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
758 'Expected string of form "WOA09INT:N" but got "'//trim(string)//'".')
759 ke = extract_integer(string(10:len_trim(string)),'',1)
760 if (ke>39 .or. ke<1) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
761 'For "WOA05INT:N" N must 0<N<40 but got "'//trim(string)//'".')
762 endif
763 allocate(dz(ke))
764 do k=1,min(ke, size(woa09_dzi))
765 dz(k) = woa09_dzi(k)
766 enddo
767 if (ke > size(woa09_dzi)) dz(ke) = dz_extra
768 elseif (index(trim(string),'WOA23INT')==1) then
769 if (len_trim(string)==8) then ! string=='WOA23INT'
770 tmpreal = 0. ; ke = 0 ; dz_extra = 0.
771 do while (tmpreal<maximum_depth)
772 ke = ke + 1
773 if (ke > size(woa23_dzi)) then
774 dz_extra = maximum_depth - tmpreal
775 exit
776 endif
777 tmpreal = tmpreal + woa23_dzi(ke)
778 enddo
779 elseif (index(trim(string),'WOA23INT:')==1) then ! string starts with 'WOA23INT:'
780 if (len_trim(string)==9) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
781 'Expected string of form "WOA23INT:N" but got "'//trim(string)//'".')
782 ke = extract_integer(string(10:len_trim(string)),'',1)
783 if (ke>39 .or. ke<1) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
784 'For "WOA05INT:N" N must 0<N<40 but got "'//trim(string)//'".')
785 endif
786 allocate(dz(ke))
787 do k=1,min(ke, size(woa23_dzi))
788 dz(k) = woa23_dzi(k)
789 enddo
790 if (ke > size(woa23_dzi)) dz(ke) = dz_extra
791 elseif (index(trim(string),'WOA09')==1) then
792 if (len_trim(string)==5) then ! string=='WOA09'
793 tmpreal = 0. ; ke = 0 ; dz_extra = 0.
794 do while (tmpreal<maximum_depth)
795 ke = ke + 1
796 if (ke > size(woa09_dz_approx)) then
797 dz_extra = maximum_depth - tmpreal
798 exit
799 endif
800 tmpreal = tmpreal + woa09_dz_approx(ke)
801 enddo
802 elseif (index(trim(string),'WOA09:')==1) then ! string starts with 'WOA09:'
803 if (len_trim(string)==6) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
804 'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
805 ke = extract_integer(string(7:len_trim(string)),'',1)
806 if (ke>40 .or. ke<1) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
807 'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
808 endif
809 allocate(dz(ke))
810 do k=1,min(ke, size(woa09_dz_approx))
811 dz(k) = woa09_dz_approx(k)
812 enddo
813 if (ke > size(woa09_dz_approx)) dz(ke) = dz_extra
814 else
815 call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
816 "Unrecognized coordinate configuration"//trim(string))
817 endif
818
819 if (main_parameters) then
820 ! This is a work around to apparently needed to work with the from_Z initialization... ???
821 if (coordinatemode(coord_mode) == regridding_zstar .or. &
822 coordinatemode(coord_mode) == regridding_hycom1 .or. &
823 coordinatemode(coord_mode) == regridding_hybgen .or. &
824 coordinatemode(coord_mode) == regridding_adaptive) then
825 if (allocated(dz)) then
826 ! Adjust target grid to be consistent with maximum_depth
827 tmpreal = sum( dz(:) )
828 if (tmpreal < maximum_depth) then
829 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
830 elseif (tmpreal > maximum_depth) then
831 if ( dz(ke) + ( maximum_depth - tmpreal ) > 0. ) then
832 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
833 else
834 call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
835 "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
836 endif
837 endif
838 endif !allocated(dz)
839 endif
840 endif
841
842 if (coordinatemode(coord_mode) == regridding_hycom1) then
843 allocate(dz_shallow(ke))
844 call get_param(param_file, mdl, "SHALLOW_"//trim(coord_res_param), dz_shallow, &
845 "HYBGEN-style Z-sigma-Z near surface fixed coordinate. "//&
846 "The default of all zeros turns this option off. "//&
847 "Let N_SIGMA be the number of consecutive non-zero entries, typically < NK. "//&
848 "Use SHALLOW_"//trim(coord_res_param)//" when rest depth is shallower than "//&
849 "SUM(SHALLOW_"//trim(coord_res_param)//"(1:N_SIGMA)). "//&
850 "Use "//trim(coord_res_param)//" when rest depth is deeper than "//&
851 "SUM("//trim(coord_res_param)//"(1:N_SIGMA)). "//&
852 "Otherwise use a linear sum of the two weighted by rest depth.",&
853 units="m", default=0.0)
854 n_sigma = ke
855 depth_s = 0.0
856 do k= 1,ke
857 depth_s = depth_s + dz_shallow(k)
858 if (dz_shallow(k) == 0.0) then
859 n_sigma = k-1
860 exit
861 endif
862 enddo
863 if (n_sigma > 0) then
864 if (main_parameters) call log_param(param_file, mdl, "!N_SIGMA", n_sigma, &
865 "Number of consecutive non-zero entries in SHALLOW_"//&
866 trim(coord_res_param)//".")
867 if (.not.allocated(dz_3d)) then
868 allocate(dz_3d(szi_(g),szj_(g),ke), source=0.0)
869 allocate(rho_target_3d(szi_(g),szj_(g),ke+1), source=0.0)
870 do i=g%isc-1,g%iec+1 ; do j=g%jsc-1,g%jec+1
871 if (g%mask2dT(i,j)>0.) then
872 do k=1,ke
873 dz_3d(i,j,k) = dz(k)
874 enddo
875 do k=1,ke+1
876 rho_target_3d(i,j,k) = rho_target(k)
877 enddo
878 endif !mask2dT
879 enddo ; enddo
880 endif
881 do i=g%isc-1,g%iec+1 ; do j=g%jsc-1,g%jec+1
882 if (g%mask2dT(i,j)>0.) then
883 nominaldepth = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0) * us%Z_to_m
884 if (nominaldepth <= depth_s) then
885 do k= 1,n_sigma
886 dz_3d(i,j,k) = dz_shallow(k)
887 enddo
888 do k= n_sigma+1,ke
889 dz_3d(i,j,k) = dz_shallow(n_sigma)
890 enddo
891 else ! >depth_s
892 depth_d = 0.0
893 do k= 1,n_sigma
894 depth_d = depth_d + dz_3d(i,j,k)
895 enddo
896 ! do nothing if nominalDepth >= depth_d
897 if (nominaldepth < depth_d) then
898 depth_q = (nominaldepth - depth_s) / (depth_d - depth_s)
899 do k= 1,n_sigma
900 dz_3d(i,j,k) = (1.0-depth_q)*dz_shallow(k) + depth_q*dz_3d(i,j,k)
901 enddo
902 do k= n_sigma+1,ke
903 dz_3d(i,j,k) = (1.0-depth_q)*dz_shallow(n_sigma) + depth_q*dz_3d(i,j,k)
904 enddo
905 endif !<depth_d and >depth_s
906 endif !nominalDepth
907 endif !mask2dT
908 enddo ; enddo
909 endif !n_sigma
910 deallocate(dz_shallow)
911 endif !REGRIDDING_HYCOM1
912
913 cs%nk=ke
914
915 ! Target resolution (for fixed coordinates)
916 if (allocated(dz_3d)) then
917 allocate( cs%coordinateResolution(cs%nk), source=-1.e30 )
918 allocate( cs%coordinateResolution_3d(szi_(g),szj_(g),cs%nk), source=-1.e30 )
919 allocate( cs%target_density_3d(szi_(g),szj_(g),cs%nk+1), source=-1.e30*us%kg_m3_to_R )
920 else
921 allocate( cs%coordinateResolution(cs%nk), source=-1.e30 )
922 if (state_dependent(cs%regridding_scheme)) then
923 ! Target values
924 allocate( cs%target_density(cs%nk+1), source=-1.e30*us%kg_m3_to_R )
925 endif
926 endif
927
928 if (allocated(dz_3d)) then
929 ! set both 1d and 3d fields
930 call setcoordinateresolution(dz, cs, scale=us%m_to_Z)
931 call setcoordinateresolution_3d(dz_3d, cs, scale=us%m_to_Z)
932 cs%coord_scale = us%Z_to_m
933 deallocate(dz_3d)
934 elseif (allocated(dz)) then
935 if (coordinatemode(coord_mode) == regridding_sigma) then
936 call setcoordinateresolution(dz, cs, scale=1.0)
937 elseif (coordinatemode(coord_mode) == regridding_rho) then
938 call setcoordinateresolution(dz, cs, scale=us%kg_m3_to_R)
939 elseif (coordinatemode(coord_mode) == regridding_adaptive) then
940 call setcoordinateresolution(dz, cs, scale=gv%m_to_H)
941 cs%coord_scale = gv%H_to_m
942 else
943 call setcoordinateresolution(dz, cs, scale=us%m_to_Z)
944 cs%coord_scale = us%Z_to_m
945 endif
946 endif
947
948 ! set coord_scale for RHO regridding independent of allocation status of dz
949 if (coordinatemode(coord_mode) == regridding_rho) then
950 cs%coord_scale = us%R_to_kg_m3
951 endif
952
953 ! ensure CS%ref_pressure is rescaled properly
954 cs%ref_pressure = us%Pa_to_RL2_T2 * cs%ref_pressure
955
956 if (allocated(rho_target_3d)) then
957 call set_target_densities_3d(cs, g, us%kg_m3_to_R, rho_target_3d)
958 deallocate(rho_target_3d)
959 elseif (allocated(rho_target)) then
960 call set_target_densities(cs, us%kg_m3_to_R*rho_target)
961 deallocate(rho_target)
962 elseif (coordinatemode(coord_mode) == regridding_rho) then
963 call set_target_densities_from_gv(gv, us, cs)
964 call log_param(param_file, mdl, "!TARGET_DENSITIES", us%R_to_kg_m3*cs%target_density(:), &
965 'RHO target densities for interfaces', "kg m-3")
966 endif
967
968 ! initialise coordinate-specific control structure
969 call initcoord(cs, g, gv, us, coord_mode, param_file)
970
971 if (coord_is_state_dependent) then
972 if (main_parameters) then
973 call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), &
974 p_ref, &
975 "The pressure that is used for calculating the coordinate "//&
976 "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
977 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
978 units="Pa", default=2.0e7, scale=us%Pa_to_RL2_T2)
979 else
980 call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), &
981 p_ref, &
982 "The pressure that is used for calculating the diagnostic coordinate "//&
983 "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
984 "This is only used for the RHO coordinate.", &
985 units="Pa", default=2.0e7, scale=us%Pa_to_RL2_T2)
986 endif
987 call get_param(param_file, mdl, create_coord_param(param_prefix, &
988 "REGRID_COMPRESSIBILITY_FRACTION", param_suffix), tmpreal, &
989 "When interpolating potential density profiles we can add "//&
990 "some artificial compressibility solely to make homogeneous "//&
991 "regions appear stratified.", units="nondim", default=0.)
992 call set_regrid_params(cs, compress_fraction=tmpreal, ref_pressure=p_ref)
993 endif
994
995 if (main_parameters) then
996 call get_param(param_file, mdl, "MIN_THICKNESS", tmpreal, &
997 "When regridding, this is the minimum layer "//&
998 "thickness allowed.", units="m", scale=gv%m_to_H, &
999 default=regriddingdefaultminthickness )
1000 call set_regrid_params(cs, min_thickness=tmpreal)
1001 call get_param(param_file, mdl, "USE_ADJUST_INTERFACE_MOTION", tmplogical, &
1002 "When regridding, after the primary grid generation, call a function that ensures "//&
1003 "positive layer thicknesses. Historically, this was required.", default=.true.)
1004 call set_regrid_params(cs, use_adjust_interface_motion=tmplogical)
1005 else
1006 call set_regrid_params(cs, min_thickness=0.)
1007 call set_regrid_params(cs, use_adjust_interface_motion=.true.)
1008 call set_regrid_params(cs, use_depth_based_time_filter=.true.)
1009 endif
1010
1011 if (main_parameters .and. coordinatemode(coord_mode) == regridding_hycom1) then
1012 call get_param(param_file, mdl, "HYCOM1_ONLY_IMPROVES", tmplogical, &
1013 "When regridding, an interface is only moved if this improves "//&
1014 "the fit to the target density.", default=.false.)
1015 call set_hycom_params(cs%hycom_CS, only_improves=tmplogical)
1016 endif
1017
1018 cs%use_hybgen_unmix = .false.
1019 if (coordinatemode(coord_mode) == regridding_hybgen) then
1020 call get_param(param_file, mdl, "USE_HYBGEN_UNMIX", cs%use_hybgen_unmix, &
1021 "If true, use hybgen unmixing code before regridding.", &
1022 default=.false.)
1023 endif
1024
1025 if (coordinatemode(coord_mode) == regridding_adaptive) then
1026 call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adapttimeratio, &
1027 "Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1)
1028 call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptzoom, &
1029 "Depth of near-surface zooming region.", units="m", default=200.0, scale=gv%m_to_H)
1030 call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptzoomcoeff, &
1031 "Coefficient of near-surface zooming diffusivity.", units="nondim", default=0.2)
1032 call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptbuoycoeff, &
1033 "Coefficient of buoyancy diffusivity.", units="nondim", default=0.8)
1034 call get_param(param_file, mdl, "ADAPT_ALPHA", adaptalpha, &
1035 "Scaling on optimization tendency.", units="nondim", default=1.0)
1036 call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmplogical, &
1037 "If true, make a HyCOM-like mixed layer by preventing interfaces "//&
1038 "from being shallower than the depths specified by the regridding coordinate.", &
1039 default=.false.)
1040 call get_param(param_file, mdl, "ADAPT_DRHO0", adaptdrho0, &
1041 "Reference density difference for stratification-dependent diffusion.", &
1042 units="kg m-3", default=0.5, scale=us%kg_m3_to_R)
1043
1044 call set_regrid_params(cs, adapttimeratio=adapttimeratio, adaptzoom=adaptzoom, &
1045 adaptzoomcoeff=adaptzoomcoeff, adaptbuoycoeff=adaptbuoycoeff, adaptalpha=adaptalpha, &
1046 adaptdomin=tmplogical, adaptdrho0=adaptdrho0)
1047 endif
1048
1049 if (main_parameters .and. coord_is_state_dependent) then
1050 call get_param(param_file, mdl, "MAXIMUM_INT_DEPTH_CONFIG", string, &
1051 "Determines how to specify the maximum interface depths.\n"//&
1052 "Valid options are:\n"//&
1053 " NONE - there are no maximum interface depths\n"//&
1054 " PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
1055 " FILE:string - read from a file. The string specifies\n"//&
1056 " the filename and variable name, separated\n"//&
1057 " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
1058 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
1059 default='NONE')
1060 message = "The list of maximum depths for each interface."
1061 allocate(z_max(ke+1))
1062 allocate(dz_max(ke))
1063 if ( trim(string) == "NONE") then
1064 ! Do nothing.
1065 elseif ( trim(string) == "PARAM") then
1066 call get_param(param_file, mdl, "MAXIMUM_INTERFACE_DEPTHS", z_max, &
1067 trim(message), units="m", scale=gv%m_to_H, fail_if_missing=.true.)
1068 call set_regrid_max_depths(cs, z_max)
1069 elseif (index(trim(string),'FILE:')==1) then
1070 if (string(6:6)=='.' .or. string(6:6)=='/') then
1071 ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
1072 filename = trim( extractword(trim(string(6:80)), 1) )
1073 else
1074 ! Otherwise assume we should look for the file in INPUTDIR
1075 filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
1076 endif
1077 if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)// &
1078 ", initialize_regridding: "// &
1079 "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
1080
1081 do_sum = .false.
1082 varname = trim( extractword(trim(string(6:)), 2) )
1083 if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)// &
1084 ", initialize_regridding: "// &
1085 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
1086 if (len_trim(varname)==0) then
1087 if (field_exists(filename,'z_max')) then ; varname = 'z_max'
1088 elseif (field_exists(filename,'dz')) then ; varname = 'dz' ; do_sum = .true.
1089 elseif (field_exists(filename,'dz_max')) then ; varname = 'dz_max' ; do_sum = .true.
1090 else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
1091 "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
1092 endif
1093 endif
1094 if (do_sum) then
1095 call mom_read_data(trim(filename), trim(varname), dz_max)
1096 z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
1097 else
1098 call mom_read_data(trim(filename), trim(varname), z_max)
1099 endif
1100 call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
1101 trim(message), units=coordinateunits(coord_mode))
1102 call set_regrid_max_depths(cs, z_max, gv%m_to_H)
1103 elseif (index(trim(string),'FNC1:')==1) then
1104 call dz_function1( trim(string(6:)), dz_max )
1105 z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
1106 call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
1107 trim(message), units=coordinateunits(coord_mode))
1108 call set_regrid_max_depths(cs, z_max, gv%m_to_H)
1109 else
1110 call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
1111 "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
1112 endif
1113 deallocate(z_max)
1114 deallocate(dz_max)
1115
1116 ! Optionally specify maximum thicknesses for each layer, enforced by moving
1117 ! the interface below a layer downward.
1118 call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", longstring, &
1119 "Determines how to specify the maximum layer thicknesses.\n"//&
1120 "Valid options are:\n"//&
1121 " NONE - there are no maximum layer thicknesses\n"//&
1122 " PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
1123 " FILE:string - read from a file. The string specifies\n"//&
1124 " the filename and variable name, separated\n"//&
1125 " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
1126 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
1127 default='NONE')
1128 message = "The list of maximum thickness for each layer."
1129 allocate(h_max(ke))
1130 if ( trim(longstring) == "NONE") then
1131 ! Do nothing.
1132 elseif ( trim(longstring) == "PARAM") then
1133 call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, &
1134 trim(message), units="m", fail_if_missing=.true., scale=gv%m_to_H)
1135 call set_regrid_max_thickness(cs, h_max)
1136 elseif (index(trim(longstring),'FILE:')==1) then
1137 if (longstring(6:6)=='.' .or. longstring(6:6)=='/') then
1138 ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
1139 filename = trim( extractword(trim(longstring(6:200)), 1) )
1140 else
1141 ! Otherwise assume we should look for the file in INPUTDIR
1142 filename = trim(inputdir) // trim( extractword(trim(longstring(6:200)), 1) )
1143 endif
1144 if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)// &
1145 ", initialize_regridding: "// &
1146 "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(longstring)//")")
1147
1148 varname = trim( extractword(trim(longstring(6:)), 2) )
1149 if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)// &
1150 ", initialize_regridding: "// &
1151 "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(longstring)//")")
1152 if (len_trim(varname)==0) then
1153 if (field_exists(filename,'h_max')) then ; varname = 'h_max'
1154 elseif (field_exists(filename,'dz_max')) then ; varname = 'dz_max'
1155 else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
1156 "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
1157 endif
1158 endif
1159 call mom_read_data(trim(filename), trim(varname), h_max)
1160 call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
1161 trim(message), units=coordinateunits(coord_mode))
1162 call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
1163 elseif (index(trim(longstring),'FNC1:')==1) then
1164 call dz_function1( trim(longstring(6:)), h_max )
1165 call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
1166 trim(message), units=coordinateunits(coord_mode))
1167 call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
1168 else
1169 call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
1170 "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(longstring))
1171 endif
1172 deallocate(h_max)
1173 endif
1174
1175 if (allocated(dz)) deallocate(dz)
1176end subroutine initialize_regridding
1177
1178
1179
1180!> Deallocation of regridding memory
1181subroutine end_regridding(CS)
1182 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
1183
1184 if (associated(cs%zlike_CS)) call end_coord_zlike(cs%zlike_CS)
1185 if (associated(cs%sigma_CS)) call end_coord_sigma(cs%sigma_CS)
1186 if (associated(cs%rho_CS)) call end_coord_rho(cs%rho_CS)
1187 if (associated(cs%hycom_CS)) call end_coord_hycom(cs%hycom_CS)
1188 if (associated(cs%adapt_CS)) call end_coord_adapt(cs%adapt_CS)
1189 if (associated(cs%hybgen_CS)) call end_hybgen_regrid(cs%hybgen_CS)
1190
1191 deallocate( cs%coordinateResolution )
1192 if (allocated(cs%coordinateResolution_3d)) deallocate( cs%coordinateResolution_3d )
1193 if (allocated(cs%target_density_3d)) deallocate( cs%target_density_3d )
1194 if (allocated(cs%target_density)) deallocate( cs%target_density )
1195 if (allocated(cs%max_interface_depths) ) deallocate( cs%max_interface_depths )
1196 if (allocated(cs%max_layer_thickness) ) deallocate( cs%max_layer_thickness )
1197
1198end subroutine end_regridding
1199
1200!------------------------------------------------------------------------------
1201!> Dispatching regridding routine for orchestrating regridding & remapping
1202subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, &
1203 frac_shelf_h, PCM_cell)
1204!------------------------------------------------------------------------------
1205! This routine takes care of (1) building a new grid and (2) remapping between
1206! the old grid and the new grid. The creation of the new grid can be based
1207! on z coordinates, target interface densities, sigma coordinates or any
1208! arbitrary coordinate system.
1209! The MOM6 interface positions are always calculated from the bottom up by
1210! accumulating the layer thicknesses starting at z=-G%bathyT. z increases
1211! upwards (decreasing k-index).
1212! The new grid is defined by the change in position of those interfaces in z
1213! dzInterface = zNew - zOld.
1214! Thus, if the regridding inflates the top layer, hNew(1) > hOld(1), then the
1215! second interface moves downward, zNew(2) < zOld(2), and dzInterface(2) < 0.
1216! hNew(k) = hOld(k) - dzInterface(k+1) + dzInterface(k)
1217! IMPORTANT NOTE:
1218! This is the converse of the sign convention used in the remapping code!
1219!------------------------------------------------------------------------------
1220
1221 ! Arguments
1222 type(remapping_cs), intent(in) :: remapcs !< Remapping parameters and options
1223 type(regridding_cs), intent(in) :: cs !< Regridding control structure
1224 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1225 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1226 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1227 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Current 3D grid obtained after
1228 !! the last time step [H ~> m or kg m-2]
1229 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamical variables (T, S, ...)
1230 real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target
1231 !! coordinate [H ~> m or kg m-2]
1232 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzinterface !< The change in position of each
1233 !! interface [H ~> m or kg m-2]
1234 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]
1235 logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1236 optional, intent(out ) :: pcm_cell !< Use PCM remapping in cells where true
1237
1238 ! Local variables
1239 real :: nom_depth_h(szi_(g),szj_(g)) !< The nominal ocean depth at each point in thickness units [H ~> m or kg m-2]
1240 real :: tot_h(szi_(g),szj_(g)) !< The total thickness of the water column [H ~> m or kg m-2]
1241 real :: tot_dz(szi_(g),szj_(g)) !< The total distance between the top and bottom of the water column [Z ~> m]
1242 real :: z_to_h ! A conversion factor used by some routines to convert coordinate
1243 ! parameters to depth units [H Z-1 ~> nondim or kg m-3]
1244 character(len=128) :: mesg ! A string for error messages
1245 integer :: i, j, k
1246
1247 if (present(pcm_cell)) pcm_cell(:,:,:) = .false.
1248
1249 z_to_h = us%Z_to_m * gv%m_to_H ! Often this is equivalent to GV%Z_to_H.
1250
1251 if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then
1252 if (tv%valid_SpV_halo < 0) then
1253 mesg = "invalid values of SpV_avg."
1254 else
1255 mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed."
1256 endif
1257 call mom_error(fatal, "Regridding_main called in fully non-Boussinesq mode with "//trim(mesg))
1258 endif
1259
1260 if (allocated(tv%SpV_avg)) then ! This is the fully non-Boussinesq case
1261 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
1262 tot_h(i,j) = 0.0 ; tot_dz(i,j) = 0.0
1263 enddo ; enddo
1264 do k=1,gv%ke ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
1265 tot_h(i,j) = tot_h(i,j) + h(i,j,k)
1266 tot_dz(i,j) = tot_dz(i,j) + gv%H_to_RZ * tv%SpV_avg(i,j,k) * h(i,j,k)
1267 enddo ; enddo ; enddo
1268 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
1269 if (tot_dz(i,j) > 0.0) then
1270 nom_depth_h(i,j) = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0) * (tot_h(i,j) / tot_dz(i,j))
1271 else
1272 nom_depth_h(i,j) = 0.0
1273 endif
1274 enddo ; enddo
1275 else
1276 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
1277 nom_depth_h(i,j) = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0) * z_to_h
1278 enddo ; enddo
1279 endif
1280
1281 select case ( cs%regridding_scheme )
1282
1283 case ( regridding_zstar )
1284 call build_zstar_grid( cs, g, gv, h, nom_depth_h, dzinterface, frac_shelf_h, zscale=z_to_h )
1285 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1286 case ( regridding_sigma_shelf_zstar)
1287 call build_zstar_grid( cs, g, gv, h, nom_depth_h, dzinterface, zscale=z_to_h )
1288 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1289 case ( regridding_sigma )
1290 call build_sigma_grid( cs, g, gv, h, nom_depth_h, dzinterface )
1291 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1292 case ( regridding_rho )
1293 call build_rho_grid( g, gv, g%US, h, nom_depth_h, tv, dzinterface, remapcs, cs, frac_shelf_h )
1294 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1295 case ( regridding_hycom1 )
1296 call build_grid_hycom1( g, gv, g%US, h, nom_depth_h, tv, h_new, dzinterface, remapcs, cs, &
1297 frac_shelf_h, zscale=z_to_h )
1298 case ( regridding_hybgen )
1299 call hybgen_regrid(g, gv, g%US, h, nom_depth_h, tv, cs%hybgen_CS, dzinterface, pcm_cell)
1300 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1301 case ( regridding_adaptive )
1302 call build_grid_adaptive(g, gv, g%US, h, nom_depth_h, tv, dzinterface, remapcs, cs)
1303 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1304
1305 case ( regridding_arbitrary )
1306 call mom_error(fatal,'MOM_regridding, regridding_main: '//&
1307 'Regridding mode "ARB" is not implemented.')
1308 case default
1309 call mom_error(fatal,'MOM_regridding, regridding_main: '//&
1310 'Unknown regridding scheme selected!')
1311
1312 end select ! type of grid
1313
1314#ifdef __DO_SAFETY_CHECKS__
1315 if (cs%nk == gv%ke) then
1316 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1 ; if (g%mask2dT(i,j)>0.) then
1317 call check_grid_column( gv%ke, h(i,j,:), dzinterface(i,j,:), 'in regridding_main')
1318 endif ; enddo ; enddo
1319 endif
1320#endif
1321 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (g%mask2dT(i,j) > 0.) then
1322 if (minval(h(i,j,:)) < 0.0) then
1323 write(0,*) 'regridding_main check_grid: i,j=', i, j, 'h_new(i,j,:)=', h_new(i,j,:)
1324 call mom_error(fatal, "regridding_main: negative thickness encountered.")
1325 endif
1326 endif ; enddo ; enddo
1327
1328end subroutine regridding_main
1329
1330!------------------------------------------------------------------------------
1331!> This routine returns flags indicating which pre-remapping state adjustments
1332!! are needed depending on the coordinate mode in use.
1333subroutine regridding_preadjust_reqs(CS, do_conv_adj, do_hybgen_unmix, hybgen_CS)
1334
1335 ! Arguments
1336 type(regridding_cs), intent(in) :: cs !< Regridding control structure
1337 logical, intent(out) :: do_conv_adj !< Convective adjustment should be done
1338 logical, intent(out) :: do_hybgen_unmix !< Hybgen unmixing should be done
1339 type(hybgen_regrid_cs), pointer, &
1340 optional, intent(out) :: hybgen_cs !< Control structure for hybgen regridding for sharing parameters.
1341
1342
1343 do_conv_adj = .false. ; do_hybgen_unmix = .false.
1344 select case ( cs%regridding_scheme )
1345
1346 case ( regridding_zstar, regridding_sigma_shelf_zstar, regridding_sigma, regridding_arbitrary, &
1347 regridding_hycom1, regridding_adaptive )
1348 do_conv_adj = .false. ; do_hybgen_unmix = .false.
1349 case ( regridding_rho )
1350 do_conv_adj = .true. ; do_hybgen_unmix = .false.
1351 case ( regridding_hybgen )
1352 do_conv_adj = .false. ; do_hybgen_unmix = cs%use_hybgen_unmix
1353 case default
1354 call mom_error(fatal,'MOM_regridding, regridding_preadjust_reqs: '//&
1355 'Unknown regridding scheme selected!')
1356 end select ! type of grid
1357
1358 if (present(hybgen_cs) .and. do_hybgen_unmix) hybgen_cs => cs%hybgen_CS
1359
1360end subroutine regridding_preadjust_reqs
1361
1362
1363!> Calculates h_new from h + delta_k dzInterface
1364subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1365 type(regridding_cs), intent(in) :: CS !< Regridding control structure
1366 type(ocean_grid_type), intent(in) :: G !< Grid structure
1367 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1368 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Old layer thicknesses [H ~> m or kg m-2]
1369 !! or other units
1370 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: dzInterface !< Change in interface positions
1371 !! in the same units as h [H ~> m or kg m-2]
1372 real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses in the same
1373 !! units as h [H ~> m or kg m-2]
1374 ! Local variables
1375 integer :: i, j, k, nki
1376
1377 nki = min(cs%nk, gv%ke)
1378
1379 !$OMP parallel do default(shared)
1380 do j = g%jsc-1,g%jec+1
1381 do i = g%isc-1,g%iec+1
1382 if (g%mask2dT(i,j)>0.) then
1383 do k=1,nki
1384 h_new(i,j,k) = max( 0., h(i,j,k) + ( dzinterface(i,j,k) - dzinterface(i,j,k+1) ) )
1385 enddo
1386 if (cs%nk > gv%ke) then
1387 do k=nki+1, cs%nk
1388 h_new(i,j,k) = max( 0., dzinterface(i,j,k) - dzinterface(i,j,k+1) )
1389 enddo
1390 endif
1391 else
1392 h_new(i,j,1:nki) = h(i,j,1:nki)
1393 if (cs%nk > gv%ke) h_new(i,j,nki+1:cs%nk) = 0.
1394 ! On land points, why are we keeping the original h rather than setting to zero? -AJA
1395 endif
1396 enddo
1397 enddo
1398
1399end subroutine calc_h_new_by_dz
1400
1401
1402!> Check that the total thickness of new and old grids are consistent
1403subroutine check_grid_column( nk, h, dzInterface, msg )
1404 integer, intent(in) :: nk !< Number of cells
1405 real, dimension(nk), intent(in) :: h !< Cell thicknesses [Z ~> m] or arbitrary units
1406 real, dimension(nk+1), intent(in) :: dzinterface !< Change in interface positions (same units as h), often [Z ~> m]
1407 character(len=*), intent(in) :: msg !< Message to append to errors
1408 ! Local variables
1409 integer :: k
1410 real :: eps ! A tiny relative thickness [nondim]
1411 real :: total_h_old ! The total thickness in the old column, in [Z ~> m] or arbitrary units
1412 real :: total_h_new ! The total thickness in the updated column, in [Z ~> m] or arbitrary units
1413 real :: h_new ! A thickness in the updated column, in [Z ~> m] or arbitrary units
1414
1415 eps =1. ; eps = epsilon(eps)
1416
1417 ! Total thickness of grid h
1418 total_h_old = 0.
1419 do k = 1,nk
1420 total_h_old = total_h_old + h(k)
1421 enddo
1422
1423 total_h_new = 0.
1424 do k = nk,1,-1
1425 h_new = h(k) + ( dzinterface(k) - dzinterface(k+1) ) ! New thickness
1426 if (h_new<0.) then
1427 write(0,*) 'k,h,hnew=',k,h(k),h_new
1428 write(0,*) 'dzI(k+1),dzI(k)=',dzinterface(k+1),dzinterface(k)
1429 call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
1430 'Negative layer thickness implied by re-gridding, '//trim(msg))
1431 endif
1432 total_h_new = total_h_new + h_new
1433
1434 enddo
1435
1436 ! Conservation by implied h_new
1437 if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps) then
1438 write(0,*) 'nk=',nk
1439 do k = 1,nk
1440 write(0,*) 'k,h,hnew=',k,h(k),h(k)+(dzinterface(k)-dzinterface(k+1))
1441 enddo
1442 write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
1443 write(0,*) 'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
1444 call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
1445 'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
1446 endif
1447
1448 ! Check that the top and bottom are intentionally moving
1449 if (dzinterface(1) /= 0.) call mom_error( fatal, &
1450 'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
1451 if (dzinterface(nk+1) /= 0.) call mom_error( fatal, &
1452 'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
1453
1454end subroutine check_grid_column
1455
1456!> Returns the change in interface position motion after filtering and
1457!! assuming the top and bottom interfaces do not move. The filtering is
1458!! a function of depth, and is applied as the integrated average filtering
1459!! over the trajectory of the interface. By design, this code can not give
1460!! tangled interfaces provided that z_old and z_new are not already tangled.
1461subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
1462 type(regridding_cs), intent(in) :: CS !< Regridding control structure
1463 integer, intent(in) :: nk !< Number of cells in source grid
1464 real, dimension(nk+1), intent(in) :: z_old !< Old grid position [H ~> m or kg m-2]
1465 real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position before filtering [H ~> m or kg m-2]
1466 real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions including
1467 !! the effects of filtering [H ~> m or kg m-2]
1468 ! Local variables
1469 real :: sgn ! The sign convention for downward [nondim].
1470 real :: dz_tgt ! The target grid movement of the unfiltered grid [H ~> m or kg m-2]
1471 real :: zr1 ! The old grid position of an interface relative to the surface [H ~> m or kg m-2]
1472 real :: z_old_k ! The corrected position of the old grid [H ~> m or kg m-2]
1473 real :: Aq ! A temporary variable related to the grid weights [nondim]
1474 real :: Bq ! A temporary variable used in the linear term in the quadratic expression for the
1475 ! filtered grid movement [H ~> m or kg m-2]
1476 real :: z0, dz0 ! Together these give the position of an interface relative to a reference hieght
1477 ! that may be adjusted for numerical accuracy in a solver [H ~> m or kg m-2]
1478 real :: F0 ! An estimated grid movement [H ~> m or kg m-2]
1479 real :: zs ! The depth at which the shallow filtering timescale applies [H ~> m or kg m-2]
1480 real :: zd ! The depth at which the deep filtering timescale applies [H ~> m or kg m-2]
1481 real :: dzwt ! The depth range over which the transition in the filtering timescale occurs [H ~> m or kg m-2]
1482 real :: Idzwt ! The Adcroft reciprocal of dzwt [H-1 ~> m-1 or m2 kg-1]
1483 real :: wtd ! The weight given to the new grid when time filtering [nondim]
1484 real :: Iwtd ! The inverse of wtd [nondim]
1485 real :: Int_zs ! A depth integral of the weights in [H ~> m or kg m-2]
1486 real :: Int_zd ! A depth integral of the weights in [H ~> m or kg m-2]
1487 real :: dInt_zs_zd ! The depth integral of the weights between the deep and shallow depths in [H ~> m or kg m-2]
1488! For debugging:
1489 real, dimension(nk+1) :: z_act ! The final grid positions after the filtered movement [H ~> m or kg m-2]
1490! real, dimension(nk+1) :: ddz_g_s, ddz_g_d
1491 logical :: debug = .false.
1492 integer :: k
1493
1494 if ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) < 0.0) then
1495 call mom_error(fatal, "filtered_grid_motion: z_old and z_new use different sign conventions.")
1496 elseif ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) == 0.0) then
1497 ! This is a massless column, so do nothing and return.
1498 do k=1,cs%nk+1 ; dz_g(k) = 0.0 ; enddo ; return
1499 elseif ((z_old(nk+1) - z_old(1)) + (z_new(cs%nk+1) - z_new(1)) > 0.0) then
1500 sgn = 1.0
1501 else
1502 sgn = -1.0
1503 endif
1504
1505 if (debug) then
1506 do k=2,cs%nk+1
1507 if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
1508 call mom_error(fatal, "filtered_grid_motion: z_new is tangled.")
1509 enddo
1510 do k=2,nk+1
1511 if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
1512 call mom_error(fatal, "filtered_grid_motion: z_old is tangled.")
1513 enddo
1514 ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0
1515 endif
1516
1517 zs = cs%depth_of_time_filter_shallow
1518 zd = cs%depth_of_time_filter_deep
1519 wtd = 1.0 - cs%old_grid_weight
1520 iwtd = 1.0 / wtd
1521
1522 dzwt = (zd - zs)
1523 idzwt = 0.0 ; if (abs(zd - zs) > 0.0) idzwt = 1.0 / (zd - zs)
1524 dint_zs_zd = 0.5*(1.0 + iwtd) * (zd - zs)
1525 aq = 0.5*(iwtd - 1.0)
1526
1527 dz_g(1) = 0.0
1528 z_old_k = z_old(1)
1529 do k = 2,cs%nk+1
1530 if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1531 ! zr1 is positive and increases with depth, and dz_tgt is positive downward.
1532 dz_tgt = sgn*(z_new(k) - z_old_k)
1533 zr1 = sgn*(z_old_k - z_old(1))
1534
1535 ! First, handle the two simple and common cases that do not pass through
1536 ! the adjustment rate transition zone.
1537 if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then
1538 dz_g(k) = sgn * wtd * dz_tgt
1539 elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then
1540 dz_g(k) = sgn * dz_tgt
1541 else
1542 ! Find the new value by inverting the equation
1543 ! integral(0 to dz_new) Iwt(z) dz = dz_tgt
1544 ! This is trivial where Iwt is a constant, and agrees with the two limits above.
1545
1546 ! Take test values at the transition points to figure out which segment
1547 ! the new value will be found in.
1548 if (zr1 >= zd) then
1549 int_zd = iwtd*(zd - zr1)
1550 int_zs = int_zd - dint_zs_zd
1551 elseif (zr1 <= zs) then
1552 int_zs = (zs - zr1)
1553 int_zd = dint_zs_zd + (zs - zr1)
1554 else
1555! Int_zd = (zd - zr1) * (Iwtd + 0.5*(1.0 - Iwtd) * (zd - zr1) / (zd - zs))
1556 int_zd = (zd - zr1) * (iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * idzwt
1557 int_zs = (zs - zr1) * (0.5*iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * idzwt
1558 ! It has been verified that Int_zs = Int_zd - dInt_zs_zd to within roundoff.
1559 endif
1560
1561 if (dz_tgt >= int_zd) then ! The new location is in the deep, slow region.
1562 dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - int_zd))
1563 elseif (dz_tgt <= int_zs) then ! The new location is in the shallow region.
1564 dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - int_zs))
1565 else ! We need to solve a quadratic equation for z_new.
1566 ! For accuracy, do the integral from the starting depth or the nearest
1567 ! edge of the transition region. The results with each choice are
1568 ! mathematically equivalent, but differ in roundoff, and this choice
1569 ! should minimize the likelihood of inadvertently overlapping interfaces.
1570 if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; f0 = dz_tgt - int_zs
1571 elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; f0 = dz_tgt - int_zd
1572 else ; dz0 = 0.0 ; z0 = zr1 ; f0 = dz_tgt ; endif
1573
1574 bq = (dzwt + 2.0*aq*(z0-zs))
1575 ! Solve the quadratic: Aq*(zn-z0)**2 + Bq*(zn-z0) - F0*dzwt = 0
1576 ! Note that b>=0, and the two terms in the standard form cancel for the right root.
1577 dz_g(k) = sgn * (dz0 + 2.0*f0*dzwt / (bq + sqrt(bq**2 + 4.0*aq*f0*dzwt) ))
1578
1579! if (debug) then
1580! dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0*Aq*(z0-zs))
1581! ddz_g_s(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
1582! dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0*Aq*(z0-zs))
1583! ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
1584!
1585! if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
1586! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).")
1587! 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)))) &
1588! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.")
1589! endif
1590 endif
1591
1592 endif
1593 enddo
1594 !dz_g(CS%nk+1) = 0.0
1595
1596 if (debug) then
1597 z_old_k = z_old(1)
1598 do k=1,cs%nk+1
1599 if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1600 z_act(k) = z_old_k + dz_g(k)
1601 enddo
1602 do k=2,cs%nk+1
1603 if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) &
1604 call mom_error(fatal, "filtered_grid_motion: z_output is tangled.")
1605 enddo
1606 endif
1607
1608end subroutine filtered_grid_motion
1609
1610!> Builds a z*-coordinate grid with partial steps (Adcroft and Campin, 2004).
1611!! z* is defined as
1612!! z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H .
1613subroutine build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, frac_shelf_h, zScale)
1614
1615 ! Arguments
1616 type(regridding_cs), intent(in) :: CS !< Regridding control structure
1617 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1618 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
1619 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1620 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
1621 !! relative to mean sea level or another locally
1622 !! valid reference height, converted to thickness
1623 !! units [H ~> m or kg m-2]
1624 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1625 !! [H ~> m or kg m-2].
1626 real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: frac_shelf_h !< Fractional
1627 !! ice shelf coverage [nondim].
1628 real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate
1629 !! resolution in Z to desired units for zInterface,
1630 !! usually Z_to_H in which case it is in
1631 !! units of [H Z-1 ~> nondim or kg m-3]
1632 ! Local variables
1633 real :: nominalDepth, minThickness, totalThickness ! Depths and thicknesses [H ~> m or kg m-2]
1634#ifdef __DO_SAFETY_CHECKS__
1635 real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2]
1636#endif
1637 real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
1638 real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
1639 integer :: i, j, k, nz
1640 logical :: ice_shelf
1641
1642 nz = gv%ke
1643 minthickness = cs%min_thickness
1644 ice_shelf = present(frac_shelf_h)
1645
1646 !$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, &
1647 !$OMP ice_shelf,minThickness,zScale,nom_depth_H) &
1648 !$OMP private(nominalDepth,totalThickness, &
1649#ifdef __DO_SAFETY_CHECKS__
1650 !$OMP dh, &
1651#endif
1652 !$OMP zNew,zOld)
1653 do j = g%jsc-1,g%jec+1
1654 do i = g%isc-1,g%iec+1
1655
1656 if (g%mask2dT(i,j)==0.) then
1657 dzinterface(i,j,:) = 0.
1658 cycle
1659 endif
1660
1661 ! Local depth (positive downward)
1662 nominaldepth = nom_depth_h(i,j)
1663
1664 ! Determine water column thickness
1665 totalthickness = 0.0
1666 do k = 1,nz
1667 totalthickness = totalthickness + h(i,j,k)
1668 enddo
1669
1670 ! if (GV%Boussinesq) then
1671 zold(nz+1) = - nominaldepth
1672 do k = nz,1,-1
1673 zold(k) = zold(k+1) + h(i,j,k)
1674 enddo
1675 ! else ! Work downward?
1676 ! endif
1677
1678 if (ice_shelf) then
1679 if (frac_shelf_h(i,j) > 0.) then ! under ice shelf
1680 call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, znew, &
1681 z_rigid_top=totalthickness-nominaldepth, &
1682 eta_orig=zold(1), zscale=zscale)
1683 else
1684 call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1685 znew, zscale=zscale)
1686 endif
1687 else
1688 call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1689 znew, zscale=zscale)
1690 endif
1691
1692 ! Calculate the final change in grid position after blending new and old grids
1693 if (cs%use_depth_based_time_filter .or. cs%old_grid_weight>0.) &
1694 call filtered_grid_motion(cs, nz, zold, znew, dzinterface(i,j,:))
1695
1696#ifdef __DO_SAFETY_CHECKS__
1697 dh = max(nominaldepth,totalthickness)
1698 if (abs(znew(1)-zold(1)) > (nz-1)*0.5*epsilon(dh)*dh) then
1699 write(0,*) 'min_thickness=',cs%min_thickness
1700 write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1701 write(0,*) 'dzInterface(1) = ', dzinterface(i,j,1), epsilon(dh), nz, cs%nk
1702 do k=1,min(nz,cs%nk)+1
1703 write(0,*) k,zold(k),znew(k)
1704 enddo
1705 do k=min(nz,cs%nk)+2,cs%nk+1
1706 write(0,*) k,zold(nz+1),znew(k)
1707 enddo
1708 do k=1,min(nz,cs%nk)
1709 write(0,*) k,h(i,j,k),znew(k)-znew(k+1),cs%coordinateResolution(k)
1710 enddo
1711 do k=min(nz,cs%nk)+1,cs%nk
1712 write(0,*) k, 0.0, znew(k)-znew(k+1), cs%coordinateResolution(k)
1713 enddo
1714 call mom_error( fatal, &
1715 'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
1716 endif
1717#endif
1718
1719 if (cs%use_adjust_interface_motion) call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
1720
1721 enddo
1722 enddo
1723
1724end subroutine build_zstar_grid
1725
1726!------------------------------------------------------------------------------
1727! Build sigma grid
1728!> This routine builds a grid based on terrain-following coordinates.
1729subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface )
1730!------------------------------------------------------------------------------
1731! This routine builds a grid based on terrain-following coordinates.
1732! The module parameter coordinateResolution(:) determines the resolution in
1733! sigma coordinate, dSigma(:). sigma-coordinates are defined by
1734! sigma = (eta-z)/(H+eta) s.t. sigma=0 at z=eta and sigma=1 at z=-H .
1735!------------------------------------------------------------------------------
1736
1737 ! Arguments
1738 type(regridding_cs), intent(in) :: CS !< Regridding control structure
1739 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1740 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
1741 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1742 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
1743 !! relative to mean sea level or another locally
1744 !! valid reference height, converted to thickness
1745 !! units [H ~> m or kg m-2]
1746 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1747 !! [H ~> m or kg m-2]
1748
1749 ! Local variables
1750 real :: nominalDepth ! The nominal depth of the sea-floor in thickness units [H ~> m or kg m-2]
1751 real :: totalThickness ! The total thickness of the water column [H ~> m or kg m-2]
1752#ifdef __DO_SAFETY_CHECKS__
1753 real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2]
1754#endif
1755 real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
1756 real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
1757 integer :: i, j, k, nz
1758
1759 nz = gv%ke
1760
1761 do i = g%isc-1,g%iec+1
1762 do j = g%jsc-1,g%jec+1
1763
1764 if (g%mask2dT(i,j)==0.) then
1765 dzinterface(i,j,:) = 0.
1766 cycle
1767 endif
1768
1769 ! Determine water column height
1770 totalthickness = 0.0
1771 do k = 1,nz
1772 totalthickness = totalthickness + h(i,j,k)
1773 enddo
1774
1775 ! In sigma coordinates, the bathymetric depth is only used as an arbitrary offset that
1776 ! cancels out when determining coordinate motion, so referencing the column postions to
1777 ! the surface is perfectly acceptable, but for preservation of previous answers the
1778 ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
1779 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1780 nominaldepth = nom_depth_h(i,j)
1781 else
1782 nominaldepth = totalthickness
1783 endif
1784
1785 call build_sigma_column(cs%sigma_CS, nominaldepth, totalthickness, znew)
1786
1787 ! Calculate the final change in grid position after blending new and old grids
1788 zold(nz+1) = -nominaldepth
1789 do k = nz,1,-1
1790 zold(k) = zold(k+1) + h(i,j,k)
1791 enddo
1792
1793 if (cs%use_depth_based_time_filter .or. cs%old_grid_weight>0.) &
1794 call filtered_grid_motion(cs, nz, zold, znew, dzinterface(i,j,:))
1795
1796#ifdef __DO_SAFETY_CHECKS__
1797 dh = max(nominaldepth,totalthickness)
1798 if (abs(znew(1)-zold(1)) > (cs%nk-1)*0.5*epsilon(dh)*dh) then
1799 write(0,*) 'min_thickness=',cs%min_thickness
1800 write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1801 write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz,cs%nk
1802 do k=1,min(nz,cs%nk)+1
1803 write(0,*) k,zold(k),znew(k)
1804 enddo
1805 do k=min(nz,cs%nk)+2,cs%nk+1
1806 write(0,*) k,zold(nz+1),znew(k)
1807 enddo
1808 do k=1,min(nz,cs%nk)
1809 write(0,*) k,h(i,j,k),znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k), &
1810 cs%coordinateResolution(k)
1811 enddo
1812 do k=min(nz,cs%nk)+1,cs%nk
1813 write(0,*) k,0.0,znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k), &
1814 cs%coordinateResolution(k)
1815 enddo
1816 call mom_error( fatal, &
1817 'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
1818 endif
1819 dzinterface(i,j,1) = 0.
1820 dzinterface(i,j,cs%nk+1) = 0.
1821#endif
1822
1823 enddo
1824 enddo
1825
1826end subroutine build_sigma_grid
1827
1828!------------------------------------------------------------------------------
1829! Build grid based on target interface densities
1830!------------------------------------------------------------------------------
1831!> This routine builds a new grid based on a given set of target interface densities.
1832subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, CS, frac_shelf_h )
1833!------------------------------------------------------------------------------
1834! This routine builds a new grid based on a given set of target interface
1835! densities (these target densities are computed by taking the mean value
1836! of given layer densities). The algorithm operates as follows within each
1837! column:
1838! 1. Given T & S within each layer, the layer densities are computed.
1839! 2. Based on these layer densities, a global density profile is reconstructed
1840! (this profile is monotonically increasing and may be discontinuous)
1841! 3. The new grid interfaces are determined based on the target interface
1842! densities.
1843! 4. T & S are remapped onto the new grid.
1844! 5. Return to step 1 until convergence or until the maximum number of
1845! iterations is reached, whichever comes first.
1846!------------------------------------------------------------------------------
1847
1848 ! Arguments
1849 type(regridding_cs), intent(in) :: CS !< Regridding control structure
1850 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1851 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1852 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1853 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1854 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
1855 !! relative to mean sea level or another locally
1856 !! valid reference height, converted to thickness
1857 !! units [H ~> m or kg m-2]
1858 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1859 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1860 !! [H ~> m or kg m-2]
1861 type(remapping_cs), intent(in) :: remapCS !< The remapping control structure
1862 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice
1863 !! shelf coverage [nondim]
1864 ! Local variables
1865 integer :: nz ! The number of layers in the input grid
1866 integer :: i, j, k
1867 real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2]
1868 real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
1869 real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
1870 real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
1871 real :: totalThickness ! Total thicknesses [H ~> m or kg m-2]
1872#ifdef __DO_SAFETY_CHECKS__
1873 real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2]
1874#endif
1875 logical :: ice_shelf
1876
1877 h_neglect = set_h_neglect(gv, cs%remap_answer_date, h_neglect_edge)
1878
1879 nz = gv%ke
1880 ice_shelf = present(frac_shelf_h)
1881
1882 if (.not.cs%target_density_set) call mom_error(fatal, "build_rho_grid: "//&
1883 "Target densities must be set before build_rho_grid is called.")
1884
1885 ! Build grid based on target interface densities
1886 do j = g%jsc-1,g%jec+1
1887 do i = g%isc-1,g%iec+1
1888
1889 if (g%mask2dT(i,j)==0.) then
1890 dzinterface(i,j,:) = 0.
1891 cycle
1892 endif
1893
1894 ! Determine total water column thickness
1895 totalthickness = 0.0
1896 do k=1,nz
1897 totalthickness = totalthickness + h(i,j,k)
1898 enddo
1899
1900 ! In rho coordinates, the bathymetric depth is only used as an arbitrary offset that
1901 ! cancels out when determining coordinate motion, so referencing the column postions to
1902 ! the surface is perfectly acceptable, but for preservation of previous answers the
1903 ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
1904 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1905 nominaldepth = nom_depth_h(i,j)
1906 else
1907 nominaldepth = totalthickness
1908 endif
1909
1910 ! Determine absolute interface positions
1911 zold(nz+1) = - nominaldepth
1912 do k = nz,1,-1
1913 zold(k) = zold(k+1) + h(i,j,k)
1914 enddo
1915
1916 if (ice_shelf) then
1917 call build_rho_column(cs%rho_CS, nz, nominaldepth, h(i,j,:), &
1918 tv%T(i,j,:), tv%S(i,j,:), tv%eqn_of_state, znew, &
1919 z_rigid_top=totalthickness - nominaldepth, eta_orig = zold(1), &
1920 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1921 else
1922 call build_rho_column(cs%rho_CS, nz, nominaldepth, h(i,j,:), &
1923 tv%T(i,j,:), tv%S(i,j,:), tv%eqn_of_state, znew, &
1924 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1925 endif
1926
1927 if (cs%integrate_downward_for_e) then
1928 zold(1) = 0.
1929 do k = 1,nz
1930 zold(k+1) = zold(k) - h(i,j,k)
1931 enddo
1932 else
1933 ! The rest of the model defines grids integrating up from the bottom
1934 zold(nz+1) = - nominaldepth
1935 do k = nz,1,-1
1936 zold(k) = zold(k+1) + h(i,j,k)
1937 enddo
1938 endif
1939
1940 ! Calculate the final change in grid position after blending new and old grids
1941 if (cs%use_depth_based_time_filter .or. cs%old_grid_weight>0.) &
1942 call filtered_grid_motion(cs, nz, zold, znew, dzinterface(i,j,:))
1943
1944#ifdef __DO_SAFETY_CHECKS__
1945 do k=2,cs%nk
1946 if (znew(k) > zold(1)) then
1947 write(0,*) 'zOld=',zold
1948 write(0,*) 'zNew=',znew
1949 call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
1950 'interior interface above surface!' )
1951 endif
1952 if (znew(k) > znew(k-1)) then
1953 write(0,*) 'zOld=',zold
1954 write(0,*) 'zNew=',znew
1955 call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
1956 'interior interfaces cross!' )
1957 endif
1958 enddo
1959
1960 totalthickness = 0.0
1961 do k = 1,nz
1962 totalthickness = totalthickness + h(i,j,k)
1963 enddo
1964
1965 dh = max(nominaldepth, totalthickness)
1966 if (abs(znew(1)-zold(1)) > (nz-1)*0.5*epsilon(dh)*dh) then
1967 write(0,*) 'min_thickness=',cs%min_thickness
1968 write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1969 write(0,*) 'zNew(1)-zOld(1) = ',znew(1)-zold(1),epsilon(dh),nz
1970 do k=1,min(nz,cs%nk)+1
1971 write(0,*) k,zold(k),znew(k)
1972 enddo
1973 do k=min(nz,cs%nk)+2,cs%nk+1
1974 write(0,*) k,zold(nz+1),znew(k)
1975 enddo
1976 do k=1,nz
1977 write(0,*) k,h(i,j,k),znew(k)-znew(k+1)
1978 enddo
1979 do k=min(nz,cs%nk)+1,cs%nk
1980 write(0,*) k, 0.0, znew(k)-znew(k+1), cs%coordinateResolution(k)
1981 enddo
1982 call mom_error( fatal, &
1983 'MOM_regridding, build_rho_grid: top surface has moved!!!' )
1984 endif
1985#endif
1986
1987 enddo ! end loop on i
1988 enddo ! end loop on j
1989
1990end subroutine build_rho_grid
1991
1992!> Builds a simple HyCOM-like grid with the deepest location of potential
1993!! density interpolated from the column profile and a clipping of depth for
1994!! each interface to a fixed z* or p* grid. This should probably be (optionally?)
1995!! changed to find the nearest location of the target density.
1996!! \remark { Based on Bleck, 2002: An ocean-ice general circulation model framed in
1997!! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88.
1998!! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 }
1999subroutine build_grid_hycom1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, remapCS, CS, &
2000 frac_shelf_h, zScale )
2001 type(ocean_grid_type), intent(in) :: G !< Grid structure
2002 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
2003 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2004 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2]
2005 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
2006 !! relative to mean sea level or another locally
2007 !! valid reference height, converted to thickness
2008 !! units [H ~> m or kg m-2]
2009 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
2010 type(remapping_cs), intent(in) :: remapCS !< The remapping control structure
2011 type(regridding_cs), intent(in) :: CS !< Regridding control structure
2012 real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
2013 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position
2014 !! in thickness units [H ~> m or kg m-2]
2015 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf
2016 !! coverage [nondim]
2017 real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate
2018 !! resolution in Z to desired units for zInterface,
2019 !! usually Z_to_H in which case it is in
2020 !! units of [H Z-1 ~> nondim or kg m-3]
2021
2022 ! Local variables
2023 real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
2024 real, dimension(SZK_(GV)) :: p_col ! Layer center pressure in the input column [R L2 T-2 ~> Pa]
2025 real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2]
2026 real, dimension(CS%nk+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
2027 real :: nominalDepth ! The nominal depth of the seafloor in thickness units [H ~> m or kg m-2]
2028 real :: h_neglect, h_neglect_edge ! Negligible thicknesses used for remapping [H ~> m or kg m-2]
2029 real :: z_top_col ! The nominal height of the sea surface or ice-ocean interface
2030 ! in thickness units [H ~> m or kg m-2]
2031 real :: totalThickness ! The total thickness of the water column [H ~> m or kg m-2]
2032 logical :: ice_shelf
2033 integer :: i, j, k, nki
2034
2035 h_neglect = set_h_neglect(gv, cs%remap_answer_date, h_neglect_edge)
2036
2037 if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_HyCOM1 : "//&
2038 "Target densities must be set before build_grid_HyCOM1 is called.")
2039
2040 nki = min(gv%ke, cs%nk)
2041 ice_shelf = present(frac_shelf_h)
2042
2043 ! Build grid based on target interface densities
2044 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
2045 if (g%mask2dT(i,j)>0.) then
2046
2047 nominaldepth = nom_depth_h(i,j)
2048
2049 if (ice_shelf) then
2050 totalthickness = 0.0
2051 do k=1,gv%ke
2052 totalthickness = totalthickness + h(i,j,k)
2053 enddo
2054 z_top_col = max(nominaldepth-totalthickness,0.0)
2055 else
2056 z_top_col = 0.0
2057 endif
2058
2059 z_col(1) = z_top_col ! Work downward rather than bottom up
2060 do k = 1, gv%ke
2061 z_col(k+1) = z_col(k) + h(i,j,k)
2062 p_col(k) = tv%P_Ref + cs%compressibility_fraction * &
2063 ( 0.5 * ( z_col(k) + z_col(k+1) ) * (gv%H_to_RZ*gv%g_Earth) - tv%P_Ref )
2064 enddo
2065
2066 call build_hycom1_column(cs%hycom_CS, remapcs, tv%eqn_of_state, gv%ke, i, j, nominaldepth, &
2067 h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, &
2068 z_col, z_col_new, zscale=zscale, &
2069 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
2070
2071 ! Calculate the final change in grid position after blending new and old grids
2072 if (cs%use_depth_based_time_filter .or. cs%old_grid_weight>0.) &
2073 call filtered_grid_motion( cs, gv%ke, z_col, z_col_new, dz_col )
2074
2075 ! This adjusts things robust to round-off errors
2076 dz_col(:) = -dz_col(:)
2077 if (cs%use_adjust_interface_motion) call adjust_interface_motion( cs, gv%ke, h(i,j,:), dz_col(:) )
2078
2079 dzinterface(i,j,1:nki+1) = dz_col(1:nki+1)
2080 if (nki<cs%nk) dzinterface(i,j,nki+2:cs%nk+1) = 0.
2081
2082 else ! on land
2083 dzinterface(i,j,:) = 0.
2084 endif ! mask2dT
2085 enddo ; enddo ! i,j
2086
2087 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
2088
2089end subroutine build_grid_hycom1
2090
2091!> This subroutine builds an adaptive grid that follows density surfaces where
2092!! possible, subject to constraints on the smoothness of interface heights.
2093subroutine build_grid_adaptive(G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, CS)
2094 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2095 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
2096 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2097 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
2098 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
2099 !! relative to mean sea level or another locally
2100 !! valid reference height, converted to thickness
2101 !! units [H ~> m or kg m-2]
2102 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
2103 !! thermodynamic variables
2104 type(regridding_cs), intent(in) :: CS !< Regridding control structure
2105 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
2106 !! [H ~> m or kg m-2]
2107 type(remapping_cs), intent(in) :: remapCS !< The remapping control structure
2108
2109 ! local variables
2110 integer :: i, j, k, nz ! indices and dimension lengths
2111 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt ! Temperature on interfaces [C ~> degC]
2112 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: sInt ! Salinity on interfaces [S ~> ppt]
2113 ! current interface positions and after tendency term is applied
2114 ! positive downward
2115 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt ! Interface depths [H ~> m or kg m-2]
2116 real, dimension(SZK_(GV)+1) :: zNext ! New interface depths [H ~> m or kg m-2]
2117
2118 nz = gv%ke
2119
2120 call assert((gv%ke == cs%nk), "build_grid_adaptive is only written to work "//&
2121 "with the same number of input and target layers.")
2122
2123 ! position surface at z = 0.
2124 zint(:,:,1) = 0.
2125
2126 ! work on interior interfaces
2127 do k = 2, nz ; do j = g%jsc-2,g%jec+2 ; do i = g%isc-2,g%iec+2
2128 tint(i,j,k) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
2129 sint(i,j,k) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
2130 zint(i,j,k) = zint(i,j,k-1) + h(i,j,k-1) ! zInt in [H]
2131 enddo ; enddo ; enddo
2132
2133 ! top and bottom temp/salt interfaces are just the layer
2134 ! average values
2135 tint(:,:,1) = tv%T(:,:,1) ; tint(:,:,nz+1) = tv%T(:,:,nz)
2136 sint(:,:,1) = tv%S(:,:,1) ; sint(:,:,nz+1) = tv%S(:,:,nz)
2137
2138 ! set the bottom interface depth
2139 zint(:,:,nz+1) = zint(:,:,nz) + h(:,:,nz)
2140
2141 ! calculate horizontal density derivatives (alpha/beta)
2142 ! between cells in a 5-point stencil, columnwise
2143 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
2144 if (g%mask2dT(i,j) < 0.5) then
2145 dzinterface(i,j,:) = 0. ! land point, don't move interfaces, and skip
2146 cycle
2147 endif
2148
2149 call build_adapt_column(cs%adapt_CS, g, gv, us, tv, i, j, zint, tint, sint, h, &
2150 nom_depth_h, znext)
2151
2152 if (cs%use_depth_based_time_filter .or. cs%old_grid_weight>0.) &
2153 call filtered_grid_motion(cs, nz, zint(i,j,:), znext, dzinterface(i,j,:))
2154 ! convert from depth to z
2155 do k = 1, nz+1 ; dzinterface(i,j,k) = -dzinterface(i,j,k) ; enddo
2156 if (cs%use_adjust_interface_motion) call adjust_interface_motion(cs, nz, h(i,j,:), dzinterface(i,j,:))
2157 enddo ; enddo
2158end subroutine build_grid_adaptive
2159
2160!> Adjust dz_Interface to ensure non-negative future thicknesses
2161subroutine adjust_interface_motion( CS, nk, h_old, dz_int )
2162 type(regridding_cs), intent(in) :: CS !< Regridding control structure
2163 integer, intent(in) :: nk !< Number of layers in h_old
2164 real, dimension(nk), intent(in) :: h_old !< Layer thicknesses on the old grid [H ~> m or kg m-2]
2165 real, dimension(CS%nk+1), intent(inout) :: dz_int !< Interface movements, adjusted to keep the thicknesses
2166 !! thicker than their minimum value [H ~> m or kg m-2]
2167 ! Local variables
2168 real :: h_new ! A layer thickness on the new grid [H ~> m or kg m-2]
2169 real :: eps ! A tiny relative thickness [nondim]
2170 real :: h_total ! The total thickness of the old grid [H ~> m or kg m-2]
2171 real :: h_err ! An error tolerance that use used to flag unacceptably large negative layer thicknesses
2172 ! that can not be explained by roundoff errors [H ~> m or kg m-2]
2173 integer :: k
2174
2175 eps = 1. ; eps = epsilon(eps)
2176
2177 h_total = 0. ; h_err = 0.
2178 do k = 1, min(cs%nk,nk)
2179 h_total = h_total + h_old(k)
2180 h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
2181 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2182 if (h_new < -3.0*h_err) then
2183 write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
2184 'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
2185 'h_new=',h_new,'h_err=',h_err
2186 call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
2187 'implied h<0 is larger than roundoff!')
2188 endif
2189 enddo
2190 if (cs%nk>nk) then
2191 do k = nk+1, cs%nk
2192 h_err = h_err + max( abs(dz_int(k)), abs(dz_int(k+1)) )*eps
2193 h_new = ( dz_int(k) - dz_int(k+1) )
2194 if (h_new < -3.0*h_err) then
2195 write(0,*) 'h<0 at k=',k,'h_old was empty',&
2196 'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
2197 'h_new=',h_new,'h_err=',h_err
2198 call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
2199 'implied h<0 is larger than roundoff!')
2200 endif
2201 enddo
2202 endif
2203 do k = min(cs%nk,nk),2,-1
2204 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2205 if (h_new<cs%min_thickness) &
2206 dz_int(k) = ( dz_int(k+1) - h_old(k) ) + cs%min_thickness ! Implies next h_new = min_thickness
2207 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2208 if (h_new<0.) &
2209 dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) ) ! Backup in case min_thickness==0
2210 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2211 if (h_new<0.) then
2212 write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
2213 'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
2214 'h_new=',h_new
2215 stop 'Still did not work!'
2216 call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
2217 'Repeated adjustment for roundoff h<0 failed!')
2218 endif
2219 enddo
2220 !if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'
2221
2222end subroutine adjust_interface_motion
2223
2224
2225!------------------------------------------------------------------------------
2226!> make sure all layers are at least as thick as the minimum thickness allowed
2227!! for regridding purposes by inflating thin layers. This breaks mass conservation
2228!! and adds mass to the model when there are excessively thin layers.
2229subroutine inflate_vanished_layers_old( CS, G, GV, h )
2230!------------------------------------------------------------------------------
2231! This routine is called when initializing the regridding options. The
2232! objective is to make sure all layers are at least as thick as the minimum
2233! thickness allowed for regridding purposes (this parameter is set in the
2234! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they
2235! are inflated up to the minimum thickness.
2236!------------------------------------------------------------------------------
2237
2238 ! Arguments
2239 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2240 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
2241 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
2242 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
2243
2244 ! Local variables
2245 integer :: i, j, k
2246 real :: htmp(gv%ke) ! A copy of a 1-d column of h [H ~> m or kg m-2]
2247
2248 do i = g%isc-1,g%iec+1
2249 do j = g%jsc-1,g%jec+1
2250
2251 ! Build grid for current column
2252 do k = 1,gv%ke
2253 htmp(k) = h(i,j,k)
2254 enddo
2255
2256 call old_inflate_layers_1d( cs%min_thickness, gv%ke, htmp )
2257
2258 ! Save modified grid
2259 do k = 1,gv%ke
2260 h(i,j,k) = htmp(k)
2261 enddo
2262
2263 enddo
2264 enddo
2265
2266end subroutine inflate_vanished_layers_old
2267
2268!------------------------------------------------------------------------------
2269!> Achieve convective adjustment by swapping layers
2270subroutine convective_adjustment(G, GV, h, tv)
2271 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
2272 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
2273 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2274 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
2275 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
2276!------------------------------------------------------------------------------
2277! Check each water column to see whether it is stratified. If not, sort the
2278! layers by successive swappings of water masses (bubble sort algorithm)
2279!------------------------------------------------------------------------------
2280
2281 ! Local variables
2282 real :: t0, t1 ! temperatures of two layers [C ~> degC]
2283 real :: s0, s1 ! salinities of two layers [S ~> ppt]
2284 real :: r0, r1 ! densities of two layers [R ~> kg m-3]
2285 real :: h0, h1 ! Layer thicknesses [H ~> m or kg m-2]
2286 real, dimension(GV%ke) :: p_col ! A column of zero pressures [R L2 T-2 ~> Pa]
2287 real, dimension(GV%ke) :: densities ! Densities in the column [R ~> kg m-3]
2288 logical :: stratified
2289 integer :: i, j, k
2290
2291 !### Doing convective adjustment based on potential densities with zero pressure seems
2292 ! questionable, although it does avoid ambiguous sorting. -RWH
2293 p_col(:) = 0.
2294
2295 ! Loop on columns
2296 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
2297
2298 ! Compute densities within current water column
2299 call calculate_density(tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state)
2300
2301 ! Repeat restratification until complete
2302 do
2303 stratified = .true.
2304 do k = 1,gv%ke-1
2305 ! Gather information of current and next cells
2306 t0 = tv%T(i,j,k) ; t1 = tv%T(i,j,k+1)
2307 s0 = tv%S(i,j,k) ; s1 = tv%S(i,j,k+1)
2308 r0 = densities(k) ; r1 = densities(k+1)
2309 h0 = h(i,j,k) ; h1 = h(i,j,k+1)
2310 ! If the density of the current cell is larger than the density
2311 ! below it, we swap the cells and recalculate the densitiies
2312 ! within the swapped cells
2313 if ( r0 > r1 ) then
2314 tv%T(i,j,k) = t1 ; tv%T(i,j,k+1) = t0
2315 tv%S(i,j,k) = s1 ; tv%S(i,j,k+1) = s0
2316 h(i,j,k) = h1 ; h(i,j,k+1) = h0
2317 ! Recompute densities at levels k and k+1
2318 call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_col(k), &
2319 densities(k), tv%eqn_of_state)
2320 call calculate_density(tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
2321 densities(k+1), tv%eqn_of_state )
2322 ! Because p_col is has uniform values, these calculate_density calls are equivalent to
2323 ! densities(k) = r1 ; densities(k+1) = r0
2324 stratified = .false.
2325 endif
2326 enddo ! k
2327
2328 if ( stratified ) exit
2329 enddo
2330
2331 enddo ; enddo ! i & j
2332
2333end subroutine convective_adjustment
2334
2335
2336!------------------------------------------------------------------------------
2337!> Return a uniform resolution vector in the units of the coordinate
2338function uniformresolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy)
2339!------------------------------------------------------------------------------
2340! Calculate a vector of uniform resolution in the units of the coordinate
2341!------------------------------------------------------------------------------
2342 ! Arguments
2343 integer, intent(in) :: nk !< Number of cells in source grid
2344 character(len=*), intent(in) :: coordmode !< A string indicating the coordinate mode.
2345 !! See the documentation for regrid_consts
2346 !! for the recognized values.
2347 real, intent(in) :: maxdepth !< The range of the grid values in some modes, in coordinate
2348 !! dependent units that might be [m] or [kg m-3] or [nondim]
2349 !! or something else.
2350 real, intent(in) :: rholight !< The minimum value of the grid in RHO mode [kg m-3]
2351 real, intent(in) :: rhoheavy !< The maximum value of the grid in RHO mode [kg m-3]
2352
2353 real :: uniformresolution(nk) !< The returned uniform resolution grid, in
2354 !! coordinate dependent units that might be [m] or
2355 !! [kg m-3] or [nondim] or something else.
2356
2357 ! Local variables
2358 integer :: scheme
2359
2360 scheme = coordinatemode(coordmode)
2361 select case ( scheme )
2362
2363 case ( regridding_zstar, regridding_hycom1, regridding_hybgen, &
2364 regridding_sigma_shelf_zstar, regridding_adaptive )
2365 uniformresolution(:) = maxdepth / real(nk)
2366
2367 case ( regridding_rho )
2368 uniformresolution(:) = (rhoheavy - rholight) / real(nk)
2369
2370 case ( regridding_sigma )
2371 uniformresolution(:) = 1. / real(nk)
2372
2373 case default
2374 call mom_error(fatal, "MOM_regridding, uniformResolution: "//&
2375 "Unrecognized choice for coordinate mode ("//trim(coordmode)//").")
2376
2377 end select ! type of grid
2378
2379end function uniformresolution
2380
2381!> Initialize the coordinate resolutions by calling the appropriate initialization
2382!! routine for the specified coordinate mode.
2383subroutine initcoord(CS, G, GV, US, coord_mode, param_file)
2384 type(regridding_cs), intent(inout) :: CS !< Regridding control structure
2385 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2386 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
2387 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2388 character(len=*), intent(in) :: coord_mode !< A string indicating the coordinate mode.
2389 !! See the documentation for regrid_consts
2390 !! for the recognized values.
2391 type(param_file_type), intent(in) :: param_file !< Parameter file
2392
2393 select case (coordinatemode(coord_mode))
2394 case (regridding_zstar)
2395 call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
2396 case (regridding_sigma_shelf_zstar)
2397 call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
2398 case (regridding_sigma)
2399 call init_coord_sigma(cs%sigma_CS, cs%nk, cs%coordinateResolution)
2400 case (regridding_rho)
2401 call init_coord_rho(cs%rho_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS)
2402 case (regridding_hycom1)
2403 if (allocated(cs%coordinateResolution_3d)) then
2404 call init_3d_coord_hycom(cs%hycom_CS, g, cs%nk, &
2405 cs%coordinateResolution_3d, cs%target_density_3d, &
2406 cs%interp_CS)
2407 else
2408 call init_coord_hycom(cs%hycom_CS, cs%nk, cs%coordinateResolution, cs%target_density, &
2409 cs%interp_CS)
2410 endif
2411 case (regridding_hybgen)
2412 call init_hybgen_regrid(cs%hybgen_CS, gv, us, param_file)
2413 case (regridding_adaptive)
2414 call init_coord_adapt(cs%adapt_CS, cs%nk, cs%coordinateResolution, gv%m_to_H, us%kg_m3_to_R)
2415 end select
2416end subroutine initcoord
2417
2418!------------------------------------------------------------------------------
2419!> Set the fixed resolution data
2420subroutine setcoordinateresolution( dz, CS, scale )
2421 real, dimension(:), intent(in) :: dz !< A vector of vertical grid spacings, in arbitrary coordinate
2422 !! dependent units, such as [m] for a z-coordinate or [kg m-3]
2423 !! for a density coordinate.
2424 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2425 real, optional, intent(in) :: scale !< A scaling factor converting dz to the internal represetation
2426 !! of coordRes, in various units that depend on the coordinate,
2427 !! such as [Z m-1 ~> 1] for a z-coordinate or [R m3 kg-1 ~> 1] for
2428 !! a density coordinate.
2429
2430 if (size(dz)/=cs%nk) call mom_error( fatal, &
2431 'setCoordinateResolution: inconsistent number of levels' )
2432
2433 if (present(scale)) then
2434 cs%coordinateResolution(:) = scale*dz(:)
2435 else
2436 cs%coordinateResolution(:) = dz(:)
2437 endif
2438
2439end subroutine setcoordinateresolution
2440
2441!> Set the 3d fixed resolution data
2442subroutine setcoordinateresolution_3d( dz_3d, CS, scale )
2443 real, dimension(:,:,:), intent(in) :: dz_3d !< A vector of vertical grid spacings, in arbitrary coordinate
2444 !! dependent units, such as [m] for a z-coordinate or [kg m-3]
2445 !! for a density coordinate.
2446 type(regridding_cs), intent(inout) :: CS !< Regridding control structure
2447 real, optional, intent(in) :: scale !< A scaling factor converting dz to coordRes [Z m-1 ~> 1]
2448
2449 if (.not.allocated(cs%coordinateResolution_3d)) &
2450 call mom_error(fatal,'setCoordinateResolution_3d: '//&
2451 'CS%coordinateResolution_3d not allocated.')
2452
2453 if (present(scale)) then
2454 cs%coordinateResolution_3d(:,:,:) = scale*dz_3d(:,:,:)
2455 else
2456 cs%coordinateResolution_3d(:,:,:) = dz_3d(:,:,:)
2457 endif
2458
2459end subroutine setcoordinateresolution_3d
2460
2461!> Set target densities based on the old Rlay variable
2462subroutine set_target_densities_from_gv( GV, US, CS )
2463 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
2464 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2465 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2466 ! Local variables
2467 integer :: k, nz
2468
2469 nz = cs%nk
2470 if (nz == 1) then ! Set a broad range of bounds. Regridding may not be meaningful in this case.
2471 cs%target_density(1) = 0.0
2472 cs%target_density(2) = 2.0*gv%Rlay(1)
2473 else
2474 cs%target_density(1) = (gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)))
2475 cs%target_density(nz+1) = (gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)))
2476 do k=2,nz
2477 cs%target_density(k) = cs%target_density(k-1) + cs%coordinateResolution(k)
2478 enddo
2479 endif
2480 cs%target_density_set = .true.
2481
2482end subroutine set_target_densities_from_gv
2483
2484!> Set target densities based on vector of interface values
2485subroutine set_target_densities_3d( CS, G, scale, rho_int_3d )
2486 type(regridding_cs), intent(inout) :: CS !< Regridding control structure
2487 type(ocean_grid_type),intent(in) :: G !< Ocean grid structure
2488 real, intent(in) :: scale !< A scaling factor converting densities [R m3 kg-1 ~> 1]
2489 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: rho_int_3d !< Interface densities [kg m-3]
2490
2491 if (.not.allocated(cs%target_density_3d)) &
2492 call mom_error(fatal,'set_target_densities_3d: '//&
2493 'CS%target_density_3d not allocated.')
2494
2495 cs%target_density_3d(:,:,:) = scale * rho_int_3d(:,:,:)
2496 cs%target_density_set = .true.
2497
2498end subroutine set_target_densities_3d
2499
2500!> Set target densities based on vector of interface values
2501subroutine set_target_densities( CS, rho_int )
2502 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2503 real, dimension(CS%nk+1), intent(in) :: rho_int !< Interface densities [R ~> kg m-3]
2504
2505 if (size(cs%target_density)/=size(rho_int)) then
2506 call mom_error(fatal, "set_target_densities inconsistent args!")
2507 endif
2508
2509 cs%target_density(:) = rho_int(:)
2510 cs%target_density_set = .true.
2511
2512end subroutine set_target_densities
2513
2514!> Set maximum interface depths based on a vector of input values.
2515subroutine set_regrid_max_depths( CS, max_depths, units_to_H )
2516 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2517 real, dimension(CS%nk+1), intent(in) :: max_depths !< Maximum interface depths, in arbitrary units, often [m]
2518 real, optional, intent(in) :: units_to_h !< A conversion factor for max_depths into H units,
2519 !! often in [H m-1 ~> 1 or kg m-3]
2520 ! Local variables
2521 real :: val_to_h ! A conversion factor from the units for max_depths into H units, often [H m-1 ~> 1 or kg m-3]
2522 ! if units_to_H is present, or [nondim] if it is absent.
2523 integer :: k
2524
2525 if (.not.allocated(cs%max_interface_depths)) allocate(cs%max_interface_depths(1:cs%nk+1))
2526
2527 val_to_h = 1.0 ; if (present(units_to_h)) val_to_h = units_to_h
2528 if (max_depths(cs%nk+1) < max_depths(1)) val_to_h = -1.0*val_to_h
2529
2530 ! Check for sign reversals in the depths.
2531 if (max_depths(cs%nk+1) < max_depths(1)) then
2532 do k=1,cs%nk
2533 if (max_depths(k+1) > max_depths(k)) &
2534 call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths!")
2535 enddo
2536 else
2537 do k=1,cs%nk
2538 if (max_depths(k+1) < max_depths(k)) &
2539 call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths.")
2540 enddo
2541 endif
2542
2543 do k=1,cs%nk+1
2544 cs%max_interface_depths(k) = val_to_h * max_depths(k)
2545 enddo
2546
2547 ! set max depths for coordinate
2548 select case (cs%regridding_scheme)
2549 case (regridding_hycom1)
2550 call set_hycom_params(cs%hycom_CS, max_interface_depths=cs%max_interface_depths)
2551 end select
2552end subroutine set_regrid_max_depths
2553
2554!> Set maximum layer thicknesses based on a vector of input values.
2555subroutine set_regrid_max_thickness( CS, max_h, units_to_H )
2556 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2557 real, dimension(CS%nk+1), intent(in) :: max_h !< Maximum layer thicknesses, in arbitrary units, often [m]
2558 real, optional, intent(in) :: units_to_h !< A conversion factor for max_h into H units,
2559 !! often [H m-1 ~> 1 or kg m-3]
2560 ! Local variables
2561 real :: val_to_h ! A conversion factor from the units for max_h into H units, often [H m-1 ~> 1 or kg m-3]
2562 ! if units_to_H is present, or [nondim] if it is absent.
2563 integer :: k
2564
2565 if (.not.allocated(cs%max_layer_thickness)) allocate(cs%max_layer_thickness(1:cs%nk))
2566
2567 val_to_h = 1.0 ; if (present( units_to_h)) val_to_h = units_to_h
2568
2569 do k=1,cs%nk
2570 cs%max_layer_thickness(k) = val_to_h * max_h(k)
2571 enddo
2572
2573 ! set max thickness for coordinate
2574 select case (cs%regridding_scheme)
2575 case (regridding_hycom1)
2576 call set_hycom_params(cs%hycom_CS, max_layer_thickness=cs%max_layer_thickness)
2577 end select
2578end subroutine set_regrid_max_thickness
2579
2580
2581!> Write the vertical coordinate information into a file.
2582!! This subroutine writes out a file containing any available data related
2583!! to the vertical grid used by the MOM ocean model when in ALE mode.
2584subroutine write_regrid_file( CS, GV, filepath )
2585 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2586 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
2587 character(len=*), intent(in) :: filepath !< The full path to the file to write
2588
2589 type(vardesc) :: vars(2)
2590 type(mom_field) :: fields(2)
2591 type(mom_netcdf_file) :: io_handle ! The I/O handle of the fileset
2592 real :: ds(gv%ke), dsi(gv%ke+1) ! The labeling layer and interface coordinates for output
2593 ! in axes in files, in coordinate-dependent units that can
2594 ! be obtained from getCoordinateUnits [various]
2595
2596 if (cs%regridding_scheme == regridding_hybgen) then
2597 call write_hybgen_coord_file(gv, cs%hybgen_CS, filepath)
2598 return
2599 endif
2600
2601 ds(:) = cs%coord_scale * cs%coordinateResolution(:)
2602 dsi(1) = 0.5*ds(1)
2603 dsi(2:gv%ke) = 0.5*( ds(1:gv%ke-1) + ds(2:gv%ke) )
2604 dsi(gv%ke+1) = 0.5*ds(gv%ke)
2605
2606 vars(1) = var_desc('ds', getcoordinateunits( cs ), &
2607 'Layer Coordinate Thickness', '1', 'L', '1')
2608 vars(2) = var_desc('ds_interface', getcoordinateunits( cs ), &
2609 'Layer Center Coordinate Separation', '1', 'i', '1')
2610
2611 call create_mom_file(io_handle, trim(filepath), vars, 2, fields, &
2612 single_file, gv=gv)
2613 call mom_write_field(io_handle, fields(1), ds)
2614 call mom_write_field(io_handle, fields(2), dsi)
2615 call io_handle%close()
2616
2617end subroutine write_regrid_file
2618
2619!> Set appropriate values for the negligible thicknesses used for remapping based on an answer date.
2620function set_h_neglect(GV, remap_answer_date, h_neglect_edge) result(h_neglect)
2621 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
2622 integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use
2623 !! for remapping. Values below 20190101 recover the
2624 !! remapping answers from 2018. Higher values use more
2625 !! robust forms of the same remapping algorithms.
2626 real, intent(out) :: h_neglect_edge !< A negligibly small thickness used in
2627 !! remapping edge value calculations [H ~> m or kg m-2]
2628 real :: h_neglect !< A negligibly small thickness used in
2629 !! remapping cell reconstructions [H ~> m or kg m-2]
2630
2631 if (remap_answer_date >= 20190101) then
2632 h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
2633 elseif (gv%Boussinesq) then
2634 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
2635 else
2636 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
2637 endif
2638end function set_h_neglect
2639
2640!> Set appropriate values for the negligible vertical layer extents used for remapping based on an answer date.
2641function set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) result(dz_neglect)
2642 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
2643 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2644 integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use
2645 !! for remapping. Values below 20190101 recover the
2646 !! remapping answers from 2018. Higher values use more
2647 !! robust forms of the same remapping algorithms.
2648 real, intent(out) :: dz_neglect_edge !< A negligibly small vertical layer extent
2649 !! used in remapping edge value calculations [Z ~> m]
2650 real :: dz_neglect !< A negligibly small vertical layer extent
2651 !! used in remapping cell reconstructions [Z ~> m]
2652
2653 if (remap_answer_date >= 20190101) then
2654 dz_neglect = gv%dZ_subroundoff ; dz_neglect_edge = gv%dZ_subroundoff
2655 elseif (gv%Boussinesq) then
2656 dz_neglect = us%m_to_Z*1.0e-30 ; dz_neglect_edge = us%m_to_Z*1.0e-10
2657 else
2658 dz_neglect = gv%kg_m2_to_H * (gv%H_to_m*us%m_to_Z) * 1.0e-30
2659 dz_neglect_edge = gv%kg_m2_to_H * (gv%H_to_m*us%m_to_Z) * 1.0e-10
2660 endif
2661end function set_dz_neglect
2662
2663!------------------------------------------------------------------------------
2664!> Query the fixed resolution data
2665function getcoordinateresolution( CS, undo_scaling )
2666 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2667 logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2668 !! rescaling of the resolution data.
2669 real, dimension(CS%nk) :: getcoordinateresolution !< The resolution or delta of the target coordinate,
2670 !! in units that depend on the coordinate [various]
2671
2672 logical :: unscale
2673 unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2674
2675 if (unscale) then
2676 getcoordinateresolution(:) = cs%coord_scale * cs%coordinateResolution(:)
2677 else
2678 getcoordinateresolution(:) = cs%coordinateResolution(:)
2679 endif
2680
2681end function getcoordinateresolution
2682
2683!> Query the target coordinate interface positions
2684function getcoordinateinterfaces( CS, undo_scaling )
2685 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2686 logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2687 !! rescaling of the resolution data.
2688 real, dimension(CS%nk+1) :: getcoordinateinterfaces !< Interface positions in target coordinate,
2689 !! in units that depend on the coordinate [various]
2690
2691 integer :: k
2692 logical :: unscale
2693 unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2694
2695 ! When using a coordinate with target densities, we need to get the actual
2696 ! densities, rather than computing the interfaces based on resolution
2697 if (cs%regridding_scheme == regridding_rho) then
2698 if (.not. cs%target_density_set) &
2699 call mom_error(fatal, 'MOM_regridding, getCoordinateInterfaces: '//&
2700 'target densities not set!')
2701
2702 if (unscale) then
2703 getcoordinateinterfaces(:) = cs%coord_scale * cs%target_density(:)
2704 else
2705 getcoordinateinterfaces(:) = cs%target_density(:)
2706 endif
2707 else
2708 if (unscale) then
2710 do k = 1, cs%nk
2712 cs%coord_scale * cs%coordinateResolution(k)
2713 enddo
2714 else
2716 do k = 1, cs%nk
2718 cs%coordinateResolution(k)
2719 enddo
2720 endif
2721 ! The following line has an "abs()" to allow ferret users to reference
2722 ! data by index. It is a temporary work around... :( -AJA
2724 endif
2725
2726end function getcoordinateinterfaces
2727
2728!------------------------------------------------------------------------------
2729!> Query the target coordinate units
2730function getcoordinateunits( CS )
2731 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2732 character(len=20) :: getcoordinateunits
2733
2734 select case ( cs%regridding_scheme )
2735 case ( regridding_zstar, regridding_hycom1, regridding_hybgen, &
2736 regridding_adaptive )
2737 getcoordinateunits = 'meter'
2738 case ( regridding_sigma_shelf_zstar )
2739 getcoordinateunits = 'meter/fraction'
2740 case ( regridding_sigma )
2741 getcoordinateunits = 'fraction'
2742 case ( regridding_rho )
2743 getcoordinateunits = 'kg/m3'
2744 case ( regridding_arbitrary )
2745 getcoordinateunits = 'unknown'
2746 case default
2747 call mom_error(fatal,'MOM_regridding, getCoordinateUnits: '//&
2748 'Unknown regridding scheme selected!')
2749 end select ! type of grid
2750
2751end function getcoordinateunits
2752
2753!------------------------------------------------------------------------------
2754!> Query the short name of the coordinate
2755function getcoordinateshortname( CS )
2756 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2757 character(len=20) :: getcoordinateshortname
2758
2759 select case ( cs%regridding_scheme )
2760 case ( regridding_zstar )
2761 !getCoordinateShortName = 'z*'
2762 ! The following line is a temporary work around... :( -AJA
2763 getcoordinateshortname = 'pseudo-depth, -z*'
2764 case ( regridding_sigma_shelf_zstar )
2765 getcoordinateshortname = 'pseudo-depth, -z*/sigma'
2766 case ( regridding_sigma )
2767 getcoordinateshortname = 'sigma'
2768 case ( regridding_rho )
2770 case ( regridding_arbitrary )
2771 getcoordinateshortname = 'coordinate'
2772 case ( regridding_hycom1 )
2773 getcoordinateshortname = 'z-rho'
2774 case ( regridding_hybgen )
2775 getcoordinateshortname = 'hybrid'
2776 case ( regridding_adaptive )
2777 getcoordinateshortname = 'adaptive'
2778 case default
2779 call mom_error(fatal,'MOM_regridding, getCoordinateShortName: '//&
2780 'Unknown regridding scheme selected!')
2781 end select ! type of grid
2782
2783end function getcoordinateshortname
2784
2785!> Can be used to set any of the parameters for MOM_regridding.
2786subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, &
2787 use_depth_based_time_filter, depth_of_time_filter_shallow, depth_of_time_filter_deep, &
2788 interp_scheme, use_adjust_interface_motion, compress_fraction, ref_pressure, &
2789 integrate_downward_for_e, remap_answers_2018, remap_answer_date, regrid_answer_date, &
2790 adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, &
2791 adaptAlpha, adaptDoMin, adaptDrho0)
2792 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2793 logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
2794 real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the
2795 !! new grid [H ~> m or kg m-2]
2796 real, optional, intent(in) :: old_grid_weight !< Weight given to old coordinate when time-filtering grid [nondim]
2797 logical, optional, intent(in) :: use_depth_based_time_filter !< Allow depth-based time filtering
2798 real, optional, intent(in) :: depth_of_time_filter_shallow !< Depth to start cubic [H ~> m or kg m-2]
2799 real, optional, intent(in) :: depth_of_time_filter_deep !< Depth to end cubic [H ~> m or kg m-2]
2800 character(len=*), optional, intent(in) :: interp_scheme !< Interpolation method for state-dependent coordinates
2801 logical, optional, intent(in) :: use_adjust_interface_motion !< Call adjust_interface_motion()
2802 real, optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density [nondim]
2803 real, optional, intent(in) :: ref_pressure !< The reference pressure for density-dependent
2804 !! coordinates [R L2 T-2 ~> Pa]
2805 logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward
2806 !! from the top.
2807 logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
2808 !! that recover the remapping answers from 2018. Otherwise
2809 !! use more robust but mathematically equivalent expressions.
2810 integer, optional, intent(in) :: remap_answer_date !< The vintage of the expressions to use for remapping
2811 integer, optional, intent(in) :: regrid_answer_date !< The vintage of the expressions to use for regridding
2812 real, optional, intent(in) :: adapttimeratio !< Ratio of the ALE timestep to the grid timescale [nondim].
2813 real, optional, intent(in) :: adaptzoom !< Depth of near-surface zooming region [H ~> m or kg m-2].
2814 real, optional, intent(in) :: adaptzoomcoeff !< Coefficient of near-surface zooming diffusivity [nondim].
2815 real, optional, intent(in) :: adaptbuoycoeff !< Coefficient of buoyancy diffusivity [nondim].
2816 real, optional, intent(in) :: adaptalpha !< Scaling factor on optimization tendency [nondim].
2817 logical, optional, intent(in) :: adaptdomin !< If true, make a HyCOM-like mixed layer by
2818 !! preventing interfaces from being shallower than
2819 !! the depths specified by the regridding coordinate.
2820 real, optional, intent(in) :: adaptdrho0 !< Reference density difference for stratification-dependent
2821 !! diffusion. [R ~> kg m-3]
2822
2823 if (present(interp_scheme)) call set_interp_scheme(cs%interp_CS, interp_scheme)
2824 if (present(boundary_extrapolation)) call set_interp_extrap(cs%interp_CS, boundary_extrapolation)
2825 if (present(regrid_answer_date)) call set_interp_answer_date(cs%interp_CS, regrid_answer_date)
2826
2827 if (present(old_grid_weight)) then
2828 if (old_grid_weight<0. .or. old_grid_weight>1.) &
2829 call mom_error(fatal,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
2830 cs%old_grid_weight = old_grid_weight
2831 endif
2832 if (present(use_depth_based_time_filter)) cs%use_depth_based_time_filter = &
2833 use_depth_based_time_filter
2834 if (present(depth_of_time_filter_shallow)) cs%depth_of_time_filter_shallow = &
2835 depth_of_time_filter_shallow
2836 if (present(depth_of_time_filter_deep)) cs%depth_of_time_filter_deep = &
2837 depth_of_time_filter_deep
2838 if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then
2839 if (cs%depth_of_time_filter_deep<cs%depth_of_time_filter_shallow) call mom_error(fatal, &
2840 'MOM_regridding, '//&
2841 'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
2842 endif
2843
2844 if (present(min_thickness)) cs%min_thickness = min_thickness
2845 if (present(use_adjust_interface_motion)) cs%use_adjust_interface_motion = use_adjust_interface_motion
2846 if (present(compress_fraction)) cs%compressibility_fraction = compress_fraction
2847 if (present(ref_pressure)) cs%ref_pressure = ref_pressure
2848 if (present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
2849 if (present(remap_answers_2018)) then
2850 if (remap_answers_2018) then
2851 cs%remap_answer_date = 20181231
2852 else
2853 cs%remap_answer_date = 20190101
2854 endif
2855 endif
2856 if (present(remap_answer_date)) cs%remap_answer_date = remap_answer_date
2857
2858 select case (cs%regridding_scheme)
2859 case (regridding_zstar)
2860 if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2861 case (regridding_sigma_shelf_zstar)
2862 if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2863 case (regridding_sigma)
2864 if (present(min_thickness)) call set_sigma_params(cs%sigma_CS, min_thickness=min_thickness)
2865 case (regridding_rho)
2866 if (present(min_thickness)) call set_rho_params(cs%rho_CS, min_thickness=min_thickness)
2867 if (present(ref_pressure)) call set_rho_params(cs%rho_CS, ref_pressure=ref_pressure)
2868 if (present(integrate_downward_for_e)) &
2869 call set_rho_params(cs%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
2870 if (associated(cs%rho_CS) .and. (present(interp_scheme) .or. &
2871 present(boundary_extrapolation))) &
2872 call set_rho_params(cs%rho_CS, interp_cs=cs%interp_CS)
2873 case (regridding_hycom1)
2874 if (associated(cs%hycom_CS) .and. (present(interp_scheme) .or. &
2875 present(boundary_extrapolation))) &
2876 call set_hycom_params(cs%hycom_CS, interp_cs=cs%interp_CS)
2877 case (regridding_hybgen)
2878 ! Do nothing for now.
2879 case (regridding_adaptive)
2880 if (present(adapttimeratio)) call set_adapt_params(cs%adapt_CS, adapttimeratio=adapttimeratio)
2881 if (present(adaptzoom)) call set_adapt_params(cs%adapt_CS, adaptzoom=adaptzoom)
2882 if (present(adaptzoomcoeff)) call set_adapt_params(cs%adapt_CS, adaptzoomcoeff=adaptzoomcoeff)
2883 if (present(adaptbuoycoeff)) call set_adapt_params(cs%adapt_CS, adaptbuoycoeff=adaptbuoycoeff)
2884 if (present(adaptalpha)) call set_adapt_params(cs%adapt_CS, adaptalpha=adaptalpha)
2885 if (present(adaptdomin)) call set_adapt_params(cs%adapt_CS, adaptdomin=adaptdomin)
2886 if (present(adaptdrho0)) call set_adapt_params(cs%adapt_CS, adaptdrho0=adaptdrho0)
2887 end select
2888
2889end subroutine set_regrid_params
2890
2891!> Returns the number of levels/layers in the regridding control structure
2892integer function get_regrid_size(CS)
2893 type(regridding_cs), intent(inout) :: cs !< Regridding control structure
2894
2895 get_regrid_size = cs%nk
2896
2897end function get_regrid_size
2898
2899!> This returns a copy of the zlike_CS stored in the regridding control structure.
2900function get_zlike_cs(CS)
2901 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2902 type(zlike_cs) :: get_zlike_cs
2903
2904 get_zlike_cs = cs%zlike_CS
2905end function get_zlike_cs
2906
2907!> This returns a copy of the sigma_CS stored in the regridding control structure.
2908function get_sigma_cs(CS)
2909 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2910 type(sigma_cs) :: get_sigma_cs
2911
2912 get_sigma_cs = cs%sigma_CS
2913end function get_sigma_cs
2914
2915!> This returns a copy of the rho_CS stored in the regridding control structure.
2916function get_rho_cs(CS)
2917 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2918 type(rho_cs) :: get_rho_cs
2919
2920 get_rho_cs = cs%rho_CS
2921end function get_rho_cs
2922
2923!------------------------------------------------------------------------------
2924!> Return coordinate-derived thicknesses for fixed coordinate systems
2925function getstaticthickness( CS, SSH, depth )
2926 type(regridding_cs), intent(in) :: cs !< Regridding control structure
2927 real, intent(in) :: ssh !< The sea surface height, in the same units as depth, often [Z ~> m]
2928 real, intent(in) :: depth !< The maximum depth of the grid, often [Z ~> m]
2929 real, dimension(CS%nk) :: getstaticthickness !< The returned thicknesses in the units of
2930 !! depth, often [Z ~> m]
2931 ! Local
2932 integer :: k
2933 real :: z, dz ! Vertical positions and grid spacing [Z ~> m]
2934
2935 select case ( cs%regridding_scheme )
2936 case ( regridding_zstar, regridding_sigma_shelf_zstar, regridding_hycom1, &
2937 regridding_hybgen, regridding_adaptive )
2938 if (depth>0.) then
2939 z = ssh
2940 do k = 1, cs%nk
2941 dz = cs%coordinateResolution(k) * ( 1. + ssh/depth ) ! Nominal dz*
2942 dz = max(dz, 0.) ! Avoid negative incase ssh=-depth
2943 dz = min(dz, depth - z) ! Clip if below topography
2944 z = z + dz ! Bottom of layer
2945 getstaticthickness(k) = dz
2946 enddo
2947 else
2948 getstaticthickness(:) = 0. ! On land ...
2949 endif
2950 case ( regridding_sigma )
2951 getstaticthickness(:) = cs%coordinateResolution(:) * ( depth + ssh )
2952 case ( regridding_rho )
2953 getstaticthickness(:) = 0. ! Not applicable
2954 case ( regridding_arbitrary )
2955 getstaticthickness(:) = 0. ! Not applicable
2956 case default
2957 call mom_error(fatal,'MOM_regridding, getStaticThickness: '//&
2958 'Unknown regridding scheme selected!')
2959 end select ! type of grid
2960
2961end function getstaticthickness
2962
2963!> Parses a string and generates a dz(:) profile that goes like k**power.
2964subroutine dz_function1( string, dz )
2965 character(len=*), intent(in) :: string !< String with list of parameters in form
2966 !! dz_min, H_total, power, precision
2967 real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses [m] or other units
2968 ! Local variables
2969 integer :: nk, k
2970 real :: dz_min ! minimum grid spacing [m] or other units
2971 real :: power ! A power to raise the relative position in index space [nondim]
2972 real :: prec ! The precision with which positions are returned [m] or other units
2973 real :: H_total ! The sum of the nominal thicknesses [m] or other units
2974
2975 nk = size(dz) ! Number of cells
2976 prec = -1024.
2977 read( string, *) dz_min, h_total, power, prec
2978 if (prec == -1024.) call mom_error(fatal,"dz_function1: "// &
2979 "Problem reading FNC1: string ="//trim(string))
2980 ! Create profile of ( dz - dz_min )
2981 do k = 1, nk
2982 dz(k) = (real(k-1)/real(nk-1))**power
2983 enddo
2984 dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2985 dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2986 dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2987 dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2988 dz(nk) = dz(nk) + ( h_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer
2989 dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2990 dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min
2991
2992end subroutine dz_function1
2993
2994!> Construct the name of a parameter for a specific coordinate based on param_prefix and param_suffix. For the main,
2995!! prognostic coordinate this will simply return the parameter name (e.g. P_REF)
2996function create_coord_param(param_prefix, param_name, param_suffix) result(coord_param)
2997 character(len=*) :: param_name !< The base name of the parameter (e.g. the one used for the main coordinate)
2998 character(len=*) :: param_prefix !< String to prefix to parameter names.
2999 character(len=*) :: param_suffix !< String to append to parameter names.
3000 character(len=MAX_PARAM_LENGTH) :: coord_param !< Parameter name prepended by param_prefix
3001 !! and appended with param_suffix
3002 integer :: out_length
3003
3004 if (len_trim(param_prefix) + len_trim(param_suffix) == 0) then
3005 coord_param = param_name
3006 else
3007 ! Note the +2 is because of two underscores
3008 out_length = len_trim(param_name)+len_trim(param_prefix)+len_trim(param_suffix)+2
3009 if (out_length > max_param_length) then
3010 call mom_error(fatal,"Coordinate parameter is too long; increase MAX_PARAM_LENGTH")
3011 endif
3012 coord_param = trim(param_prefix)//"_"//trim(param_name)//"_"//trim(param_suffix)
3013 endif
3014
3015end function create_coord_param
3016
3017!> Parses a string and generates a rho_target(:) profile with refined resolution downward
3018!! and returns the number of levels
3019integer function rho_function1( string, rho_target )
3020 character(len=*), intent(in) :: string !< String with list of parameters in form
3021 !! dz_min, H_total, power, precision
3022 real, dimension(:), allocatable, intent(inout) :: rho_target !< Profile of interface densities [kg m-3]
3023 ! Local variables
3024 integer :: nki, k, nk
3025 real :: dx ! Fractional distance from interface nki [nondim]
3026 real :: ddx ! Change in dx between interfaces [nondim]
3027 real :: rho_1, rho_2 ! Density of the top two layers in a profile [kg m-3]
3028 real :: rho_3 ! Density in the third layer, below which the density increase linearly
3029 ! in subsequent layers [kg m-3]
3030 real :: drho ! Change in density over the linear region [kg m-3]
3031 real :: rho_4 ! The densest density in this profile [kg m-3], which might be very large.
3032 real :: drho_min ! A minimal fractional density difference [nondim]?
3033
3034 read( string, *) nk, rho_1, rho_2, rho_3, drho, rho_4, drho_min
3035 allocate(rho_target(nk+1))
3036 nki = nk + 1 - 4 ! Number of interfaces minus 4 specified values
3037 rho_target(1) = rho_1
3038 rho_target(2) = rho_2
3039 dx = 0.
3040 do k = 0, nki
3041 ddx = max( drho_min, real(nki-k)/real(nki*nki) )
3042 dx = dx + ddx
3043 rho_target(3+k) = rho_3 + (2. * drho) * dx
3044 enddo
3045 rho_target(nki+4) = rho_4
3046
3047 rho_function1 = nk
3048
3049end function rho_function1
3050
3051!> \namespace mom_regridding
3052!!
3053!! A vertical grid is defined solely by the cell thicknesses, \f$h\f$.
3054!! Most calculations in this module start with the coordinate at the bottom
3055!! of the column set to -depth, and use a increasing value of coordinate with
3056!! decreasing k. This is consistent with the rest of MOM6 that uses position,
3057!! \f$z\f$ which is a negative quantity for most of the ocean.
3058!!
3059!! A change in grid is define through a change in position of the interfaces:
3060!! \f[
3061!! z^n_{k+1/2} = z^{n-1}_{k+1/2} + \Delta z_{k+1/2}
3062!! \f]
3063!! with the positive upward coordinate convention
3064!! \f[
3065!! z_{k-1/2} = z_{k+1/2} + h_k
3066!! \f]
3067!! so that
3068!! \f[
3069!! h^n_k = h^{n-1}_k + ( \Delta z_{k-1/2} - \Delta z_{k+1/2} )
3070!! \f]
3071!!
3072!! Original date of creation: 2008.06.09 by L. White
3073
3074end module mom_regridding