MOM6
coord_rho Module Reference

Detailed Description

Regrid columns for the continuous isopycnal (rho) coordinate.

Data Types

type  rho_cs
 Control structure containing required parameters for the rho coordinate. More...
 

Functions/Subroutines

subroutine, public init_coord_rho (CS, nk, ref_pressure, target_density, interp_CS, rho_scale)
 Initialise a rho_CS with pointers to parameters. More...
 
subroutine, public end_coord_rho (CS)
 This subroutine deallocates memory in the control structure for the coord_rho module. More...
 
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. More...
 
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. More...
 
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. More...
 
subroutine copy_finite_thicknesses (nk, h_in, threshold, nout, h_out, mapping)
 Copy column thicknesses with vanished layers removed. More...
 
subroutine, public old_inflate_layers_1d (min_thickness, nk, h)
 Inflate vanished layers to finite (nonzero) width. More...
 

Variables

integer, parameter nb_regridding_iterations = 1
 Maximum number of regridding iterations. More...
 
real, parameter deviation_tolerance = 1e-10
 Deviation tolerance between succesive grids in regridding iterations. More...
 

Function/Subroutine Documentation

◆ build_rho_column()

subroutine, public coord_rho::build_rho_column ( type(rho_cs), intent(in)  CS,
integer, intent(in)  nz,
real, intent(in)  depth,
real, dimension(nz), intent(in)  h,
real, dimension(nz), intent(in)  T,
real, dimension(nz), intent(in)  S,
type(eos_type), pointer  eqn_of_state,
real, dimension(cs%nk+1), intent(inout)  z_interface,
real, intent(in), optional  h_neglect,
real, intent(in), optional  h_neglect_edge 
)

