MOM_hybgen_regrid.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> This module contains the hybgen regridding routines from HYCOM, with minor
6!! modifications to follow the MOM6 coding conventions
8
9use mom_eos, only : eos_type, calculate_density
10use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, assert
11use mom_file_parser, only : get_param, param_file_type, log_param
12use mom_io, only : create_mom_file, file_exists
13use mom_io, only : mom_infra_file, mom_field
14use mom_io, only : mom_read_data, mom_write_field, vardesc, var_desc, single_file
17use mom_variables, only : ocean_grid_type, thermo_var_ptrs
19
20implicit none ; private
21
22#include <MOM_memory.h>
23
24!> Control structure containing required parameters for the hybgen coordinate generator
25type, public :: hybgen_regrid_cs ; private
26
27 real :: min_thickness !< Minimum thickness allowed for layers [H ~> m or kg m-2]
28
29 integer :: nk !< Number of layers on the target grid
30
31 !> Reference pressure for density calculations [R L2 T-2 ~> Pa]
32 real :: ref_pressure
33
34 !> Hybgen uses PCM if layer is within hybiso of target density [R ~> kg m-3]
35 real :: hybiso
36 !> Number of sigma levels used by HYBGEN
37 integer :: nsigma
38
39 real :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2]
40 real :: qhybrlx !< Fractional relaxation within a regridding step [nondim]
41
42 real, allocatable, dimension(:) :: &
43 dp0k, & !< minimum deep z-layer separation [H ~> m or kg m-2]
44 ds0k !< minimum shallow z-layer separation [H ~> m or kg m-2]
45
46 real :: coord_scale = 1.0 !< A scaling factor to restores the depth coordinates to
47 !! values in m [m H-1 ~> 1 or m3 kg-1]
48 real :: rho_coord_scale = 1.0 !< A scaling factor to restores the denesity coordinates to
49 !! values in kg m-3 [kg m-3 R-1 ~> 1]
50
51 real :: dpns !< depth to start terrain following [H ~> m or kg m-2]
52 real :: dsns !< depth to stop terrain following [H ~> m or kg m-2]
53 real :: min_dilate !< The minimum amount of dilation that is permitted when converting target
54 !! coordinates from z to z* [nondim]. This limit applies when wetting occurs.
55 real :: max_dilate !< The maximum amount of dilation that is permitted when converting target
56 !! coordinates from z to z* [nondim]. This limit applies when drying occurs.
57
58 real :: thkbot !< Thickness of a bottom boundary layer, within which hybgen does
59 !! something different. [H ~> m or kg m-2]
60
61 !> Shallowest depth for isopycnal layers [H ~> m or kg m-2]
62 real :: topiso_const
63 ! real, dimension(:,:), allocatable :: topiso
64
65 !> Nominal density of interfaces [R ~> kg m-3]
66 real, allocatable, dimension(:) :: target_density
67
68 real :: dp_far_from_sfc !< A distance that determines when an interface is suffiently far from
69 !! the surface that certain adjustments can be made in the Hybgen regridding
70 !! code [H ~> m or kg m-2]. In Hycom, this is set to tenm (nominally 10 m).
71 real :: dp_far_from_bot !< A distance that determines when an interface is suffiently far from
72 !! the bottom that certain adjustments can be made in the Hybgen regridding
73 !! code [H ~> m or kg m-2]. In Hycom, this is set to onem (nominally 1 m).
74 real :: h_thin !< A layer thickness below which a layer is considered to be too thin for
75 !! certain adjustments to be made in the Hybgen regridding code [H ~> m or kg m-2].
76 !! In Hycom, this is set to onemm (nominally 0.001 m).
77
78 real :: rho_eps !< A small nonzero density that is used to prevent division by zero
79 !! in several expressions in the Hybgen regridding code [R ~> kg m-3].
80
81 real :: onem !< Nominally one m in thickness units [H ~> m or kg m-2], used only in
82 !! certain debugging tests.
83
84end type hybgen_regrid_cs
85
86
89
90contains
91
92!> Initialise a hybgen_regrid_CS control structure and store its parameters
93subroutine init_hybgen_regrid(CS, GV, US, param_file)
94 type(hybgen_regrid_cs), pointer :: cs !< Unassociated pointer to hold the control structure
95 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
96 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
97 type(param_file_type), intent(in) :: param_file !< Parameter file
98
99 character(len=40) :: mdl = "MOM_hybgen_regrid" ! This module's name.
100 real :: hybrlx ! The number of remappings over which to move toward the target coordinate [timesteps]
101 character(len=40) :: dp0_coord_var, ds0_coord_var, rho_coord_var
102 character(len=200) :: filename, coord_file, inputdir ! Strings for file/path
103 logical :: use_coord_file
104 integer :: k
105
106 if (associated(cs)) call mom_error(fatal, "init_hybgen_regrid: CS already associated!")
107 allocate(cs)
108
109 cs%nk = gv%ke
110
111 allocate(cs%target_density(cs%nk))
112 allocate(cs%dp0k(cs%nk), source=0.0) ! minimum deep z-layer separation
113 allocate(cs%ds0k(cs%nk), source=0.0) ! minimum shallow z-layer separation
114
115 do k=1,cs%nk ; cs%target_density(k) = gv%Rlay(k) ; enddo
116
117 call get_param(param_file, mdl, "P_REF", cs%ref_pressure, &
118 "The pressure that is used for calculating the coordinate "//&
119 "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
120 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
121 units="Pa", default=2.0e7, scale=us%Pa_to_RL2_T2)
122
123 call get_param(param_file, mdl, "HYBGEN_MIN_THICKNESS", cs%min_thickness, &
124 "The minimum layer thickness allowed when regridding with Hybgen.", &
125 units="m", default=0.0, scale=gv%m_to_H )
126
127 call get_param(param_file, mdl, "HYBGEN_N_SIGMA", cs%nsigma, &
128 "The number of sigma-coordinate (terrain-following) layers with Hybgen regridding.", &
129 default=0)
130 call get_param(param_file, mdl, "HYBGEN_COORD_FILE", coord_file, &
131 "The file from which the Hybgen profile is read, or blank to use a list of "//&
132 "real input parameters from the MOM_input file.", default="")
133
134 use_coord_file = (len_trim(coord_file) > 0)
135 call get_param(param_file, mdl, "HYBGEN_DEEP_DZ_PR0FILE", cs%dp0k, &
136 "The layerwise list of deep z-level minimum thicknesses for Hybgen (dp0k in Hycom).", &
137 units="m", default=0.0, scale=gv%m_to_H, do_not_log=use_coord_file)
138 call get_param(param_file, mdl, "HYBGEN_SHALLOW_DZ_PR0FILE", cs%ds0k, &
139 "The layerwise list of shallow z-level minimum thicknesses for Hybgen (ds0k in Hycom).", &
140 units="m", default=0.0, scale=gv%m_to_H, do_not_log=use_coord_file)
141
142 if (use_coord_file) then
143 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
144 inputdir = slasher(inputdir)
145 filename = trim(inputdir)//trim(coord_file)
146 call log_param(param_file, mdl, "INPUTDIR/HYBGEN_COORD_FILE", filename)
147 if (.not.file_exists(filename)) call mom_error(fatal, &
148 " set_coord_from_file: Unable to open "//trim(filename))
149
150 call get_param(param_file, mdl, "HYBGEN_DEEP_DZ_VAR", dp0_coord_var, &
151 "The variable in HYBGEN_COORD_FILE that is to be used for the "//&
152 "deep z-level minimum thicknesses for Hybgen (dp0k in Hycom).", &
153 default="dp0")
154 call get_param(param_file, mdl, "HYBGEN_SHALLOW_DZ_VAR", ds0_coord_var, &
155 "The variable in HYBGEN_COORD_FILE that is to be used for the "//&
156 "shallow z-level minimum thicknesses for Hybgen (ds0k in Hycom).", &
157 default="ds0")
158 call get_param(param_file, mdl, "HYBGEN_TGT_DENSITY_VAR", rho_coord_var, &
159 "The variable in HYBGEN_COORD_FILE that is to be used for the Hybgen "//&
160 "target layer densities, or blank to reuse the values in GV%Rlay.", &
161 default="")
162
163 call mom_read_data(filename, dp0_coord_var, cs%dp0k, scale=gv%m_to_H)
164
165 call mom_read_data(filename, ds0_coord_var, cs%ds0k, scale=gv%m_to_H)
166
167 if (len_trim(rho_coord_var) > 0) &
168 call mom_read_data(filename, rho_coord_var, cs%target_density, scale=us%kg_m3_to_R)
169 endif
170
171 call get_param(param_file, mdl, "HYBGEN_ISOPYCNAL_DZ_MIN", cs%dp00i, &
172 "The Hybgen deep isopycnal spacing minimum thickness (dp00i in Hycom)", &
173 units="m", default=0.0, scale=gv%m_to_H)
174 call get_param(param_file, mdl, "HYBGEN_MIN_ISO_DEPTH", cs%topiso_const, &
175 "The Hybgen shallowest depth for isopycnal layers (isotop in Hycom)", &
176 units="m", default=0.0, scale=gv%m_to_H)
177 call get_param(param_file, mdl, "HYBGEN_RELAX_PERIOD", hybrlx, &
178 "The Hybgen coordinate relaxation period in timesteps, or 1 to move to "//&
179 "the new target coordinates in a single step. This must be >= 1.", &
180 units="timesteps", default=1.0)
181 if (hybrlx < 1.0) call mom_error(fatal, "init_hybgen_regrid: HYBGEN_RELAX_PERIOD must be at least 1.")
182 cs%qhybrlx = 1.0 / hybrlx
183 call get_param(param_file, mdl, "HYBGEN_BBL_THICKNESS", cs%thkbot, &
184 "A bottom boundary layer thickness within which Hybgen is able to move "//&
185 "overlying layers upward to match a target density.", &
186 units="m", default=0.0, scale=gv%m_to_H)
187 call get_param(param_file, mdl, "HYBGEN_FAR_FROM_SURFACE", cs%dp_far_from_sfc, &
188 "A distance that determines when an interface is suffiently far "//&
189 "from the surface that certain adjustments can be made in the Hybgen "//&
190 "regridding code. In Hycom, this is set to tenm (nominally 10 m).", &
191 units="m", default=10.0, scale=gv%m_to_H)
192 call get_param(param_file, mdl, "HYBGEN_FAR_FROM_BOTTOM", cs%dp_far_from_bot, &
193 "A distance that determines when an interface is suffiently far "//&
194 "from the bottom that certain adjustments can be made in the Hybgen "//&
195 "regridding code. In Hycom, this is set to onem (nominally 1 m).", &
196 units="m", default=1.0, scale=gv%m_to_H)
197 call get_param(param_file, mdl, "HYBGEN_H_THIN", cs%h_thin, &
198 "A layer thickness below which a layer is considered to be too thin for "//&
199 "certain adjustments to be made in the Hybgen regridding code. "//&
200 "In Hycom, this is set to onemm (nominally 0.001 m).", &
201 units="m", default=0.001, scale=gv%m_to_H)
202
203 call get_param(param_file, mdl, "HYBGEN_DENSITY_EPSILON", cs%rho_eps, &
204 "A small nonzero density that is used to prevent division by zero "//&
205 "in several expressions in the Hybgen regridding code.", &
206 units="kg m-3", default=1e-11, scale=us%kg_m3_to_R)
207
208
209 call get_param(param_file, mdl, "HYBGEN_REMAP_DENSITY_MATCH", cs%hybiso, &
210 "A tolerance between the layer densities and their target, within which "//&
211 "Hybgen determines that remapping uses PCM for a layer.", &
212 units="kg m-3", default=0.0, scale=us%kg_m3_to_R)
213 call get_param(param_file, mdl, "HYBGEN_REMAP_MIN_ZSTAR_DILATE", cs%min_dilate, &
214 "The maximum amount of dilation that is permitted when converting target "//&
215 "coordinates from z to z* [nondim]. This limit applies when drying occurs.", &
216 units="nondim", default=0.5)
217 call get_param(param_file, mdl, "HYBGEN_REMAP_MAX_ZSTAR_DILATE", cs%max_dilate, &
218 "The maximum amount of dilation that is permitted when converting target "//&
219 "coordinates from z to z* [nondim]. This limit applies when drying occurs.", &
220 units="nondim", default=2.0)
221
222 cs%onem = 1.0 * gv%m_to_H
223
224 do k=1,cs%nk ; cs%dp0k(k) = max(cs%dp0k(k), cs%min_thickness) ; enddo
225 cs%dp00i = max(cs%dp00i, cs%min_thickness)
226
227 ! Determine the depth range over which to use a sigma (terrain-following) coordinate.
228 ! --- terrain following starts at depth dpns and ends at depth dsns
229 if (cs%nsigma == 0) then
230 cs%dpns = cs%dp0k(1)
231 cs%dsns = 0.0
232 else
233 cs%dpns = 0.0
234 cs%dsns = 0.0
235 do k=1,cs%nsigma
236 cs%dpns = cs%dpns + cs%dp0k(k)
237 cs%dsns = cs%dsns + cs%ds0k(k)
238 enddo !k
239 endif !nsigma
240
241 cs%coord_scale = gv%H_to_m
242 cs%Rho_coord_scale = us%R_to_kg_m3
243
244end subroutine init_hybgen_regrid
245
246!> Writes out a file containing any available data related
247!! to the vertical grid used by the MOM ocean model.
248subroutine write_hybgen_coord_file(GV, CS, filepath)
249 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
250 type(hybgen_regrid_cs), intent(in) :: cs !< Control structure for this module
251 character(len=*), intent(in) :: filepath !< The full path to the file to write
252 ! Local variables
253 type(vardesc) :: vars(3)
254 type(mom_field) :: fields(3)
255 type(mom_infra_file) :: io_handle ! The I/O handle of the fileset
256
257 vars(1) = var_desc("dp0", "meter", "Deep z-level minimum thicknesses for Hybgen", '1', 'L', '1')
258 vars(2) = var_desc("ds0", "meter", "Shallow z-level minimum thicknesses for Hybgen", '1', 'L', '1')
259 vars(3) = var_desc("Rho_tgt", "kg m-3", "Target coordinate potential densities for Hybgen", '1', 'L', '1')
260 call create_mom_file(io_handle, trim(filepath), vars, 3, fields, &
261 single_file, gv=gv)
262
263 call mom_write_field(io_handle, fields(1), cs%dp0k, unscale=cs%coord_scale)
264 call mom_write_field(io_handle, fields(2), cs%ds0k, unscale=cs%coord_scale)
265 call mom_write_field(io_handle, fields(3), cs%target_density, unscale=cs%Rho_coord_scale)
266
267 call io_handle%close()
268end subroutine write_hybgen_coord_file
269
270!> This subroutine deallocates memory in the control structure for the hybgen module
271subroutine end_hybgen_regrid(CS)
272 type(hybgen_regrid_cs), pointer :: cs !< Coordinate control structure
273
274 ! nothing to do
275 if (.not. associated(cs)) return
276
277 deallocate(cs%target_density)
278 deallocate(cs%dp0k, cs%ds0k)
279 deallocate(cs)
280end subroutine end_hybgen_regrid
281
282!> This subroutine can be used to retrieve the parameters for the hybgen regrid module
283subroutine get_hybgen_regrid_params(CS, nk, ref_pressure, hybiso, nsigma, dp00i, qhybrlx, &
284 dp0k, ds0k, dpns, dsns, min_dilate, max_dilate, &
285 thkbot, topiso_const, target_density)
286 type(hybgen_regrid_cs), pointer :: cs !< Coordinate regridding control structure
287 integer, optional, intent(out) :: nk !< Number of layers on the target grid
288 real, optional, intent(out) :: ref_pressure !< Reference pressure for density calculations [R L2 T-2 ~> Pa]
289 real, optional, intent(out) :: hybiso !< Hybgen uses PCM if layer is within hybiso of target density [R ~> kg m-3]
290 integer, optional, intent(out) :: nsigma !< Number of sigma levels used by HYBGEN
291 real, optional, intent(out) :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2]
292 real, optional, intent(out) :: qhybrlx !< Fractional relaxation amount per timestep, 0 < qyhbrlx <= 1 [nondim]
293 real, optional, intent(out) :: dp0k(:) !< minimum deep z-layer separation [H ~> m or kg m-2]
294 real, optional, intent(out) :: ds0k(:) !< minimum shallow z-layer separation [H ~> m or kg m-2]
295 real, optional, intent(out) :: dpns !< depth to start terrain following [H ~> m or kg m-2]
296 real, optional, intent(out) :: dsns !< depth to stop terrain following [H ~> m or kg m-2]
297 real, optional, intent(out) :: min_dilate !< The minimum amount of dilation that is permitted when
298 !! converting target coordinates from z to z* [nondim].
299 !! This limit applies when wetting occurs.
300 real, optional, intent(out) :: max_dilate !< The maximum amount of dilation that is permitted when
301 !! converting target coordinates from z to z* [nondim].
302 !! This limit applies when drying occurs.
303 real, optional, intent(out) :: thkbot !< Thickness of a bottom boundary layer, within which
304 !! hybgen does something different. [H ~> m or kg m-2]
305 real, optional, intent(out) :: topiso_const !< Shallowest depth for isopycnal layers [H ~> m or kg m-2]
306 ! real, dimension(:,:), allocatable :: topiso
307 real, optional, intent(out) :: target_density(:) !< Nominal density of interfaces [R ~> kg m-3]
308
309 if (.not. associated(cs)) call mom_error(fatal, "get_hybgen_params: CS not associated")
310
311 if (present(nk)) nk = cs%nk
312 if (present(ref_pressure)) ref_pressure = cs%ref_pressure
313 if (present(hybiso)) hybiso = cs%hybiso
314 if (present(nsigma)) nsigma = cs%nsigma
315 if (present(dp00i)) dp00i = cs%dp00i
316 if (present(qhybrlx)) qhybrlx = cs%qhybrlx
317 if (present(dp0k)) then
318 if (size(dp0k) < cs%nk) call mom_error(fatal, "get_hybgen_regrid_params: "//&
319 "The dp0k argument is not allocated with enough space.")
320 dp0k(1:cs%nk) = cs%dp0k(1:cs%nk)
321 endif
322 if (present(ds0k)) then
323 if (size(ds0k) < cs%nk) call mom_error(fatal, "get_hybgen_regrid_params: "//&
324 "The ds0k argument is not allocated with enough space.")
325 ds0k(1:cs%nk) = cs%ds0k(1:cs%nk)
326 endif
327 if (present(dpns)) dpns = cs%dpns
328 if (present(dsns)) dsns = cs%dsns
329 if (present(min_dilate)) min_dilate = cs%min_dilate
330 if (present(max_dilate)) max_dilate = cs%max_dilate
331 if (present(thkbot)) thkbot = cs%thkbot
332 if (present(topiso_const)) topiso_const = cs%topiso_const
333 if (present(target_density)) then
334 if (size(target_density) < cs%nk) call mom_error(fatal, "get_hybgen_regrid_params: "//&
335 "The target_density argument is not allocated with enough space.")
336 target_density(1:cs%nk) = cs%target_density(1:cs%nk)
337 endif
338
339end subroutine get_hybgen_regrid_params
340
341
342!> Modify the input grid to give a new vertical grid based on the HYCOM hybgen code.
343subroutine hybgen_regrid(G, GV, US, dp, nom_depth_H, tv, CS, dzInterface, PCM_cell)
344 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
345 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
346 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
347 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
348 intent(in) :: dp !< Source grid layer thicknesses [H ~> m or kg m-2]
349 real, dimension(SZI_(G),SZJ_(G)), &
350 intent(in) :: nom_depth_h !< The bathymetric depth of this column
351 !! relative to mean sea level or another locally
352 !! valid reference height, converted to thickness
353 !! units [H ~> m or kg m-2]
354 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
355 type(hybgen_regrid_cs), intent(in) :: cs !< hybgen control structure
356 real, dimension(SZI_(G),SZJ_(G),CS%nk+1), &
357 intent(inout) :: dzinterface !< The change in height of each interface,
358 !! using a sign convention opposite to the change
359 !! in pressure [H ~> m or kg m-2]
360 logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
361 optional, intent(inout) :: pcm_cell !< If true, PCM remapping should be used in a cell.
362 !! This is effectively intent out, but values in wide
363 !! halo regions and land points are reused.
364
365 ! --- -------------------------------------
366 ! --- hybrid grid generator from HYCOM
367 ! --- -------------------------------------
368
369 ! These notes on the parameters for the hybrid grid generator are inhereted from the
370 ! Hycom source code for these algorithms.
371 !
372 ! From blkdat.input (units may have changed from m to pressure):
373 !
374 ! --- 'nsigma' = number of sigma levels
375 ! --- 'dp0k ' = layer k deep z-level spacing minimum thickness (m)
376 ! --- k=1,nk
377 ! --- 'ds0k ' = layer k shallow z-level spacing minimum thickness (m)
378 ! --- k=1,nsigma
379 ! --- 'dp00i' = deep isopycnal spacing minimum thickness (m)
380 ! --- 'isotop' = shallowest depth for isopycnal layers (m)
381 ! now in topiso(:,:)
382 ! --- 'sigma ' = isopycnal layer target densities (sigma units)
383 ! --- now in theta(:,:,1:nk)
384 !
385 ! --- the above specifies a vertical coord. that is isopycnal or:
386 ! --- near surface z in deep water, based on dp0k
387 ! --- near surface z in shallow water, based on ds0k and nsigma
388 ! --- terrain-following between them, based on ds0k and nsigma
389 !
390 ! --- terrain following starts at depth dpns=sum(dp0k(k),k=1,nsigma) and
391 ! --- ends at depth dsns=sum(ds0k(k),k=1,nsigma), and the depth of the
392 ! --- k-th layer interface varies linearly with total depth between
393 ! --- these two reference depths, i.e. a z-sigma-z fixed coordinate.
394 !
395 ! --- near the surface (i.e. shallower than isotop), layers are always
396 ! --- fixed depth (z or sigma).
397 ! -- layer 1 is always fixed, so isotop=0.0 is not realizable.
398 ! --- near surface layers can also be forced to be fixed depth
399 ! --- by setting target densities (sigma(k)) very small.
400 !
401 ! --- away from the surface, the minimum layer thickness is dp00i.
402 !
403 ! --- for fixed depth targets to be:
404 ! --- z-only set nsigma=0,
405 ! --- sigma-z (shallow-deep) use a very small ds0k(:),
406 ! --- sigma-only set nsigma=nk, dp0k large, and ds0k small.
407
408 ! These arrays work with the input column
409 real :: p_col(gv%ke) ! A column of reference pressures [R L2 T-2 ~> Pa]
410 real :: temp_in(gv%ke) ! A column of input potential temperatures [C ~> degC]
411 real :: saln_in(gv%ke) ! A column of input layer salinities [S ~> ppt]
412 real :: rcv_in(gv%ke) ! An input column of coordinate potential density [R ~> kg m-3]
413 real :: dp_in(gv%ke) ! The input column of layer thicknesses [H ~> m or kg m-2]
414 logical :: pcm_lay(gv%ke) ! If true for a layer, use PCM remapping for that layer
415
416 ! These arrays are on the target grid.
417 real :: rcv_tgt(cs%nk) ! Target potential density [R ~> kg m-3]
418 real :: rcv(cs%nk) ! Initial values of coordinate potential density on the target grid [R ~> kg m-3]
419 real :: h_col(cs%nk) ! A column of layer thicknesses [H ~> m or kg m-2]
420 real :: dz_int(cs%nk+1) ! The change in interface height due to remapping [H ~> m or kg m-2]
421 real :: rcv_integral ! Integrated coordinate potential density in a layer [R H ~> kg m-2 or kg2 m-5]
422
423 real :: qhrlx(cs%nk+1) ! Fractional relaxation within a timestep (between 0 and 1) [nondim]
424 real :: dp0ij(cs%nk) ! minimum layer thickness [H ~> m or kg m-2]
425 real :: dp0cum(cs%nk+1) ! minimum interface depth [H ~> m or kg m-2]
426
427 real :: h_tot ! Total thickness of the water column [H ~> m or kg m-2]
428 real :: nominaldepth ! Depth of ocean bottom (positive downward) [H ~> m or kg m-2]
429 real :: dilate ! A factor by which to dilate the target positions from z to z* [nondim]
430 integer :: fixlay ! Deepest fixed coordinate layer
431 integer, dimension(0:CS%nk) :: k_end ! The index of the deepest source layer that contributes to
432 ! each target layer, in the unusual case where the input grid is
433 ! larger than the new grid. This situation only occurs during certain
434 ! types of initialization or when generating output diagnostics.
435 integer :: i, j, k, nk, k2, nk_in
436
437 nk = cs%nk
438
439 p_col(:) = cs%ref_pressure
440
441 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1 ; if (g%mask2dT(i,j)>0.) then
442
443 ! Store one-dimensional arrays of thicknesses for the 'old' vertical grid before regridding
444 h_tot = 0.0
445 do k=1,gv%ke
446 temp_in(k) = tv%T(i,j,k)
447 saln_in(k) = tv%S(i,j,k)
448 dp_in(k) = dp(i,j,k)
449 h_tot = h_tot + dp_in(k)
450 enddo
451
452 ! This sets the input column's coordinate potential density from T and S.
453 call calculate_density(temp_in, saln_in, p_col, rcv_in, tv%eqn_of_state)
454
455 ! Set the initial properties on the new grid from the old grid.
456 nk_in = gv%ke
457 if (gv%ke > cs%nk) then ; do k=gv%ke,cs%nk+1,-1
458 ! Remove any excess massless layers from the bottom of the input column.
459 if (dp_in(k) > 0.0) exit
460 nk_in = k-1
461 enddo ; endif
462
463 if (cs%nk >= nk_in) then
464 ! Simply copy over the common layers. This is the usual case.
465 do k=1,min(cs%nk,gv%ke)
466 h_col(k) = dp_in(k)
467 rcv(k) = rcv_in(k)
468 enddo
469 if (cs%nk > gv%ke) then
470 ! Pad out the input column with additional massless layers with the bottom properties.
471 ! This case only occurs during initialization or perhaps when writing diagnostics.
472 do k=gv%ke+1,cs%nk
473 rcv(k) = rcv_in(gv%ke)
474 h_col(k) = 0.0
475 enddo
476 endif
477 else ! (CS%nk < nk_in)
478 ! The input column has more data than the output. For now, combine layers to
479 ! make them the same size, but there may be better approaches that should be taken.
480 ! This case only occurs during initialization or perhaps when writing diagnostics.
481 ! This case was not handled by the original Hycom code in hybgen.F90.
482 do k=0,cs%nk ; k_end(k) = (k * nk_in) / cs%nk ; enddo
483 do k=1,cs%nk
484 h_col(k) = 0.0 ; rcv_integral = 0.0
485 do k2=k_end(k-1) + 1,k_end(k)
486 h_col(k) = h_col(k) + dp_in(k2)
487 rcv_integral = rcv_integral + dp_in(k2)*rcv_in(k2)
488 enddo
489 if (h_col(k) > gv%H_subroundoff) then
490 ! Take the volume-weighted average properties.
491 rcv(k) = rcv_integral / h_col(k)
492 else ! Take the properties of the topmost source layer that contributes.
493 rcv(k) = rcv_in(k_end(k-1) + 1)
494 endif
495 enddo
496 endif
497
498 ! Set the target densities for the new layers.
499 do k=1,cs%nk
500 ! Rcv_tgt(k) = theta(i,j,k) ! If a 3-d target density were set up in theta, use that here.
501 rcv_tgt(k) = cs%target_density(k) ! MOM6 does not yet support 3-d target densities.
502 enddo
503
504 ! The following block of code is used to trigger z* stretching of the targets heights.
505 nominaldepth = nom_depth_h(i,j)
506 if (h_tot <= cs%min_dilate*nominaldepth) then
507 dilate = cs%min_dilate
508 elseif (h_tot >= cs%max_dilate*nominaldepth) then
509 dilate = cs%max_dilate
510 else
511 dilate = h_tot / nominaldepth
512 endif
513
514 ! Convert the regridding parameters into specific constraints for this column.
515 call hybgen_column_init(nk, cs%nsigma, cs%dp0k, cs%ds0k, cs%dp00i, &
516 cs%topiso_const, cs%qhybrlx, cs%dpns, cs%dsns, h_tot, dilate, &
517 h_col, fixlay, qhrlx, dp0ij, dp0cum)
518
519 ! Determine whether to require the use of PCM remapping from each source layer.
520 do k=1,gv%ke
521 if (cs%hybiso > 0.0) then
522 ! --- thin or isopycnal source layers are remapped with PCM.
523 pcm_lay(k) = (k > fixlay) .and. (abs(rcv(k) - rcv_tgt(k)) < cs%hybiso)
524 else ! hybiso==0.0, so purely isopycnal layers use PCM
525 pcm_lay(k) = .false.
526 endif ! hybiso
527 enddo !k
528
529 ! Determine the new layer thicknesses.
530 call hybgen_column_regrid(cs, nk, cs%thkbot, rcv_tgt, fixlay, qhrlx, dp0ij, &
531 dp0cum, rcv, h_col, dz_int)
532
533 ! Store the output from hybgenaij_regrid in 3-d arrays.
534 if (present(pcm_cell)) then ; do k=1,gv%ke
535 pcm_cell(i,j,k) = pcm_lay(k)
536 enddo ; endif
537
538 do k=1,nk+1
539 ! Note that dzInterface uses the opposite sign convention from the change in p.
540 dzinterface(i,j,k) = -dz_int(k)
541 enddo
542
543 else
544 if (present(pcm_cell)) then ; do k=1,gv%ke
545 pcm_cell(i,j,k) = .false.
546 enddo ; endif
547 do k=1,cs%nk+1 ; dzinterface(i,j,k) = 0.0 ; enddo
548 endif ; enddo ; enddo !i & j.
549
550end subroutine hybgen_regrid
551
552!> Initialize some of the variables that are used for regridding or unmixing, including the
553!! stretched contraits on where the new interfaces can be.
554subroutine hybgen_column_init(nk, nsigma, dp0k, ds0k, dp00i, topiso_i_j, &
555 qhybrlx, dpns, dsns, h_tot, dilate, h_col, &
556 fixlay, qhrlx, dp0ij, dp0cum)
557 integer, intent(in) :: nk !< The number of layers in the new grid
558 integer, intent(in) :: nsigma !< The number of sigma levels
559 real, intent(in) :: dp0k(nk) !< Layer deep z-level spacing minimum thicknesses [H ~> m or kg m-2]
560 real, intent(in) :: ds0k(nsigma) !< Layer shallow z-level spacing minimum thicknesses [H ~> m or kg m-2]
561 real, intent(in) :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2]
562 real, intent(in) :: topiso_i_j !< Shallowest depth for isopycnal layers [H ~> m or kg m-2]
563 real, intent(in) :: qhybrlx !< Fractional relaxation amount per timestep, 0 < qyhbrlx <= 1 [nondim]
564 real, intent(in) :: h_tot !< The sum of the initial layer thicknesses [H ~> m or kg m-2]
565 real, intent(in) :: dilate !< A factor by which to dilate the target positions
566 !! from z to z* [nondim]
567 real, intent(in) :: h_col(nk) !< Initial layer thicknesses [H ~> m or kg m-2]
568 real, intent(in) :: dpns !< Vertical sum of dp0k [H ~> m or kg m-2]
569 real, intent(in) :: dsns !< Vertical sum of ds0k [H ~> m or kg m-2]
570 integer, intent(out) :: fixlay !< Deepest fixed coordinate layer
571 real, intent(out) :: qhrlx(nk+1) !< Fractional relaxation within a timestep (between 0 and 1) [nondim]
572 real, intent(out) :: dp0ij(nk) !< minimum layer thickness [H ~> m or kg m-2]
573 real, intent(out) :: dp0cum(nk+1) !< minimum interface depth [H ~> m or kg m-2]
574
575 ! --- --------------------------------------------------------------
576 ! --- hybrid grid generator, single column - initialization.
577 ! --- --------------------------------------------------------------
578
579 ! Local variables
580 real :: qdep ! Total water column thickness as a fraction of dp0k (vs ds0k) [nondim]
581 real :: q ! A portion of the thickness that contributes to the new cell [H ~> m or kg m-2]
582 real :: p_int(nk+1) ! Interface depths [H ~> m or kg m-2]
583 integer :: k, fixall
584
585 ! --- dpns = sum(dp0k(k),k=1,nsigma)
586 ! --- dsns = sum(ds0k(k),k=1,nsigma)
587 ! --- terrain following starts (on the deep side) at depth dpns and ends (on the
588 ! --- shallow side) at depth dsns and the depth of the k-th layer interface varies
589 ! --- linearly with total depth between these two reference depths.
590 if ((h_tot >= dilate * dpns) .or. (dpns <= dsns)) then
591 qdep = 1.0 ! Not terrain following - this column is too thick or terrain following is disabled.
592 elseif (h_tot <= dilate * dsns) then
593 qdep = 0.0 ! Not terrain following - this column is too thin
594 else
595 qdep = (h_tot - dilate * dsns) / (dilate * (dpns - dsns))
596 endif
597
598 if (qdep < 1.0) then
599 ! Terrain following or shallow fixed coordinates, qhrlx=1 and ignore dp00
600 p_int( 1) = 0.0
601 dp0cum(1) = 0.0
602 qhrlx( 1) = 1.0
603 dp0ij( 1) = dilate * (qdep*dp0k(1) + (1.0-qdep)*ds0k(1))
604
605 dp0cum(2) = dp0cum(1) + dp0ij(1)
606 qhrlx( 2) = 1.0
607 p_int( 2) = p_int(1) + h_col(1)
608 do k=2,nk
609 qhrlx( k+1) = 1.0
610 dp0ij( k) = dilate * (qdep*dp0k(k) + (1.0-qdep)*ds0k(k))
611 dp0cum(k+1) = dp0cum(k) + dp0ij(k)
612 p_int( k+1) = p_int(k) + h_col(k)
613 enddo !k
614 else
615 ! Not terrain following
616 p_int( 1) = 0.0
617 dp0cum(1) = 0.0
618 qhrlx( 1) = 1.0 !no relaxation in top layer
619 dp0ij( 1) = dilate * dp0k(1)
620
621 dp0cum(2) = dp0cum(1) + dp0ij(1)
622 qhrlx( 2) = 1.0 !no relaxation in top layer
623 p_int( 2) = p_int(1) + h_col(1)
624 do k=2,nk
625 if ((dp0k(k) <= dp00i) .or. (dilate * dp0k(k) >= p_int(k) - dp0cum(k))) then
626 ! This layer is in fixed surface coordinates.
627 dp0ij(k) = dp0k(k)
628 qhrlx(k+1) = 1.0
629 else
630 q = dp0k(k) * (dilate * dp0k(k) / ( p_int(k) - dp0cum(k)) ) ! A fraction between 0 and 1 of dp0 to use here.
631 if (dp00i >= q) then
632 ! This layer is much deeper than the fixed surface coordinates.
633 dp0ij(k) = dp00i
634 qhrlx(k+1) = qhybrlx
635 else
636 ! This layer spans the margins of the fixed surface coordinates.
637 ! In this case dp00i < q < dp0k.
638 dp0ij(k) = dilate * q
639 qhrlx(k+1) = qhybrlx * (dp0k(k) - dp00i) / &
640 ((dp0k(k) - q) + (q - dp00i)*qhybrlx) ! 1 at dp0k, qhybrlx at dp00i
641 endif
642
643 ! The old equivalent code is:
644 ! hybrlx = 1.0 / qhybrlx
645 ! q = max( dp00i, dp0k(k) * (dp0k(k) / max(dp0k( k), p_int(k) - dp0cum(k)) ) )
646 ! qts = 1.0 - (q-dp00i) / (dp0k(k) - dp00i) !0 at q = dp0k, 1 at q=dp00i
647 ! qhrlx( k+1) = 1.0 / (1.0 + qts*(hybrlx-1.0)) !1 at dp0k, qhybrlx at dp00i
648 endif
649 dp0cum(k+1) = dp0cum(k) + dp0ij(k)
650 p_int(k+1) = p_int(k) + h_col(k)
651 enddo !k
652 endif !qdep<1:else
653
654 ! Identify the current fixed coordinate layers
655 fixlay = 1 !layer 1 always fixed
656 do k=2,nk
657 if (dp0cum(k) >= dilate * topiso_i_j) then
658 exit !layers k to nk might be isopycnal
659 endif
660 ! Top of layer is above topiso, i.e. always fixed coordinate layer
661 qhrlx(k+1) = 1.0 !no relaxation in fixed layers
662 fixlay = fixlay+1
663 enddo !k
664
665 fixall = fixlay
666 do k=fixall+1,nk
667 if (p_int(k+1) > dp0cum(k+1) + 0.1*dp0ij(k)) then
668 if ( (fixlay > fixall) .and. (p_int(k) > dp0cum(k)) ) then
669 ! --- The previous layer should remain fixed.
670 fixlay = fixlay-1
671 endif
672 exit !layers k to nk might be isopycnal
673 endif
674 ! Sometimes fixed coordinate layer
675 qhrlx(k) = 1.0 !no relaxation in fixed layers
676 fixlay = fixlay+1
677 enddo !k
678
679end subroutine hybgen_column_init
680
681!> The cushion function from Bleck & Benjamin, 1992, which returns a smoothly varying
682!! but limited value that goes between dp0 and delp
683real function cushn(delp, dp0)
684 real, intent(in) :: delp ! A thickness change [H ~> m or kg m-2]
685 real, intent(in) :: dp0 ! A non-negative reference thickness [H ~> m or kg m-2]
686
687 ! These are the nondimensional parameters that define the cushion function.
688 real, parameter :: qqmn=-4.0, qqmx=2.0 ! shifted range for cushn [nondim]
689! real, parameter :: qqmn=-2.0, qqmx=4.0 ! traditional range for cushn [nondim]
690! real, parameter :: qqmn=-4.0, qqmx=6.0 ! somewhat wider range for cushn [nondim]
691 ! These are derivative nondimensional parameters.
692 ! real, parameter :: cusha = qqmn**2 * (qqmx-1.0) / (qqmx-qqmn)**2
693 ! real, parameter :: I_qqmn = 1.0 / qqmn
694 real, parameter :: qq_scale = (qqmx-1.0) / (qqmx-qqmn)**2 ! A scaling factor based on qqmn and qqmx [nondim]
695 real, parameter :: i_qqmx = 1.0 / qqmx ! The inverse of qqmx [nondim]
696
697 ! --- if delp >= qqmx*dp0 >> dp0, cushn returns delp.
698 ! --- if delp <= qqmn*dp0 << -dp0, cushn returns dp0.
699
700 ! This is the original version from Hycom.
701 ! qq = max(qqmn, min(qqmx, delp/dp0))
702 ! cushn = dp0 * (1.0 + cusha * (1.0-I_qqmn*qq)**2) * max(1.0, delp/(dp0*qqmx))
703
704 ! This is mathematically equivalent, has one fewer divide, and works as intended even if dp0 = 0.
705 if (delp >= qqmx*dp0) then
706 cushn = delp
707 elseif (delp < qqmn*dp0) then
708 cushn = max(dp0, delp * i_qqmx)
709 else
710 cushn = max(dp0, delp * i_qqmx) * (1.0 + qq_scale * ((delp / dp0) - qqmn)**2)
711 endif
712
713end function cushn
714
715!> Create a new grid for a column of water using the Hybgen algorithm.
716subroutine hybgen_column_regrid(CS, nk, thkbot, Rcv_tgt, &
717 fixlay, qhrlx, dp0ij, dp0cum, Rcv, h_in, dp_int)
718 type(hybgen_regrid_cs), intent(in) :: CS !< hybgen regridding control structure
719 integer, intent(in) :: nk !< number of layers
720 real, intent(in) :: thkbot !< thickness of bottom boundary layer [H ~> m or kg m-2]
721 real, intent(in) :: Rcv_tgt(nk) !< Target potential density [R ~> kg m-3]
722 integer, intent(in) :: fixlay !< deepest fixed coordinate layer
723 real, intent(in) :: qhrlx( nk+1) !< relaxation coefficient per timestep [nondim]
724 real, intent(in) :: dp0ij( nk) !< minimum layer thickness [H ~> m or kg m-2]
725 real, intent(in) :: dp0cum(nk+1) !< minimum interface depth [H ~> m or kg m-2]
726 real, intent(in) :: Rcv(nk) !< Coordinate potential density [R ~> kg m-3]
727 real, intent(in) :: h_in(nk) !< Layer thicknesses [H ~> m or kg m-2]
728 real, intent(out) :: dp_int(nk+1) !< The change in interface positions [H ~> m or kg m-2]
729
730 ! --- ------------------------------------------------------
731 ! --- hybrid grid generator, single column - regrid.
732 ! --- ------------------------------------------------------
733
734 ! Local variables
735 real :: p_new ! A new interface position [H ~> m or kg m-2]
736 real :: pres_in(nk+1) ! layer interface positions [H ~> m or kg m-2]
737 real :: p_int(nk+1) ! layer interface positions [H ~> m or kg m-2]
738 real :: h_col(nk) ! Updated layer thicknesses [H ~> m or kg m-2]
739 real :: q_frac ! A fraction of a layer to entrain [nondim]
740 real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]
741 real :: h_hat3 ! Thickness movement upward across the interface between layers k-2 and k-3 [H ~> m or kg m-2]
742 real :: h_hat2 ! Thickness movement upward across the interface between layers k-1 and k-2 [H ~> m or kg m-2]
743 real :: h_hat ! Thickness movement upward across the interface between layers k and k-1 [H ~> m or kg m-2]
744 real :: h_hat0 ! A first guess at thickness movement upward across the interface
745 ! between layers k and k-1 [H ~> m or kg m-2]
746 real :: dh_cor ! Thickness changes [H ~> m or kg m-2]
747 logical :: trap_errors
748 integer :: k
749 character(len=256) :: mesg ! A string for output messages
750
751 ! This line needs to be consistent with the parameters set in cushn().
752 real, parameter :: qqmn=-4.0, qqmx=2.0 ! shifted range for cushn [nondim]
753! real, parameter :: qqmn=-2.0, qqmx=4.0 ! traditional range for cushn [nondim]
754! real, parameter :: qqmn=-4.0, qqmx=6.0 ! somewhat wider range for cushn [nondim]
755
756 trap_errors = .true.
757
758 do k=1,nk+1 ; dp_int(k) = 0.0 ; enddo
759
760 p_int(1) = 0.0
761 do k=1,nk
762 h_col(k) = max(h_in(k), 0.0)
763 p_int(k+1) = p_int(k) + h_col(k)
764 enddo
765 h_min = min( cs%min_thickness, p_int(nk+1)/real(cs%nk) )
766
767 if (trap_errors) then
768 do k=1,nk+1 ; pres_in(k) = p_int(k) ; enddo
769 endif
770
771 ! Try to restore isopycnic conditions by moving layer interfaces
772 ! qhrlx(k) are relaxation amounts per timestep.
773
774 ! Maintain prescribed thickness in layer k <= fixlay
775 ! There may be massless layers at the bottom, so work upwards.
776 do k=min(nk-1,fixlay),1,-1
777 p_new = min(dp0cum(k+1), p_int(nk+1) - (nk-k)*h_min) ! This could be positive or negative.
778 dh_cor = p_new - p_int(k+1)
779 if (k<fixlay) dh_cor = min(dh_cor, h_col(k+1) - h_min)
780 h_col(k) = h_col(k) + dh_cor
781 h_col(k+1) = h_col(k+1) - dh_cor
782 dp_int(k+1) = dp_int(k+1) + dh_cor
783 p_int(k+1) = p_new
784 enddo
785
786 ! Eliminate negative thicknesses below the fixed layers, entraining from below as necessary.
787 do k=fixlay+1,nk-1
788 if (h_col(k) >= h_min) exit ! usually get here quickly
789 dh_cor = h_min - h_col(k) ! This is positive.
790 h_col(k) = h_min ! = h_col(k) + dh_cor
791 h_col(k+1) = h_col(k+1) - dh_cor
792 dp_int(k+1) = dp_int(k+1) + dh_cor
793 p_int(k+1) = p_int(fixlay+1)
794 enddo
795 if (h_col(nk) < h_min) then ! This should be uncommon, and should only arise at the level of roundoff.
796 do k=nk,2,-1
797 if (h_col(k) >= h_min) exit
798 dh_cor = h_col(k) - h_min ! dh_cor is negative.
799 h_col(k-1) = h_col(k-1) + dh_cor
800 h_col(k) = h_min ! = h_col(k) - dh_cor
801 dp_int(k) = dp_int(k) + dh_cor
802 p_int(k) = p_int(k) + dh_cor
803 enddo
804 endif
805
806 ! Remap the non-fixed layers.
807
808 ! In the Hycom version, this loop was fused with the loop correcting water that is
809 ! too light, and it ran down the water column, but if there are a set of layers
810 ! that are very dense, that structure can lead to all of the water being remapped
811 ! into a single thick layer. Splitting the loops and running the loop upwards
812 ! (as is done here) avoids that catastrophic problem for layers that are far from
813 ! their targets. However, this code is still prone to a thin-thick-thin null mode.
814 do k=nk,fixlay+2,-1
815 ! This is how the Hycom code would do this loop: do k=fixlay+1,nk ; if (k>fixlay+1) then
816
817 if ((rcv(k) > rcv_tgt(k) + cs%rho_eps)) then
818 ! Water in layer k is too dense, so try to dilute with water from layer k-1
819 ! Do not move interface if k = fixlay + 1
820
821 if ((rcv(k-1) >= rcv_tgt(k-1)) .or. &
822 (p_int(k) <= dp0cum(k) + cs%dp_far_from_bot) .or. &
823 (h_col(k) <= h_col(k-1))) then
824 ! If layer k-1 is too light, there is a conflict in the direction the
825 ! inteface between them should move, so thicken the thinner of the two.
826
827 if ((rcv_tgt(k) - rcv(k-1)) <= cs%rho_eps) then
828 ! layer k-1 is far too dense, take the entire layer
829 ! If this code is working downward and this branch is repeated in a series
830 ! of successive layers, it can accumulate into a very thick homogenous layers.
831 h_hat0 = 0.0 ! This line was not in the Hycom version of hybgen.F90.
832 h_hat = dp0ij(k-1) - h_col(k-1)
833 else
834 ! Entrain enough from the layer above to bring layer k to its target density.
835 q_frac = (rcv_tgt(k) - rcv(k)) / (rcv_tgt(k) - rcv(k-1)) ! -1 <= q_frac < 0
836 h_hat0 = q_frac*h_col(k) ! -h_col(k-1) <= h_hat0 < 0
837 if (k == fixlay+2) then
838 ! Treat layer k-1 as fixed.
839 h_hat = max(h_hat0, dp0ij(k-1) - h_col(k-1))
840 else
841 ! Maintain the minimum thickess of layer k-1.
842 h_hat = cushn(h_hat0 + h_col(k-1), dp0ij(k-1)) - h_col(k-1)
843 endif !fixlay+2:else
844 endif
845 ! h_hat is usually negative, so this check may be unnecessary if the values of
846 ! dp0ij are limited to not be below the seafloor?
847 h_hat = min(h_hat, p_int(nk+1) - p_int(k))
848
849 ! If isopycnic conditions cannot be achieved because of a blocking
850 ! layer (thinner than its minimum thickness) in the interior ocean,
851 ! move interface k-1 (and k-2 if necessary) upward
852 ! Only work on layers that are sufficiently far from the fixed near-surface layers.
853 if ((h_hat >= 0.0) .and. (k > fixlay+2) .and. (p_int(k-1) > dp0cum(k-1) + cs%dp_far_from_sfc)) then
854
855 ! Only act if interface k-1 is near the bottom or layer k-2 could donate water.
856 if ( (p_int(nk+1) - p_int(k-1) < thkbot) .or. &
857 (h_col(k-2) > qqmx*dp0ij(k-2)) ) then
858 ! Determine how much water layer k-2 could supply without becoming too thin.
859 if (k == fixlay+3) then
860 ! Treat layer k-2 as fixed.
861 h_hat2 = max(h_hat0 - h_hat, dp0ij(k-2) - h_col(k-2))
862 else
863 ! Maintain minimum thickess of layer k-2.
864 h_hat2 = cushn(h_col(k-2) + (h_hat0 - h_hat), dp0ij(k-2)) - h_col(k-2)
865 endif !fixlay+3:else
866
867 if (h_hat2 < -cs%h_thin) then
868 dh_cor = qhrlx(k-1) * max(h_hat2, -h_hat - h_col(k-1))
869 h_col(k-2) = h_col(k-2) + dh_cor
870 h_col(k-1) = h_col(k-1) - dh_cor
871 dp_int(k-1) = dp_int(k-1) + dh_cor
872 p_int(k-1) = p_int(k-1) + dh_cor
873 ! Recalculate how much layer k-1 could donate to layer k.
874 h_hat = cushn(h_hat0 + h_col(k-1), dp0ij(k-1)) - h_col(k-1)
875 elseif (k <= fixlay+3) then
876 ! Do nothing.
877 elseif ( (p_int(k-2) > dp0cum(k-2) + cs%dp_far_from_sfc) .and. &
878 ( (p_int(nk+1) - p_int(k-2) < thkbot) .or. &
879 (h_col(k-3) > qqmx*dp0ij(k-3)) ) ) then
880
881 ! Determine how much water layer k-3 could supply without becoming too thin.
882 if (k == fixlay+4) then
883 ! Treat layer k-3 as fixed.
884 h_hat3 = max(h_hat0 - h_hat, dp0ij(k-3) - h_col(k-3))
885 else
886 ! Maintain minimum thickess of layer k-3.
887 h_hat3 = cushn(h_col(k-3) + (h_hat0 - h_hat), dp0ij(k-3)) - h_col(k-3)
888 endif !fixlay+4:else
889 if (h_hat3 < -cs%h_thin) then
890 ! Water is moved from layer k-3 to k-2, but do not dilute layer k-2 too much.
891 dh_cor = qhrlx(k-2) * max(h_hat3, -h_col(k-2))
892 h_col(k-3) = h_col(k-3) + dh_cor
893 h_col(k-2) = h_col(k-2) - dh_cor
894 dp_int(k-2) = dp_int(k-2) + dh_cor
895 p_int(k-2) = p_int(k-2) + dh_cor
896
897 ! Now layer k-2 might be able donate to layer k-1.
898 h_hat2 = cushn(h_col(k-2) + (h_hat0 - h_hat), dp0ij(k-2)) - h_col(k-2)
899 if (h_hat2 < -cs%h_thin) then
900 dh_cor = qhrlx(k-1) * (max(h_hat2, -h_hat - h_col(k-1)) )
901 h_col(k-2) = h_col(k-2) + dh_cor
902 h_col(k-1) = h_col(k-1) - dh_cor
903 dp_int(k-1) = dp_int(k-1) + dh_cor
904 p_int(k-1) = p_int(k-1) + dh_cor
905 ! Recalculate how much layer k-1 could donate to layer k.
906 h_hat = cushn(h_hat0 + h_col(k-1), dp0ij(k-1)) - h_col(k-1)
907 endif !h_hat2
908 endif !h_hat3
909 endif !h_hat2:blocking
910 endif ! Layer k-2 could move.
911 endif ! blocking, i.e., h_hat >= 0, and far enough from the fixed layers to permit a change.
912
913 if (h_hat < 0.0) then
914 ! entrain layer k-1 water into layer k, move interface up.
915 dh_cor = qhrlx(k) * h_hat
916 h_col(k-1) = h_col(k-1) + dh_cor
917 h_col(k) = h_col(k) - dh_cor
918 dp_int(k) = dp_int(k) + dh_cor
919 p_int(k) = p_int(k) + dh_cor
920 endif !entrain
921
922 endif !too-dense adjustment
923 endif
924
925 ! In the original Hycom version, there is not a break between these two loops.
926 enddo
927
928 do k=fixlay+1,nk
929 if (rcv(k) < rcv_tgt(k) - cs%rho_eps) then ! layer too light
930 ! Water in layer k is too light, so try to dilute with water from layer k+1.
931 ! Entrainment is not possible if layer k touches the bottom.
932 if (p_int(k+1) < p_int(nk+1)) then ! k<nk
933 if ((rcv(k+1) <= rcv_tgt(k+1)) .or. &
934 (p_int(k+1) <= dp0cum(k+1) + cs%dp_far_from_bot) .or. &
935 (h_col(k) < h_col(k+1))) then
936 ! If layer k+1 is too dense, there is a conflict in the direction the
937 ! inteface between them should move, so thicken the thinner of the two.
938
939 if ((rcv(k+1) - rcv_tgt(k)) <= cs%rho_eps) then
940 ! layer k+1 is far too light, so take the entire layer
941 ! Because this code is working downward, this flux does not accumulate across
942 ! successive layers.
943 h_hat = h_col(k+1)
944 else
945 q_frac = (rcv_tgt(k) - rcv(k)) / (rcv(k+1) - rcv_tgt(k)) ! 0 < q_frac <= 1
946 h_hat = q_frac*h_col(k)
947 endif
948
949 ! If layer k+1, or layer k+2, does not touch the bottom, maintain minimum
950 ! thicknesses of layers k and k+1 as much as possible. otherwise, permit
951 ! layers to collapse to zero thickness at the bottom.
952 if (p_int(min(k+3,nk+1)) < p_int(nk+1)) then
953 if (p_int(nk+1) - p_int(k) > dp0ij(k) + dp0ij(k+1)) then
954 h_hat = h_col(k+1) - cushn(h_col(k+1) - h_hat, dp0ij(k+1))
955 endif
956 ! Try to bring layer layer k up to its minimum thickness.
957 h_hat = max(h_hat, dp0ij(k) - h_col(k))
958 ! Do not drive layer k+1 below its minimum thickness or take more than half of it.
959 h_hat = min(h_hat, max(0.5*h_col(k+1), h_col(k+1) - dp0ij(k+1)) )
960 else
961 ! Layers that touch the bottom can lose their entire contents.
962 h_hat = min(h_col(k+1), h_hat)
963 endif !p.k+2<p.nk+1
964
965 if (h_hat > 0.0) then
966 ! Entrain layer k+1 water into layer k.
967 dh_cor = qhrlx(k+1) * h_hat
968 h_col(k) = h_col(k) + dh_cor
969 h_col(k+1) = h_col(k+1) - dh_cor
970 dp_int(k+1) = dp_int(k+1) + dh_cor
971 p_int(k+1) = p_int(k+1) + dh_cor
972 endif !entrain
973
974 endif !too-light adjustment
975 endif !above bottom
976 endif !too light
977
978 ! If layer above is still too thin, move interface down.
979 dh_cor = min(qhrlx(k-1) * min(dp0ij(k-1) - h_col(k-1), p_int(nk+1) - p_int(k)), h_col(k))
980 if (dh_cor > 0.0) then
981 h_col(k-1) = h_col(k-1) + dh_cor
982 h_col(k) = h_col(k) - dh_cor
983 dp_int(k) = dp_int(k) + dh_cor
984 p_int(k) = p_int(k) + dh_cor
985 endif
986
987 enddo !k Hybrid vertical coordinate relocation moving interface downward
988
989 if (trap_errors) then
990 ! Verify that everything is consistent.
991 do k=1,nk
992 if (abs((h_col(k) - h_in(k)) + (dp_int(k) - dp_int(k+1))) > 1.0e-13*max(p_int(nk+1), cs%onem)) then
993 write(mesg, '("k ",I0," h ",es13.4," h_in ",es13.4, " dp ",2es13.4," err ",es13.4)') &
994 k, h_col(k), h_in(k), dp_int(k), dp_int(k+1), (h_col(k) - h_in(k)) + (dp_int(k) - dp_int(k+1))
995 call mom_error(fatal, "Mismatched thickness changes in hybgen_regrid: "//trim(mesg))
996 endif
997 if (h_col(k) < 0.0) then ! Could instead do: -1.0e-15*max(p_int(nk+1), CS%onem)) then
998 write(mesg, '("k ",I0," h ",es13.4," h_in ",es13.4, " dp ",2es13.4, " fixlay ",I0)') &
999 k, h_col(k), h_in(k), dp_int(k), dp_int(k+1), fixlay
1000 call mom_error(fatal, "Significantly negative final thickness in hybgen_regrid: "//trim(mesg))
1001 endif
1002 enddo
1003 do k=1,nk+1
1004 if (abs(dp_int(k) - (p_int(k) - pres_in(k))) > 1.0e-13*max(p_int(nk+1), cs%onem)) then
1005 call mom_error(fatal, "Mismatched interface height changes in hybgen_regrid.")
1006 endif
1007 enddo
1008 endif
1009
1010end subroutine hybgen_column_regrid
1011
1012end module mom_hybgen_regrid
1013
1014! This code was translated in 2022 from the HYCOM hybgen code, which was primarily developed
1015! between 2000 and 2015, with some minor subsequent changes and bug fixes.