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