Recon1d_type.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!> A generic type for vertical 1D reconstructions
6module recon1d_type
7
9
10implicit none ; private
11
12public recon1d
13public testing
14
15!> The base class for implementations of 1D reconstructions
16type, abstract :: recon1d
17
18 integer :: n = 0 !< Number of cells in column
19 real, allocatable, dimension(:) :: u_mean !< Cell mean [A]
20 real :: h_neglect = 0. !< A negligibly small width used in cell reconstructions in the same units as h [H]
21 real :: x_tolerance = 1. * epsilon(1.) !< Solver tolerance for x in element (0,1) [nondim]
22 logical :: check = .false. !< If true, enable some consistency checking
23
24 logical :: debug = .false. !< If true, dump info as calculations are made (do not enable)
25contains
26
27 ! The following functions/subroutines are deferred and must be provided specifically by each scheme
28
29 !> Deferred implementation of initialization
30 procedure(i_init), deferred :: init
31 !> Deferred implementation of reconstruction function
32 procedure(i_reconstruct), deferred :: reconstruct
33 !> Deferred implementation of the average over an interval
34 procedure(i_average), deferred :: average
35 !> Deferred implementation of evaluating the reconstruction at a point
36 procedure(i_f), deferred :: f
37 !> Deferred implementation of the derivative of the reconstruction at a point
38 procedure(i_dfdx), deferred :: dfdx
39 !> Deferred implementation of check_reconstruction
40 !!
41 !! Returns True if a check fails. Returns False if all checks pass.
42 !! Checks are about internal, or inferred, state for arbitrary inputs.
43 !! Checks should cover all the expected properties of a reconstruction.
44 procedure(i_check_reconstruction), deferred :: check_reconstruction
45 !> Deferred implementation of unit tests for the reconstruction
46 !!
47 !! Returns True if a test fails. Returns False if all tests pass.
48 !! Tests in unit_tests() are usually checks against known (e.g. analytic) solutions.
49 procedure(i_unit_tests), deferred :: unit_tests
50 !> Deferred implementation of deallocation
51 procedure(i_destroy), deferred :: destroy
52
53 ! The following functions/subroutines are shared across all reconstructions and provided by this module
54 ! unless replaced for the purpose of optimization
55
56 !> Solves for x such that f(x)=t
57 procedure :: x => x
58 !> Remaps the column to subgrid h_sub
60 !> Set debugging
61 procedure :: set_debug => a_set_debug
62
63 ! The following functions usually point to the same implementation as above but
64 ! for derived secondary children these allow invocation of the parent class function.
65
66 !> Second interface to init(), used to reach the primary class if derived from a primary implementation
67 procedure(i_init_parent), deferred :: init_parent
68 !> Second interface to reconstruct(), used to reach the primary class if derived from a primary implementation
69 procedure(i_reconstruct_parent), deferred :: reconstruct_parent
70
71end type recon1d
72
73interface
74
75 !> Initialize a 1D reconstruction for n cells
76 subroutine i_init(this, n, h_neglect, check)
77 import :: recon1d
78 class(recon1d), intent(out) :: this !< This reconstruction
79 integer, intent(in) :: n !< Number of cells in this column
80 real, optional, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H]
81 logical, optional, intent(in) :: check !< If true, enable some consistency checking
82 end subroutine i_init
83
84 !> Calculate a 1D reconstructions based on h(:) and u(:)
85 subroutine i_reconstruct(this, h, u)
86 import :: recon1d
87 class(recon1d), intent(inout) :: this !< This reconstruction
88 real, intent(in) :: h(*) !< Grid spacing (thickness), typically in [H]
89 real, intent(in) :: u(*) !< Cell mean values [A]
90 end subroutine i_reconstruct
91
92 !> Average between xa and xb for cell k of a 1D reconstruction [A]
93 !!
94 !! It is assumed that 0<=xa<=1, 0<=xb<=1, and xa<=xb
95 real function i_average(this, k, xa, xb)
96 import :: recon1d
97 class(recon1d), intent(in) :: this !< This reconstruction
98 integer, intent(in) :: k !< Cell number
99 real, intent(in) :: xa !< Start of averaging interval on element (0 to 1)
100 real, intent(in) :: xb !< End of averaging interval on element (0 to 1)
101 end function i_average
102
103 !> Point-wise value of reconstruction [A]
104 !!
105 !! The function is only valid for 0 <= x <= 1. x is effectively clipped to this range.
106 real function i_f(this, k, x)
107 import :: recon1d
108 class(recon1d), intent(in) :: this !< This reconstruction
109 integer, intent(in) :: k !< Cell number
110 real, intent(in) :: x !< Non-dimensional position within element [nondim]
111 end function i_f
112
113 !> Point-wise value of derivative reconstruction [A]
114 !!
115 !! The function is only valid for 0 <= x <= 1. x is effectively clipped to this range.
116 real function i_dfdx(this, k, x)
117 import :: recon1d
118 class(recon1d), intent(in) :: this !< This reconstruction
119 integer, intent(in) :: k !< Cell number
120 real, intent(in) :: x !< Non-dimensional position within element [nondim]
121 end function i_dfdx
122
123 !> Point-wise solver for x: f(x)=t [nondim]
124 !!
125 !! The function solves for the non-dimensional position x within the cell where
126 !! the reconstruction f(x)=t. The solver returns x=0 or x=1 if the target, t,
127 !! is outside of the cell.
128 real function i_x(this, k, t)
129 import :: recon1d
130 class(recon1d), intent(in) :: this !< This reconstruction
131 integer, intent(in) :: k !< Cell number
132 real, intent(in) :: t !< Value to solve for [A]
133 end function i_x
134
135 !> Returns true if some inconsistency is detected, false otherwise
136 !!
137 !! The nature of "consistency" is defined by the implementations
138 !! and might be no-ops.
139 logical function i_check_reconstruction(this, h, u)
140 import :: recon1d
141 class(recon1d), intent(in) :: this !< This reconstruction
142 real, intent(in) :: h(*) !< Grid spacing (thickness), typically in [H]
143 real, intent(in) :: u(*) !< Cell mean values [A]
144 end function i_check_reconstruction
145
146 !> Deallocate a 1D reconstruction
147 subroutine i_destroy(this)
148 import :: recon1d
149 class(recon1d), intent(inout) :: this !< This reconstruction
150 end subroutine i_destroy
151
152 !> Second interface to init(), or to parent init()
153 subroutine i_init_parent(this, n, h_neglect, check)
154 import :: recon1d
155 class(recon1d), intent(out) :: this !< This reconstruction
156 integer, intent(in) :: n !< Number of cells in this column
157 real, optional, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H]
158 logical, optional, intent(in) :: check !< If true, enable some consistency checking
159 end subroutine i_init_parent
160
161 !> Second interface to reconstruct(), or to parent reconstruct()
162 subroutine i_reconstruct_parent(this, h, u)
163 import :: recon1d
164 class(recon1d), intent(inout) :: this !< This reconstruction
165 real, intent(in) :: h(*) !< Grid spacing (thickness), typically in [H]
166 real, intent(in) :: u(*) !< Cell mean values [A]
167 end subroutine i_reconstruct_parent
168
169 !> Runs reconstruction unit tests and returns True for any fails, False otherwise
170 !!
171 !! Assumes single process/thread context
172 logical function i_unit_tests(this, verbose, stdout, stderr)
173 import :: recon1d
174 class(recon1d), intent(inout) :: this !< This reconstruction
175 logical, intent(in) :: verbose !< True, if verbose
176 integer, intent(in) :: stdout !< I/O channel for stdout
177 integer, intent(in) :: stderr !< I/O channel for stderr
178 end function i_unit_tests
179
180end interface
181
182contains
183
184!> Solve for x such that f(x)=t
185!!
186!! This solver uses bounded Newton-Raphson method with a fixed
187!! number of iterations
188real function x(this, k, t)
189 class(recon1d), intent(in) :: this !< This reconstruction
190 integer, intent(in) :: k !< Cell number
191 real, intent(in) :: t !< Value to solve for [A]
192 real :: xl, xr, xo ! Left/right bounds and guess [nondim]
193 real :: fl, fr ! Left right values [A]
194 real :: slp ! Difference across cell or derivative wrt nondim x [A]
195 real :: f_at_x ! Value at current x [A]
196 integer :: iter
197
198 x = 0.5 ! Fall back for special conditions
199 fl = this%f(k, 0.)
200 fr = this%f(k, 1.)
201 slp = fr - fl
202 if ( ( fl - t ) * ( t - fr ) > 0. ) then
203 ! t is inside the range fl..fr
204 xl = 0.
205 xr = 1.
206 xo = ( t - this%f(k, 0.) ) / slp ! First guess by regula falsi
207 f_at_x = this%f(k, xo)
208 do iter = 1,10
209 slp = this%dfdx(k, xo)
210 x = xo - ( f_at_x - t ) / slp ! Newton-Raphson step
211 if ( x < xl ) x = 0.5 * ( xl + xo ) ! Replace with bi-section
212 if ( x > xr ) x = 0.5 * ( xr + xo ) ! Replace with bi-section
213 f_at_x = this%f(k, x)
214 if ( abs(f_at_x - t) <= 0. .or. abs(x - xo) < this%x_tolerance ) return
215 if ( f_at_x < t ) xl = x ! Replace left bound
216 if ( f_at_x > t ) xr = x ! Replace right bound
217 xo = x
218 enddo
219 elseif ( abs(slp) > 0. ) then
220 slp = sign(1., slp)
221 ! if t>u_mean & slp=1 then x=1
222 ! if t<u_mean & slp=1 then x=0
223 ! if t>u_mean & slp=-1 then x=0
224 ! if t<u_mean & slp=-1 then x=1
225 x = 0.5 + slp * sign(0.5, t - this%u_mean(k))
226 else
227 ! slp=0 so estimate "direction" from neighbors
228 slp = this%f(min(k+1,this%n), 0.) - this%f(max(k-1,1), 1.)
229 if ( abs(slp) > 0. ) slp = sign(1., slp)
230 ! if t>u_mean & slp=1 then x=1
231 ! if t<u_mean & slp=1 then x=0
232 ! if t>u_mean & slp=-1 then x=0
233 ! if t<u_mean & slp=-1 then x=1
234 ! if t=u_mean then x=0.5
235 ! if slp=0 then x=0.5
236 if ( abs(t - this%u_mean(k)) > 0. ) x = 0.5 + slp * sign(0.5, t - this%u_mean(k))
237 endif
238end function x
239
240!> Remaps the column to subgrid h_sub
241!!
242!! It is assumed that h_sub is a perfect sub-grid of h0, meaning each h0 cell
243!! can be constructed by joining a contiguous set of h_sub cells. The integer
244!! indices isrc_start, isrc_end, isub_src provide this mapping, and are
245!! calculated in MOM_remapping
246subroutine remap_to_sub_grid(this, h0, u0, n1, h_sub, &
247 isrc_start, isrc_end, isrc_max, isub_src, &
248 u_sub, uh_sub, u02_err)
249 class(recon1d), intent(in) :: this !< 1-D reconstruction type
250 real, intent(in) :: h0(*) !< Source grid widths (size n0) [H]
251 real, intent(in) :: u0(*) !< Source grid widths (size n0) [H]
252 integer, intent(in) :: n1 !< Number of cells in target grid
253 real, intent(in) :: h_sub(*) !< Overlapping sub-cell thicknesses, h_sub [H]
254 integer, intent(in) :: isrc_start(*) !< Index of first sub-cell within each source cell
255 integer, intent(in) :: isrc_end(*) !< Index of last sub-cell within each source cell
256 integer, intent(in) :: isrc_max(*) !< Index of thickest sub-cell within each source cell
257 integer, intent(in) :: isub_src(*) !< Index of source cell for each sub-cell
258 real, intent(out) :: u_sub(*) !< Sub-cell cell averages (size n1) [A]
259 real, intent(out) :: uh_sub(*) !< Sub-cell cell integrals (size n1) [A H]
260 real, intent(out) :: u02_err !< Integrated reconstruction error estimates [A H]
261 ! Local variables
262 integer :: i_sub ! Index of sub-cell
263 integer :: i0 ! Index into h0(1:n0), source column
264 integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell
265 real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H]
266 real :: xa, xb ! Non-dimensional position within a source cell (0..1) [nondim]
267 real :: dh ! The width of the sub-cell [H]
268 real :: duh ! The total amount of accumulated stuff (u*h) [A H]
269 real :: dh0_eff ! Running sum of source cell thickness [H]
270 integer :: i0_last_thick_cell, n0
271! real :: u0_min(this%n), u0_max(this%n) ! Min/max of u0 for each source cell [A]
272! real :: ul,ur ! Left/right edge values [A]
273
274 n0 = this%n
275
276 i0_last_thick_cell = 0
277 do i0 = 1, n0
278! ul = this%f(i0, 0.)
279! ur = this%f(i0, 1.)
280! u0_min(i0) = min(ul, ur)
281! u0_max(i0) = max(ul, ur)
282 if (h0(i0)>0.) i0_last_thick_cell = i0
283 enddo
284
285 ! Loop over each sub-cell to calculate average/integral values within each sub-cell.
286 ! Uses: h_sub, isub_src, h0_eff
287 ! Sets: u_sub, uh_sub
288 xa = 0.
289 dh0_eff = 0.
290 u02_err = 0.
291 do i_sub = 1, n0+n1
292
293 ! Sub-cell thickness from loop above
294 dh = h_sub(i_sub)
295
296 ! Source cell
297 i0 = isub_src(i_sub)
298
299 ! Evaluate average and integral for sub-cell i_sub.
300 ! Integral is over distance dh but expressed in terms of non-dimensional
301 ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).
302 dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
303 if (h0(i0)>0.) then
304 xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0
305 xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
306 u_sub(i_sub) = this%average( i0, xa, xb )
307 else ! Vanished cell
308 xb = 1.
309 u_sub(i_sub) = u0(i0)
310 endif
311! u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )
312! u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )
313 uh_sub(i_sub) = dh * u_sub(i_sub)
314
315 if (isub_src(i_sub+1) /= i0) then
316 ! If the next sub-cell is in a different source cell, reset the position counters
317 dh0_eff = 0.
318 xa = 0.
319 else
320 xa = xb ! Next integral will start at end of last
321 endif
322
323 enddo
324 i_sub = n0+n1+1
325 ! Sub-cell thickness from loop above
326 dh = h_sub(i_sub)
327 ! Source cell
328 i0 = isub_src(i_sub)
329
330 ! Evaluate average and integral for sub-cell i_sub.
331 ! Integral is over distance dh but expressed in terms of non-dimensional
332 ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).
333 dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
334 if (h0(i0)>0.) then
335 xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0
336 xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
337 u_sub(i_sub) = this%average( i0, xa, xb )
338 else ! Vanished cell
339 xb = 1.
340 u_sub(i_sub) = u0(i0)
341 endif
342! u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )
343! u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )
344 uh_sub(i_sub) = dh * u_sub(i_sub)
345
346 ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within
347 ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals
348 ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (\@Hallberg-NOAA).
349 ! Uses: i0_last_thick_cell, isrc_max, h_sub, isrc_start, isrc_end, uh_sub, u0, h0
350 ! Updates: uh_sub
351 do i0 = 1, i0_last_thick_cell
352 i_max = isrc_max(i0)
353 dh_max = h_sub(i_max)
354 if (dh_max > 0.) then
355 ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell.
356 duh = 0.
357 do i_sub = isrc_start(i0), isrc_end(i0)
358 if (i_sub /= i_max) duh = duh + uh_sub(i_sub)
359 enddo
360 uh_sub(i_max) = u0(i0)*h0(i0) - duh
361 u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) )
362 endif
363 enddo
364
365 ! This should not generally be used
366 if (this%check) then
367 if ( this%check_reconstruction(h0, u0) ) stop 910 ! A debugger is required to understand why this failed
368 endif
369
370end subroutine remap_to_sub_grid
371
372!> Turns on debugging
373subroutine a_set_debug(this)
374 class(recon1d), intent(inout) :: this !< 1-D reconstruction type
375
376 this%debug = .true.
377
378end subroutine a_set_debug
379
380!> \namespace recon1d_type
381!!
382!! \section section_recon1d_type Generic vertical reconstruction type
383!!
384!! A class to describe generic reconstruction in 1-D. This module has no implementations
385!! but defines the interfaces for members that implement a reconstruction.
386!!
387!! e.g. a chain of derived reconstructions might look like
388!! Recon1d_type <- Recond1d_XYZ <- Recon1d_XYZ_v2
389!! where
390!! Recon1d_type - defines the interfaces (this module)
391!! Recon1d_XYZ - extends Recon1d_type, implements the XYZ reconstruction in reconstruct(),
392!! and reconstruc_parent() -> reconstruct() of the same Recon1d_XYZ module
393!! Recon1d_XYZ_v2 - implements a slight variant of Recon1d_XYZ via reconstruct()
394!! but reconstruc_parent() is not redefined so that it still is defined by Recon1d_XYZ
395!!
396!! The schemes that use this structure are described in \ref Vertical_Reconstruction
397end module recon1d_type