MOM_hybgen_unmix.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 unmixing routines from HYCOM, with
6!! modifications to follow the MOM6 coding conventions and several bugs fixed
8
9use mom_eos, only : eos_type, calculate_density, calculate_density_derivs
10use mom_error_handler, only : mom_mesg, mom_error, fatal, warning
11use mom_file_parser, only : get_param, param_file_type, log_param
15use mom_tracer_registry, only : tracer_registry_type, tracer_type, mom_tracer_chkinv
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_unmix_cs ; private
26
27 integer :: nsigma !< Number of sigma levels used by HYBGEN
28 real :: hybiso !< Hybgen uses PCM if layer is within hybiso of target density [R ~> kg m-3]
29
30 real :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2]
31 real :: qhybrlx !< Hybgen relaxation amount per thermodynamic time steps [nondim]
32
33 real, allocatable, dimension(:) :: &
34 dp0k, & !< minimum deep z-layer separation [H ~> m or kg m-2]
35 ds0k !< minimum shallow z-layer separation [H ~> m or kg m-2]
36
37 real :: dpns !< depth to start terrain following [H ~> m or kg m-2]
38 real :: dsns !< depth to stop terrain following [H ~> m or kg m-2]
39 real :: min_dilate !< The minimum amount of dilation that is permitted when converting target
40 !! coordinates from z to z* [nondim]. This limit applies when wetting occurs.
41 real :: max_dilate !< The maximum amount of dilation that is permitted when converting target
42 !! coordinates from z to z* [nondim]. This limit applies when drying occurs.
43
44 real :: topiso_const !< Shallowest depth for isopycnal layers [H ~> m or kg m-2]
45 ! real, dimension(:,:), allocatable :: topiso
46
47 real :: ref_pressure !< Reference pressure for density calculations [R L2 T-2 ~> Pa]
48 real, allocatable, dimension(:) :: target_density !< Nominal density of interfaces [R ~> kg m-3]
49
50end type hybgen_unmix_cs
51
54
55contains
56
57!> Initialise a hybgen_unmix_CS control structure and store its parameters
58subroutine init_hybgen_unmix(CS, GV, US, param_file, hybgen_regridCS)
59 type(hybgen_unmix_cs), pointer :: cs !< Unassociated pointer to hold the control structure
60 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
61 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
62 type(param_file_type), intent(in) :: param_file !< Parameter file
63 type(hybgen_regrid_cs), pointer :: hybgen_regridcs !< Control structure for hybgen
64 !! regridding for sharing parameters.
65 integer :: k
66
67 if (associated(cs)) call mom_error(fatal, "init_hybgen_unmix: CS already associated!")
68 allocate(cs)
69 allocate(cs%target_density(gv%ke))
70
71 allocate(cs%dp0k(gv%ke), source=0.0) ! minimum deep z-layer separation
72 allocate(cs%ds0k(gv%ke), source=0.0) ! minimum shallow z-layer separation
73
74 ! Set the parameters for the hybgen unmixing from a hybgen regridding control structure.
75 call get_hybgen_regrid_params(hybgen_regridcs, ref_pressure=cs%ref_pressure, &
76 nsigma=cs%nsigma, dp0k=cs%dp0k, ds0k=cs%ds0k, &
77 dp00i=cs%dp00i, topiso_const=cs%topiso_const, qhybrlx=cs%qhybrlx, &
78 hybiso=cs%hybiso, min_dilate=cs%min_dilate, max_dilate=cs%max_dilate, &
79 target_density=cs%target_density)
80
81 ! Determine the depth range over which to use a sigma (terrain-following) coordinate.
82 ! --- terrain following starts at depth dpns and ends at depth dsns
83 if (cs%nsigma == 0) then
84 cs%dpns = cs%dp0k(1)
85 cs%dsns = 0.0
86 else
87 cs%dpns = 0.0
88 cs%dsns = 0.0
89 do k=1,cs%nsigma
90 cs%dpns = cs%dpns + cs%dp0k(k)
91 cs%dsns = cs%dsns + cs%ds0k(k)
92 enddo !k
93 endif !nsigma
94
95end subroutine init_hybgen_unmix
96
97!> This subroutine deallocates memory in the control structure for the hybgen unmixing module
98subroutine end_hybgen_unmix(CS)
99 type(hybgen_unmix_cs), pointer :: cs !< Coordinate control structure
100
101 ! nothing to do
102 if (.not. associated(cs)) return
103
104 deallocate(cs%target_density)
105 deallocate(cs%dp0k, cs%ds0k)
106 deallocate(cs)
107end subroutine end_hybgen_unmix
108
109!> This subroutine can be used to set the parameters for the hybgen module
110subroutine set_hybgen_unmix_params(CS, min_thickness)
111 type(hybgen_unmix_cs), pointer :: cs !< Coordinate unmixing control structure
112 real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
113
114 if (.not. associated(cs)) call mom_error(fatal, "set_hybgen_params: CS not associated")
115
116! if (present(min_thickness)) CS%min_thickness = min_thickness
117end subroutine set_hybgen_unmix_params
118
119
120!> Unmix the properties in the lowest layer with mass if it is too light, and make
121!! any other changes to the water column to prepare for regridding.
122subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
123 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
124 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
125 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
126 type(hybgen_unmix_cs), intent(in) :: cs !< hybgen control structure
127 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
128 type(tracer_registry_type), pointer :: reg !< Tracer registry structure
129 integer, intent(in) :: ntr !< The number of tracers in the registry, or
130 !! 0 if the registry is not in use.
131 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
132 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
133
134! --- --------------------------------------------
135! --- hybrid grid generator, single j-row (part A).
136! --- --------------------------------------------
137
138 character(len=256) :: mesg ! A string for output messages
139 integer :: fixlay ! deepest fixed coordinate layer
140 real :: qhrlx( gv%ke+1) ! relaxation coefficient per timestep [nondim]
141 real :: dp0ij( gv%ke) ! minimum layer thickness [H ~> m or kg m-2]
142 real :: dp0cum(gv%ke+1) ! minimum interface depth [H ~> m or kg m-2]
143
144 real :: rcv_tgt(gv%ke) ! Target potential density [R ~> kg m-3]
145 real :: temp(gv%ke) ! A column of potential temperature [C ~> degC]
146 real :: saln(gv%ke) ! A column of salinity [S ~> ppt]
147 real :: rcv(gv%ke) ! A column of coordinate potential density [R ~> kg m-3]
148 real :: h_col(gv%ke) ! A column of layer thicknesses [H ~> m or kg m-2]
149 real :: p_col(gv%ke) ! A column of reference pressures [R L2 T-2 ~> Pa]
150 real :: tracer(gv%ke,max(ntr,1)) ! Columns of each tracer [Conc]
151 real :: h_tot ! Total thickness of the water column [H ~> m or kg m-2]
152 real :: dz_tot ! Vertical distance between the top and bottom of the water column [Z ~> m]
153 real :: nominaldepth ! Depth of ocean bottom in thickness units (positive downward) [H ~> m or kg m-2]
154 real :: h_thin ! A negligibly small thickness to identify essentially
155 ! vanished layers [H ~> m or kg m-2]
156 real :: dilate ! A factor by which to dilate the target positions from z to z* [nondim]
157
158 real :: th_tot_in, th_tot_out ! Column integrated temperature [C H ~> degC m or degC kg m-2]
159 real :: sh_tot_in, sh_tot_out ! Column integrated salinity [S H ~> ppt m or ppt kg m-2]
160 real :: trh_tot_in(max(ntr,1)) ! Initial column integrated tracer amounts [conc H ~> conc m or conc kg m-2]
161 real :: trh_tot_out(max(ntr,1)) ! Final column integrated tracer amounts [conc H ~> conc m or conc kg m-2]
162
163 logical :: debug_conservation ! If true, test for non-conservation.
164 logical :: terrain_following ! True if this column is terrain following.
165 integer :: trcflg(max(ntr,1)) ! Hycom tracer type flag for each tracer
166 integer :: i, j, k, nk, m
167
168 nk = gv%ke
169
170 ! Set all tracers to be passive. Setting this to 2 treats a tracer like temperature.
171 trcflg(:) = 3
172
173 h_thin = 1e-6*gv%m_to_H
174 debug_conservation = .false. ! Set this to true for debugging
175
176 if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then
177 if (tv%valid_SpV_halo < 0) then
178 mesg = "invalid values of SpV_avg."
179 else
180 mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed."
181 endif
182 call mom_error(fatal, "hybgen_unmix called in fully non-Boussinesq mode with "//trim(mesg))
183 endif
184
185 p_col(:) = cs%ref_pressure
186
187 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1 ; if (g%mask2dT(i,j)>0.) then
188
189 h_tot = 0.0
190 do k=1,nk
191 ! Rcv_tgt(k) = theta(i,j,k) ! If a 3-d target density were set up in theta, use that here.
192 rcv_tgt(k) = cs%target_density(k) ! MOM6 does not yet support 3-d target densities.
193 h_col(k) = h(i,j,k)
194 h_tot = h_tot + h_col(k)
195 temp(k) = tv%T(i,j,k)
196 saln(k) = tv%S(i,j,k)
197 enddo
198
199 ! This sets the potential density from T and S.
200 call calculate_density(temp, saln, p_col, rcv, tv%eqn_of_state)
201
202 do m=1,ntr ; do k=1,nk
203 tracer(k,m) = reg%Tr(m)%t(i,j,k)
204 enddo ; enddo
205
206 ! Store original amounts to test for conservation of temperature, salinity, and tracers.
207 if (debug_conservation) then
208 th_tot_in = 0.0 ; sh_tot_in = 0.0 ; trh_tot_in(:) = 0.0
209 do k=1,nk
210 sh_tot_in = sh_tot_in + h_col(k)*saln(k)
211 th_tot_in = th_tot_in + h_col(k)*temp(k)
212 enddo
213 do m=1,ntr ; do k=1,nk
214 trh_tot_in(m) = trh_tot_in(m) + h_col(k)*tracer(k,m)
215 enddo ; enddo
216 endif
217
218 ! The following block of code is used to trigger z* stretching of the targets heights.
219 if (allocated(tv%SpV_avg)) then ! This is the fully non-Boussinesq version
220 dz_tot = 0.0
221 do k=1,nk
222 dz_tot = dz_tot + gv%H_to_RZ * tv%SpV_avg(i,j,k) * h_col(k)
223 enddo
224 if (dz_tot <= cs%min_dilate * (g%meanSL(i,j) + g%bathyT(i,j))) then
225 dilate = cs%min_dilate
226 elseif (dz_tot >= cs%max_dilate * (g%meanSL(i,j) + g%bathyT(i,j))) then
227 dilate = cs%max_dilate
228 else
229 dilate = dz_tot / (g%meanSL(i,j) + g%bathyT(i,j))
230 endif
231 else
232 nominaldepth = (g%meanSL(i,j) + g%bathyT(i,j)) * gv%Z_to_H
233 if (h_tot <= cs%min_dilate * nominaldepth) then
234 dilate = cs%min_dilate
235 elseif (h_tot >= cs%max_dilate * nominaldepth) then
236 dilate = cs%max_dilate
237 else
238 dilate = h_tot / nominaldepth
239 endif
240 endif
241
242 terrain_following = (h_tot < dilate*cs%dpns) .and. (cs%dpns >= cs%dsns)
243
244 ! Convert the regridding parameters into specific constraints for this column.
245 call hybgen_column_init(nk, cs%nsigma, cs%dp0k, cs%ds0k, cs%dp00i, &
246 cs%topiso_const, cs%qhybrlx, cs%dpns, cs%dsns, h_tot, dilate, &
247 h_col, fixlay, qhrlx, dp0ij, dp0cum)
248
249 ! Do any unmixing of the column that is needed to move the layer properties toward their targets.
250 call hybgen_column_unmix(cs, nk, rcv_tgt, temp, saln, rcv, tv%eqn_of_state, &
251 ntr, tracer, trcflg, fixlay, qhrlx, h_col, &
252 terrain_following, h_thin)
253
254 ! Store the output from hybgen_unmix in the 3-d arrays.
255 do k=1,nk
256 h(i,j,k) = h_col(k)
257 enddo
258 ! Note that temperature and salinity are among the tracers unmixed here.
259 do m=1,ntr ; do k=1,nk
260 reg%Tr(m)%t(i,j,k) = tracer(k,m)
261 enddo ; enddo
262 ! However, temperature and salinity may have been treated differently from other tracers.
263 do k=1,nk
264 tv%T(i,j,k) = temp(k)
265 tv%S(i,j,k) = saln(k)
266 enddo
267
268 ! Test for conservation of temperature, salinity, and tracers.
269 if (debug_conservation) then
270 th_tot_out = 0.0 ; sh_tot_out = 0.0 ; trh_tot_out(:) = 0.0
271 do k=1,nk
272 sh_tot_out = sh_tot_out + h_col(k)*saln(k)
273 th_tot_out = th_tot_out + h_col(k)*temp(k)
274 enddo
275 do m=1,ntr ; do k=1,nk
276 trh_tot_out(m) = trh_tot_out(m) + h_col(k)*tracer(k,m)
277 enddo ; enddo
278 if (abs(sh_tot_in - sh_tot_out) > 1.e-15*(abs(sh_tot_in) + abs(sh_tot_out))) then
279 write(mesg, '("i,j=",I0,",",I0," Sh_tot = ",2es17.8," err = ",es13.4)') &
280 i, j, sh_tot_in, sh_tot_out, (sh_tot_in - sh_tot_out)
281 call mom_error(fatal, "Mismatched column salinity in hybgen_unmix: "//trim(mesg))
282 endif
283 if (abs(th_tot_in - th_tot_out) > 1.e-10*(abs(th_tot_in) + abs(th_tot_out))) then
284 write(mesg, '("i,j=",I0,",",I0," Th_tot = ",2es17.8," err = ",es13.4)') &
285 i, j, th_tot_in, th_tot_out, (th_tot_in - th_tot_out)
286 call mom_error(fatal, "Mismatched column temperature in hybgen_unmix: "//trim(mesg))
287 endif
288 do m=1,ntr
289 if (abs(trh_tot_in(m) - trh_tot_out(m)) > 1.e-10*(abs(trh_tot_in(m)) + abs(trh_tot_out(m)))) then
290 write(mesg, '("i,j=",I0,",",I0," Trh_tot(",i0,") = ",2es17.8," err = ",es13.4)') &
291 i, j, m, trh_tot_in(m), trh_tot_out(m), (trh_tot_in(m) - trh_tot_out(m))
292 call mom_error(fatal, "Mismatched column tracer in hybgen_unmix: "//trim(mesg))
293 endif
294 enddo
295 endif
296 endif ; enddo ; enddo !i & j.
297
298 ! Update the layer properties
299 if (allocated(tv%SpV_avg)) call calc_derived_thermo(tv, h, g, gv, us, halo=1)
300
301end subroutine hybgen_unmix
302
303
304!> Unmix the properties in the lowest layer if it is too light.
305subroutine hybgen_column_unmix(CS, nk, Rcv_tgt, temp, saln, Rcv, eqn_of_state, &
306 ntr, tracer, trcflg, fixlay, qhrlx, h_col, &
307 terrain_following, h_thin)
308 type(hybgen_unmix_cs), intent(in) :: CS !< hybgen unmixing control structure
309 integer, intent(in) :: nk !< The number of layers
310 integer, intent(in) :: fixlay !< deepest fixed coordinate layer
311 real, intent(in) :: qhrlx(nk+1) !< Relaxation fraction per timestep [nondim], < 1.
312 real, intent(in) :: Rcv_tgt(nk) !< Target potential density [R ~> kg m-3]
313 real, intent(inout) :: temp(nk) !< A column of potential temperature [C ~> degC]
314 real, intent(inout) :: saln(nk) !< A column of salinity [S ~> ppt]
315 real, intent(inout) :: Rcv(nk) !< Coordinate potential density [R ~> kg m-3]
316 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
317 integer, intent(in) :: ntr !< The number of registered passive tracers
318 real, intent(inout) :: tracer(nk, max(ntr,1)) !< Columns of the passive tracers [Conc]
319 integer, intent(in) :: trcflg(max(ntr,1)) !< Hycom tracer type flag for each tracer
320 real, intent(inout) :: h_col(nk+1) !< Layer thicknesses [H ~> m or kg m-2]
321 logical, intent(in) :: terrain_following !< True if this column is terrain following
322 real, intent(in) :: h_thin !< A negligibly small thickness to identify
323 !! essentially vanished layers [H ~> m or kg m-2]
324
325!
326! --- ------------------------------------------------------------------
327! --- hybrid grid generator, single column - ummix lowest massive layer.
328! --- ------------------------------------------------------------------
329!
330 ! Local variables
331 real :: h_hat ! A portion of a layer to move across an interface [H ~> m or kg m-2]
332 real :: delt, deltm ! Temperature differences between successive layers [C ~> degC]
333 real :: dels, delsm ! Salinity differences between successive layers [S ~> ppt]
334 real :: abs_dRdT ! The absolute value of the derivative of the coordinate density
335 ! with temperature [R C-1 ~> kg m-3 degC-1]
336 real :: abs_dRdS ! The absolute value of the derivative of the coordinate density
337 ! with salinity [R S-1 ~> kg m-3 ppt-1]
338 real :: q, qts ! Nondimensional fractions in the range of 0 to 1 [nondim]
339 real :: frac_dts ! The fraction of the temperature or salinity difference between successive
340 ! layers by which the source layer's property changes by the loss of water
341 ! that matches the destination layers properties via unmixing [nondim].
342 real :: qtr ! The fraction of the water that will come from the layer below,
343 ! used for updating the concentration of passive tracers [nondim]
344 real :: swap_T ! A swap variable for temperature [C ~> degC]
345 real :: swap_S ! A swap variable for salinity [S ~> ppt]
346 real :: swap_tr ! A temporary swap variable for the tracers [conc]
347 logical, parameter :: lunmix=.true. ! unmix a too light deepest layer
348 integer :: k, ka, kp, kt, m
349
350 ! --- identify the deepest layer kp with significant thickness (> h_thin)
351 kp = 2 !minimum allowed value
352 do k=nk,3,-1
353 if (h_col(k) >= h_thin) then
354 kp = k
355 exit
356 endif
357 enddo !k
358
359 k = kp !at least 2
360 ka = max(k-2,1) !k might be 2
361!
362 if ( ((k > fixlay+1) .and. (.not.terrain_following)) .and. & ! layer not fixed depth
363 (h_col(k-1) >= h_thin) .and. & ! layer above not too thin
364 (rcv_tgt(k) > rcv(k)) .and. & ! layer is lighter than its target
365 ((rcv(k-1) > rcv(k)) .and. (rcv(ka) > rcv(k))) ) then
366!
367! --- water in the deepest inflated layer with significant thickness
368! --- (kp) is too light, and it is lighter than the two layers above.
369! ---
370! --- this should only occur when relaxing or nudging layer thickness
371! --- and is a bug (bad interaction with tsadvc) even in those cases
372! ---
373! --- entrain the entire layer into the one above
374!--- note the double negative in T=T-q*(T-T'), equiv. to T=T+q*(T'-T)
375 q = h_col(k) / (h_col(k) + h_col(k-1))
376 temp(k-1) = temp(k-1) - q*(temp(k-1) - temp(k))
377 saln(k-1) = saln(k-1) - q*(saln(k-1) - saln(k))
378 call calculate_density(temp(k-1), saln(k-1), cs%ref_pressure, rcv(k-1), eqn_of_state)
379
380 do m=1,ntr
381 tracer(k-1,m) = tracer(k-1,m) - q*(tracer(k-1,m) - tracer(k,m) )
382 enddo !m
383! --- entrained the entire layer into the one above, so now kp=kp-1
384 h_col(k-1) = h_col(k-1) + h_col(k)
385 h_col(k) = 0.0
386 kp = k-1
387 elseif ( ((k > fixlay+1) .and. (.not.terrain_following)) .and. & ! layer not fixed depth
388 (h_col(k-1) >= h_thin) .and. & ! layer above not too thin
389 (rcv_tgt(k) > rcv(k)) .and. & ! layer is lighter than its target
390 (rcv(k-1) > rcv(k)) ) then
391! --- water in the deepest inflated layer with significant thickness
392! --- (kp) is too light, and it is lighter than the layer above, but not the layer two above.
393! ---
394! --- swap the entire layer with the one above.
395 if (h_col(k) <= h_col(k-1)) then
396 ! The bottom layer is thinner; swap the entire bottom layer with a portion of the layer above.
397 q = h_col(k) / h_col(k-1) !<=1.0
398
399 swap_t = temp(k-1)
400 temp(k-1) = temp(k-1) + q*(temp(k) - temp(k-1))
401 temp(k) = swap_t
402
403 swap_s = saln(k-1)
404 saln(k-1) = saln(k-1) + q*(saln(k) - saln(k-1))
405 saln(k) = swap_s
406
407 rcv(k) = rcv(k-1)
408 call calculate_density(temp(k-1), saln(k-1), cs%ref_pressure, rcv(k-1), eqn_of_state)
409
410 do m=1,ntr
411 swap_tr = tracer(k-1,m)
412 tracer(k-1,m) = tracer(k-1,m) - q * (tracer(k-1,m) - tracer(k,m))
413 tracer(k,m) = swap_tr
414 enddo !m
415 else
416 ! The bottom layer is thicker; swap the entire layer above with a portion of the bottom layer.
417 q = h_col(k-1) / h_col(k) !<1.0
418
419 swap_t = temp(k)
420 temp(k) = temp(k) + q*(temp(k-1) - temp(k))
421 temp(k-1) = swap_t
422
423 swap_s = saln(k)
424 saln(k) = saln(k) + q*(saln(k-1) - saln(k))
425 saln(k-1) = swap_s
426
427 rcv(k-1) = rcv(k)
428 call calculate_density(temp(k), saln(k), cs%ref_pressure, rcv(k), eqn_of_state)
429
430 do m=1,ntr
431 swap_tr = tracer(k,m)
432 tracer(k,m) = tracer(k,m) + q * (tracer(k-1,m) - tracer(k,m))
433 tracer(k-1,m) = swap_tr
434 enddo !m
435 endif !bottom too light
436 endif
437
438 k = kp !at least 2
439 ka = max(k-2,1) !k might be 2
440
441 if ( lunmix .and. & ! usually .true.
442 ((k > fixlay+1) .and. (.not.terrain_following)) .and. & ! layer not fixed depth
443 (h_col(k-1) >= h_thin) .and. & ! layer above not too thin
444 (rcv(k) < rcv_tgt(k)) .and. & ! layer is lighter than its target
445 (rcv(k) > rcv_tgt(k-1)) .and. & ! layer is denser than the target above
446 (abs(rcv_tgt(k-1) - rcv(k-1)) < cs%hybiso) .and. & ! layer above is near its target
447 (rcv(k) - rcv(k-1) > 0.001*(rcv_tgt(k) - rcv_tgt(k-1))) ) then
448!
449! --- water in the deepest inflated layer with significant thickness (kp) is too
450! --- light but denser than the layer above, with the layer above near-isopycnal
451! ---
452! --- split layer into 2 sublayers, one near the desired density
453! --- and one exactly matching the T&S properties of layer k-1.
454! --- To prevent "runaway" T or S, the result satisfies either
455! --- abs(T.k - T.k-1) <= abs(T.k-N - T.k-1) or
456! --- abs(S.k - S.k-1) <= abs(S.k-N - S.k-1) where
457! --- Rcv.k-1 - Rcv.k-N is at least Rcv_tgt(k-1) - Rcv_tgt(k-2)
458! --- It is also limited to a 50% change in layer thickness.
459
460 ka = 1
461 do kt=k-2,2,-1
462 if ( rcv(k-1) - rcv(kt) >= rcv_tgt(k-1) - rcv_tgt(k-2) ) then
463 ka = kt !usually k-2
464 exit
465 endif
466 enddo
467
468 delsm = abs(saln(ka) - saln(k-1))
469 dels = abs(saln(k-1) - saln(k))
470 deltm = abs(temp(ka) - temp(k-1))
471 delt = abs(temp(k-1) - temp(k))
472
473 call calculate_density_derivs(temp(k-1), saln(k-1), cs%ref_pressure, abs_drdt, abs_drds, eqn_of_state)
474 ! Bound deltm and delsm based on the equation of state and density differences between layers.
475 abs_drdt = abs(abs_drdt) ; abs_drds = abs(abs_drds)
476 if (abs_drdt * deltm > rcv_tgt(k)-rcv_tgt(k-1)) deltm = (rcv_tgt(k)-rcv_tgt(k-1)) / abs_drdt
477 if (abs_drds * delsm > rcv_tgt(k)-rcv_tgt(k-1)) delsm = (rcv_tgt(k)-rcv_tgt(k-1)) / abs_drds
478
479 qts = 0.0
480 if (qts*dels < min(delsm-dels, dels)) qts = min(delsm-dels, dels) / dels
481 if (qts*delt < min(deltm-delt, delt)) qts = min(deltm-delt, delt) / delt
482
483 ! Note that Rcv_tgt(k) > Rcv(k) > Rcv(k-1), and 0 <= qts <= 1.
484 ! qhrlx is relaxation coefficient (inverse baroclinic time steps), 0 <= qhrlx <= 1.
485 ! This takes the minimum of the two estimates.
486 if ((1.0+qts) * (rcv_tgt(k)-rcv(k)) < qts * (rcv_tgt(k)-rcv(k-1))) then
487 q = qhrlx(k) * ((rcv_tgt(k)-rcv(k)) / (rcv_tgt(k)-rcv(k-1)))
488 else
489 q = qhrlx(k) * (qts / (1.0+qts)) ! upper sublayer <= 50% of total
490 endif
491 frac_dts = q / (1.0-q) ! 0 <= q <= 0.5, so 0 <= frac_dts <= 1
492
493 h_hat = q * h_col(k)
494 h_col(k-1) = h_col(k-1) + h_hat
495 h_col(k) = h_col(k) - h_hat
496
497 temp(k) = temp(k) + frac_dts * (temp(k) - temp(k-1))
498 saln(k) = saln(k) + frac_dts * (saln(k) - saln(k-1))
499 call calculate_density(temp(k), saln(k), cs%ref_pressure, rcv(k), eqn_of_state)
500
501 if ((ntr > 0) .and. (h_hat /= 0.0)) then
502 ! qtr is the fraction of the new upper layer from the old lower layer.
503 ! The nonconservative original from Hycom: qtr = h_hat / max(h_hat, h_col(k)) !between 0 and 1
504 qtr = h_hat / h_col(k-1) ! Between 0 and 1, noting the h_col(k-1) = h_col(k-1) + h_hat above.
505 do m=1,ntr
506 if (trcflg(m) == 2) then !temperature tracer
507 tracer(k,m) = tracer(k,m) + frac_dts * (tracer(k,m) - tracer(k-1,m))
508 else !standard tracer - not split into two sub-layers
509 tracer(k-1,m) = tracer(k-1,m) + qtr * (tracer(k,m) - tracer(k-1,m))
510 endif !trcflg
511 enddo !m
512 endif !tracers
513 endif !too light
514
515! ! Fill properties of massless or near-massless (thickness < h_thin) layers
516! ! This was in the Hycom verion, but it appears to be unnecessary in MOM6.
517! do k=kp+1,nk
518! ! --- fill thin and massless layers on sea floor with fluid from above
519! Rcv(k) = Rcv(k-1)
520! do m=1,ntr
521! tracer(k,m) = tracer(k-1,m)
522! enddo !m
523! saln(k) = saln(k-1)
524! temp(k) = temp(k-1)
525! enddo !k
526
527end subroutine hybgen_column_unmix
528
529end module mom_hybgen_unmix