Build a rho coordinate column.

  1. Density profiles are calculated on the source grid.
  2. Positions of target densities (for interfaces) are found by interpolation.
    Parameters
    [in]cscoord_rho control structure
    [in]nzNumber of levels on source grid (i.e. length of h, T, S)
    [in]depthDepth of ocean bottom (positive in m)
    [in]hLayer thicknesses [H ~> m or kg m-2]
    [in]tT for source column
    [in]sS for source column
    eqn_of_stateEquation of state structure
    [in,out]z_interfaceAbsolute positions of interfaces
    [in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h0.
    [in]h_neglect_edgeA negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 102 of file coord_rho.F90.

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 

References regrid_interp::build_and_interpolate_grid(), copy_finite_thicknesses(), and old_inflate_layers_1d().

Referenced by mom_diag_remap::diag_remap_update().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ build_rho_column_iteratively()

subroutine coord_rho::build_rho_column_iteratively ( type(rho_cs), intent(in)  CS,
type(remapping_cs), intent(in)  remapCS,
integer, intent(in)  nz,
real, intent(in)  depth,
real, dimension(nz), intent(in)  h,
real, dimension(nz), intent(in)  T,
real, dimension(nz), intent(in)  S,
type(eos_type), pointer  eqn_of_state,
real, dimension(nz+1), intent(inout)  zInterface,
real, intent(in), optional  h_neglect,
real, intent(in), optional  h_neglect_edge 
)
private

Iteratively build a rho coordinate column.

The algorithm operates as follows within each column:

  1. Given T & S within each layer, the layer densities are computed.
  2. Based on these layer densities, a global density profile is reconstructed (this profile is monotonically increasing and may be discontinuous)
  3. The new grid interfaces are determined based on the target interface densities.
  4. T & S are remapped onto the new grid.
  5. Return to step 1 until convergence or until the maximum number of iterations is reached, whichever comes first.
    Parameters
    [in]csRegridding control structure
    [in]remapcsRemapping parameters and options
    [in]nzNumber of levels
    [in]depthDepth of ocean bottom [Z ~> m]
    [in]hLayer thicknesses in Z coordinates [Z ~> m]
    [in]tT for column [degC]
    [in]sS for column [ppt]
    eqn_of_stateEquation of state structure
    [in,out]zinterfaceAbsolute positions of interfaces
    [in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h [Z ~> m]
    [in]h_neglect_edgeA negligibly small width for the purpose of edge value calculations in the same units as h [Z ~> m]

Definition at line 196 of file coord_rho.F90.

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 

References regrid_interp::build_and_interpolate_grid(), copy_finite_thicknesses(), deviation_tolerance, nb_regridding_iterations, old_inflate_layers_1d(), and mom_remapping::remapping_core_h().

Here is the call graph for this function:

◆ copy_finite_thicknesses()

subroutine coord_rho::copy_finite_thicknesses ( integer, intent(in)  nk,
real, dimension(nk), intent(in)  h_in,
real, intent(in)  threshold,
integer, intent(out)  nout,
real, dimension(nk), intent(out)  h_out,
integer, dimension(nk), intent(out)  mapping 
)
private

Copy column thicknesses with vanished layers removed.

Parameters
[in]nkNumber of layer for h_in, T_in, S_in
[in]h_inThickness of input column
[in]thresholdThickness threshold defining vanished layers
[out]noutNumber of non-vanished layers
[out]h_outThickness of output column
[out]mappingIndex of k-out corresponding to k-in

Definition at line 313 of file coord_rho.F90.

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 

Referenced by build_rho_column(), and build_rho_column_iteratively().

Here is the caller graph for this function:

◆ end_coord_rho()

subroutine, public coord_rho::end_coord_rho ( type(rho_cs), pointer  CS)

This subroutine deallocates memory in the control structure for the coord_rho module.

Parameters
csCoordinate control structure

Definition at line 71 of file coord_rho.F90.

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)

◆ init_coord_rho()

subroutine, public coord_rho::init_coord_rho ( type(rho_cs), pointer  CS,
integer, intent(in)  nk,
real, intent(in)  ref_pressure,
real, dimension(:), intent(in)  target_density,
type(interp_cs_type), intent(in)  interp_CS,
real, intent(in), optional  rho_scale 
)

Initialise a rho_CS with pointers to parameters.

Parameters
csUnassociated pointer to hold the control structure
[in]nkNumber of layers in the grid
[in]ref_pressureCoordinate reference pressure [Pa]
[in]target_densityNominal density of interfaces [kg m-3 or R ~> kg m-3]
[in]interp_csControls for interpolation
[in]rho_scaleA dimensional scaling factor for target_density

Definition at line 50 of file coord_rho.F90.

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 

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ old_inflate_layers_1d()

subroutine, public coord_rho::old_inflate_layers_1d ( real, intent(in)  min_thickness,
integer, intent(in)  nk,
real, dimension(:), intent(inout)  h 
)

Inflate vanished layers to finite (nonzero) width.

Parameters
[in]min_thicknessMinimum allowed thickness [H ~> m or kg m-2]
[in]nkNumber of layers in the grid
[in,out]hLayer thicknesses [H ~> m or kg m-2]

Definition at line 357 of file coord_rho.F90.

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 

Referenced by build_rho_column(), build_rho_column_iteratively(), and mom_regridding::inflate_vanished_layers_old().

Here is the caller graph for this function:

◆ set_rho_params()

subroutine, public coord_rho::set_rho_params ( type(rho_cs), pointer  CS,
real, intent(in), optional  min_thickness,
logical, intent(in), optional  integrate_downward_for_e,
type(interp_cs_type), intent(in), optional  interp_CS 
)

This subroutine can be used to set the parameters for the coord_rho module.

Parameters
csCoordinate control structure
[in]min_thicknessMinimum allowed thickness [H ~> m or kg m-2]
[in]integrate_downward_for_eIf true, integrate for interface positions from the top downward. If false, integrate from the bottom upward, as does the rest of the model.
[in]interp_csControls for interpolation

Definition at line 81 of file coord_rho.F90.

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

References mom_error_handler::mom_error().

Here is the call graph for this function:

Variable Documentation

◆ deviation_tolerance

real, parameter coord_rho::deviation_tolerance = 1e-10
private

Deviation tolerance between succesive grids in regridding iterations.

Definition at line 42 of file coord_rho.F90.

42 real, parameter :: DEVIATION_TOLERANCE = 1e-10

Referenced by build_rho_column_iteratively().

◆ nb_regridding_iterations

integer, parameter coord_rho::nb_regridding_iterations = 1

Maximum number of regridding iterations.

Definition at line 40 of file coord_rho.F90.

40 integer, parameter :: NB_REGRIDDING_ITERATIONS = 1

Referenced by build_rho_column_iteratively().