MOM6
coord_rho.F90
Go to the documentation of this file.
1 !> Regrid columns for the continuous isopycnal (rho) coordinate
2 module coord_rho
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_error_handler, only : mom_error, fatal
10 
11 implicit none ; private
12 
13 !> Control structure containing required parameters for the rho coordinate
14 type, public :: rho_cs ; private
15 
16  !> Number of layers
17  integer :: nk
18 
19  !> Minimum thickness allowed for layers, often in [H ~> m or kg m-2]
20  real :: min_thickness = 0.
21 
22  !> Reference pressure for density calculations [Pa]
23  real :: ref_pressure
24 
25  !> If true, integrate for interface positions from the top downward.
26  !! If false, integrate from the bottom upward, as does the rest of the model.
27  logical :: integrate_downward_for_e = .false.
28 
29  !> Nominal density of interfaces [R ~> kg m-3]
30  real, allocatable, dimension(:) :: target_density
31 
32  !> Density scaling factor [R m3 kg-1 ~> 1]
33  real :: kg_m3_to_r
34 
35  !> Interpolation control structure
36  type(interp_cs_type) :: interp_cs
37 end type rho_cs
38 
39 !> Maximum number of regridding iterations
40 integer, parameter :: nb_regridding_iterations = 1
41 !> Deviation tolerance between succesive grids in regridding iterations
42 real, parameter :: deviation_tolerance = 1e-10
43 
45 
46 contains
47 
48 !> Initialise a rho_CS with pointers to parameters
49 subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS, rho_scale)
50  type(rho_cs), pointer :: cs !< Unassociated pointer to hold the control structure
51  integer, intent(in) :: nk !< Number of layers in the grid
52  real, intent(in) :: ref_pressure !< Coordinate reference pressure [Pa]
53  real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [kg m-3 or R ~> kg m-3]
54  type(interp_cs_type), intent(in) :: interp_cs !< Controls for interpolation
55  real, optional, intent(in) :: rho_scale !< A dimensional scaling factor for target_density
56 
57  if (associated(cs)) call mom_error(fatal, "init_coord_rho: CS already associated!")
58  allocate(cs)
59  allocate(cs%target_density(nk+1))
60 
61  cs%nk = nk
62  cs%ref_pressure = ref_pressure
63  cs%target_density(:) = target_density(:)
64  cs%interp_CS = interp_cs
65  cs%kg_m3_to_R = 1.0 ; if (present(rho_scale)) cs%kg_m3_to_R = rho_scale
66 
67 end subroutine init_coord_rho
68 
69 !> This subroutine deallocates memory in the control structure for the coord_rho module
70 subroutine end_coord_rho(CS)
71  type(rho_cs), pointer :: cs !< Coordinate control structure
72 
73  ! nothing to do
74  if (.not. associated(cs)) return
75  deallocate(cs%target_density)
76  deallocate(cs)
77 end subroutine end_coord_rho
78 
79 !> This subroutine can be used to set the parameters for the coord_rho module
80 subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS)
81  type(rho_cs), pointer :: cs !< Coordinate control structure
82  real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
83  logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface
84  !! positions from the top downward. If false, integrate
85  !! from the bottom upward, as does the rest of the model.
86 
87  type(interp_cs_type), optional, intent(in) :: interp_cs !< Controls for interpolation
88 
89  if (.not. associated(cs)) call mom_error(fatal, "set_rho_params: CS not associated")
90 
91  if (present(min_thickness)) cs%min_thickness = min_thickness
92  if (present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
93  if (present(interp_cs)) cs%interp_CS = interp_cs
94 end subroutine set_rho_params
95 
96 !> Build a rho coordinate column
97 !!
98 !! 1. Density profiles are calculated on the source grid.
99 !! 2. Positions of target densities (for interfaces) are found by interpolation.
100 subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, &
101  h_neglect, h_neglect_edge)
102  type(rho_cs), intent(in) :: cs !< coord_rho control structure
103  integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S)
104  real, intent(in) :: depth !< Depth of ocean bottom (positive in m)
105  real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
106  real, dimension(nz), intent(in) :: t !< T for source column
107  real, dimension(nz), intent(in) :: s !< S for source column
108  type(eos_type), pointer :: eqn_of_state !< Equation of state structure
109  real, dimension(CS%nk+1), &
110  intent(inout) :: z_interface !< Absolute positions of interfaces
111  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
112  !! purpose of cell reconstructions
113  !! in the same units as h0.
114  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
115  !! for the purpose of edge value calculations
116  !! in the same units as h0.
117  ! Local variables
118  integer :: k, count_nonzero_layers
119  integer, dimension(nz) :: mapping
120  real, dimension(nz) :: p, h_nv
121  real, dimension(nz) :: densities ! Layer density [R ~> kg m-3]
122  real, dimension(nz+1) :: xtmp
123  real, dimension(CS%nk) :: h_new ! New thicknesses
124  real, dimension(CS%nk+1) :: x1
125 
126  ! Construct source column with vanished layers removed (stored in h_nv)
127  call copy_finite_thicknesses(nz, h, cs%min_thickness, count_nonzero_layers, h_nv, mapping)
128 
129  if (count_nonzero_layers > 1) then
130  xtmp(1) = 0.0
131  do k = 1,count_nonzero_layers
132  xtmp(k+1) = xtmp(k) + h_nv(k)
133  enddo
134 
135  ! Compute densities on source column
136  p(:) = cs%ref_pressure
137  call calculate_density(t, s, p, densities, 1, nz, eqn_of_state, scale=cs%kg_m3_to_R)
138  do k = 1,count_nonzero_layers
139  densities(k) = densities(mapping(k))
140  enddo
141 
142  ! Based on source column density profile, interpolate to generate a new grid
143  call build_and_interpolate_grid(cs%interp_CS, densities, count_nonzero_layers, &
144  h_nv, xtmp, cs%target_density, cs%nk, h_new, &
145  x1, h_neglect, h_neglect_edge)
146 
147  ! Inflate vanished layers
148  call old_inflate_layers_1d(cs%min_thickness, cs%nk, h_new)
149 
150  ! Comment: The following adjustment of h_new, and re-calculation of h_new via x1 needs to be removed
151  x1(1) = 0.0 ; do k = 1,cs%nk ; x1(k+1) = x1(k) + h_new(k) ; enddo
152  do k = 1,cs%nk
153  h_new(k) = x1(k+1) - x1(k)
154  enddo
155 
156  else ! count_nonzero_layers <= 1
157  if (nz == cs%nk) then
158  h_new(:) = h(:) ! This keeps old behavior
159  else
160  h_new(:) = 0.
161  h_new(1) = h(1)
162  endif
163  endif
164 
165  ! Return interface positions
166  if (cs%integrate_downward_for_e) then
167  ! Remapping is defined integrating from zero
168  z_interface(1) = 0.
169  do k = 1,cs%nk
170  z_interface(k+1) = z_interface(k) - h_new(k)
171  enddo
172  else
173  ! The rest of the model defines grids integrating up from the bottom
174  z_interface(cs%nk+1) = -depth
175  do k = cs%nk,1,-1
176  z_interface(k) = z_interface(k+1) + h_new(k)
177  enddo
178  endif
179 
180 end subroutine build_rho_column
181 
182 !> Iteratively build a rho coordinate column
183 !!
184 !! The algorithm operates as follows within each column:
185 !!
186 !! 1. Given T & S within each layer, the layer densities are computed.
187 !! 2. Based on these layer densities, a global density profile is reconstructed
188 !! (this profile is monotonically increasing and may be discontinuous)
189 !! 3. The new grid interfaces are determined based on the target interface
190 !! densities.
191 !! 4. T & S are remapped onto the new grid.
192 !! 5. Return to step 1 until convergence or until the maximum number of
193 !! iterations is reached, whichever comes first.
194 subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_state, &
195  zInterface, h_neglect, h_neglect_edge)
196  type(rho_cs), intent(in) :: CS !< Regridding control structure
197  type(remapping_cs), intent(in) :: remapCS !< Remapping parameters and options
198  integer, intent(in) :: nz !< Number of levels
199  real, intent(in) :: depth !< Depth of ocean bottom [Z ~> m]
200  real, dimension(nz), intent(in) :: h !< Layer thicknesses in Z coordinates [Z ~> m]
201  real, dimension(nz), intent(in) :: T !< T for column [degC]
202  real, dimension(nz), intent(in) :: S !< S for column [ppt]
203  type(eos_type), pointer :: eqn_of_state !< Equation of state structure
204  real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces
205  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
206  !! purpose of cell reconstructions
207  !! in the same units as h [Z ~> m]
208  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
209  !! for the purpose of edge value calculations
210  !! in the same units as h [Z ~> m]
211  ! Local variables
212  integer :: k, m
213  integer :: count_nonzero_layers
214  real :: deviation ! When iterating to determine the final
215  ! grid, this is the deviation between two
216  ! successive grids.
217  real :: threshold
218  real, dimension(nz) :: p, densities, T_tmp, S_tmp, Tmp
219  integer, dimension(nz) :: mapping
220  real, dimension(nz) :: h0, h1, hTmp
221  real, dimension(nz+1) :: x0, x1, xTmp
222 
223  threshold = cs%min_thickness
224  p(:) = cs%ref_pressure
225  t_tmp(:) = t(:)
226  s_tmp(:) = s(:)
227  h0(:) = h(:)
228 
229  ! Start iterations to build grid
230  m = 1
231  deviation = 1e10
232  do while ( ( m <= nb_regridding_iterations ) .and. &
233  ( deviation > deviation_tolerance ) )
234 
235  ! Construct column with vanished layers removed
236  call copy_finite_thicknesses(nz, h0, threshold, count_nonzero_layers, htmp, mapping)
237  if ( count_nonzero_layers <= 1 ) then
238  h1(:) = h0(:)
239  exit ! stop iterations here
240  endif
241 
242  xtmp(1) = 0.0
243  do k = 1,count_nonzero_layers
244  xtmp(k+1) = xtmp(k) + htmp(k)
245  enddo
246 
247  ! Compute densities within current water column
248  call calculate_density( t_tmp, s_tmp, p, densities, &
249  1, nz, eqn_of_state, scale=cs%kg_m3_to_R)
250 
251  do k = 1,count_nonzero_layers
252  densities(k) = densities(mapping(k))
253  enddo
254 
255  ! One regridding iteration
256  ! Based on global density profile, interpolate to generate a new grid
257  call build_and_interpolate_grid(cs%interp_CS, densities, count_nonzero_layers, &
258  htmp, xtmp, cs%target_density, nz, h1, x1, h_neglect, h_neglect_edge)
259 
260  call old_inflate_layers_1d( cs%min_thickness, nz, h1 )
261  x1(1) = 0.0 ; do k = 1,nz ; x1(k+1) = x1(k) + h1(k) ; enddo
262 
263  ! Remap T and S from previous grid to new grid
264  do k = 1,nz
265  h1(k) = x1(k+1) - x1(k)
266  enddo
267 
268  call remapping_core_h(remapcs, nz, h0, s, nz, h1, tmp, h_neglect, h_neglect_edge)
269  s_tmp(:) = tmp(:)
270 
271  call remapping_core_h(remapcs, nz, h0, t, nz, h1, tmp, h_neglect, h_neglect_edge)
272  t_tmp(:) = tmp(:)
273 
274  ! Compute the deviation between two successive grids
275  deviation = 0.0
276  x0(1) = 0.0
277  x1(1) = 0.0
278  do k = 2,nz
279  x0(k) = x0(k-1) + h0(k-1)
280  x1(k) = x1(k-1) + h1(k-1)
281  deviation = deviation + (x0(k)-x1(k))**2
282  enddo
283  deviation = sqrt( deviation / (nz-1) )
284 
285  m = m + 1
286 
287  ! Copy final grid onto start grid for next iteration
288  h0(:) = h1(:)
289 
290  enddo ! end regridding iterations
291 
292  if (cs%integrate_downward_for_e) then
293  zinterface(1) = 0.
294  do k = 1,nz
295  zinterface(k+1) = zinterface(k) - h1(k)
296  ! Adjust interface position to accommodate inflating layers
297  ! without disturbing the interface above
298  enddo
299  else
300  ! The rest of the model defines grids integrating up from the bottom
301  zinterface(nz+1) = -depth
302  do k = nz,1,-1
303  zinterface(k) = zinterface(k+1) + h1(k)
304  ! Adjust interface position to accommodate inflating layers
305  ! without disturbing the interface above
306  enddo
307  endif
308 
309 end subroutine build_rho_column_iteratively
310 
311 !> Copy column thicknesses with vanished layers removed
312 subroutine copy_finite_thicknesses(nk, h_in, threshold, nout, h_out, mapping)
313  integer, intent(in) :: nk !< Number of layer for h_in, T_in, S_in
314  real, dimension(nk), intent(in) :: h_in !< Thickness of input column
315  real, intent(in) :: threshold !< Thickness threshold defining vanished layers
316  integer, intent(out) :: nout !< Number of non-vanished layers
317  real, dimension(nk), intent(out) :: h_out !< Thickness of output column
318  integer, dimension(nk), intent(out) :: mapping !< Index of k-out corresponding to k-in
319  ! Local variables
320  integer :: k, k_thickest
321  real :: thickness_in_vanished, thickest_h_out
322 
323  ! Build up new grid
324  nout = 0
325  thickness_in_vanished = 0.0
326  thickest_h_out = h_in(1)
327  k_thickest = 1
328  do k = 1, nk
329  mapping(k) = nout ! Note k>=nout always
330  h_out(k) = 0. ! Make sure h_out is set everywhere
331  if (h_in(k) > threshold) then
332  ! For non-vanished layers
333  nout = nout + 1
334  mapping(nout) = k
335  h_out(nout) = h_in(k)
336  if (h_out(nout) > thickest_h_out) then
337  thickest_h_out = h_out(nout)
338  k_thickest = nout
339  endif
340  else
341  ! Add up mass in vanished layers
342  thickness_in_vanished = thickness_in_vanished + h_in(k)
343  endif
344  enddo
345 
346  ! No finite layers
347  if (nout <= 1) return
348 
349  ! Adjust for any lost volume in vanished layers
350  h_out(k_thickest) = h_out(k_thickest) + thickness_in_vanished
351 
352 end subroutine copy_finite_thicknesses
353 
354 !------------------------------------------------------------------------------
355 !> Inflate vanished layers to finite (nonzero) width
356 subroutine old_inflate_layers_1d( min_thickness, nk, h )
357 
358  ! Argument
359  real, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
360  integer, intent(in) :: nk !< Number of layers in the grid
361  real, dimension(:), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
362 
363  ! Local variable
364  integer :: k
365  integer :: k_found
366  integer :: count_nonzero_layers
367  real :: delta
368  real :: correction
369  real :: maxthickness
370 
371  ! Count number of nonzero layers
372  count_nonzero_layers = 0
373  do k = 1,nk
374  if ( h(k) > min_thickness ) then
375  count_nonzero_layers = count_nonzero_layers + 1
376  endif
377  enddo
378 
379  ! If all layer thicknesses are greater than the threshold, exit routine
380  if ( count_nonzero_layers == nk ) return
381 
382  ! If all thicknesses are zero, inflate them all and exit
383  if ( count_nonzero_layers == 0 ) then
384  do k = 1,nk
385  h(k) = min_thickness
386  enddo
387  return
388  endif
389 
390  ! Inflate zero layers
391  correction = 0.0
392  do k = 1,nk
393  if ( h(k) <= min_thickness ) then
394  delta = min_thickness - h(k)
395  correction = correction + delta
396  h(k) = h(k) + delta
397  endif
398  enddo
399 
400  ! Modify thicknesses of nonzero layers to ensure volume conservation
401  maxthickness = h(1)
402  k_found = 1
403  do k = 1,nk
404  if ( h(k) > maxthickness ) then
405  maxthickness = h(k)
406  k_found = k
407  endif
408  enddo
409 
410  h(k_found) = h(k_found) - correction
411 
412 end subroutine old_inflate_layers_1d
413 
414 end module coord_rho
coord_rho::deviation_tolerance
real, parameter deviation_tolerance
Deviation tolerance between succesive grids in regridding iterations.
Definition: coord_rho.F90:42
coord_rho::nb_regridding_iterations
integer, parameter nb_regridding_iterations
Maximum number of regridding iterations.
Definition: coord_rho.F90:40
coord_rho::set_rho_params
subroutine, public set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS)
This subroutine can be used to set the parameters for the coord_rho module.
Definition: coord_rho.F90:81
coord_rho
Regrid columns for the continuous isopycnal (rho) coordinate.
Definition: coord_rho.F90:2
mom_remapping::remapping_core_h
subroutine, public remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned.
Definition: MOM_remapping.F90:188
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
coord_rho::rho_cs
Control structure containing required parameters for the rho coordinate.
Definition: coord_rho.F90:14
mom_remapping::remapping_cs
Container for remapping parameters.
Definition: MOM_remapping.F90:22
coord_rho::init_coord_rho
subroutine, public init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS, rho_scale)
Initialise a rho_CS with pointers to parameters.
Definition: coord_rho.F90:50
mom_remapping
Provides column-wise vertical remapping functions.
Definition: MOM_remapping.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::degree_max
integer, parameter, public degree_max
Interpolant degrees.
Definition: regrid_interp.F90:56
regrid_interp::interp_cs_type
Control structure for regrid_interp module.
Definition: regrid_interp.F90:23
coord_rho::build_rho_column
subroutine, public build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, h_neglect, h_neglect_edge)
Build a rho coordinate column.
Definition: coord_rho.F90:102
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
coord_rho::end_coord_rho
subroutine, public end_coord_rho(CS)
This subroutine deallocates memory in the control structure for the coord_rho module.
Definition: coord_rho.F90:71
regrid_interp::build_and_interpolate_grid
subroutine, public build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, n1, h1, x1, h_neglect, h_neglect_edge)
Build a grid by interpolating for target values.
Definition: regrid_interp.F90:307
coord_rho::copy_finite_thicknesses
subroutine copy_finite_thicknesses(nk, h_in, threshold, nout, h_out, mapping)
Copy column thicknesses with vanished layers removed.
Definition: coord_rho.F90:313
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
coord_rho::old_inflate_layers_1d
subroutine, public old_inflate_layers_1d(min_thickness, nk, h)
Inflate vanished layers to finite (nonzero) width.
Definition: coord_rho.F90:357
coord_rho::build_rho_column_iteratively
subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_state, zInterface, h_neglect, h_neglect_edge)
Iteratively build a rho coordinate column.
Definition: coord_rho.F90:196
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60