MOM_hybgen_remap.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!> This module contains the hybgen remapping routines from HYCOM, with minor
6!! modifications to follow the MOM6 coding conventions
8
9implicit none ; private
10
12
13contains
14
15!> Set up the coefficients for PLM remapping of a set of scalars
16subroutine hybgen_plm_coefs(si, dpi, slope, nk, ns, thin, PCM_lay)
17 integer, intent(in) :: nk !< The number of input layers
18 integer, intent(in) :: ns !< The number of scalar fields to work on
19 real, intent(in) :: si(nk,ns) !< The cell-averaged input scalar fields [A]
20 real, intent(in) :: dpi(nk) !< The input grid layer thicknesses [H ~> m or kg m-2]
21 real, intent(out) :: slope(nk,ns) !< The PLM slope times cell width [A]
22 real, intent(in) :: thin !< A negligible layer thickness that can be ignored [H ~> m or kg m-2]
23 logical, optional, intent(in) :: pcm_lay(nk) !< If true for a layer, use PCM remapping for that layer
24
25!-----------------------------------------------------------------------
26! 1) coefficients for remapping from one set of vertical cells to another.
27! method: piecewise linear across each input cell with
28! monotonized central-difference limiter.
29!
30! van Leer, B., 1977, J. Comp. Phys., 23 276-299.
31!
32! 2) input arguments:
33! si - initial scalar fields in pi-layer space
34! dpi - initial layer thicknesses (dpi(k) = pi(k+1)-pi(k))
35! nk - number of layers
36! ns - number of fields
37! thin - layer thickness (>0) that can be ignored
38! PCM_lay - use PCM for selected layers (optional)
39!
40! 3) output arguments:
41! slope - coefficients for hybgen_plm_remap
42! profile(y) = si+slope*(y-1), -0.5 <= y <= 0.5
43!
44! 4) Tim Campbell, Mississippi State University, October 2002.
45! Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007.
46!-----------------------------------------------------------------------
47!
48 real :: qcen ! A layer's thickness divided by the distance between the centers
49 ! of the adjacent cells, usually ~0.5, but always <= 1 [nondim]
50 real :: zbot, zcen, ztop ! Tracer slopes times the layer thickness [A]
51 integer :: i, k
52
53 do i=1,ns
54 slope(1, i) = 0.0
55 slope(nk,i) = 0.0
56 enddo !i
57 do k= 2,nk-1
58 if (dpi(k) <= thin) then !use PCM
59 do i=1,ns ; slope(k,i) = 0.0 ; enddo
60 else
61! --- use qcen in place of 0.5 to allow for non-uniform grid
62 qcen = dpi(k) / (dpi(k)+0.5*(dpi(k-1)+dpi(k+1))) !dpi(k)>thin
63 do i=1,ns
64! --- PLM (non-zero slope, but no new extrema)
65! --- layer value is si-0.5*slope at top interface,
66! --- and si+0.5*slope at bottom interface.
67!
68! --- monotonized central-difference limiter (van Leer, 1977,
69! --- JCP 23 pp 276-299). For a discussion of PLM limiters, see
70! --- Finite Volume Methods for Hyperbolic Problems by R.J. Leveque.
71 ztop = 2.0*(si(k, i)-si(k-1,i))
72 zbot = 2.0*(si(k+1,i)-si(k, i))
73 zcen = qcen*(si(k+1,i)-si(k-1,i))
74 if (ztop*zbot > 0.0) then !ztop,zbot are the same sign
75 slope(k,i) = sign(min(abs(zcen),abs(zbot),abs(ztop)), zbot)
76 else
77 slope(k,i) = 0.0 !local extrema, so no slope
78 endif
79 enddo !i
80 endif !PCM:PLM
81 enddo !k
82
83 if (present(pcm_lay)) then
84 do k=1,nk ; if (pcm_lay(k)) then
85 do i=1,ns ; slope(k,i) = 0.0 ; enddo
86 endif ; enddo
87 endif
88
89end subroutine hybgen_plm_coefs
90
91
92!> Set up the coefficients for PPM remapping of a set of scalars
93subroutine hybgen_ppm_coefs(s, h_src, edges, nk, ns, thin, PCM_lay)
94 integer, intent(in) :: nk !< The number of input layers
95 integer, intent(in) :: ns !< The scalar fields to work on
96 real, intent(in) :: s(nk,ns) !< The input scalar fields [A]
97 real, intent(in) :: h_src(nk) !< The input grid layer thicknesses [H ~> m or kg m-2]
98 real, intent(out) :: edges(nk,2,ns) !< The PPM interpolation edge values of the scalar fields [A]
99 real, intent(in) :: thin !< A negligible layer thickness that can be ignored [H ~> m or kg m-2]
100 logical, optional, intent(in) :: pcm_lay(nk) !< If true for a layer, use PCM remapping for that layer
101
102!-----------------------------------------------------------------------
103! 1) coefficients for remapping from one set of vertical cells to another.
104! method: monotonic piecewise parabolic across each input cell
105!
106! Colella, P. & P.R. Woodward, 1984, J. Comp. Phys., 54, 174-201.
107!
108! 2) input arguments:
109! s - initial scalar fields in pi-layer space
110! h_src - initial layer thicknesses (>=0)
111! nk - number of layers
112! ns - number of fields
113! thin - layer thickness (>0) that can be ignored
114! PCM_lay - use PCM for selected layers (optional)
115!
116! 3) output arguments:
117! edges - cell edge scalar values for the PPM reconstruction
118! edges.1 is value at interface above
119! edges.2 is value at interface below
120!
121! 4) Tim Campbell, Mississippi State University, October 2002.
122! Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007.
123!-----------------------------------------------------------------------
124!
125 real :: dp(nk) ! Input grid layer thicknesses, but with a minimum thickness given by thin [H ~> m or kg m-2]
126 logical :: pcm_layer(nk) ! True for layers that should use PCM remapping, either because they are
127 ! very thin, or because this is specified by PCM_lay.
128 real :: da ! Difference between the unlimited scalar edge value estimates [A]
129 real :: a6 ! Scalar field differences that are proportional to the curvature [A]
130 real :: slk, srk ! Differences between adjacent cell averages of scalars [A]
131 real :: sck ! Scalar differences across a cell [A]
132 real :: as(nk) ! Scalar field difference across each cell [A]
133 real :: al(nk), ar(nk) ! Scalar field at the left and right edges of a cell [A]
134 real :: h112(nk+1), h122(nk+1) ! Combinations of thicknesses [H ~> m or kg m-2]
135 real :: i_h12(nk+1) ! Inverses of combinations of thickesses [H-1 ~> m-1 or m2 kg-1]
136 real :: h2_h123(nk) ! A ratio of a layer thickness of the sum of 3 adjacent thicknesses [nondim]
137 real :: i_h0123(nk) ! Inverse of the sum of 4 adjacent thicknesses [H-1 ~> m-1 or m2 kg-1]
138 real :: h01_h112(nk+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses.
139 real :: h23_h122(nk+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses.
140 integer :: k, i
141
142 ! This PPM remapper is not currently written to work with massless layers, so set
143 ! the thicknesses for very thin layers to some minimum value.
144 do k=1,nk ; dp(k) = max(h_src(k), thin) ; enddo
145
146 ! Specify the layers that will use PCM remapping.
147 if (present(pcm_lay)) then
148 do k=1,nk ; pcm_layer(k) = (pcm_lay(k) .or. dp(k) <= thin) ; enddo
149 else
150 do k=1,nk ; pcm_layer(k) = (dp(k) <= thin) ; enddo
151 endif
152
153 !compute grid metrics
154 do k=2,nk
155 h112(k) = 2.*dp(k-1) + dp(k)
156 h122(k) = dp(k-1) + 2.*dp(k)
157 i_h12(k) = 1.0 / (dp(k-1) + dp(k))
158 enddo !k
159 do k=2,nk-1
160 h2_h123(k) = dp(k) / (dp(k) + (dp(k-1)+dp(k+1)))
161 enddo
162 do k=3,nk-1
163 i_h0123(k) = 1.0 / ((dp(k-2) + dp(k-1)) + (dp(k) + dp(k+1)))
164
165 h01_h112(k) = (dp(k-2) + dp(k-1)) / (2.0*dp(k-1) + dp(k))
166 h23_h122(k) = (dp(k) + dp(k+1)) / (dp(k-1) + 2.0*dp(k))
167 enddo
168
169 do i=1,ns
170 !Compute average slopes: Colella, Eq. (1.8)
171 as(1) = 0.
172 do k=2,nk-1
173 if (pcm_layer(k)) then !use PCM
174 as(k) = 0.0
175 else
176 slk = s(k, i)-s(k-1,i)
177 srk = s(k+1,i)-s(k, i)
178 if (slk*srk > 0.) then
179 sck = h2_h123(k)*( h112(k)*srk*i_h12(k+1) + h122(k+1)*slk*i_h12(k) )
180 as(k) = sign(min(abs(2.0*slk), abs(sck), abs(2.0*srk)), sck)
181 else
182 as(k) = 0.
183 endif
184 endif !PCM:PPM
185 enddo !k
186 as(nk) = 0.
187 !Compute "first guess" edge values: Colella, Eq. (1.6)
188 al(1) = s(1,i) ! 1st layer PCM
189 ar(1) = s(1,i) ! 1st layer PCM
190 al(2) = s(1,i) ! 1st layer PCM
191 do k=3,nk-1
192 ! This is a 4th order explicit edge value estimate.
193 al(k) = (dp(k)*s(k-1,i) + dp(k-1)*s(k,i)) * i_h12(k) &
194 + i_h0123(k)*( 2.*dp(k)*dp(k-1)*i_h12(k)*(s(k,i)-s(k-1,i)) * &
195 ( h01_h112(k) - h23_h122(k) ) &
196 + (dp(k)*as(k-1)*h23_h122(k) - dp(k-1)*as(k)*h01_h112(k)) )
197 ar(k-1) = al(k)
198 enddo !k
199 ar(nk-1) = s(nk,i) ! last layer PCM
200 al(nk) = s(nk,i) ! last layer PCM
201 ar(nk) = s(nk,i) ! last layer PCM
202 !Impose monotonicity: Colella, Eq. (1.10)
203 do k=2,nk-1
204 if ((pcm_layer(k)) .or. ((s(k+1,i)-s(k,i))*(s(k,i)-s(k-1,i)) <= 0.)) then !local extremum
205 al(k) = s(k,i)
206 ar(k) = s(k,i)
207 else
208 da = ar(k)-al(k)
209 a6 = 6.0*s(k,i) - 3.0*(al(k)+ar(k))
210 if (da*a6 > da*da) then !peak in right half of zone
211 al(k) = 3.0*s(k,i) - 2.0*ar(k)
212 elseif (da*a6 < -da*da) then !peak in left half of zone
213 ar(k) = 3.0*s(k,i) - 2.0*al(k)
214 endif
215 endif
216 enddo !k
217 !Set coefficients
218 do k=1,nk
219 edges(k,1,i) = al(k)
220 edges(k,2,i) = ar(k)
221 enddo !k
222 enddo !i
223
224end subroutine hybgen_ppm_coefs
225
226
227!> Set up the coefficients for PPM remapping of a set of scalars
228subroutine hybgen_weno_coefs(s, h_src, edges, nk, ns, thin, PCM_lay)
229 integer, intent(in) :: nk !< The number of input layers
230 integer, intent(in) :: ns !< The number of scalar fields to work on
231 real, intent(in) :: s(nk,ns) !< The input scalar fields [A]
232 real, intent(in) :: h_src(nk) !< The input grid layer thicknesses [H ~> m or kg m-2]
233 real, intent(out) :: edges(nk,2,ns) !< The WENO interpolation edge values of the scalar fields [A]
234 real, intent(in) :: thin !< A negligible layer thickness that can be ignored [H ~> m or kg m-2]
235 logical, optional, intent(in) :: pcm_lay(nk) !< If true for a layer, use PCM remapping for that layer
236
237!-----------------------------------------------------------------------
238! 1) coefficients for remapping from one set of vertical cells to another.
239! method: monotonic WENO-like alternative to PPM across each input cell
240! a second order polynomial approximation of the profiles
241! using a WENO reconciliation of the slopes to compute the
242! interfacial values
243!
244! This scheme might have ben developed by Shchepetkin. A.F., personal communication.
245! See also Engwirda, D., and M. Kelley, A WENO-type slope-limiter for a family of piecewise
246! polynomial methods, arXive:1606.08188v1, 27 June 2016.
247!
248! 2) input arguments:
249! s - initial scalar fields in pi-layer space
250! h_src - initial layer thicknesses (>=0)
251! nk - number of layers
252! ns - number of fields
253! thin - layer thickness (>0) that can be ignored
254! PCM_lay - use PCM for selected layers (optional)
255!
256! 3) output arguments:
257! edges - cell edge scalar values for the WENO reconstruction
258! edges.1 is value at interface above
259! edges.2 is value at interface below
260!
261! 4) Laurent Debreu, Grenoble.
262! Alan J. Wallcraft, Naval Research Laboratory, July 2008.
263!-----------------------------------------------------------------------
264!
265! real, parameter :: dsmll=1.0e-8 ! This has units of [A2], and hence can not be a parameter.
266!
267 real :: curv_cell ! An estimate of the tracer curvature centered on a cell times the grid
268 ! spacing [A H-1 ~> A m-1 or A m2 kg-1]
269 real :: seh1, seh2 ! Tracer slopes at the cell edges times the cell grid spacing [A]
270 real :: q01, q02 ! Various tracer differences between a cell average and the edge values [A]
271 real :: q001, q002 ! Tracer slopes at the cell edges times the cell grid spacing [A]
272 logical :: pcm_layer(nk) ! True for layers that should use PCM remapping, either because they are
273 ! very thin, or because this is specified by PCM_lay.
274 real :: dp(nk) ! Input grid layer thicknesses, but with a minimum thickness given by thin [H ~> m or kg m-2]
275 real :: qdpkm(nk) ! Inverse of the sum of two adjacent thicknesses [H-1 ~> m-1 or m2 kg-1]
276 real :: qdpkmkp(nk) ! Inverse of the sum of three adjacent thicknesses [H-1 ~> m-1 or m2 kg-1]
277 real :: dpkm2kp(nk) ! Twice the distance between the centers of the layers two apart [H ~> m or kg m-2]
278 real :: zw(nk,2) ! Squared combinations of the differences between the cell average tracer
279 ! concentrations and the left and right edges [A2]
280 real :: min_ratio ! The minimum ratio of the values of zw used to interpolate the edge values [nondim]
281 real :: wt1 ! The weight of the upper layer in the interpolated shared edge value [nondim]
282 real :: slope_edge(nk+1) ! Tracer slopes at the edges [A H-1 ~> A m-1 or A m2 kg-1]
283 real :: val_edge(nk+1) ! A weighted average edge concentration [A]
284 integer :: i, k
285
286 min_ratio = 1.0e-8
287
288 ! The WENO remapper is not currently written to work with massless layers, so set
289 ! the thicknesses for very thin layers to some minimum value.
290 do k=1,nk ; dp(k) = max(h_src(k), thin) ; enddo
291
292 ! Specify the layers that will use PCM remapping.
293 if (present(pcm_lay)) then
294 do k=1,nk ; pcm_layer(k) = (pcm_lay(k) .or. dp(k) <= thin) ; enddo
295 else
296 do k=1,nk ; pcm_layer(k) = (dp(k) <= thin) ; enddo
297 endif
298
299 !compute grid metrics
300 do k=2,nk-1
301 qdpkm( k) = 1.0 / (dp(k-1) + dp(k))
302 qdpkmkp(k) = 1.0 / (dp(k-1) + dp(k) + dp(k+1))
303 dpkm2kp(k) = dp(k-1) + 2.0*dp(k) + dp(k+1)
304 enddo !k
305 qdpkm(nk) = 1.0 / (dp(nk-1) + dp(nk))
306
307 do i=1,ns
308 do k=2,nk
309 slope_edge(k) = qdpkm(k) * (s(k,i)-s(k-1,i))
310 enddo !k
311 k = 1 !PCM first layer
312 edges(k,1,i) = s(k,i)
313 edges(k,2,i) = s(k,i)
314 zw(k,1) = 0.0
315 zw(k,2) = 0.0
316 do k=2,nk-1
317 if ((slope_edge(k)*slope_edge(k+1) < 0.0) .or. pcm_layer(k)) then !use PCM
318 edges(k,1,i) = s(k,i)
319 edges(k,2,i) = s(k,i)
320 zw(k,1) = 0.0
321 zw(k,2) = 0.0
322 else
323 seh1 = dp(k)*slope_edge(k+1)
324 seh2 = dp(k)*slope_edge(k)
325 q01 = dpkm2kp(k)*slope_edge(k+1)
326 q02 = dpkm2kp(k)*slope_edge(k)
327 if (abs(seh1) > abs(q02)) then
328 seh1 = q02
329 endif
330 if (abs(seh2) > abs(q01)) then
331 seh2 = q01
332 endif
333 curv_cell = (seh1 - seh2) * qdpkmkp(k)
334 q001 = seh1 - curv_cell*dp(k+1)
335 q002 = seh2 + curv_cell*dp(k-1)
336 ! q001 = (seh1 * (dp(k-1) + dp(k)) + seh2 * dp(k+1)) * qdpkmkp(k)
337 ! q002 = (seh2 * (dp(k+1) + dp(k)) + seh1 * dp(k-1)) * qdpkmkp(k)
338
339 edges(k,2,i) = s(k,i) + q001
340 edges(k,1,i) = s(k,i) - q002
341 zw(k,1) = (2.0*q001 - q002)**2
342 zw(k,2) = (2.0*q002 - q001)**2
343 endif !PCM:WENO
344 enddo !k
345 k = nk !PCM last layer
346 edges(k,1,i) = s(k,i)
347 edges(k,2,i) = s(k,i)
348 zw(k, 1) = 0.0
349 zw(k, 2) = 0.0
350
351 do k=2,nk
352 ! This was the original code based on that in Hycom, but because zw has
353 ! dimensions of [A2], it can not use a constant (hard coded) value of dsmll.
354 ! ds2a = max(zw(k-1,2), dsmll)
355 ! ds2b = max(zw(k, 1), dsmll)
356 ! val_edge(K) = (ds2b*edges(k-1,2,i)+ds2a*edges(k,1,i)) / (ds2b+ds2a)
357 ! Use a weighted average of the two layers' estimated edge values as the actual edge value.
358 if (zw(k,1) + zw(k-1,2) <= 0.0) then
359 wt1 = 0.5
360 elseif (zw(k,1) <= min_ratio * (zw(k,1) + zw(k-1,2))) then
361 wt1 = min_ratio
362 elseif (zw(k-1,2) <= min_ratio * (zw(k,1) + zw(k-1,2))) then
363 wt1 = (1.0 - min_ratio)
364 else
365 wt1 = zw(k,1) / (zw(k,1) + zw(k-1,2))
366 endif
367 val_edge(k) = wt1*edges(k-1,2,i) + (1.0-wt1)*edges(k,1,i)
368 enddo !k
369 val_edge( 1) = 2.0*s( 1,i)-val_edge( 2) !not used?
370 val_edge(nk+1) = 2.0*s(nk,i)-val_edge(nk) !not used?
371
372 do k=2,nk-1
373 if (.not.pcm_layer(k)) then !don't use PCM
374 q01 = val_edge(k+1) - s(k,i)
375 q02 = s(k,i) - val_edge(k)
376 if (q01*q02 < 0.0) then
377 q01 = 0.0
378 q02 = 0.0
379 elseif (abs(q01) > abs(2.0*q02)) then
380 q01 = 2.0*q02
381 elseif (abs(q02) > abs(2.0*q01)) then
382 q02 = 2.0*q01
383 endif
384 edges(k,1,i) = s(k,i) - q02
385 edges(k,2,i) = s(k,i) + q01
386 endif ! PCM:WENO
387 enddo !k
388 enddo !i
389
390end subroutine hybgen_weno_coefs
391
392end module mom_hybgen_remap