Recon1d_EMPLM_WA_poly.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!> Extrapolated-Monotonized Piecewise Linear Method 1D reconstruction
6!!
7!! This extends MPLM_poly, following White and Adcroft, 2008 \cite white2008, by extraplating for the slopes of the
8!! first and last cells. This extrapolation is used by White et al., 2009, during grid-generation.
9!!
10!! This stores and evaluates the reconstruction using a polynomial representation which is not preferred
11!! but was the form used in OM4.
13
14use recon1d_mplm_wa_poly, only : mplm_wa_poly, testing
15
16implicit none ; private
17
18public emplm_wa_poly
19
20!> Extrapolation Limited Monotonic PLM reconstruction following White and Adcroft, 2008
21!!
22!! The source for the methods ultimately used by this class are:
23!! - init() -> recon1d_mplm_wa_poly.init()
24!! - reconstruct() -> recon1d_mplm_wa_poly.reconstruct()
25!! - average() -> recon1d_mplm_wa_poly.average()
26!! - f() -> recon1d_mplm_wa_poly -> recon1d_mplm_wa -> recon1d_plm_cw.f()
27!! - dfdx() -> recon1d_mplm_wa_poly -> recon1d_mplm_wa -> recon1d_plm_cw.dfdx()
28!! - check_reconstruction() *locally defined
29!! - unit_tests() *locally defined
30!! - destroy() -> recon1d_mplm_wa_poly -> recon1d_mplm_wa -> recon1d_plm_cw.destroy()
31!! - remap_to_sub_grid() *locally defined
32!! - init_parent() -> recon1d_mplm_wa_poly -> recon1d_mplm_wa.init()
33!! - reconstruct_parent() -> recon1d_mplm_wa_poly -> recon1d_mplm_wa.reconstruct()
34type, extends (mplm_wa_poly) :: emplm_wa_poly
35
36contains
37 !> Implementation of the EMPLM_WA_poly reconstruction with boundary extrapolation
38 procedure :: reconstruct => reconstruct
39 !> Implementation of check reconstruction for the EMPLM_WA_poly reconstruction
41 !> Implementation of unit tests for the EMPLM_WA_poly reconstruction
42 procedure :: unit_tests => unit_tests
43
44end type emplm_wa_poly
45
46contains
47
48!> Calculate a 1D PLM reconstruction based on h(:) and u(:)
49subroutine reconstruct(this, h, u)
50 class(emplm_wa_poly), intent(inout) :: this !< This reconstruction
51 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
52 real, intent(in) :: u(*) !< Cell mean values [A]
53 ! Local variables
54 integer :: n
55 real :: slope ! Difference of u across cell [A]
56
57 ! Use parent (MPLM_WA) reconstruction
58 call this%reconstruct_parent(h, u)
59
60 n = this%n
61
62 ! Fix reconstruction for first cell
63 slope = - plm_extrapolate_slope( h(2), h(1), this%h_neglect, u(2), u(1) )
64 this%ul(1) = u(1) - 0.5 * slope
65 this%ur(1) = u(1) + 0.5 * slope
66 this%poly_coef(1,1) = this%ul(1)
67 this%poly_coef(1,2) = this%ur(1) - this%ul(1)
68
69 ! Fix reconstruction for last cell
70 slope = plm_extrapolate_slope( h(n-1), h(n), this%h_neglect, u(n-1), u(n) )
71 this%ul(n) = u(n) - 0.5 * slope
72 this%ur(n) = u(n) + 0.5 * slope
73 this%poly_coef(n,1) = this%ul(n)
74 this%poly_coef(n,2) = this%ur(n) - this%ul(n)
75
76end subroutine reconstruct
77
78!> Returns a PLM slope using h2 extrapolation from a cell to the left, in the same
79!! arbitrary units as the input values [A].
80!! Use the negative to extrapolate from the cell to the right.
81real elemental pure function plm_extrapolate_slope(h_l, h_c, h_neglect, u_l, u_c)
82 real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H]
83 real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H]
84 real, intent(in) :: h_neglect !< A negligible thickness [H]
85 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
86 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
87 ! Local variables
88 real :: left_edge ! Left edge value [A]
89 real :: hl, hc ! Left and central cell thicknesses [H]
90
91 ! Avoid division by zero for vanished cells
92 hl = h_l + h_neglect
93 hc = h_c + h_neglect
94
95 ! The h2 scheme is used to compute the left edge value
96 left_edge = (u_l*hc + u_c*hl) / (hl + hc)
97
98 plm_extrapolate_slope = 2.0 * ( u_c - left_edge )
99
100end function plm_extrapolate_slope
101
102!> Checks the EMPLM_WA_poly reconstruction for consistency
103logical function check_reconstruction(this, h, u)
104 class(emplm_wa_poly), intent(in) :: this !< This reconstruction
105 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
106 real, intent(in) :: u(*) !< Cell mean values [A]
107 ! Local variables
108 integer :: k
109
110 check_reconstruction = .false.
111
112 do k = 1, this%n
113 if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
114 enddo
115
116 ! Check implied curvature
117 do k = 1, this%n
118 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ur(k) - this%u_mean(k) ) < 0. ) check_reconstruction = .true.
119 enddo
120
121 ! These two checks fail MOM_remapping:test_recon_consistency in the presence of vanished layers
122 ! e.g. intel/2023.2.0 on gaea at iter=26
123
124! ! Check bounding of right edges, w.r.t. the cell means
125! do K = 1, this%n-1
126! if ( ( this%ur(k) - this%u_mean(k) ) * ( this%u_mean(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
127! enddo
128
129! ! Check bounding of left edges, w.r.t. the cell means
130! do K = 2, this%n
131! if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%u_mean(k-1) ) < 0. ) check_reconstruction = .true.
132! enddo
133
134 ! Check order of u, ur, ul
135 ! Note that in the OM4-era implementation, we were not consistent for top and bottom layers due
136 ! extrapolation using cell means rather than edge values, hence reduced range for K
137 do k = 2, this%n-2
138 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%ul(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
139 enddo
140
141 ! Check bounding of left edges, w.r.t. this cell mean and the previous cell right edge
142 do k = 3, this%n-1
143 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%ur(k-1) ) < 0. ) check_reconstruction = .true.
144 enddo
145
146end function check_reconstruction
147
148
149!> Runs PLM reconstruction unit tests and returns True for any fails, False otherwise
150logical function unit_tests(this, verbose, stdout, stderr)
151 class(emplm_wa_poly), intent(inout) :: this !< This reconstruction
152 logical, intent(in) :: verbose !< True, if verbose
153 integer, intent(in) :: stdout !< I/O channel for stdout
154 integer, intent(in) :: stderr !< I/O channel for stderr
155 ! Local variables
156 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
157 real, allocatable :: ull(:), urr(:) ! test values [A]
158 type(testing) :: test ! convenience functions
159 integer :: k
160
161 call test%set( stdout=stdout ) ! Sets the stdout channel in test
162 call test%set( stderr=stderr ) ! Sets the stderr channel in test
163 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
164
165 call this%init(3)
166 call test%test( this%n /= 3, 'Setting number of levels')
167 allocate( um(3), ul(3), ur(3), ull(3), urr(3) )
168
169 call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) )
170 call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values')
171
172 do k = 1, 3
173 ul(k) = this%f(k, 0.)
174 um(k) = this%f(k, 0.5)
175 ur(k) = this%f(k, 1.)
176 enddo
177 call test%real_arr(3, ul, (/0.,2.,4./), 'Evaluation on left edge')
178 call test%real_arr(3, um, (/1.,3.,5./), 'Evaluation in center')
179 call test%real_arr(3, ur, (/2.,4.,6./), 'Evaluation on right edge')
180
181 do k = 1, 3
182 ul(k) = this%dfdx(k, 0.)
183 um(k) = this%dfdx(k, 0.5)
184 ur(k) = this%dfdx(k, 1.)
185 enddo
186 call test%real_arr(3, ul, (/2.,2.,2./), 'dfdx on left edge')
187 call test%real_arr(3, um, (/2.,2.,2./), 'dfdx in center')
188 call test%real_arr(3, ur, (/2.,2.,2./), 'dfdx on right edge')
189
190 do k = 1, 3
191 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
192 enddo
193 call test%real_arr(3, um, (/1.25,3.25,5.25/), 'Return interval average')
194
195 unit_tests = test%summarize('EMPLM_WA_poly:unit_tests')
196
197end function unit_tests
198
199!> \namespace recon1d_emplm_wa_poly
200!!
201
202end module recon1d_emplm_wa_poly