MOM6
coord_slight.F90
Go to the documentation of this file.
1 !> Regrid columns for the SLight coordinate
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_error_handler, only : mom_error, fatal
11 
12 implicit none ; private
13 
14 !> Control structure containing required parameters for the SLight coordinate
15 type, public :: slight_cs ; private
16 
17  !> Number of layers/levels
18  integer :: nk
19 
20  !> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2]
21  real :: min_thickness
22 
23  !> Reference pressure for potential density calculations [Pa]
24  real :: ref_pressure
25 
26  !> Fraction (between 0 and 1) of compressibility to add to potential density
27  !! profiles when interpolating for target grid positions. [nondim]
28  real :: compressibility_fraction
29 
30  ! The following 4 parameters were introduced for use with the SLight coordinate:
31  !> Depth over which to average to determine the mixed layer potential density [H ~> m or kg m-2]
32  real :: rho_ml_avg_depth
33 
34  !> Number of layers to offset the mixed layer density to find resolved stratification [nondim]
35  real :: nlay_ml_offset
36 
37  !> The number of fixed-thickness layers at the top of the model
38  integer :: nz_fixed_surface = 2
39 
40  !> The fixed resolution in the topmost SLight_nkml_min layers [H ~> m or kg m-2]
41  real :: dz_ml_min
42 
43  !> If true, detect regions with much weaker stratification in the coordinate
44  !! than based on in-situ density, and use a stretched coordinate there.
45  logical :: fix_haloclines = .false.
46 
47  !> A length scale over which to filter T & S when looking for spuriously
48  !! unstable water mass profiles [H ~> m or kg m-2].
49  real :: halocline_filter_length
50 
51  !> A value of the stratification ratio that defines a problematic halocline region [nondim].
52  real :: halocline_strat_tol
53 
54  !> Nominal density of interfaces [R ~> kg m-3].
55  real, allocatable, dimension(:) :: target_density
56 
57  !> Density scaling factor [R m3 kg-1 ~> 1]
58  real :: kg_m3_to_r
59 
60  !> Maximum depths of interfaces [H ~> m or kg m-2].
61  real, allocatable, dimension(:) :: max_interface_depths
62 
63  !> Maximum thicknesses of layers [H ~> m or kg m-2].
64  real, allocatable, dimension(:) :: max_layer_thickness
65 
66  !> Interpolation control structure
67  type(interp_cs_type) :: interp_cs
68 end type slight_cs
69 
71 
72 contains
73 
74 !> Initialise a slight_CS with pointers to parameters
75 subroutine init_coord_slight(CS, nk, ref_pressure, target_density, interp_CS, m_to_H, rho_scale)
76  type(slight_cs), pointer :: cs !< Unassociated pointer to hold the control structure
77  integer, intent(in) :: nk !< Number of layers in the grid
78  real, intent(in) :: ref_pressure !< Coordinate reference pressure [Pa]
79  real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [kg m-3]
80  type(interp_cs_type), intent(in) :: interp_cs !< Controls for interpolation
81  real, optional, intent(in) :: m_to_h !< A conversion factor from m to the units of thicknesses
82  real, optional, intent(in) :: rho_scale !< A dimensional scaling factor for target_density
83 
84  real :: m_to_h_rescale ! A unit conversion factor.
85 
86  if (associated(cs)) call mom_error(fatal, "init_coord_slight: CS already associated!")
87  allocate(cs)
88  allocate(cs%target_density(nk+1))
89 
90  m_to_h_rescale = 1.0 ; if (present(m_to_h)) m_to_h_rescale = m_to_h
91 
92  cs%nk = nk
93  cs%ref_pressure = ref_pressure
94  cs%target_density(:) = target_density(:)
95  cs%interp_CS = interp_cs
96 
97  ! Set real parameter default values
98  cs%compressibility_fraction = 0. ! Nondim.
99  cs%Rho_ML_avg_depth = 1.0 * m_to_h_rescale
100  cs%nlay_ml_offset = 2.0 ! Nondim.
101  cs%dz_ml_min = 1.0 * m_to_h_rescale
102  cs%halocline_filter_length = 2.0 * m_to_h_rescale
103  cs%halocline_strat_tol = 0.25 ! Nondim.
104  cs%kg_m3_to_R = 1.0 ; if (present(rho_scale)) cs%kg_m3_to_R = rho_scale
105 
106 end subroutine init_coord_slight
107 
108 !> This subroutine deallocates memory in the control structure for the coord_slight module
109 subroutine end_coord_slight(CS)
110  type(slight_cs), pointer :: cs !< Coordinate control structure
111 
112  ! nothing to do
113  if (.not. associated(cs)) return
114  deallocate(cs%target_density)
115  deallocate(cs)
116 end subroutine end_coord_slight
117 
118 !> This subroutine can be used to set the parameters for the coord_slight module
119 subroutine set_slight_params(CS, max_interface_depths, max_layer_thickness, &
120  min_thickness, compressibility_fraction, dz_ml_min, &
121  nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_offset, fix_haloclines, &
122  halocline_filter_length, halocline_strat_tol, interp_CS)
123  type(slight_cs), pointer :: cs !< Coordinate control structure
124  real, dimension(:), &
125  optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces [H ~> m or kg m-2]
126  real, dimension(:), &
127  optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers [H ~> m or kg m-2]
128  real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the
129  !! new grid through regridding [H ~> m or kg m-2]
130  real, optional, intent(in) :: compressibility_fraction !< Fraction (between 0 and 1) of
131  !! compressibility to add to potential density profiles when
132  !! interpolating for target grid positions. [nondim]
133  real, optional, intent(in) :: dz_ml_min !< The fixed resolution in the topmost
134  !! SLight_nkml_min layers [H ~> m or kg m-2]
135  integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickness layers at the
136  !! top of the model
137  real, optional, intent(in) :: rho_ml_avg_depth !< Depth over which to average to determine
138  !! the mixed layer potential density [H ~> m or kg m-2]
139  real, optional, intent(in) :: nlay_ml_offset !< Number of layers to offset the mixed layer
140  !! density to find resolved stratification [nondim]
141  logical, optional, intent(in) :: fix_haloclines !< If true, detect regions with much weaker than
142  !! based on in-situ density, and use a stretched coordinate there.
143  real, optional, intent(in) :: halocline_filter_length !< A length scale over which to filter T & S
144  !! when looking for spuriously unstable water mass profiles [H ~> m or kg m-2].
145  real, optional, intent(in) :: halocline_strat_tol !< A value of the stratification ratio that
146  !! defines a problematic halocline region [nondim].
147  type(interp_cs_type), &
148  optional, intent(in) :: interp_cs !< Controls for interpolation
149 
150  if (.not. associated(cs)) call mom_error(fatal, "set_slight_params: CS not associated")
151 
152  if (present(max_interface_depths)) then
153  if (size(max_interface_depths) /= cs%nk+1) &
154  call mom_error(fatal, "set_slight_params: max_interface_depths inconsistent size")
155  allocate(cs%max_interface_depths(cs%nk+1))
156  cs%max_interface_depths(:) = max_interface_depths(:)
157  endif
158 
159  if (present(max_layer_thickness)) then
160  if (size(max_layer_thickness) /= cs%nk) &
161  call mom_error(fatal, "set_slight_params: max_layer_thickness inconsistent size")
162  allocate(cs%max_layer_thickness(cs%nk))
163  cs%max_layer_thickness(:) = max_layer_thickness(:)
164  endif
165 
166  if (present(min_thickness)) cs%min_thickness = min_thickness
167  if (present(compressibility_fraction)) cs%compressibility_fraction = compressibility_fraction
168 
169  if (present(dz_ml_min)) cs%dz_ml_min = dz_ml_min
170  if (present(nz_fixed_surface)) cs%nz_fixed_surface = nz_fixed_surface
171  if (present(rho_ml_avg_depth)) cs%Rho_ML_avg_depth = rho_ml_avg_depth
172  if (present(nlay_ml_offset)) cs%nlay_ML_offset = nlay_ml_offset
173  if (present(fix_haloclines)) cs%fix_haloclines = fix_haloclines
174  if (present(halocline_filter_length)) cs%halocline_filter_length = halocline_filter_length
175  if (present(halocline_strat_tol)) then
176  if (halocline_strat_tol > 1.0) call mom_error(fatal, "set_slight_params: "//&
177  "HALOCLINE_STRAT_TOL must not exceed 1.0.")
178  cs%halocline_strat_tol = halocline_strat_tol
179  endif
180 
181  if (present(interp_cs)) cs%interp_CS = interp_cs
182 end subroutine set_slight_params
183 
184 !> Build a SLight coordinate column
185 subroutine build_slight_column(CS, eqn_of_state, H_to_Pa, H_subroundoff, &
186  nz, depth, h_col, T_col, S_col, p_col, z_col, z_col_new, &
187  h_neglect, h_neglect_edge)
188  type(slight_cs), intent(in) :: cs !< Coordinate control structure
189  type(eos_type), pointer :: eqn_of_state !< Equation of state structure
190  real, intent(in) :: h_to_pa !< GV%H_to_Pa
191  real, intent(in) :: h_subroundoff !< GV%H_subroundoff
192  integer, intent(in) :: nz !< Number of levels
193  real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2])
194  real, dimension(nz), intent(in) :: t_col !< T for column
195  real, dimension(nz), intent(in) :: s_col !< S for column
196  real, dimension(nz), intent(in) :: h_col !< Layer thicknesses [H ~> m or kg m-2]
197  real, dimension(nz), intent(in) :: p_col !< Layer quantities
198  real, dimension(nz+1), intent(in) :: z_col !< Interface positions relative to the surface [H ~> m or kg m-2]
199  real, dimension(nz+1), intent(inout) :: z_col_new !< Absolute positions of interfaces [H ~> m or kg m-2]
200  real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of
201  !! cell reconstructions [H ~> m or kg m-2].
202  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose
203  !! of edge value calculations [H ~> m or kg m-2].
204  ! Local variables
205  real, dimension(nz) :: rho_col ! Layer densities [R ~> kg m-3]
206  real, dimension(nz) :: t_f, s_f ! Filtered ayer quantities
207  logical, dimension(nz+1) :: reliable ! If true, this interface is in a reliable position.
208  real, dimension(nz+1) :: t_int, s_int ! Temperature and salinity interpolated to interfaces.
209  real, dimension(nz+1) :: rho_tmp ! A temporary density [R ~> kg m-3]
210  real, dimension(nz+1) :: drho_dp ! The partial derivative of density with pressure [kg m-3 Pa-1]
211  real, dimension(nz+1) :: p_is, p_r
212  real, dimension(nz+1) :: drhois_dt ! The partial derivative of in situ density with temperature
213  ! in [R degC-1 ~> kg m-3 degC-1]
214  real, dimension(nz+1) :: drhois_ds ! The partial derivative of in situ density with salinity
215  ! in [R ppt-1 ~> kg m-3 ppt-1]
216  real, dimension(nz+1) :: drhor_dt ! The partial derivative of reference density with temperature
217  ! in [R degC-1 ~> kg m-3 degC-1]
218  real, dimension(nz+1) :: drhor_ds ! The partial derivative of reference density with salinity
219  ! in [R ppt-1 ~> kg m-3 ppt-1]
220  real, dimension(nz+1) :: strat_rat
221  real :: h_to_cpa
222  real :: dris, drr ! In situ and reference density differences [R ~> kg m-3]
223  real :: fn_now, i_hstol, fn_zero_val
224  real :: z_int_unst
225  real :: dz ! A uniform layer thickness in very shallow water [H ~> m or kg m-2].
226  real :: dz_ur ! The total thickness of an unstable region [H ~> m or kg m-2].
227  real :: wgt, cowgt ! A weight and its complement, nondim.
228  real :: rho_ml_av ! The average potential density in a near-surface region [R ~> kg m-3].
229  real :: h_ml_av ! A thickness to try to use in taking the near-surface average [H ~> m or kg m-2].
230  real :: rho_x_z ! A cumulative integral of a density [R H ~> kg m-2 or kg2 m-5].
231  real :: z_wt ! The thickness actually used in taking the near-surface average [H ~> m or kg m-2].
232  real :: k_interior ! The (real) value of k where the interior grid starts.
233  real :: k_int2 ! The (real) value of k where the interior grid starts.
234  real :: z_interior ! The depth where the interior grid starts [H ~> m or kg m-2].
235  real :: z_ml_fix ! The depth at which the fixed-thickness near-surface layers end [H ~> m or kg m-2].
236  real :: dz_dk ! The thickness of layers between the fixed-thickness
237  ! near-surface layars and the interior [H ~> m or kg m-2].
238  real :: lfilt ! A filtering lengthscale [H ~> m or kg m-2].
239  logical :: maximum_depths_set ! If true, the maximum depths of interface have been set.
240  logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set.
241  real :: k2_used, k2here, dz_sum, z_max
242  integer :: k2
243  real :: h_tr, b_denom_1, b1, d1 ! Temporary variables used by the tridiagonal solver.
244  real, dimension(nz) :: c1 ! Temporary variables used by the tridiagonal solver.
245  integer :: kur1, kur2 ! The indicies at the top and bottom of an unreliable region.
246  integer :: kur_ss ! The index to start with in the search for the next unstable region.
247  integer :: i, j, k, nkml
248 
249  maximum_depths_set = allocated(cs%max_interface_depths)
250  maximum_h_set = allocated(cs%max_layer_thickness)
251 
252  if (z_col(nz+1) - z_col(1) < nz*cs%min_thickness) then
253  ! This is a nearly massless total depth, so distribute the water evenly.
254  dz = (z_col(nz+1) - z_col(1)) / real(nz)
255  do k=2,nz ; z_col_new(k) = z_col(1) + dz*real(k-1) ; enddo
256  else
257  call calculate_density(t_col, s_col, p_col, rho_col, 1, nz, &
258  eqn_of_state, scale=cs%kg_m3_to_R)
259 
260  ! Find the locations of the target potential densities, flagging
261  ! locations in apparently unstable regions as not reliable.
262  call rho_interfaces_col(rho_col, h_col, z_col, cs%target_density, nz, &
263  z_col_new, cs, reliable, debug=.true., &
264  h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
265 
266  ! Ensure that the interfaces are at least CS%min_thickness apart.
267  if (cs%min_thickness > 0.0) then
268  ! Move down interfaces below overly thin layers.
269  do k=2,nz ; if (z_col_new(k) < z_col_new(k-1) + cs%min_thickness) then
270  z_col_new(k) = z_col_new(k-1) + cs%min_thickness
271  endif ; enddo
272  ! Now move up any interfaces that are too close to the bottom.
273  do k=nz,2,-1 ; if (z_col_new(k) > z_col_new(k+1) - cs%min_thickness) then
274  z_col_new(k) = z_col_new(k+1) - cs%min_thickness
275  else
276  exit ! No more interfaces can be too close to the bottom.
277  endif ; enddo
278  endif
279 
280  ! Fix up the unreliable regions.
281  kur_ss = 2 ! reliable(1) and reliable(nz+1) must always be true.
282  do
283  ! Search for the uppermost unreliable interface postion.
284  kur1 = nz+2
285  do k=kur_ss,nz ; if (.not.reliable(k)) then
286  kur1 = k ; exit
287  endif ; enddo
288  if (kur1 > nz) exit ! Everything is now reliable.
289 
290  kur2 = kur1-1 ! For error checking.
291  do k=kur1+1,nz+1 ; if (reliable(k)) then
292  kur2 = k-1 ; kur_ss = k ; exit
293  endif ; enddo
294  if (kur2 < kur1) call mom_error(fatal, "Bad unreliable range.")
295 
296  dz_ur = z_col_new(kur2+1) - z_col_new(kur1-1)
297  ! drho = CS%target_density(kur2+1) - CS%target_density(kur1-1)
298  ! Perhaps reset the wgt and cowgt depending on how bad the old interface
299  ! locations were.
300  wgt = 1.0 ; cowgt = 0.0 ! = 1.0-wgt
301  do k=kur1,kur2
302  z_col_new(k) = cowgt*z_col_new(k) + &
303  wgt * (z_col_new(kur1-1) + dz_ur*(k - (kur1-1)) / ((kur2 - kur1) + 2))
304  enddo
305  enddo
306 
307  ! Determine which interfaces are in the s-space region and the depth extent
308  ! of this region.
309  z_wt = 0.0 ; rho_x_z = 0.0
310  h_ml_av = cs%Rho_ml_avg_depth
311  do k=1,nz
312  if (z_wt + h_col(k) >= h_ml_av) then
313  rho_x_z = rho_x_z + rho_col(k) * (h_ml_av - z_wt)
314  z_wt = h_ml_av
315  exit
316  else
317  rho_x_z = rho_x_z + rho_col(k) * h_col(k)
318  z_wt = z_wt + h_col(k)
319  endif
320  enddo
321  if (z_wt > 0.0) rho_ml_av = rho_x_z / z_wt
322 
323  nkml = cs%nz_fixed_surface
324  ! Find the interface that matches rho_ml_av.
325  if (rho_ml_av <= cs%target_density(nkml)) then
326  k_interior = cs%nlay_ml_offset + real(nkml)
327  elseif (rho_ml_av > cs%target_density(nz+1)) then
328  k_interior = real(nz+1)
329  else ; do k=nkml,nz
330  if ((rho_ml_av >= cs%target_density(k)) .and. &
331  (rho_ml_av < cs%target_density(k+1))) then
332  k_interior = (cs%nlay_ml_offset + k) + &
333  (rho_ml_av - cs%target_density(k)) / &
334  (cs%target_density(k+1) - cs%target_density(k))
335  exit
336  endif
337  enddo ; endif
338  if (k_interior > real(nz+1)) k_interior = real(nz+1)
339 
340  ! Linearly interpolate to find z_interior. This could be made more sophisticated.
341  k = int(ceiling(k_interior))
342  z_interior = (k-k_interior)*z_col_new(k-1) + (1.0+(k_interior-k))*z_col_new(k)
343 
344  if (cs%fix_haloclines) then
345  ! ! Identify regions above the reference pressure where the chosen
346  ! ! potential density significantly underestimates the actual
347  ! ! stratification, and use these to find a second estimate of
348  ! ! z_int_unst and k_interior.
349 
350  if (cs%halocline_filter_length > 0.0) then
351  lfilt = cs%halocline_filter_length
352 
353  ! Filter the temperature and salnity with a fixed lengthscale.
354  h_tr = h_col(1) + h_subroundoff
355  b1 = 1.0 / (h_tr + lfilt) ; d1 = h_tr * b1
356  t_f(1) = (b1*h_tr)*t_col(1) ; s_f(1) = (b1*h_tr)*s_col(1)
357  do k=2,nz
358  c1(k) = lfilt * b1
359  h_tr = h_col(k) + h_subroundoff ; b_denom_1 = h_tr + d1*lfilt
360  b1 = 1.0 / (b_denom_1 + lfilt) ; d1 = b_denom_1 * b1
361  t_f(k) = b1 * (h_tr*t_col(k) + lfilt*t_f(k-1))
362  s_f(k) = b1 * (h_tr*s_col(k) + lfilt*s_f(k-1))
363  enddo
364  do k=nz-1,1,-1
365  t_f(k) = t_f(k) + c1(k+1)*t_f(k+1) ; s_f(k) = s_f(k) + c1(k+1)*s_f(k+1)
366  enddo
367  else
368  do k=1,nz ; t_f(k) = t_col(k) ; s_f(k) = s_col(k) ; enddo
369  endif
370 
371  t_int(1) = t_f(1) ; s_int(1) = s_f(1)
372  do k=2,nz
373  t_int(k) = 0.5*(t_f(k-1) + t_f(k)) ; s_int(k) = 0.5*(s_f(k-1) + s_f(k))
374  p_is(k) = z_col(k) * h_to_pa
375  p_r(k) = cs%ref_pressure + cs%compressibility_fraction * ( p_is(k) - cs%ref_pressure )
376  enddo
377  t_int(nz+1) = t_f(nz) ; s_int(nz+1) = s_f(nz)
378  p_is(nz+1) = z_col(nz+1) * h_to_pa
379  call calculate_density_derivs(t_int, s_int, p_is, drhois_dt, drhois_ds, 2, nz-1, &
380  eqn_of_state, scale=cs%kg_m3_to_R)
381  call calculate_density_derivs(t_int, s_int, p_r, drhor_dt, drhor_ds, 2, nz-1, &
382  eqn_of_state, scale=cs%kg_m3_to_R)
383  if (cs%compressibility_fraction > 0.0) then
384  call calculate_compress(t_int, s_int, p_r, rho_tmp, drho_dp, 2, nz-1, &
385  eqn_of_state)
386  else
387  do k=2,nz ; drho_dp(k) = 0.0 ; enddo
388  endif
389 
390  h_to_cpa = cs%compressibility_fraction*cs%kg_m3_to_R*h_to_pa
391  strat_rat(1) = 1.0
392  do k=2,nz
393  dris = drhois_dt(k) * (t_f(k) - t_f(k-1)) + &
394  drhois_ds(k) * (s_f(k) - s_f(k-1))
395  drr = (drhor_dt(k) * (t_f(k) - t_f(k-1)) + &
396  drhor_ds(k) * (s_f(k) - s_f(k-1))) + &
397  drho_dp(k) * (h_to_cpa*0.5*(h_col(k) + h_col(k-1)))
398 
399  if (dris <= 0.0) then
400  strat_rat(k) = 2.0 ! Maybe do this? => ; if (drR < 0.0) strat_rat(K) = -2.0
401  else
402  strat_rat(k) = 2.0*max(drr,0.0) / (dris + abs(drr))
403  endif
404  enddo
405  strat_rat(nz+1) = 1.0
406 
407  z_int_unst = 0.0 ; fn_now = 0.0
408  fn_zero_val = min(2.0*cs%halocline_strat_tol, &
409  0.5*(1.0 + cs%halocline_strat_tol))
410  if (cs%halocline_strat_tol > 0.0) then
411  ! Use Adcroft's reciprocal rule.
412  i_hstol = 0.0 ; if (fn_zero_val - cs%halocline_strat_tol > 0.0) &
413  i_hstol = 1.0 / (fn_zero_val - cs%halocline_strat_tol)
414  do k=nz,1,-1 ; if (cs%ref_pressure > p_is(k+1)) then
415  z_int_unst = z_int_unst + fn_now * h_col(k)
416  if (strat_rat(k) <= fn_zero_val) then
417  if (strat_rat(k) <= cs%halocline_strat_tol) then ; fn_now = 1.0
418  else
419  fn_now = max(fn_now, (fn_zero_val - strat_rat(k)) * i_hstol)
420  endif
421  endif
422  endif ; enddo
423  else
424  do k=nz,1,-1 ; if (cs%ref_pressure > p_is(k+1)) then
425  z_int_unst = z_int_unst + fn_now * h_col(k)
426  if (strat_rat(k) <= cs%halocline_strat_tol) fn_now = 1.0
427  endif ; enddo
428  endif
429 
430  if (z_interior < z_int_unst) then
431  ! Find a second estimate of the extent of the s-coordinate region.
432  kur1 = max(int(ceiling(k_interior)),2)
433  if (z_col_new(kur1-1) < z_interior) then
434  k_int2 = kur1
435  do k = kur1,nz+1 ; if (z_col_new(k) >= z_int_unst) then
436  ! This is linear interpolation again.
437  if (z_col_new(k-1) >= z_int_unst) &
438  call mom_error(fatal,"build_grid_SLight, bad halocline structure.")
439  k_int2 = real(k-1) + (z_int_unst - z_col_new(k-1)) / &
440  (z_col_new(k) - z_col_new(k-1))
441  exit
442  endif ; enddo
443  if (z_col_new(nz+1) < z_int_unst) then
444  ! This should be unnecessary.
445  z_int_unst = z_col_new(nz+1) ; k_int2 = real(nz+1)
446  endif
447 
448  ! Now take the larger values.
449  if (k_int2 > k_interior) then
450  k_interior = k_int2 ; z_interior = z_int_unst
451  endif
452  endif
453  endif
454  endif ! fix_haloclines
455 
456  z_col_new(1) = 0.0
457  do k=2,nkml+1
458  z_col_new(k) = min((k-1)*cs%dz_ml_min, &
459  z_col_new(nz+1) - cs%min_thickness*(nz+1-k))
460  enddo
461  z_ml_fix = z_col_new(nkml+1)
462  if (z_interior > z_ml_fix) then
463  dz_dk = (z_interior - z_ml_fix) / (k_interior - (nkml+1))
464  do k=nkml+2,int(floor(k_interior))
465  z_col_new(k) = z_ml_fix + dz_dk * (k - (nkml+1))
466  enddo
467  else ! The fixed-thickness z-region penetrates into the interior.
468  do k=nkml+2,nz
469  if (z_col_new(k) <= z_col_new(cs%nz_fixed_surface+1)) then
470  z_col_new(k) = z_col_new(cs%nz_fixed_surface+1)
471  else ; exit ; endif
472  enddo
473  endif
474 
475  if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz
476  ! The loop bounds are 2 & nz so the top and bottom interfaces do not move.
477  ! Recall that z_col_new is positive downward.
478  z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k), &
479  z_col_new(k-1) + cs%max_layer_thickness(k-1))
480  enddo ; elseif (maximum_depths_set) then ; do k=2,nz
481  z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k))
482  enddo ; elseif (maximum_h_set) then ; do k=2,nz
483  z_col_new(k) = min(z_col_new(k), z_col_new(k-1) + cs%max_layer_thickness(k-1))
484  enddo ; endif
485 
486  endif ! Total thickness exceeds nz*CS%min_thickness.
487 
488 end subroutine build_slight_column
489 
490 !> Finds the new interface locations in a column of water that match the
491 !! prescribed target densities.
492 subroutine rho_interfaces_col(rho_col, h_col, z_col, rho_tgt, nz, z_col_new, &
493  CS, reliable, debug, h_neglect, h_neglect_edge)
494  integer, intent(in) :: nz !< Number of layers
495  real, dimension(nz), intent(in) :: rho_col !< Initial layer reference densities.
496  real, dimension(nz), intent(in) :: h_col !< Initial layer thicknesses.
497  real, dimension(nz+1), intent(in) :: z_col !< Initial interface heights.
498  real, dimension(nz+1), intent(in) :: rho_tgt !< Interface target densities.
499  real, dimension(nz+1), intent(inout) :: z_col_new !< New interface heights.
500  type(slight_cs), intent(in) :: CS !< Coordinate control structure
501  logical, dimension(nz+1), intent(inout) :: reliable !< If true, the interface positions
502  !! are well defined from a stable region.
503  logical, optional, intent(in) :: debug !< If present and true, do debugging checks.
504  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
505  !! purpose of cell reconstructions
506  !! in the same units as h_col.
507  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
508  !! for the purpose of edge value calculations
509  !! in the same units as h_col.
510 
511  real, dimension(nz+1) :: ru_max_int ! The maximum and minimum densities in
512  real, dimension(nz+1) :: ru_min_int ! an unstable region around an interface.
513  real, dimension(nz) :: ru_max_lay ! The maximum and minimum densities in
514  real, dimension(nz) :: ru_min_lay ! an unstable region containing a layer.
515  real, dimension(nz,2) :: ppoly_i_E ! Edge value of polynomial
516  real, dimension(nz,2) :: ppoly_i_S ! Edge slope of polynomial
517  real, dimension(nz,DEGREE_MAX+1) :: ppoly_i_coefficients ! Coefficients of polynomial
518  logical, dimension(nz) :: unstable_lay ! If true, this layer is in an unstable region.
519  logical, dimension(nz+1) :: unstable_int ! If true, this interface is in an unstable region.
520  real :: rt ! The current target density [kg m-3].
521  real :: zf ! The fractional z-position within a layer of the target density.
522  real :: rfn
523  real :: a(5) ! Coefficients of a local polynomial minus the target density.
524  real :: zf1, zf2, rfn1, rfn2
525  real :: drfn_dzf, sgn, delta_zf, zf_prev
526  real :: tol
527  logical :: k_found ! If true, the position has been found.
528  integer :: k_layer ! The index of the stable layer containing an interface.
529  integer :: ppoly_degree
530  integer :: k, k1, k1_min, itt, max_itt, m
531 
532  real :: z_sgn ! 1 or -1, depending on whether z increases with increasing K.
533  logical :: debugging
534 
535  debugging = .false. ; if (present(debug)) debugging = debug
536  max_itt = nr_iterations
537  tol = nr_tolerance
538 
539  z_sgn = 1.0 ; if ( z_col(1) > z_col(nz+1) ) z_sgn = -1.0
540  if (debugging) then
541  do k=1,nz
542  if (abs((z_col(k+1) - z_col(k)) - z_sgn*h_col(k)) > &
543  1.0e-14*(abs(z_col(k+1)) + abs(z_col(k)) + abs(h_col(k))) ) &
544  call mom_error(fatal, "rho_interfaces_col: Inconsistent z_col and h_col")
545  enddo
546  endif
547 
548  if ( z_col(1) == z_col(nz+1) ) then
549  ! This is a massless column!
550  do k=1,nz+1 ; z_col_new(k) = z_col(1) ; reliable(k) = .true. ; enddo
551  return
552  endif
553 
554  ! This sets up the piecewise polynomials based on the rho_col profile.
555  call regridding_set_ppolys(cs%interp_CS, rho_col, nz, h_col, ppoly_i_e, ppoly_i_s, &
556  ppoly_i_coefficients, ppoly_degree, h_neglect, h_neglect_edge)
557 
558  ! Determine the density ranges of unstably stratified segments.
559  ! Interfaces that start out in an unstably stratified segment can
560  ! only escape if they are outside of the bounds of that segment, and no
561  ! interfaces are ever mapped into an unstable segment.
562  unstable_int(1) = .false.
563  ru_max_int(1) = ppoly_i_e(1,1)
564 
565  unstable_lay(1) = (ppoly_i_e(1,1) > ppoly_i_e(1,2))
566  ru_max_lay(1) = max(ppoly_i_e(1,1), ppoly_i_e(1,2))
567 
568  do k=2,nz
569  unstable_int(k) = (ppoly_i_e(k-1,2) > ppoly_i_e(k,1))
570  ru_max_int(k) = max(ppoly_i_e(k-1,2), ppoly_i_e(k,1))
571  ru_min_int(k) = min(ppoly_i_e(k-1,2), ppoly_i_e(k,1))
572  if (unstable_int(k) .and. unstable_lay(k-1)) &
573  ru_max_int(k) = max(ru_max_lay(k-1), ru_max_int(k))
574 
575  unstable_lay(k) = (ppoly_i_e(k,1) > ppoly_i_e(k,2))
576  ru_max_lay(k) = max(ppoly_i_e(k,1), ppoly_i_e(k,2))
577  ru_min_lay(k) = min(ppoly_i_e(k,1), ppoly_i_e(k,2))
578  if (unstable_lay(k) .and. unstable_int(k)) &
579  ru_max_lay(k) = max(ru_max_int(k), ru_max_lay(k))
580  enddo
581  unstable_int(nz+1) = .false.
582  ru_min_int(nz+1) = ppoly_i_e(nz,2)
583 
584  do k=nz,1,-1
585  if (unstable_lay(k) .and. unstable_int(k+1)) &
586  ru_min_lay(k) = min(ru_min_int(k+1), ru_min_lay(k))
587 
588  if (unstable_int(k) .and. unstable_lay(k)) &
589  ru_min_int(k) = min(ru_min_lay(k), ru_min_int(k))
590  enddo
591 
592  z_col_new(1) = z_col(1) ; reliable(1) = .true.
593  k1_min = 1
594  do k=2,nz ! Find the locations of the various target densities for the interfaces.
595  rt = rho_tgt(k)
596  k_layer = -1
597  k_found = .false.
598 
599  ! Many light layers are found at the top, so start there.
600  if (rt <= ppoly_i_e(k1_min,1)) then
601  z_col_new(k) = z_col(k1_min)
602  k_found = .true.
603  ! Do not change k1_min for the next layer.
604  elseif (k1_min == nz+1) then
605  z_col_new(k) = z_col(nz+1)
606  else
607  ! Start with the previous location and search outward.
608  if (unstable_int(k) .and. (rt >= ru_min_int(k)) .and. (rt <= ru_max_int(k))) then
609  ! This interface started in an unstable region and should not move due to remapping.
610  z_col_new(k) = z_col(k) ; reliable(k) = .false.
611  k1_min = k ; k_found = .true.
612  elseif ((rt >= ppoly_i_e(k-1,2)) .and. (rt <= ppoly_i_e(k,1))) then
613  ! This interface is already in the right place and does not move.
614  z_col_new(k) = z_col(k) ; reliable(k) = .true.
615  k1_min = k ; k_found = .true.
616  elseif (rt < ppoly_i_e(k-1,2)) then ! Search upward
617  do k1=k-1,k1_min,-1
618  ! Check whether rt is in layer k.
619  if ((rt < ppoly_i_e(k1,2)) .and. (rt > ppoly_i_e(k1,1))) then
620  ! rt is in layer k.
621  k_layer = k1
622  k1_min = k1 ; k_found = .true. ; exit
623  elseif (unstable_lay(k1) .and. (rt >= ru_min_lay(k1)) .and. (rt <= ru_max_lay(k1))) then
624  ! rt would be found at unstable layer that it can not penetrate.
625  ! It is possible that this can never happen?
626  z_col_new(k) = z_col(k1+1) ; reliable(k) = .false.
627  k1_min = k1 ; k_found = .true. ; exit
628  endif
629  ! Check whether rt is at interface K.
630  if (k1 > 1) then ; if ((rt <= ppoly_i_e(k1,1)) .and. (rt >= ppoly_i_e(k1-1,2))) then
631  ! rt is at interface K1
632  z_col_new(k) = z_col(k1) ; reliable(k) = .true.
633  k1_min = k1 ; k_found = .true. ; exit
634  elseif (unstable_int(k1) .and. (rt >= ru_min_int(k1)) .and. (rt <= ru_max_int(k1))) then
635  ! rt would be found at an unstable interface that it can not pass.
636  ! It is possible that this can never happen?
637  z_col_new(k) = z_col(k1) ; reliable(k) = .false.
638  k1_min = k1 ; k_found = .true. ; exit
639  endif ; endif
640  enddo
641 
642  if (.not.k_found) then
643  ! This should not happen unless k1_min = 1.
644  if (k1_min < 2) then
645  z_col_new(k) = z_col(k1_min)
646  else
647  z_col_new(k) = z_col(k1_min)
648  endif
649  endif
650 
651  else ! Search downward
652  do k1=k,nz
653  if ((rt < ppoly_i_e(k1,2)) .and. (rt > ppoly_i_e(k1,1))) then
654  ! rt is in layer k.
655  k_layer = k1
656  k1_min = k1 ; k_found = .true. ; exit
657  elseif (unstable_lay(k1) .and. (rt >= ru_min_lay(k1)) .and. (rt <= ru_max_lay(k1))) then
658  ! rt would be found at unstable layer that it can not penetrate.
659  ! It is possible that this can never happen?
660  z_col_new(k) = z_col(k1)
661  reliable(k) = .false.
662  k1_min = k1 ; k_found = .true. ; exit
663  endif
664  if (k1 < nz) then ; if ((rt <= ppoly_i_e(k1+1,1)) .and. (rt >= ppoly_i_e(k1,2))) then
665  ! rt is at interface K1+1
666 
667  z_col_new(k) = z_col(k1+1) ; reliable(k) = .true.
668  k1_min = k1+1 ; k_found = .true. ; exit
669  elseif (unstable_int(k1+1) .and. (rt >= ru_min_int(k1+1)) .and. (rt <= ru_max_int(k1+1))) then
670  ! rt would be found at an unstable interface that it can not pass.
671  ! It is possible that this can never happen?
672  z_col_new(k) = z_col(k1+1)
673  reliable(k) = .false.
674  k1_min = k1+1 ; k_found = .true. ; exit
675  endif ; endif
676  enddo
677  if (.not.k_found) then
678  z_col_new(k) = z_col(nz+1)
679  if (rt >= ppoly_i_e(nz,2)) then
680  reliable(k) = .true.
681  else
682  reliable(k) = .false.
683  endif
684  endif
685  endif
686 
687  if (k_layer > 0) then ! The new location is inside of layer k_layer.
688  ! Note that this is coded assuming that this layer is stably stratified.
689  if (.not.(ppoly_i_e(k1,2) > ppoly_i_e(k1,1))) call mom_error(fatal, &
690  "build_grid_SLight: Erroneously searching for an interface in an unstratified layer.") !### COMMENT OUT LATER?
691 
692  ! Use the false position method to find the location (degree <= 1) or the first guess.
693  zf = (rt - ppoly_i_e(k1,1)) / (ppoly_i_e(k1,2) - ppoly_i_e(k1,1))
694 
695  if (ppoly_degree > 1) then ! Iterate to find the solution.
696  a(:) = 0.0 ; a(1) = ppoly_i_coefficients(k_layer,1) - rt
697  do m=2,ppoly_degree+1 ; a(m) = ppoly_i_coefficients(k_layer,m) ; enddo
698  ! Bracket the root.
699  zf1 = 0.0 ; rfn1 = a(1)
700  zf2 = 1.0 ; rfn2 = a(1) + (a(2) + (a(3) + (a(4) + a(5))))
701  if (rfn1 * rfn2 > 0.0) call mom_error(fatal, "build_grid_SLight: Bad bracketing.") !### COMMENT OUT LATER?
702 
703  do itt=1,max_itt
704  rfn = a(1) + zf*(a(2) + zf*(a(3) + zf*(a(4) + zf*a(5))))
705  ! Reset one of the ends of the bracket.
706  if (rfn * rfn1 > 0.0) then
707  zf1 = zf ; rfn1 = rfn
708  else
709  zf2 = zf ; rfn2 = rfn
710  endif
711  if (rfn1 == rfn2) exit
712 
713  drfn_dzf = (a(2) + zf*(2.0*a(3) + zf*(3.0*a(4) + zf*4.0*a(5))))
714  sgn = 1.0 ; if (drfn_dzf < 0.0) sgn = -1.0
715 
716  if ((sgn*(zf - rfn) >= zf1 * abs(drfn_dzf)) .and. &
717  (sgn*(zf - rfn) <= zf2 * abs(drfn_dzf))) then
718  delta_zf = -rfn / drfn_dzf
719  zf = zf + delta_zf
720  else ! Newton's method goes out of bounds, so use a false position method estimate
721  zf_prev = zf
722  zf = ( rfn2 * zf1 - rfn1 * zf2 ) / (rfn2 - rfn1)
723  delta_zf = zf - zf_prev
724  endif
725 
726  if (abs(delta_zf) < tol) exit
727  enddo
728  endif
729  z_col_new(k) = z_col(k_layer) + zf * z_sgn * h_col(k_layer)
730  reliable(k) = .true.
731  endif
732 
733  endif
734 
735  enddo
736  z_col_new(nz+1) = z_col(nz+1) ; reliable(nz+1) = .true.
737 
738 end subroutine rho_interfaces_col
739 
740 end module coord_slight
coord_slight::slight_cs
Control structure containing required parameters for the SLight coordinate.
Definition: coord_slight.F90:15
coord_slight::end_coord_slight
subroutine, public end_coord_slight(CS)
This subroutine deallocates memory in the control structure for the coord_slight module.
Definition: coord_slight.F90:110
mom_eos::calculate_compress
Calculates the compressibility of water from T, S, and P.
Definition: MOM_EOS.F90:86
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
regrid_interp
Vertical interpolation for regridding.
Definition: regrid_interp.F90:2
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
regrid_interp::nr_iterations
integer, parameter, public nr_iterations
Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are used to build the new grid...
Definition: regrid_interp.F90:66
coord_slight::set_slight_params
subroutine, public set_slight_params(CS, max_interface_depths, max_layer_thickness, min_thickness, compressibility_fraction, dz_ml_min, nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_offset, fix_haloclines, halocline_filter_length, halocline_strat_tol, interp_CS)
This subroutine can be used to set the parameters for the coord_slight module.
Definition: coord_slight.F90:123
regrid_interp::nr_tolerance
real, parameter, public nr_tolerance
Tolerance for Newton-Raphson iterations (stop when increment falls below this)
Definition: regrid_interp.F90:68
coord_slight::build_slight_column
subroutine, public build_slight_column(CS, eqn_of_state, H_to_Pa, H_subroundoff, nz, depth, h_col, T_col, S_col, p_col, z_col, z_col_new, h_neglect, h_neglect_edge)
Build a SLight coordinate column.
Definition: coord_slight.F90:188
regrid_interp::degree_max
integer, parameter, public degree_max
Interpolant degrees.
Definition: regrid_interp.F90:56
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
coord_slight::init_coord_slight
subroutine, public init_coord_slight(CS, nk, ref_pressure, target_density, interp_CS, m_to_H, rho_scale)
Initialise a slight_CS with pointers to parameters.
Definition: coord_slight.F90:76
regrid_interp::interp_cs_type
Control structure for regrid_interp module.
Definition: regrid_interp.F90:23
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
coord_slight::rho_interfaces_col
subroutine rho_interfaces_col(rho_col, h_col, z_col, rho_tgt, nz, z_col_new, CS, reliable, debug, h_neglect, h_neglect_edge)
Finds the new interface locations in a column of water that match the prescribed target densities.
Definition: coord_slight.F90:494
regrid_interp::regridding_set_ppolys
subroutine, public regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, ppoly0_coefs, degree, h_neglect, h_neglect_edge)
Builds an interpolated profile for the densities within each grid cell.
Definition: regrid_interp.F90:80
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
coord_slight
Regrid columns for the SLight coordinate.
Definition: coord_slight.F90:2
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60