MOM6
MOM_remapping.F90
Go to the documentation of this file.
1 !> Provides column-wise vertical remapping functions
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 ! Original module written by Laurent White, 2008.06.09
6 
7 use mom_error_handler, only : mom_error, fatal
16 
17 implicit none ; private
18 
19 #include <MOM_memory.h>
20 
21 !> Container for remapping parameters
22 type, public :: remapping_cs
23  private
24  !> Determines which reconstruction to use
25  integer :: remapping_scheme = -911
26  !> Degree of polynomial reconstruction
27  integer :: degree = 0
28  !> If true, extrapolate boundaries
29  logical :: boundary_extrapolation = .true.
30  !> If true, reconstructions are checked for consistency.
31  logical :: check_reconstruction = .false.
32  !> If true, the result of remapping are checked for conservation and bounds.
33  logical :: check_remapping = .false.
34  !> If true, the intermediate values used in remapping are forced to be bounded.
35  logical :: force_bounds_in_subcell = .false.
36  !> If true use older, less acccurate expressions.
37  logical :: answers_2018 = .true.
38 end type
39 
40 ! The following routines are visible to the outside world
44 public dzfromh1h2
45 
46 ! The following are private parameter constants
47 integer, parameter :: remapping_pcm = 0 !< O(h^1) remapping scheme
48 integer, parameter :: remapping_plm = 1 !< O(h^2) remapping scheme
49 integer, parameter :: remapping_ppm_h4 = 2 !< O(h^3) remapping scheme
50 integer, parameter :: remapping_ppm_ih4 = 3 !< O(h^3) remapping scheme
51 integer, parameter :: remapping_pqm_ih4ih3 = 4 !< O(h^4) remapping scheme
52 integer, parameter :: remapping_pqm_ih6ih5 = 5 !< O(h^5) remapping scheme
53 
54 integer, parameter :: integration_pcm = 0 !< Piecewise Constant Method
55 integer, parameter :: integration_plm = 1 !< Piecewise Linear Method
56 integer, parameter :: integration_ppm = 3 !< Piecewise Parabolic Method
57 integer, parameter :: integration_pqm = 5 !< Piecewise Quartic Method
58 
59 character(len=40) :: mdl = "MOM_remapping" !< This module's name.
60 
61 !> Documentation for external callers
62 character(len=256), public :: remappingschemesdoc = &
63  "PCM (1st-order accurate)\n"//&
64  "PLM (2nd-order accurate)\n"//&
65  "PPM_H4 (3rd-order accurate)\n"//&
66  "PPM_IH4 (3rd-order accurate)\n"//&
67  "PQM_IH4IH3 (4th-order accurate)\n"//&
68  "PQM_IH6IH5 (5th-order accurate)\n"
69 character(len=3), public :: remappingdefaultscheme = "PLM" !< Default remapping method
70 
71 ! This CPP macro turns on/off bounding of integrations limits so that they are
72 ! always within the cell. Roundoff can lead to the non-dimensional bounds being
73 ! outside of the range 0 to 1.
74 #define __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
75 
76 real, parameter :: hneglect_dflt = 1.e-30 !< A thickness [H ~> m or kg m-2] that can be
77  !! added to thicknesses in a denominator without
78  !! changing the numerical result, except where
79  !! a division by zero would otherwise occur.
80 
81 logical, parameter :: old_algorithm = .false. !< Use the old "broken" algorithm.
82  !! This is a temporary measure to assist
83  !! debugging until we delete the old algorithm.
84 
85 contains
86 
87 !> Set parameters within remapping object
88 subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, &
89  check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
90  type(remapping_cs), intent(inout) :: cs !< Remapping control structure
91  character(len=*), optional, intent(in) :: remapping_scheme !< Remapping scheme to use
92  logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
93  logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions
94  logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping
95  logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded
96  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
97 
98  if (present(remapping_scheme)) then
99  call setreconstructiontype( remapping_scheme, cs )
100  endif
101  if (present(boundary_extrapolation)) then
102  cs%boundary_extrapolation = boundary_extrapolation
103  endif
104  if (present(check_reconstruction)) then
105  cs%check_reconstruction = check_reconstruction
106  endif
107  if (present(check_remapping)) then
108  cs%check_remapping = check_remapping
109  endif
110  if (present(force_bounds_in_subcell)) then
111  cs%force_bounds_in_subcell = force_bounds_in_subcell
112  endif
113  if (present(answers_2018)) then
114  cs%answers_2018 = answers_2018
115  endif
116 end subroutine remapping_set_param
117 
118 subroutine extract_member_remapping_cs(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, &
119  check_remapping, force_bounds_in_subcell)
120  type(remapping_cs), intent(in) :: cs !< Control structure for remapping module
121  integer, optional, intent(out) :: remapping_scheme !< Determines which reconstruction scheme to use
122  integer, optional, intent(out) :: degree !< Degree of polynomial reconstruction
123  logical, optional, intent(out) :: boundary_extrapolation !< If true, extrapolate boundaries
124  logical, optional, intent(out) :: check_reconstruction !< If true, reconstructions are checked for consistency.
125  logical, optional, intent(out) :: check_remapping !< If true, the result of remapping are checked
126  !! for conservation and bounds.
127  logical, optional, intent(out) :: force_bounds_in_subcell !< If true, the intermediate values used in
128  !! remapping are forced to be bounded.
129 
130  if (present(remapping_scheme)) remapping_scheme = cs%remapping_scheme
131  if (present(degree)) degree = cs%degree
132  if (present(boundary_extrapolation)) boundary_extrapolation = cs%boundary_extrapolation
133  if (present(check_reconstruction)) check_reconstruction = cs%check_reconstruction
134  if (present(check_remapping)) check_remapping = cs%check_remapping
135  if (present(force_bounds_in_subcell)) force_bounds_in_subcell = cs%force_bounds_in_subcell
136 
137 end subroutine extract_member_remapping_cs
138 !> Calculate edge coordinate x from cell width h
139 subroutine buildgridfromh(nz, h, x)
140  integer, intent(in) :: nz !< Number of cells
141  real, dimension(nz), intent(in) :: h !< Cell widths
142  real, dimension(nz+1), intent(inout) :: x !< Edge coordiantes starting at x(1)=0
143  ! Local variables
144  integer :: k
145 
146  x(1) = 0.0
147  do k = 1,nz
148  x(k+1) = x(k) + h(k)
149  enddo
150 
151 end subroutine buildgridfromh
152 
153 !> Compare two summation estimates of positive data and judge if due to more
154 !! than round-off.
155 !! When two sums are calculated from different vectors that should add up to
156 !! the same value, the results can differ by round off. The round off error
157 !! can be bounded to be proportional to the number of operations.
158 !! This function returns true if the difference between sum1 and sum2 is
159 !! larger than than the estimated round off bound.
160 !! \note This estimate/function is only valid for summation of positive data.
161 function ispossumerrsignificant(n1, sum1, n2, sum2)
162  integer, intent(in) :: n1 !< Number of values in sum1
163  integer, intent(in) :: n2 !< Number of values in sum2
164  real, intent(in) :: sum1 !< Sum of n1 values
165  real, intent(in) :: sum2 !< Sum of n2 values
166  logical :: ispossumerrsignificant !< True if difference in sums is large
167  ! Local variables
168  real :: sumerr, allowederr, eps
169 
170  if (sum1<0.) call mom_error(fatal,'isPosSumErrSignificant: sum1<0 is not allowed!')
171  if (sum2<0.) call mom_error(fatal,'isPosSumErrSignificant: sum2<0 is not allowed!')
172  sumerr = abs(sum1-sum2)
173  eps = epsilon(sum1)
174  allowederr = eps*0.5*(real(n1-1)*sum1+real(n2-1)*sum2)
175  if (sumerr>allowederr) then
176  write(0,*) 'isPosSumErrSignificant: sum1,sum2=',sum1,sum2
177  write(0,*) 'isPosSumErrSignificant: eps=',eps
178  write(0,*) 'isPosSumErrSignificant: err,n*eps=',sumerr,allowederr
179  write(0,*) 'isPosSumErrSignificant: err/eps,n1,n2,n1+n2=',sumerr/eps,n1,n2,n1+n2
180  ispossumerrsignificant = .true.
181  else
182  ispossumerrsignificant = .false.
183  endif
184 end function ispossumerrsignificant
185 
186 !> Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned.
187 subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
188  type(remapping_cs), intent(in) :: cs !< Remapping control structure
189  integer, intent(in) :: n0 !< Number of cells on source grid
190  real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid
191  real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid
192  integer, intent(in) :: n1 !< Number of cells on target grid
193  real, dimension(n1), intent(in) :: h1 !< Cell widths on target grid
194  real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid
195  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
196  !! purpose of cell reconstructions
197  !! in the same units as h0.
198  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
199  !! for the purpose of edge value
200  !! calculations in the same units as h0.
201  ! Local variables
202  integer :: imethod
203  real, dimension(n0,2) :: ppoly_r_e !Edge value of polynomial
204  real, dimension(n0,2) :: ppoly_r_s !Edge slope of polynomial
205  real, dimension(n0,CS%degree+1) :: ppoly_r_coefs !Coefficients of polynomial
206  integer :: k
207  real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
208  real :: hneglect, hneglect_edge
209 
210  hneglect = 1.0e-30 ; if (present(h_neglect)) hneglect = h_neglect
211  hneglect_edge = 1.0e-10 ; if (present(h_neglect_edge)) hneglect_edge = h_neglect_edge
212 
213  call build_reconstructions_1d( cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod, &
214  hneglect, hneglect_edge )
215 
216  if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &
217  cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
218 
219 
220  call remap_via_sub_cells( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h1, imethod, &
221  cs%force_bounds_in_subcell, u1, uh_err )
222 
223  if (cs%check_remapping) then
224  ! Check errors and bounds
225  call measure_input_bounds( n0, h0, u0, ppoly_r_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
226  call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
227  if (imethod<5) then ! We except PQM until we've debugged it
228  if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err .and. abs(h1tot-h0tot)<h0err+h1err) &
229  .or. (u1min<u0min .or. u1max>u0max) ) then
230  write(0,*) 'iMethod = ',imethod
231  write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err
232  if (abs(h1tot-h0tot)>h0err+h1err) &
233  write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!'
234  write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err
235  if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
236  write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!'
237  write(0,*) 'U: u0min=',u0min,'u1min=',u1min
238  if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!'
239  write(0,*) 'U: u0max=',u0max,'u1max=',u1max
240  if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!'
241  write(0,'(a3,6a24)') 'k','h0','left edge','u0','right edge','h1','u1'
242  do k = 1, max(n0,n1)
243  if (k<=min(n0,n1)) then
244  write(0,'(i3,1p6e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2),h1(k),u1(k)
245  elseif (k>n0) then
246  write(0,'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
247  else
248  write(0,'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
249  endif
250  enddo
251  write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients'
252  do k = 1, n0
253  write(0,'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
254  enddo
255  call mom_error( fatal, 'MOM_remapping, remapping_core_h: '//&
256  'Remapping result is inconsistent!' )
257  endif
258  endif ! method<5
259  endif
260 
261 end subroutine remapping_core_h
262 
263 !> Remaps column of values u0 on grid h0 to implied grid h1
264 !! where the interfaces of h1 differ from those of h0 by dx.
265 subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_edge )
266  type(remapping_cs), intent(in) :: cs !< Remapping control structure
267  integer, intent(in) :: n0 !< Number of cells on source grid
268  real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid
269  real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid
270  integer, intent(in) :: n1 !< Number of cells on target grid
271  real, dimension(n1+1), intent(in) :: dx !< Cell widths on target grid
272  real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid
273  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
274  !! purpose of cell reconstructions
275  !! in the same units as h0.
276  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
277  !! for the purpose of edge value
278  !! calculations in the same units as h0.
279  ! Local variables
280  integer :: imethod
281  real, dimension(n0,2) :: ppoly_r_e !Edge value of polynomial
282  real, dimension(n0,2) :: ppoly_r_s !Edge slope of polynomial
283  real, dimension(n0,CS%degree+1) :: ppoly_r_coefs !Coefficients of polynomial
284  integer :: k
285  real :: eps, h0tot, h0err, h1tot, h1err
286  real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
287  real, dimension(n1) :: h1 !< Cell widths on target grid
288  real :: hneglect, hneglect_edge
289 
290  hneglect = 1.0e-30 ; if (present(h_neglect)) hneglect = h_neglect
291  hneglect_edge = 1.0e-10 ; if (present(h_neglect_edge)) hneglect_edge = h_neglect_edge
292 
293  call build_reconstructions_1d( cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod,&
294  hneglect, hneglect_edge )
295 
296  if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &
297  cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
298 
299  ! This is a temporary step prior to switching to remapping_core_h()
300  do k = 1, n1
301  if (k<=n0) then
302  h1(k) = max( 0., h0(k) + ( dx(k+1) - dx(k) ) )
303  else
304  h1(k) = max( 0., dx(k+1) - dx(k) )
305  endif
306  enddo
307  call remap_via_sub_cells( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h1, imethod, &
308  cs%force_bounds_in_subcell,u1, uh_err )
309 ! call remapByDeltaZ( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, dx, iMethod, u1, hNeglect )
310 ! call remapByProjection( n0, h0, u0, CS%ppoly_r, n1, h1, iMethod, u1, hNeglect )
311 
312  if (cs%check_remapping) then
313  ! Check errors and bounds
314  call measure_input_bounds( n0, h0, u0, ppoly_r_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
315  call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
316  if (imethod<5) then ! We except PQM until we've debugged it
317  if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err .and. abs(h1tot-h0tot)<h0err+h1err) &
318  .or. (u1min<u0min .or. u1max>u0max) ) then
319  write(0,*) 'iMethod = ',imethod
320  write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err
321  if (abs(h1tot-h0tot)>h0err+h1err) &
322  write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!'
323  write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err
324  if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
325  write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!'
326  write(0,*) 'U: u0min=',u0min,'u1min=',u1min
327  if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!'
328  write(0,*) 'U: u0max=',u0max,'u1max=',u1max
329  if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!'
330  write(0,'(a3,6a24)') 'k','h0','left edge','u0','right edge','h1','u1'
331  do k = 1, max(n0,n1)
332  if (k<=min(n0,n1)) then
333  write(0,'(i3,1p6e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2),h1(k),u1(k)
334  elseif (k>n0) then
335  write(0,'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
336  else
337  write(0,'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
338  endif
339  enddo
340  write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients'
341  do k = 1, n0
342  write(0,'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
343  enddo
344  call mom_error( fatal, 'MOM_remapping, remapping_core_w: '//&
345  'Remapping result is inconsistent!' )
346  endif
347  endif ! method<5
348  endif
349 
350 end subroutine remapping_core_w
351 
352 !> Creates polynomial reconstructions of u0 on the source grid h0.
353 subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, &
354  ppoly_r_E, ppoly_r_S, iMethod, h_neglect, &
355  h_neglect_edge )
356  type(remapping_cs), intent(in) :: cs !< Remapping control structure
357  integer, intent(in) :: n0 !< Number of cells on source grid
358  real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid
359  real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid
360  real, dimension(n0,CS%degree+1), &
361  intent(out) :: ppoly_r_coefs !< Coefficients of polynomial
362  real, dimension(n0,2), intent(out) :: ppoly_r_e !< Edge value of polynomial
363  real, dimension(n0,2), intent(out) :: ppoly_r_s !< Edge slope of polynomial
364  integer, intent(out) :: imethod !< Integration method
365  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
366  !! purpose of cell reconstructions
367  !! in the same units as h0.
368  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
369  !! for the purpose of edge value
370  !! calculations in the same units as h0.
371  ! Local variables
372  integer :: local_remapping_scheme
373  integer :: remapping_scheme !< Remapping scheme
374  logical :: boundary_extrapolation !< Extrapolate at boundaries if true
375 
376  ! Reset polynomial
377  ppoly_r_e(:,:) = 0.0
378  ppoly_r_s(:,:) = 0.0
379  ppoly_r_coefs(:,:) = 0.0
380  imethod = -999
381 
382  local_remapping_scheme = cs%remapping_scheme
383  if (n0<=1) then
384  local_remapping_scheme = remapping_pcm
385  elseif (n0<=3) then
386  local_remapping_scheme = min( local_remapping_scheme, remapping_plm )
387  elseif (n0<=4) then
388  local_remapping_scheme = min( local_remapping_scheme, remapping_ppm_h4 )
389  endif
390  select case ( local_remapping_scheme )
391  case ( remapping_pcm )
392  call pcm_reconstruction( n0, u0, ppoly_r_e, ppoly_r_coefs)
393  imethod = integration_pcm
394  case ( remapping_plm )
395  call plm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
396  if ( cs%boundary_extrapolation ) then
397  call plm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect)
398  endif
399  imethod = integration_plm
400  case ( remapping_ppm_h4 )
401  call edge_values_explicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
402  call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
403  if ( cs%boundary_extrapolation ) then
404  call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
405  endif
406  imethod = integration_ppm
407  case ( remapping_ppm_ih4 )
408  call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
409  call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
410  if ( cs%boundary_extrapolation ) then
411  call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
412  endif
413  imethod = integration_ppm
414  case ( remapping_pqm_ih4ih3 )
415  call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
416  call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_s, h_neglect, answers_2018=cs%answers_2018 )
417  call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect )
418  if ( cs%boundary_extrapolation ) then
419  call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &
420  ppoly_r_coefs, h_neglect )
421  endif
422  imethod = integration_pqm
423  case ( remapping_pqm_ih6ih5 )
424  call edge_values_implicit_h6( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
425  call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_s, h_neglect, answers_2018=cs%answers_2018 )
426  call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect )
427  if ( cs%boundary_extrapolation ) then
428  call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &
429  ppoly_r_coefs, h_neglect )
430  endif
431  imethod = integration_pqm
432  case default
433  call mom_error( fatal, 'MOM_remapping, build_reconstructions_1d: '//&
434  'The selected remapping method is invalid' )
435  end select
436 
437 end subroutine build_reconstructions_1d
438 
439 !> Checks that edge values and reconstructions satisfy bounds
440 subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, &
441  ppoly_r_coefs, ppoly_r_E, ppoly_r_S)
442  integer, intent(in) :: n0 !< Number of cells on source grid
443  real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid
444  real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid
445  integer, intent(in) :: deg !< Degree of polynomial reconstruction
446  logical, intent(in) :: boundary_extrapolation !< Extrapolate at boundaries if true
447  real, dimension(n0,deg+1),intent(out) :: ppoly_r_coefs !< Coefficients of polynomial
448  real, dimension(n0,2), intent(out) :: ppoly_r_E !< Edge value of polynomial
449  real, dimension(n0,2), intent(out) :: ppoly_r_S !< Edge slope of polynomial
450  ! Local variables
451  integer :: i0, n
452  real :: u_l, u_c, u_r ! Cell averages
453  real :: u_min, u_max
454  logical :: problem_detected
455 
456  problem_detected = .false.
457  do i0 = 1, n0
458  u_l = u0(max(1,i0-1))
459  u_c = u0(i0)
460  u_r = u0(min(n0,i0+1))
461  if (i0 > 1 .or. .not. boundary_extrapolation) then
462  u_min = min(u_l, u_c)
463  u_max = max(u_l, u_c)
464  if (ppoly_r_e(i0,1) < u_min) then
465  write(0,'(a,i4,5(x,a,1pe24.16))') 'Left edge undershoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &
466  'edge=',ppoly_r_e(i0,1),'err=',ppoly_r_e(i0,1)-u_min
467  problem_detected = .true.
468  endif
469  if (ppoly_r_e(i0,1) > u_max) then
470  write(0,'(a,i4,5(x,a,1pe24.16))') 'Left edge overshoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &
471  'edge=',ppoly_r_e(i0,1),'err=',ppoly_r_e(i0,1)-u_max
472  problem_detected = .true.
473  endif
474  endif
475  if (i0 < n0 .or. .not. boundary_extrapolation) then
476  u_min = min(u_c, u_r)
477  u_max = max(u_c, u_r)
478  if (ppoly_r_e(i0,2) < u_min) then
479  write(0,'(a,i4,5(x,a,1pe24.16))') 'Right edge undershoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, &
480  'edge=',ppoly_r_e(i0,2),'err=',ppoly_r_e(i0,2)-u_min
481  problem_detected = .true.
482  endif
483  if (ppoly_r_e(i0,2) > u_max) then
484  write(0,'(a,i4,5(x,a,1pe24.16))') 'Right edge overshoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, &
485  'edge=',ppoly_r_e(i0,2),'err=',ppoly_r_e(i0,2)-u_max
486  problem_detected = .true.
487  endif
488  endif
489  if (i0 > 1) then
490  if ( (u_c-u_l)*(ppoly_r_e(i0,1)-ppoly_r_e(i0-1,2)) < 0.) then
491  write(0,'(a,i4,5(x,a,1pe24.16))') 'Non-monotonic edges at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &
492  'right edge=',ppoly_r_e(i0-1,2),'left edge=',ppoly_r_e(i0,1)
493  write(0,'(5(a,1pe24.16,x))') 'u(i0)-u(i0-1)',u_c-u_l,'edge diff=',ppoly_r_e(i0,1)-ppoly_r_e(i0-1,2)
494  problem_detected = .true.
495  endif
496  endif
497  if (problem_detected) then
498  write(0,'(a,1p9e24.16)') 'Polynomial coeffs:',ppoly_r_coefs(i0,:)
499  write(0,'(3(a,1pe24.16,x))') 'u_l=',u_l,'u_c=',u_c,'u_r=',u_r
500  write(0,'(a4,10a24)') 'i0','h0(i0)','u0(i0)','left edge','right edge','Polynomial coefficients'
501  do n = 1, n0
502  write(0,'(i4,1p10e24.16)') n,h0(n),u0(n),ppoly_r_e(n,1),ppoly_r_e(n,2),ppoly_r_coefs(n,:)
503  enddo
504  call mom_error(fatal, 'MOM_remapping, check_reconstructions_1d: '// &
505  'Edge values or polynomial coefficients were inconsistent!')
506  endif
507  enddo
508 
509 end subroutine check_reconstructions_1d
510 
511 !> Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating
512 !! the n0+n1+1 sub-integrals of the intersection of h0 and h1, and the summing the
513 !! appropriate integrals into the h1*u1 values. h0 and h1 must have the same units.
514 subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, &
515  force_bounds_in_subcell, u1, uh_err, ah_sub, aisub_src, aiss, aise )
516  integer, intent(in) :: n0 !< Number of cells in source grid
517  real, intent(in) :: h0(n0) !< Source grid widths (size n0)
518  real, intent(in) :: u0(n0) !< Source cell averages (size n0)
519  real, intent(in) :: ppoly0_E(n0,2) !< Edge value of polynomial
520  real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial
521  integer, intent(in) :: n1 !< Number of cells in target grid
522  real, intent(in) :: h1(n1) !< Target grid widths (size n1)
523  integer, intent(in) :: method !< Remapping scheme to use
524  logical, intent(in) :: force_bounds_in_subcell !< Force sub-cell values to be bounded
525  real, intent(out) :: u1(n1) !< Target cell averages (size n1)
526  real, intent(out) :: uh_err !< Estimate of bound on error in sum of u*h
527  real, optional, intent(out) :: ah_sub(n0+n1+1) !< h_sub
528  integer, optional, intent(out) :: aisub_src(n0+n1+1) !< i_sub_src
529  integer, optional, intent(out) :: aiss(n0) !< isrc_start
530  integer, optional, intent(out) :: aise(n0) !< isrc_ens
531  ! Local variables
532  integer :: i_sub ! Index of sub-cell
533  integer :: i0 ! Index into h0(1:n0), source column
534  integer :: i1 ! Index into h1(1:n1), target column
535  integer :: i_start0 ! Used to record which sub-cells map to source cells
536  integer :: i_start1 ! Used to record which sub-cells map to target cells
537  integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell
538  real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell
539  real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell
540  real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell
541  real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell
542  integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell
543  integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell
544  integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell
545  integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell
546  real, dimension(n0) :: h0_eff ! Effective thickness of source cells
547  real, dimension(n0) :: u0_min ! Minimum value of reconstructions in source cell
548  real, dimension(n0) :: u0_max ! Minimum value of reconstructions in source cell
549  integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell
550  integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell
551  real :: xa, xb ! Non-dimensional position within a source cell (0..1)
552  real :: h0_supply, h1_supply ! The amount of width available for constructing sub-cells
553  real :: dh ! The width of the sub-cell
554  real :: duh ! The total amount of accumulated stuff (u*h)
555  real :: dh0_eff ! Running sum of source cell thickness
556  ! For error checking/debugging
557  logical, parameter :: force_bounds_in_target = .true. ! To fix round-off issues
558  logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues
559  logical, parameter :: debug_bounds = .false. ! For debugging overshoots etc.
560  integer :: k, i0_last_thick_cell
561  real :: h0tot, h0err, h1tot, h1err, h2tot, h2err, u02_err
562  real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, u2tot, u2err, u2min, u2max, u_orig
563  logical :: src_has_volume !< True if h0 has not been consumed
564  logical :: tgt_has_volume !< True if h1 has not been consumed
565 
566  if (old_algorithm) isrc_max(:)=1
567 
568  i0_last_thick_cell = 0
569  do i0 = 1, n0
570  u0_min(i0) = min(ppoly0_e(i0,1), ppoly0_e(i0,2))
571  u0_max(i0) = max(ppoly0_e(i0,1), ppoly0_e(i0,2))
572  if (h0(i0)>0.) i0_last_thick_cell = i0
573  enddo
574 
575  ! Initialize algorithm
576  h0_supply = h0(1)
577  h1_supply = h1(1)
578  src_has_volume = .true.
579  tgt_has_volume = .true.
580  i0 = 1 ; i1 = 1
581  i_start0 = 1 ; i_start1 = 1
582  i_max = 1
583  dh_max = 0.
584  dh0_eff = 0.
585 
586  ! First sub-cell is always vanished
587  h_sub(1) = 0.
588  isrc_start(1) = 1
589  isrc_end(1) = 1
590  isrc_max(1) = 1
591  isub_src(1) = 1
592 
593  ! Loop over each sub-cell to calculate intersections with source and target grids
594  do i_sub = 2, n0+n1+1
595 
596  ! This is the width of the sub-cell, determined by which ever column has the least
597  ! supply available to consume.
598  dh = min(h0_supply, h1_supply)
599 
600  ! This is the running sum of the source cell thickness. After summing over each
601  ! sub-cell, the sum of sub-cell thickness might differ from the original source
602  ! cell thickness due to round off.
603  dh0_eff = dh0_eff + min(dh, h0_supply)
604 
605  ! Record the source index (i0) that this sub-cell integral belongs to. This
606  ! is needed to index the reconstruction coefficients for the source cell
607  ! used in the integrals of the sub-cell width.
608  isub_src(i_sub) = i0
609  h_sub(i_sub) = dh
610 
611  ! For recording the largest sub-cell within a source cell.
612  if (dh >= dh_max) then
613  i_max = i_sub
614  dh_max = dh
615  endif
616 
617  ! Which ever column (source or target) has the least width left to consume determined
618  ! the width, dh, of sub-cell i_sub in the expression for dh above.
619  if (h0_supply <= h1_supply .and. src_has_volume) then
620  ! h0_supply is smaller than h1_supply) so we consume h0_supply and increment the
621  ! source cell index.
622  h1_supply = h1_supply - dh ! Although this is a difference the result will
623  ! be non-negative because of the conditional.
624  ! Record the sub-cell start/end index that span the source cell i0.
625  isrc_start(i0) = i_start0
626  isrc_end(i0) = i_sub
627  i_start0 = i_sub + 1
628  ! Record the sub-cell that is the largest fraction of the source cell.
629  isrc_max(i0) = i_max
630  i_max = i_sub + 1
631  dh_max = 0.
632  ! Record the source cell thickness found by summing the sub-cell thicknesses.
633  h0_eff(i0) = dh0_eff
634  ! Move the source index.
635  if (old_algorithm) then
636  if (i0 < i0_last_thick_cell) then
637  i0 = i0 + 1
638  h0_supply = h0(i0)
639  dh0_eff = 0.
640  do while (h0_supply==0. .and. i0<i0_last_thick_cell)
641  ! This loop skips over vanished source cells
642  i0 = i0 + 1
643  h0_supply = h0(i0)
644  enddo
645  else
646  h0_supply = 1.e30
647  endif
648  else
649  if (i0 < n0) then
650  i0 = i0 + 1
651  h0_supply = h0(i0)
652  dh0_eff = 0.
653  else
654  h0_supply = 0.
655  src_has_volume = .false.
656  endif
657  endif
658  elseif (h0_supply >= h1_supply .and. tgt_has_volume) then
659  ! h1_supply is smaller than h0_supply) so we consume h1_supply and increment the
660  ! target cell index.
661  h0_supply = h0_supply - dh ! Although this is a difference the result will
662  ! be non-negative because of the conditional.
663  ! Record the sub-cell start/end index that span the target cell i1.
664  itgt_start(i1) = i_start1
665  itgt_end(i1) = i_sub
666  i_start1 = i_sub + 1
667  ! Move the target index.
668  if (i1 < n1) then
669  i1 = i1 + 1
670  h1_supply = h1(i1)
671  else
672  if (old_algorithm) then
673  h1_supply = 1.e30
674  else
675  h1_supply = 0.
676  tgt_has_volume = .false.
677  endif
678  endif
679  elseif (src_has_volume) then
680  ! We ran out of target volume but still have source cells to consume
681  h_sub(i_sub) = h0_supply
682  ! Record the sub-cell start/end index that span the source cell i0.
683  isrc_start(i0) = i_start0
684  isrc_end(i0) = i_sub
685  i_start0 = i_sub + 1
686  ! Record the sub-cell that is the largest fraction of the source cell.
687  isrc_max(i0) = i_max
688  i_max = i_sub + 1
689  dh_max = 0.
690  ! Record the source cell thickness found by summing the sub-cell thicknesses.
691  h0_eff(i0) = dh0_eff
692  if (i0 < n0) then
693  i0 = i0 + 1
694  h0_supply = h0(i0)
695  dh0_eff = 0.
696  else
697  h0_supply = 0.
698  src_has_volume = .false.
699  endif
700  elseif (tgt_has_volume) then
701  ! We ran out of source volume but still have target cells to consume
702  h_sub(i_sub) = h1_supply
703  ! Record the sub-cell start/end index that span the target cell i1.
704  itgt_start(i1) = i_start1
705  itgt_end(i1) = i_sub
706  i_start1 = i_sub + 1
707  ! Move the target index.
708  if (i1 < n1) then
709  i1 = i1 + 1
710  h1_supply = h1(i1)
711  else
712  h1_supply = 0.
713  tgt_has_volume = .false.
714  endif
715  else
716  stop 'remap_via_sub_cells: THIS SHOULD NEVER HAPPEN!'
717  endif
718 
719  enddo
720 
721  ! Loop over each sub-cell to calculate average/integral values within each sub-cell.
722  xa = 0.
723  dh0_eff = 0.
724  uh_sub(1) = 0.
725  u_sub(1) = ppoly0_e(1,1)
726  u02_err = 0.
727  do i_sub = 2, n0+n1
728 
729  ! Sub-cell thickness from loop above
730  dh = h_sub(i_sub)
731 
732  ! Source cell
733  i0 = isub_src(i_sub)
734 
735  ! Evaluate average and integral for sub-cell i_sub.
736  ! Integral is over distance dh but expressed in terms of non-dimensional
737  ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).
738  dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
739  if (h0_eff(i0)>0.) then
740  xb = dh0_eff / h0_eff(i0) ! This expression yields xa <= xb <= 1.0
741  xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
742  u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_e, ppoly0_coefs, method, i0, xa, xb)
743  else ! Vanished cell
744  xb = 1.
745  u_sub(i_sub) = u0(i0)
746  endif
747  if (debug_bounds) then
748  if (method<5 .and.(u_sub(i_sub)<u0_min(i0) .or. u_sub(i_sub)>u0_max(i0))) then
749  write(0,*) 'Sub cell average is out of bounds',i_sub,'method=',method
750  write(0,*) 'xa,xb: ',xa,xb
751  write(0,*) 'Edge values: ',ppoly0_e(i0,:),'mean',u0(i0)
752  write(0,*) 'a_c: ',(u0(i0)-ppoly0_e(i0,1))+(u0(i0)-ppoly0_e(i0,2))
753  write(0,*) 'Polynomial coeffs: ',ppoly0_coefs(i0,:)
754  write(0,*) 'Bounds min=',u0_min(i0),'max=',u0_max(i0)
755  write(0,*) 'Average: ',u_sub(i_sub),'rel to min=',u_sub(i_sub)-u0_min(i0),'rel to max=',u_sub(i_sub)-u0_max(i0)
756  call mom_error( fatal, 'MOM_remapping, remap_via_sub_cells: '//&
757  'Sub-cell average is out of bounds!' )
758  endif
759  endif
760  if (force_bounds_in_subcell) then
761  ! These next two lines should not be needed but when using PQM we found roundoff
762  ! can lead to overshoots. These lines sweep issues under the rug which need to be
763  ! properly .. later. -AJA
764  u_orig = u_sub(i_sub)
765  u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )
766  u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )
767  u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig )
768  endif
769  uh_sub(i_sub) = dh * u_sub(i_sub)
770 
771  if (isub_src(i_sub+1) /= i0) then
772  ! If the next sub-cell is in a different source cell, reset the position counters
773  dh0_eff = 0.
774  xa = 0.
775  else
776  xa = xb ! Next integral will start at end of last
777  endif
778 
779  enddo
780  u_sub(n0+n1+1) = ppoly0_e(n0,2) ! This value is only needed when total target column
781  uh_sub(n0+n1+1) = ppoly0_e(n0,2) * h_sub(n0+n1+1) ! is wider than the source column
782 
783  if (adjust_thickest_subcell) then
784  ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within
785  ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals
786  ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (@Hallberg-NOAA).
787  do i0 = 1, i0_last_thick_cell
788  i_max = isrc_max(i0)
789  dh_max = h_sub(i_max)
790  if (dh_max > 0.) then
791  ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell.
792  duh = 0.
793  do i_sub = isrc_start(i0), isrc_end(i0)
794  if (i_sub /= i_max) duh = duh + uh_sub(i_sub)
795  enddo
796  uh_sub(i_max) = u0(i0)*h0(i0) - duh
797  u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) )
798  endif
799  enddo
800  endif
801 
802  ! Loop over each target cell summing the integrals from sub-cells within the target cell.
803  uh_err = 0.
804  do i1 = 1, n1
805  if (h1(i1) > 0.) then
806  duh = 0. ; dh = 0.
807  i_sub = itgt_start(i1)
808  if (force_bounds_in_target) then
809  u1min = u_sub(i_sub)
810  u1max = u_sub(i_sub)
811  endif
812  do i_sub = itgt_start(i1), itgt_end(i1)
813  if (force_bounds_in_target) then
814  u1min = min(u1min, u_sub(i_sub))
815  u1max = max(u1max, u_sub(i_sub))
816  endif
817  dh = dh + h_sub(i_sub)
818  duh = duh + uh_sub(i_sub)
819  ! This accumulates the contribution to the error bound for the sum of u*h
820  uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh)
821  enddo
822  u1(i1) = duh / dh
823  ! This is the contribution from the division to the error bound for the sum of u*h
824  uh_err = uh_err + abs(duh)*epsilon(duh)
825  if (force_bounds_in_target) then
826  u_orig = u1(i1)
827  u1(i1) = max(u1min, min(u1max, u1(i1)))
828  ! Adjusting to be bounded contributes to the error for the sum of u*h
829  uh_err = uh_err + dh*abs( u1(i1)-u_orig )
830  endif
831  else
832  u1(i1) = u_sub(itgt_start(i1))
833  endif
834  enddo
835 
836  ! Check errors and bounds
837  if (debug_bounds) then
838  call measure_input_bounds( n0, h0, u0, ppoly0_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
839  call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
840  call measure_output_bounds( n0+n1+1, h_sub, u_sub, h2tot, h2err, u2tot, u2err, u2min, u2max )
841  if (method<5) then ! We except PQM until we've debugged it
842  if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err+u02_err .and. abs(h1tot-h0tot)<h0err+h1err) &
843  .or. (abs(u2tot-u0tot)>u0err+u2err+u02_err .and. abs(h2tot-h0tot)<h0err+h2err) &
844  .or. (u1min<u0min .or. u1max>u0max) ) then
845  write(0,*) 'method = ',method
846  write(0,*) 'Source to sub-cells:'
847  write(0,*) 'H: h0tot=',h0tot,'h2tot=',h2tot,'dh=',h2tot-h0tot,'h0err=',h0err,'h2err=',h2err
848  if (abs(h2tot-h0tot)>h0err+h2err) &
849  write(0,*) 'H non-conservation difference=',h2tot-h0tot,'allowed err=',h0err+h2err,' <-----!'
850  write(0,*) 'UH: u0tot=',u0tot,'u2tot=',u2tot,'duh=',u2tot-u0tot,'u0err=',u0err,'u2err=',u2err,&
851  'adjustment err=',u02_err
852  if (abs(u2tot-u0tot)>u0err+u2err) &
853  write(0,*) 'U non-conservation difference=',u2tot-u0tot,'allowed err=',u0err+u2err,' <-----!'
854  write(0,*) 'Sub-cells to target:'
855  write(0,*) 'H: h2tot=',h2tot,'h1tot=',h1tot,'dh=',h1tot-h2tot,'h2err=',h2err,'h1err=',h1err
856  if (abs(h1tot-h2tot)>h2err+h1err) &
857  write(0,*) 'H non-conservation difference=',h1tot-h2tot,'allowed err=',h2err+h1err,' <-----!'
858  write(0,*) 'UH: u2tot=',u2tot,'u1tot=',u1tot,'duh=',u1tot-u2tot,'u2err=',u2err,'u1err=',u1err,'uh_err=',uh_err
859  if (abs(u1tot-u2tot)>u2err+u1err) &
860  write(0,*) 'U non-conservation difference=',u1tot-u2tot,'allowed err=',u2err+u1err,' <-----!'
861  write(0,*) 'Source to target:'
862  write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err
863  if (abs(h1tot-h0tot)>h0err+h1err) &
864  write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!'
865  write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err
866  if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
867  write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!'
868  write(0,*) 'U: u0min=',u0min,'u1min=',u1min,'u2min=',u2min
869  if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!'
870  if (u2min<u0min) write(0,*) 'U2 minimum overshoot=',u2min-u0min,' <-----!'
871  write(0,*) 'U: u0max=',u0max,'u1max=',u1max,'u2max=',u2max
872  if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!'
873  if (u2max>u0max) write(0,*) 'U2 maximum overshoot=',u2max-u0max,' <-----!'
874  write(0,'(a3,6a24,2a3)') 'k','h0','left edge','u0','right edge','h1','u1','is','ie'
875  do k = 1, max(n0,n1)
876  if (k<=min(n0,n1)) then
877  write(0,'(i3,1p6e24.16,2i3)') k,h0(k),ppoly0_e(k,1),u0(k),ppoly0_e(k,2),h1(k),u1(k),itgt_start(k),itgt_end(k)
878  elseif (k>n0) then
879  write(0,'(i3,96x,1p2e24.16,2i3)') k,h1(k),u1(k),itgt_start(k),itgt_end(k)
880  else
881  write(0,'(i3,1p4e24.16)') k,h0(k),ppoly0_e(k,1),u0(k),ppoly0_e(k,2)
882  endif
883  enddo
884  write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients'
885  do k = 1, n0
886  write(0,'(i3,1p6e24.16)') k,u0(k),ppoly0_coefs(k,:)
887  enddo
888  write(0,'(a3,3a24,a3,2a24)') 'k','Sub-cell h','Sub-cell u','Sub-cell hu','i0','xa','xb'
889  xa = 0.
890  dh0_eff = 0.
891  do k = 1, n0+n1+1
892  dh = h_sub(k)
893  i0 = isub_src(k)
894  dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
895  xb = dh0_eff / h0_eff(i0) ! This expression yields xa <= xb <= 1.0
896  xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
897  write(0,'(i3,1p3e24.16,i3,1p2e24.16)') k,h_sub(k),u_sub(k),uh_sub(k),i0,xa,xb
898  if (k<=n0+n1) then
899  if (isub_src(k+1) /= i0) then
900  dh0_eff = 0.; xa = 0.
901  else
902  xa = xb
903  endif
904  endif
905  enddo
906  call mom_error( fatal, 'MOM_remapping, remap_via_sub_cells: '//&
907  'Remapping result is inconsistent!' )
908  endif
909  endif ! method<5
910  endif ! debug_bounds
911 
912  ! Include the error remapping from source to sub-cells in the estimate of total remapping error
913  uh_err = uh_err + u02_err
914 
915  if (present(ah_sub)) ah_sub(1:n0+n1+1) = h_sub(1:n0+n1+1)
916  if (present(aisub_src)) aisub_src(1:n0+n1+1) = isub_src(1:n0+n1+1)
917  if (present(aiss)) aiss(1:n0) = isrc_start(1:n0)
918  if (present(aise)) aise(1:n0) = isrc_end(1:n0)
919 
920 end subroutine remap_via_sub_cells
921 
922 !> Returns the average value of a reconstruction within a single source cell, i0,
923 !! between the non-dimensional positions xa and xb (xa<=xb) with dimensional
924 !! separation dh.
925 real function average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb)
926  integer, intent(in) :: n0 !< Number of cells in source grid
927  real, intent(in) :: u0(:) !< Cell means
928  real, intent(in) :: ppoly0_e(:,:) !< Edge value of polynomial
929  real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial
930  integer, intent(in) :: method !< Remapping scheme to use
931  integer, intent(in) :: i0 !< Source cell index
932  real, intent(in) :: xa !< Non-dimensional start position within source cell
933  real, intent(in) :: xb !< Non-dimensional end position within source cell
934  ! Local variables
935  real :: u_ave, xa_2, xb_2, xa2pxb2, xapxb
936  real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials
937 
938  real :: mx, a_l, a_r, u_c, ya, yb, my, xa2b2ab, ya2b2ab, a_c
939 
940  if (xb > xa) then
941  select case ( method )
942  case ( integration_pcm )
943  u_ave = u0(i0)
944  case ( integration_plm )
945  u_ave = ( &
946  ppoly0_coefs(i0,1) &
947  + ppoly0_coefs(i0,2) * 0.5 * ( xb + xa ) )
948  case ( integration_ppm )
949  mx = 0.5 * ( xa + xb )
950  a_l = ppoly0_e(i0, 1)
951  a_r = ppoly0_e(i0, 2)
952  u_c = u0(i0)
953  a_c = 0.5 * ( ( u_c - a_l ) + ( u_c - a_r ) ) ! a_6 / 6
954  if (mx<0.5) then
955  ! This integration of the PPM reconstruction is expressed in distances from the left edge
956  xa2b2ab = (xa*xa+xb*xb)+xa*xb
957  u_ave = a_l + ( ( a_r - a_l ) * mx &
958  + a_c * ( 3. * ( xb + xa ) - 2.*xa2b2ab ) )
959  else
960  ! This integration of the PPM reconstruction is expressed in distances from the right edge
961  ya = 1. - xa
962  yb = 1. - xb
963  my = 0.5 * ( ya + yb )
964  ya2b2ab = (ya*ya+yb*yb)+ya*yb
965  u_ave = a_r + ( ( a_l - a_r ) * my &
966  + a_c * ( 3. * ( yb + ya ) - 2.*ya2b2ab ) )
967  endif
968  case ( integration_pqm )
969  xa_2 = xa*xa
970  xb_2 = xb*xb
971  xa2pxb2 = xa_2 + xb_2
972  xapxb = xa + xb
973  u_ave = ( &
974  ppoly0_coefs(i0,1) &
975  + ( ppoly0_coefs(i0,2) * 0.5 * ( xapxb ) &
976  + ( ppoly0_coefs(i0,3) * r_3 * ( xa2pxb2 + xa*xb ) &
977  + ( ppoly0_coefs(i0,4) * 0.25* ( xa2pxb2 * xapxb ) &
978  + ppoly0_coefs(i0,5) * 0.2 * ( ( xb*xb_2 + xa*xa_2 ) * xapxb + xa_2*xb_2 ) ) ) ) )
979  case default
980  call mom_error( fatal,'The selected integration method is invalid' )
981  end select
982  else ! dh == 0.
983  select case ( method )
984  case ( integration_pcm )
985  u_ave = ppoly0_coefs(i0,1)
986  case ( integration_plm )
987  !u_ave = ppoly0_coefs(i0,1) &
988  ! + xa * ppoly0_coefs(i0,2)
989  a_l = ppoly0_e(i0, 1)
990  a_r = ppoly0_e(i0, 2)
991  ya = 1. - xa
992  if (xa < 0.5) then
993  u_ave = a_l + xa * ( a_r - a_l )
994  else
995  u_ave = a_r + ya * ( a_l - a_r )
996  endif
997  case ( integration_ppm )
998  !u_ave = ppoly0_coefs(i0,1) &
999  ! + xa * ( ppoly0_coefs(i0,2) &
1000  ! + xa * ppoly0_coefs(i0,3) )
1001  a_l = ppoly0_e(i0, 1)
1002  a_r = ppoly0_e(i0, 2)
1003  u_c = u0(i0)
1004  a_c = 3. * ( ( u_c - a_l ) + ( u_c - a_r ) ) ! a_6
1005  ya = 1. - xa
1006  if (xa < 0.5) then
1007  u_ave = a_l + xa * ( ( a_r - a_l ) + a_c * ya )
1008  else
1009  u_ave = a_r + ya * ( ( a_l - a_r ) + a_c * xa )
1010  endif
1011  case ( integration_pqm )
1012  u_ave = ppoly0_coefs(i0,1) &
1013  + xa * ( ppoly0_coefs(i0,2) &
1014  + xa * ( ppoly0_coefs(i0,3) &
1015  + xa * ( ppoly0_coefs(i0,4) &
1016  + xa * ppoly0_coefs(i0,5) ) ) )
1017  case default
1018  call mom_error( fatal,'The selected integration method is invalid' )
1019  end select
1020  endif
1021  average_value_ppoly = u_ave
1022 
1023 end function average_value_ppoly
1024 
1025 !> Measure totals and bounds on source grid
1026 subroutine measure_input_bounds( n0, h0, u0, ppoly_E, h0tot, h0err, u0tot, u0err, u0min, u0max )
1027  integer, intent(in) :: n0 !< Number of cells on source grid
1028  real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid
1029  real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid
1030  real, dimension(n0,2), intent(in) :: ppoly_E !< Cell edge values on source grid
1031  real, intent(out) :: h0tot !< Sum of cell widths
1032  real, intent(out) :: h0err !< Magnitude of round-off error in h0tot
1033  real, intent(out) :: u0tot !< Sum of cell widths times values
1034  real, intent(out) :: u0err !< Magnitude of round-off error in u0tot
1035  real, intent(out) :: u0min !< Minimum value in reconstructions of u0
1036  real, intent(out) :: u0max !< Maximum value in reconstructions of u0
1037  ! Local variables
1038  integer :: k
1039  real :: eps
1040 
1041  eps = epsilon(h0(1))
1042  h0tot = h0(1)
1043  h0err = 0.
1044  u0tot = h0(1) * u0(1)
1045  u0err = 0.
1046  u0min = min( ppoly_e(1,1), ppoly_e(1,2) )
1047  u0max = max( ppoly_e(1,1), ppoly_e(1,2) )
1048  do k = 2, n0
1049  h0tot = h0tot + h0(k)
1050  h0err = h0err + eps * max(h0tot, h0(k))
1051  u0tot = u0tot + h0(k) * u0(k)
1052  u0err = u0err + eps * max(abs(u0tot), abs(h0(k) * u0(k)))
1053  u0min = min( u0min, ppoly_e(k,1), ppoly_e(k,2) )
1054  u0max = max( u0max, ppoly_e(k,1), ppoly_e(k,2) )
1055  enddo
1056 
1057 end subroutine measure_input_bounds
1058 
1059 !> Measure totals and bounds on destination grid
1060 subroutine measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
1061  integer, intent(in) :: n1 !< Number of cells on destination grid
1062  real, dimension(n1), intent(in) :: h1 !< Cell widths on destination grid
1063  real, dimension(n1), intent(in) :: u1 !< Cell averages on destination grid
1064  real, intent(out) :: h1tot !< Sum of cell widths
1065  real, intent(out) :: h1err !< Magnitude of round-off error in h1tot
1066  real, intent(out) :: u1tot !< Sum of cell widths times values
1067  real, intent(out) :: u1err !< Magnitude of round-off error in u1tot
1068  real, intent(out) :: u1min !< Minimum value in reconstructions of u1
1069  real, intent(out) :: u1max !< Maximum value in reconstructions of u1
1070  ! Local variables
1071  integer :: k
1072  real :: eps
1073 
1074  eps = epsilon(h1(1))
1075  h1tot = h1(1)
1076  h1err = 0.
1077  u1tot = h1(1) * u1(1)
1078  u1err = 0.
1079  u1min = u1(1)
1080  u1max = u1(1)
1081  do k = 2, n1
1082  h1tot = h1tot + h1(k)
1083  h1err = h1err + eps * max(h1tot, h1(k))
1084  u1tot = u1tot + h1(k) * u1(k)
1085  u1err = u1err + eps * max(abs(u1tot), abs(h1(k) * u1(k)))
1086  u1min = min(u1min, u1(k))
1087  u1max = max(u1max, u1(k))
1088  enddo
1089 
1090 end subroutine measure_output_bounds
1091 
1092 !> Remaps column of values u0 on grid h0 to grid h1 by integrating
1093 !! over the projection of each h1 cell onto the h0 grid.
1094 subroutine remapbyprojection( n0, h0, u0, ppoly0_E, ppoly0_coefs, &
1095  n1, h1, method, u1, h_neglect )
1096  integer, intent(in) :: n0 !< Number of cells in source grid
1097  real, intent(in) :: h0(:) !< Source grid widths (size n0)
1098  real, intent(in) :: u0(:) !< Source cell averages (size n0)
1099  real, intent(in) :: ppoly0_E(:,:) !< Edge value of polynomial
1100  real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial
1101  integer, intent(in) :: n1 !< Number of cells in target grid
1102  real, intent(in) :: h1(:) !< Target grid widths (size n1)
1103  integer, intent(in) :: method !< Remapping scheme to use
1104  real, intent(out) :: u1(:) !< Target cell averages (size n1)
1105  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
1106  !! purpose of cell reconstructions
1107  !! in the same units as h.
1108  ! Local variables
1109  integer :: iTarget
1110  real :: xL, xR ! coordinates of target cell edges
1111  integer :: jStart ! Used by integrateReconOnInterval()
1112  real :: xStart ! Used by integrateReconOnInterval()
1113 
1114  ! Loop on cells in target grid (grid1). For each target cell, we need to find
1115  ! in which source cells the target cell edges lie. The associated indexes are
1116  ! noted j0 and j1.
1117  xr = 0. ! Left boundary is at x=0
1118  jstart = 1
1119  xstart = 0.
1120  do itarget = 1,n1
1121  ! Determine the coordinates of the target cell edges
1122  xl = xr
1123  xr = xl + h1(itarget)
1124 
1125  call integraterecononinterval( n0, h0, u0, ppoly0_e, ppoly0_coefs, method, &
1126  xl, xr, h1(itarget), u1(itarget), jstart, xstart, h_neglect )
1127 
1128  enddo ! end iTarget loop on target grid cells
1129 
1130 end subroutine remapbyprojection
1131 
1132 
1133 !> Remaps column of values u0 on grid h0 to implied grid h1
1134 !! where the interfaces of h1 differ from those of h0 by dx.
1135 !! The new grid is defined relative to the original grid by change
1136 !! dx1(:) = xNew(:) - xOld(:)
1137 !! and the remapping calculated so that
1138 !! hNew(k) qNew(k) = hOld(k) qOld(k) + F(k+1) - F(k)
1139 !! where
1140 !! F(k) = dx1(k) qAverage
1141 !! and where qAverage is the average qOld in the region zOld(k) to zNew(k).
1142 subroutine remapbydeltaz( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, &
1143  method, u1, h1, h_neglect )
1144  integer, intent(in) :: n0 !< Number of cells in source grid
1145  real, dimension(:), intent(in) :: h0 !< Source grid sizes (size n0)
1146  real, dimension(:), intent(in) :: u0 !< Source cell averages (size n0)
1147  real, dimension(:,:), intent(in) :: ppoly0_E !< Edge value of polynomial
1148  real, dimension(:,:), intent(in) :: ppoly0_coefs !< Coefficients of polynomial
1149  integer, intent(in) :: n1 !< Number of cells in target grid
1150  real, dimension(:), intent(in) :: dx1 !< Target grid edge positions (size n1+1)
1151  integer, intent(in) :: method !< Remapping scheme to use
1152  real, dimension(:), intent(out) :: u1 !< Target cell averages (size n1)
1153  real, dimension(:), &
1154  optional, intent(out) :: h1 !< Target grid widths (size n1)
1155  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
1156  !! purpose of cell reconstructions
1157  !! in the same units as h.
1158  ! Local variables
1159  integer :: iTarget
1160  real :: xL, xR ! coordinates of target cell edges
1161  real :: xOld, hOld, uOld
1162  real :: xNew, hNew, h_err
1163  real :: uhNew, hFlux, uAve, fluxL, fluxR
1164  integer :: jStart ! Used by integrateReconOnInterval()
1165  real :: xStart ! Used by integrateReconOnInterval()
1166 
1167  ! Loop on cells in target grid. For each cell, iTarget, the left flux is
1168  ! the right flux of the cell to the left, iTarget-1.
1169  ! The left flux is initialized by started at iTarget=0 to calculate the
1170  ! right flux which can take into account the target left boundary being
1171  ! in the interior of the source domain.
1172  fluxr = 0.
1173  h_err = 0. ! For measuring round-off error
1174  jstart = 1
1175  xstart = 0.
1176  do itarget = 0,n1
1177  fluxl = fluxr ! This does nothing for iTarget=0
1178 
1179  if (itarget == 0) then
1180  xold = 0. ! Left boundary is at x=0
1181  hold = -1.e30 ! Should not be used for iTarget = 0
1182  uold = -1.e30 ! Should not be used for iTarget = 0
1183  elseif (itarget <= n0) then
1184  xold = xold + h0(itarget) ! Position of right edge of cell
1185  hold = h0(itarget)
1186  uold = u0(itarget)
1187  h_err = h_err + epsilon(hold) * max(hold, xold)
1188  else
1189  hold = 0. ! as if for layers>n0, they were vanished
1190  uold = 1.e30 ! and the initial value should not matter
1191  endif
1192  xnew = xold + dx1(itarget+1)
1193  xl = min( xold, xnew )
1194  xr = max( xold, xnew )
1195 
1196  ! hFlux is the positive width of the remapped volume
1197  hflux = abs(dx1(itarget+1))
1198  call integraterecononinterval( n0, h0, u0, ppoly0_e, ppoly0_coefs, method, &
1199  xl, xr, hflux, uave, jstart, xstart )
1200  ! uAve is the average value of u, independent of sign of dx1
1201  fluxr = dx1(itarget+1)*uave ! Includes sign of dx1
1202 
1203  if (itarget>0) then
1204  hnew = hold + ( dx1(itarget+1) - dx1(itarget) )
1205  hnew = max( 0., hnew )
1206  uhnew = ( uold * hold ) + ( fluxr - fluxl )
1207  if (hnew>0.) then
1208  u1(itarget) = uhnew / hnew
1209  else
1210  u1(itarget) = uave
1211  endif
1212  if (present(h1)) h1(itarget) = hnew
1213  endif
1214 
1215  enddo ! end iTarget loop on target grid cells
1216 
1217 end subroutine remapbydeltaz
1218 
1219 
1220 !> Integrate the reconstructed column profile over a single cell
1221 subroutine integraterecononinterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, &
1222  xL, xR, hC, uAve, jStart, xStart, h_neglect )
1223  integer, intent(in) :: n0 !< Number of cells in source grid
1224  real, dimension(:), intent(in) :: h0 !< Source grid sizes (size n0)
1225  real, dimension(:), intent(in) :: u0 !< Source cell averages
1226  real, dimension(:,:), intent(in) :: ppoly0_E !< Edge value of polynomial
1227  real, dimension(:,:), intent(in) :: ppoly0_coefs !< Coefficients of polynomial
1228  integer, intent(in) :: method !< Remapping scheme to use
1229  real, intent(in) :: xL !< Left edges of target cell
1230  real, intent(in) :: xR !< Right edges of target cell
1231  real, intent(in) :: hC !< Cell width hC = xR - xL
1232  real, intent(out) :: uAve !< Average value on target cell
1233  integer, intent(inout) :: jStart !< The index of the cell to start searching from
1234  !< On exit, contains index of last cell used
1235  real, intent(inout) :: xStart !< The left edge position of cell jStart
1236  !< On first entry should be 0.
1237  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
1238  !! purpose of cell reconstructions
1239  !! in the same units as h.
1240  ! Local variables
1241  integer :: j, k
1242  integer :: jL, jR ! indexes of source cells containing target
1243  ! cell edges
1244  real :: q ! complete integration
1245  real :: xi0, xi1 ! interval of integration (local -- normalized
1246  ! -- coordinates)
1247  real :: x0jLl, x0jLr ! Left/right position of cell jL
1248  real :: x0jRl, x0jRr ! Left/right position of cell jR
1249  real :: hAct ! The distance actually used in the integration
1250  ! (notionally xR - xL) which differs due to roundoff.
1251  real :: x0_2, x1_2, x02px12, x0px1 ! Used in evaluation of integrated polynomials
1252  real :: hNeglect ! A negligible thicness in the same units as h.
1253  real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials
1254 
1255  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
1256 
1257  q = -1.e30
1258  x0jll = -1.e30
1259  x0jrl = -1.e30
1260 
1261  ! Find the left most cell in source grid spanned by the target cell
1262  jl = -1
1263  x0jlr = xstart
1264  do j = jstart, n0
1265  x0jll = x0jlr
1266  x0jlr = x0jll + h0(j)
1267  ! Left edge is found in cell j
1268  if ( ( xl >= x0jll ) .AND. ( xl <= x0jlr ) ) then
1269  jl = j
1270  exit ! once target grid cell is found, exit loop
1271  endif
1272  enddo
1273  jstart = jl
1274  xstart = x0jll
1275 
1276 ! ! HACK to handle round-off problems. Need only at j=n0.
1277 ! ! This moves the effective cell boundary outwards a smidgen.
1278 ! if (xL>x0jLr) x0jLr = xL
1279 
1280  ! If, at this point, jL is equal to -1, it means the vanished
1281  ! cell lies outside the source grid. In other words, it means that
1282  ! the source and target grids do not cover the same physical domain
1283  ! and there is something very wrong !
1284  if ( jl == -1 ) call mom_error(fatal, &
1285  'MOM_remapping, integrateReconOnInterval: '//&
1286  'The location of the left-most cell could not be found')
1287 
1288 
1289  ! ============================================================
1290  ! Check whether target cell is vanished. If it is, the cell
1291  ! average is simply the interpolated value at the location
1292  ! of the vanished cell. If it isn't, we need to integrate the
1293  ! quantity within the cell and divide by the cell width to
1294  ! determine the cell average.
1295  ! ============================================================
1296  ! 1. Cell is vanished
1297  !if ( abs(xR - xL) <= epsilon(xR)*max(abs(xR),abs(xL)) ) then
1298  if ( abs(xr - xl) == 0.0 ) then
1299 
1300  ! We check whether the source cell (i.e. the cell in which the
1301  ! vanished target cell lies) is vanished. If it is, the interpolated
1302  ! value is set to be mean of the edge values (which should be the same).
1303  ! If it isn't, we simply interpolate.
1304  if ( h0(jl) == 0.0 ) then
1305  uave = 0.5 * ( ppoly0_e(jl,1) + ppoly0_e(jl,2) )
1306  else
1307  ! WHY IS THIS NOT WRITTEN AS xi0 = ( xL - x0jLl ) / h0(jL) ---AJA
1308  xi0 = xl / ( h0(jl) + hneglect ) - x0jll / ( h0(jl) + hneglect )
1309 
1310  select case ( method )
1311  case ( integration_pcm )
1312  uave = ppoly0_coefs(jl,1)
1313  case ( integration_plm )
1314  uave = ppoly0_coefs(jl,1) &
1315  + xi0 * ppoly0_coefs(jl,2)
1316  case ( integration_ppm )
1317  uave = ppoly0_coefs(jl,1) &
1318  + xi0 * ( ppoly0_coefs(jl,2) &
1319  + xi0 * ppoly0_coefs(jl,3) )
1320  case ( integration_pqm )
1321  uave = ppoly0_coefs(jl,1) &
1322  + xi0 * ( ppoly0_coefs(jl,2) &
1323  + xi0 * ( ppoly0_coefs(jl,3) &
1324  + xi0 * ( ppoly0_coefs(jl,4) &
1325  + xi0 * ppoly0_coefs(jl,5) ) ) )
1326  case default
1327  call mom_error( fatal,'The selected integration method is invalid' )
1328  end select
1329 
1330  endif ! end checking whether source cell is vanished
1331 
1332  ! 2. Cell is not vanished
1333  else
1334 
1335  ! Find the right most cell in source grid spanned by the target cell
1336  jr = -1
1337  x0jrr = xstart
1338  do j = jstart,n0
1339  x0jrl = x0jrr
1340  x0jrr = x0jrl + h0(j)
1341  ! Right edge is found in cell j
1342  if ( ( xr >= x0jrl ) .AND. ( xr <= x0jrr ) ) then
1343  jr = j
1344  exit ! once target grid cell is found, exit loop
1345  endif
1346  enddo ! end loop on source grid cells
1347 
1348  ! If xR>x0jRr then the previous loop reached j=n0 and the target
1349  ! position, xR, was beyond the right edge of the source grid (h0).
1350  ! This can happen due to roundoff, in which case we set jR=n0.
1351  if (xr>x0jrr) jr = n0
1352 
1353  ! To integrate, two cases must be considered: (1) the target cell is
1354  ! entirely contained within a cell of the source grid and (2) the target
1355  ! cell spans at least two cells of the source grid.
1356 
1357  if ( jl == jr ) then
1358  ! The target cell is entirely contained within a cell of the source
1359  ! grid. This situation is represented by the following schematic, where
1360  ! the cell in which xL and xR are located has index jL=jR :
1361  !
1362  ! ----|-----o--------o----------|-------------
1363  ! xL xR
1364  !
1365  ! Determine normalized coordinates
1366 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1367  xi0 = max( 0., min( 1., ( xl - x0jll ) / ( h0(jl) + hneglect ) ) )
1368  xi1 = max( 0., min( 1., ( xr - x0jll ) / ( h0(jl) + hneglect ) ) )
1369 #else
1370  xi0 = xl / h0(jl) - x0jll / ( h0(jl) + hneglect )
1371  xi1 = xr / h0(jl) - x0jll / ( h0(jl) + hneglect )
1372 #endif
1373 
1374  hact = h0(jl) * ( xi1 - xi0 )
1375 
1376  ! Depending on which polynomial is used, integrate quantity
1377  ! between xi0 and xi1. Integration is carried out in normalized
1378  ! coordinates, hence: \int_xL^xR p(x) dx = h \int_xi0^xi1 p(xi) dxi
1379  select case ( method )
1380  case ( integration_pcm )
1381  q = ( xr - xl ) * ppoly0_coefs(jl,1)
1382  case ( integration_plm )
1383  q = ( xr - xl ) * ( &
1384  ppoly0_coefs(jl,1) &
1385  + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
1386  case ( integration_ppm )
1387  q = ( xr - xl ) * ( &
1388  ppoly0_coefs(jl,1) &
1389  + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1390  + ppoly0_coefs(jl,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
1391  case ( integration_pqm )
1392  x0_2 = xi0*xi0
1393  x1_2 = xi1*xi1
1394  x02px12 = x0_2 + x1_2
1395  x0px1 = xi1 + xi0
1396  q = ( xr - xl ) * ( &
1397  ppoly0_coefs(jl,1) &
1398  + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1399  + ( ppoly0_coefs(jl,3) * r_3 * ( x02px12 + xi0*xi1 ) &
1400  + ppoly0_coefs(jl,4) * 0.25* ( x02px12 * x0px1 ) &
1401  + ppoly0_coefs(jl,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
1402  case default
1403  call mom_error( fatal,'The selected integration method is invalid' )
1404  end select
1405 
1406  else
1407  ! The target cell spans at least two cells of the source grid.
1408  ! This situation is represented by the following schematic, where
1409  ! the cells in which xL and xR are located have indexes jL and jR,
1410  ! respectively :
1411  !
1412  ! ----|-----o---|--- ... --|---o----------|-------------
1413  ! xL xR
1414  !
1415  ! We first integrate from xL up to the right boundary of cell jL, then
1416  ! add the integrated amounts of cells located between jL and jR and then
1417  ! integrate from the left boundary of cell jR up to xR
1418 
1419  q = 0.0
1420 
1421  ! Integrate from xL up to right boundary of cell jL
1422 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1423  xi0 = max( 0., min( 1., ( xl - x0jll ) / ( h0(jl) + hneglect ) ) )
1424 #else
1425  xi0 = (xl - x0jll) / ( h0(jl) + hneglect )
1426 #endif
1427  xi1 = 1.0
1428 
1429  hact = h0(jl) * ( xi1 - xi0 )
1430 
1431  select case ( method )
1432  case ( integration_pcm )
1433  q = q + ( x0jlr - xl ) * ppoly0_coefs(jl,1)
1434  case ( integration_plm )
1435  q = q + ( x0jlr - xl ) * ( &
1436  ppoly0_coefs(jl,1) &
1437  + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
1438  case ( integration_ppm )
1439  q = q + ( x0jlr - xl ) * ( &
1440  ppoly0_coefs(jl,1) &
1441  + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1442  + ppoly0_coefs(jl,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
1443  case ( integration_pqm )
1444  x0_2 = xi0*xi0
1445  x1_2 = xi1*xi1
1446  x02px12 = x0_2 + x1_2
1447  x0px1 = xi1 + xi0
1448  q = q + ( x0jlr - xl ) * ( &
1449  ppoly0_coefs(jl,1) &
1450  + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1451  + ( ppoly0_coefs(jl,3) * r_3 * ( x02px12 + xi0*xi1 ) &
1452  + ppoly0_coefs(jl,4) * 0.25* ( x02px12 * x0px1 ) &
1453  + ppoly0_coefs(jl,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
1454  case default
1455  call mom_error( fatal, 'The selected integration method is invalid' )
1456  end select
1457 
1458  ! Integrate contents within cells strictly comprised between jL and jR
1459  if ( jr > (jl+1) ) then
1460  do k = jl+1,jr-1
1461  q = q + h0(k) * u0(k)
1462  hact = hact + h0(k)
1463  enddo
1464  endif
1465 
1466  ! Integrate from left boundary of cell jR up to xR
1467  xi0 = 0.0
1468 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1469  xi1 = max( 0., min( 1., ( xr - x0jrl ) / ( h0(jr) + hneglect ) ) )
1470 #else
1471  xi1 = (xr - x0jrl) / ( h0(jr) + hneglect )
1472 #endif
1473 
1474  hact = hact + h0(jr) * ( xi1 - xi0 )
1475 
1476  select case ( method )
1477  case ( integration_pcm )
1478  q = q + ( xr - x0jrl ) * ppoly0_coefs(jr,1)
1479  case ( integration_plm )
1480  q = q + ( xr - x0jrl ) * ( &
1481  ppoly0_coefs(jr,1) &
1482  + ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) )
1483  case ( integration_ppm )
1484  q = q + ( xr - x0jrl ) * ( &
1485  ppoly0_coefs(jr,1) &
1486  + ( ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) &
1487  + ppoly0_coefs(jr,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
1488  case ( integration_pqm )
1489  x0_2 = xi0*xi0
1490  x1_2 = xi1*xi1
1491  x02px12 = x0_2 + x1_2
1492  x0px1 = xi1 + xi0
1493  q = q + ( xr - x0jrl ) * ( &
1494  ppoly0_coefs(jr,1) &
1495  + ( ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) &
1496  + ( ppoly0_coefs(jr,3) * r_3 * ( x02px12 + xi0*xi1 ) &
1497  + ppoly0_coefs(jr,4) * 0.25* ( x02px12 * x0px1 ) &
1498  + ppoly0_coefs(jr,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
1499  case default
1500  call mom_error( fatal,'The selected integration method is invalid' )
1501  end select
1502 
1503  endif ! end integration for non-vanished cells
1504 
1505  ! The cell average is the integrated value divided by the cell width
1506 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1507 if (hact==0.) then
1508  uave = ppoly0_coefs(jl,1)
1509 else
1510  uave = q / hact
1511 endif
1512 #else
1513  uave = q / hc
1514 #endif
1515 
1516  endif ! endif clause to check if cell is vanished
1517 
1518 end subroutine integraterecononinterval
1519 
1520 !> Calculates the change in interface positions based on h1 and h2
1521 subroutine dzfromh1h2( n1, h1, n2, h2, dx )
1522  integer, intent(in) :: n1 !< Number of cells on source grid
1523  real, dimension(:), intent(in) :: h1 !< Cell widths of source grid (size n1)
1524  integer, intent(in) :: n2 !< Number of cells on target grid
1525  real, dimension(:), intent(in) :: h2 !< Cell widths of target grid (size n2)
1526  real, dimension(:), intent(out) :: dx !< Change in interface position (size n2+1)
1527  ! Local variables
1528  integer :: k
1529  real :: x1, x2
1530 
1531  x1 = 0.
1532  x2 = 0.
1533  dx(1) = 0.
1534  do k = 1, max(n1,n2)
1535  if (k <= n1) x1 = x1 + h1(k) ! Interface k+1, right of source cell k
1536  if (k <= n2) then
1537  x2 = x2 + h2(k) ! Interface k+1, right of target cell k
1538  dx(k+1) = x2 - x1 ! Change of interface k+1, target - source
1539  endif
1540  enddo
1541 
1542 end subroutine dzfromh1h2
1543 
1544 !> Constructor for remapping control structure
1545 subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, &
1546  check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
1547  ! Arguments
1548  type(remapping_cs), intent(inout) :: cs !< Remapping control structure
1549  character(len=*), intent(in) :: remapping_scheme !< Remapping scheme to use
1550  logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
1551  logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions
1552  logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping
1553  logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded
1554  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
1555 
1556  ! Note that remapping_scheme is mandatory for initialize_remapping()
1557  call remapping_set_param(cs, remapping_scheme=remapping_scheme, boundary_extrapolation=boundary_extrapolation, &
1558  check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
1559  force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018)
1560 
1561 end subroutine initialize_remapping
1562 
1563 !> Changes the method of reconstruction
1564 !! Use this routine to parse a string parameter specifying the reconstruction
1565 !! and re-allocates work arrays appropriately. It is called from
1566 !! initialize_remapping but can be called from an external module too.
1567 subroutine setreconstructiontype(string,CS)
1568  character(len=*), intent(in) :: string !< String to parse for method
1569  type(remapping_cs), intent(inout) :: CS !< Remapping control structure
1570  ! Local variables
1571  integer :: degree
1572  degree = -99
1573  select case ( uppercase(trim(string)) )
1574  case ("PCM")
1575  cs%remapping_scheme = remapping_pcm
1576  degree = 0
1577  case ("PLM")
1578  cs%remapping_scheme = remapping_plm
1579  degree = 1
1580  case ("PPM_H4")
1581  cs%remapping_scheme = remapping_ppm_h4
1582  degree = 2
1583  case ("PPM_IH4")
1584  cs%remapping_scheme = remapping_ppm_ih4
1585  degree = 2
1586  case ("PQM_IH4IH3")
1587  cs%remapping_scheme = remapping_pqm_ih4ih3
1588  degree = 4
1589  case ("PQM_IH6IH5")
1590  cs%remapping_scheme = remapping_pqm_ih6ih5
1591  degree = 4
1592  case default
1593  call mom_error(fatal, "setReconstructionType: "//&
1594  "Unrecognized choice for REMAPPING_SCHEME ("//trim(string)//").")
1595  end select
1596 
1597  cs%degree = degree
1598 
1599 end subroutine setreconstructiontype
1600 
1601 !> Destrcutor for remapping control structure
1602 subroutine end_remapping(CS)
1603  type(remapping_cs), intent(inout) :: cs !< Remapping control structure
1604 
1605  cs%degree = 0
1606 
1607 end subroutine end_remapping
1608 
1609 !> Runs unit tests on remapping functions.
1610 !! Should only be called from a single/root thread
1611 !! Returns True if a test fails, otherwise False
1612 logical function remapping_unit_tests(verbose)
1613  logical, intent(in) :: verbose !< If true, write results to stdout
1614  ! Local variables
1615  integer, parameter :: n0 = 4, n1 = 3, n2 = 6
1616  real :: h0(n0), x0(n0+1), u0(n0)
1617  real :: h1(n1), x1(n1+1), u1(n1), hn1(n1), dx1(n1+1)
1618  real :: h2(n2), x2(n2+1), u2(n2), hn2(n2), dx2(n2+1)
1619  data u0 /9., 3., -3., -9./ ! Linear profile, 4 at surface to -4 at bottom
1620  data h0 /4*0.75/ ! 4 uniform layers with total depth of 3
1621  data h1 /3*1./ ! 3 uniform layers with total depth of 3
1622  data h2 /6*0.5/ ! 6 uniform layers with total depth of 3
1623  type(remapping_cs) :: cs !< Remapping control structure
1624  real, allocatable, dimension(:,:) :: ppoly0_e, ppoly0_s, ppoly0_coefs
1625  logical :: answers_2018 ! If true use older, less acccurate expressions.
1626  integer :: i
1627  real :: err, h_neglect, h_neglect_edge
1628  logical :: thistest, v
1629 
1630  v = verbose
1631  h_neglect = hneglect_dflt
1632  h_neglect_edge = 1.0e-10
1633  answers_2018 = .true.
1634 
1635  write(*,*) '==== MOM_remapping: remapping_unit_tests ================='
1636  remapping_unit_tests = .false. ! Normally return false
1637 
1638  thistest = .false.
1639  call buildgridfromh(n0, h0, x0)
1640  do i=1,n0+1
1641  err=x0(i)-0.75*real(i-1)
1642  if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
1643  enddo
1644  if (thistest) write(*,*) 'remapping_unit_tests: Failed buildGridFromH() 1'
1646  call buildgridfromh(n1, h1, x1)
1647  do i=1,n1+1
1648  err=x1(i)-real(i-1)
1649  if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
1650  enddo
1651  if (thistest) write(*,*) 'remapping_unit_tests: Failed buildGridFromH() 2'
1653 
1654  thistest = .false.
1655  call initialize_remapping(cs, 'PPM_H4', answers_2018=answers_2018)
1656  if (verbose) write(*,*) 'h0 (test data)'
1657  if (verbose) call dumpgrid(n0,h0,x0,u0)
1658 
1659  call dzfromh1h2( n0, h0, n1, h1, dx1 )
1660  call remapping_core_w( cs, n0, h0, u0, n1, dx1, u1, h_neglect, h_neglect_edge)
1661  do i=1,n1
1662  err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1663  if (abs(err)>real(n1-1)*epsilon(err)) thistest = .true.
1664  enddo
1665  if (verbose) write(*,*) 'h1 (by projection)'
1666  if (verbose) call dumpgrid(n1,h1,x1,u1)
1667  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapping_core_w()'
1669 
1670  thistest = .false.
1671  allocate(ppoly0_e(n0,2))
1672  allocate(ppoly0_s(n0,2))
1673  allocate(ppoly0_coefs(n0,cs%degree+1))
1674 
1675  ppoly0_e(:,:) = 0.0
1676  ppoly0_s(:,:) = 0.0
1677  ppoly0_coefs(:,:) = 0.0
1678 
1679  call edge_values_explicit_h4( n0, h0, u0, ppoly0_e, h_neglect=1e-10, answers_2018=answers_2018 )
1680  call ppm_reconstruction( n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )
1681  call ppm_boundary_extrapolation( n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )
1682  u1(:) = 0.
1683  call remapbyprojection( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1684  n1, h1, integration_ppm, u1, h_neglect )
1685  do i=1,n1
1686  err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1687  if (abs(err)>2.*epsilon(err)) thistest = .true.
1688  enddo
1689  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapByProjection()'
1691 
1692  thistest = .false.
1693  u1(:) = 0.
1694  call remapbydeltaz( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1695  n1, x1-x0(1:n1+1), &
1696  integration_ppm, u1, hn1, h_neglect )
1697  if (verbose) write(*,*) 'h1 (by delta)'
1698  if (verbose) call dumpgrid(n1,h1,x1,u1)
1699  hn1=hn1-h1
1700  do i=1,n1
1701  err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1702  if (abs(err)>2.*epsilon(err)) thistest = .true.
1703  enddo
1704  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapByDeltaZ() 1'
1706 
1707  thistest = .false.
1708  call buildgridfromh(n2, h2, x2)
1709  dx2(1:n0+1) = x2(1:n0+1) - x0
1710  dx2(n0+2:n2+1) = x2(n0+2:n2+1) - x0(n0+1)
1711  call remapbydeltaz( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1712  n2, dx2, &
1713  integration_ppm, u2, hn2, h_neglect )
1714  if (verbose) write(*,*) 'h2'
1715  if (verbose) call dumpgrid(n2,h2,x2,u2)
1716  if (verbose) write(*,*) 'hn2'
1717  if (verbose) call dumpgrid(n2,hn2,x2,u2)
1718 
1719  do i=1,n2
1720  err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
1721  if (abs(err)>2.*epsilon(err)) thistest = .true.
1722  enddo
1723  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapByDeltaZ() 2'
1725 
1726  if (verbose) write(*,*) 'Via sub-cells'
1727  thistest = .false.
1728  call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1729  n2, h2, integration_ppm, .false., u2, err )
1730  if (verbose) call dumpgrid(n2,h2,x2,u2)
1731 
1732  do i=1,n2
1733  err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
1734  if (abs(err)>2.*epsilon(err)) thistest = .true.
1735  enddo
1736  if (thistest) write(*,*) 'remapping_unit_tests: Failed remap_via_sub_cells() 2'
1738 
1739  call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1740  6, (/.125,.125,.125,.125,.125,.125/), integration_ppm, .false., u2, err )
1741  if (verbose) call dumpgrid(6,h2,x2,u2)
1742 
1743  call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1744  3, (/2.25,1.5,1./), integration_ppm, .false., u2, err )
1745  if (verbose) call dumpgrid(3,h2,x2,u2)
1746 
1747  if (.not. remapping_unit_tests) write(*,*) 'Pass'
1748 
1749  write(*,*) '===== MOM_remapping: new remapping_unit_tests =================='
1750 
1751  deallocate(ppoly0_e, ppoly0_s, ppoly0_coefs)
1752  allocate(ppoly0_coefs(5,6))
1753  allocate(ppoly0_e(5,2))
1754  allocate(ppoly0_s(5,2))
1755 
1756  call pcm_reconstruction(3, (/1.,2.,4./), ppoly0_e(1:3,:), &
1757  ppoly0_coefs(1:3,:) )
1759  test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,4./), 'PCM: left edges')
1761  test_answer(v, 3, ppoly0_e(:,2), (/1.,2.,4./), 'PCM: right edges')
1763  test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,4./), 'PCM: P0')
1764 
1765  call plm_reconstruction(3, (/1.,1.,1./), (/1.,3.,5./), ppoly0_e(1:3,:), &
1766  ppoly0_coefs(1:3,:), h_neglect )
1768  test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,5./), 'Unlim PLM: left edges')
1770  test_answer(v, 3, ppoly0_e(:,2), (/1.,4.,5./), 'Unlim PLM: right edges')
1772  test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,5./), 'Unlim PLM: P0')
1774  test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Unlim PLM: P1')
1775 
1776  call plm_reconstruction(3, (/1.,1.,1./), (/1.,2.,7./), ppoly0_e(1:3,:), &
1777  ppoly0_coefs(1:3,:), h_neglect )
1779  test_answer(v, 3, ppoly0_e(:,1), (/1.,1.,7./), 'Left lim PLM: left edges')
1781  test_answer(v, 3, ppoly0_e(:,2), (/1.,3.,7./), 'Left lim PLM: right edges')
1783  test_answer(v, 3, ppoly0_coefs(:,1), (/1.,1.,7./), 'Left lim PLM: P0')
1785  test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Left lim PLM: P1')
1786 
1787  call plm_reconstruction(3, (/1.,1.,1./), (/1.,6.,7./), ppoly0_e(1:3,:), &
1788  ppoly0_coefs(1:3,:), h_neglect )
1790  test_answer(v, 3, ppoly0_e(:,1), (/1.,5.,7./), 'Right lim PLM: left edges')
1792  test_answer(v, 3, ppoly0_e(:,2), (/1.,7.,7./), 'Right lim PLM: right edges')
1794  test_answer(v, 3, ppoly0_coefs(:,1), (/1.,5.,7./), 'Right lim PLM: P0')
1796  test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Right lim PLM: P1')
1797 
1798  call plm_reconstruction(3, (/1.,2.,3./), (/1.,4.,9./), ppoly0_e(1:3,:), &
1799  ppoly0_coefs(1:3,:), h_neglect )
1801  test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,9./), 'Non-uniform line PLM: left edges')
1803  test_answer(v, 3, ppoly0_e(:,2), (/1.,6.,9./), 'Non-uniform line PLM: right edges')
1805  test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,9./), 'Non-uniform line PLM: P0')
1807  test_answer(v, 3, ppoly0_coefs(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1')
1808 
1809  call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_e, &
1810  h_neglect=1e-10, answers_2018=answers_2018 )
1811  ! The next two tests currently fail due to roundoff.
1812  thistest = test_answer(v, 5, ppoly0_e(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges')
1813  thistest = test_answer(v, 5, ppoly0_e(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges')
1814  ppoly0_e(:,1) = (/0.,2.,4.,6.,8./)
1815  ppoly0_e(:,2) = (/2.,4.,6.,8.,10./)
1816  call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_e(1:5,:), &
1817  ppoly0_coefs(1:5,:), h_neglect )
1819  test_answer(v, 5, ppoly0_coefs(:,1), (/1.,2.,4.,6.,9./), 'Line PPM: P0')
1821  test_answer(v, 5, ppoly0_coefs(:,2), (/0.,2.,2.,2.,0./), 'Line PPM: P1')
1823  test_answer(v, 5, ppoly0_coefs(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2')
1824 
1825  call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_e, &
1826  h_neglect=1e-10, answers_2018=answers_2018 )
1827  ! The next two tests currently fail due to roundoff.
1828  thistest = test_answer(v, 5, ppoly0_e(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges')
1829  thistest = test_answer(v, 5, ppoly0_e(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges')
1830  ppoly0_e(:,1) = (/0.,0.,3.,12.,27./)
1831  ppoly0_e(:,2) = (/0.,3.,12.,27.,48./)
1832  call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,1.,7.,19.,37./), ppoly0_e(1:5,:), &
1833  ppoly0_coefs(1:5,:), h_neglect )
1835  test_answer(v, 5, ppoly0_e(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: left edges')
1837  test_answer(v, 5, ppoly0_e(:,2), (/0.,3.,12.,27.,37./), 'Parabola PPM: right edges')
1839  test_answer(v, 5, ppoly0_coefs(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: P0')
1841  test_answer(v, 5, ppoly0_coefs(:,2), (/0.,0.,6.,12.,0./), 'Parabola PPM: P1')
1843  test_answer(v, 5, ppoly0_coefs(:,3), (/0.,3.,3.,3.,0./), 'Parabola PPM: P2')
1844 
1845  ppoly0_e(:,1) = (/0.,0.,6.,10.,15./)
1846  ppoly0_e(:,2) = (/0.,6.,12.,17.,15./)
1847  call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,5.,7.,16.,15./), ppoly0_e(1:5,:), &
1848  ppoly0_coefs(1:5,:), h_neglect )
1850  test_answer(v, 5, ppoly0_e(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: left edges')
1852  test_answer(v, 5, ppoly0_e(:,2), (/0.,6.,9.,16.,15./), 'Limits PPM: right edges')
1854  test_answer(v, 5, ppoly0_coefs(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: P0')
1856  test_answer(v, 5, ppoly0_coefs(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1')
1858  test_answer(v, 5, ppoly0_coefs(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2')
1859 
1860  call plm_reconstruction(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_e(1:4,:), &
1861  ppoly0_coefs(1:4,:), h_neglect )
1863  test_answer(v, 4, ppoly0_e(1:4,1), (/5.,5.,3.,1./), 'PPM: left edges h=0110')
1865  test_answer(v, 4, ppoly0_e(1:4,2), (/5.,3.,1.,1./), 'PPM: right edges h=0110')
1866  call remap_via_sub_cells( 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_e(1:4,:), &
1867  ppoly0_coefs(1:4,:), &
1868  2, (/1.,1./), integration_plm, .false., u2, err )
1870  test_answer(v, 2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11')
1871 
1872  deallocate(ppoly0_e, ppoly0_s, ppoly0_coefs)
1873 
1874  if (.not. remapping_unit_tests) write(*,*) 'Pass'
1875 
1876 end function remapping_unit_tests
1877 
1878 !> Returns true if any cell of u and u_true are not identical. Returns false otherwise.
1879 logical function test_answer(verbose, n, u, u_true, label)
1880  logical, intent(in) :: verbose !< If true, write results to stdout
1881  integer, intent(in) :: n !< Number of cells in u
1882  real, dimension(n), intent(in) :: u !< Values to test
1883  real, dimension(n), intent(in) :: u_true !< Values to test against (correct answer)
1884  character(len=*), intent(in) :: label !< Message
1885  ! Local variables
1886  integer :: k
1887 
1888  test_answer = .false.
1889  do k = 1, n
1890  if (u(k) /= u_true(k)) test_answer = .true.
1891  enddo
1892  if (test_answer .or. verbose) then
1893  write(*,'(a4,2a24,x,a)') 'k','Calculated value','Correct value',label
1894  do k = 1, n
1895  if (u(k) /= u_true(k)) then
1896  write(*,'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),' err=',u(k)-u_true(k),' < wrong'
1897  else
1898  write(*,'(i4,1p2e24.16)') k,u(k),u_true(k)
1899  endif
1900  enddo
1901  endif
1902 
1903 end function test_answer
1904 
1905 !> Convenience function for printing grid to screen
1906 subroutine dumpgrid(n,h,x,u)
1907  integer, intent(in) :: n !< Number of cells
1908  real, dimension(:), intent(in) :: h !< Cell thickness
1909  real, dimension(:), intent(in) :: x !< Interface delta
1910  real, dimension(:), intent(in) :: u !< Cell average values
1911  integer :: i
1912  write(*,'("i=",20i10)') (i,i=1,n+1)
1913  write(*,'("x=",20es10.2)') (x(i),i=1,n+1)
1914  write(*,'("i=",5x,20i10)') (i,i=1,n)
1915  write(*,'("h=",5x,20es10.2)') (h(i),i=1,n)
1916  write(*,'("u=",5x,20es10.2)') (u(i),i=1,n)
1917 end subroutine dumpgrid
1918 
1919 end module mom_remapping
ppm_functions::ppm_boundary_extrapolation
subroutine, public ppm_boundary_extrapolation(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Reconstruction by parabolas within boundary cells.
Definition: PPM_functions.F90:134
mom_remapping::measure_input_bounds
subroutine measure_input_bounds(n0, h0, u0, ppoly_E, h0tot, h0err, u0tot, u0err, u0min, u0max)
Measure totals and bounds on source grid.
Definition: MOM_remapping.F90:1027
mom_remapping::remapping_pqm_ih6ih5
integer, parameter remapping_pqm_ih6ih5
O(h^5) remapping scheme.
Definition: MOM_remapping.F90:52
plm_functions::plm_boundary_extrapolation
subroutine, public plm_boundary_extrapolation(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Reconstruction by linear polynomials within boundary cells.
Definition: PLM_functions.F90:206
regrid_edge_slopes::edge_slopes_implicit_h3
subroutine, public edge_slopes_implicit_h3(N, h, u, edge_slopes, h_neglect, answers_2018)
Compute ih4 edge slopes (implicit third order accurate) in the same units as h.
Definition: regrid_edge_slopes.F90:50
pcm_functions
Piecewise constant reconstruction functions.
Definition: PCM_functions.F90:2
ppm_functions
Provides functions used with the Piecewise-Parabolic-Method in the vertical ALE algorithm.
Definition: PPM_functions.F90:2
mom_remapping::setreconstructiontype
subroutine setreconstructiontype(string, CS)
Changes the method of reconstruction Use this routine to parse a string parameter specifying the reco...
Definition: MOM_remapping.F90:1568
mom_remapping::ispossumerrsignificant
logical function ispossumerrsignificant(n1, sum1, n2, sum2)
Compare two summation estimates of positive data and judge if due to more than round-off....
Definition: MOM_remapping.F90:162
mom_remapping::remapping_core_h
subroutine, public remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned.
Definition: MOM_remapping.F90:188
regrid_edge_slopes::edge_slopes_implicit_h5
subroutine, public edge_slopes_implicit_h5(N, h, u, edge_slopes, h_neglect, answers_2018)
Compute ih5 edge values (implicit fifth order accurate)
Definition: regrid_edge_slopes.F90:192
mom_remapping::buildgridfromh
subroutine buildgridfromh(nz, h, x)
Calculate edge coordinate x from cell width h.
Definition: MOM_remapping.F90:140
mom_remapping::remapping_ppm_h4
integer, parameter remapping_ppm_h4
O(h^3) remapping scheme.
Definition: MOM_remapping.F90:49
mom_remapping::remapbyprojection
subroutine remapbyprojection(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, u1, h_neglect)
Remaps column of values u0 on grid h0 to grid h1 by integrating over the projection of each h1 cell o...
Definition: MOM_remapping.F90:1096
mom_remapping::remapping_pqm_ih4ih3
integer, parameter remapping_pqm_ih4ih3
O(h^4) remapping scheme.
Definition: MOM_remapping.F90:51
mom_remapping::remapping_unit_tests
logical function, public remapping_unit_tests(verbose)
Runs unit tests on remapping functions. Should only be called from a single/root thread Returns True ...
Definition: MOM_remapping.F90:1613
mom_remapping::extract_member_remapping_cs
subroutine, public extract_member_remapping_cs(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell)
Definition: MOM_remapping.F90:120
mom_remapping::mdl
character(len=40) mdl
This module's name.
Definition: MOM_remapping.F90:59
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_remapping::build_reconstructions_1d
subroutine, public build_reconstructions_1d(CS, n0, h0, u0, ppoly_r_coefs, ppoly_r_E, ppoly_r_S, iMethod, h_neglect, h_neglect_edge)
Creates polynomial reconstructions of u0 on the source grid h0.
Definition: MOM_remapping.F90:356
mom_remapping::remapping_cs
Container for remapping parameters.
Definition: MOM_remapping.F90:22
mom_remapping::integration_pcm
integer, parameter integration_pcm
Piecewise Constant Method.
Definition: MOM_remapping.F90:54
regrid_edge_values
Edge value estimation for high-order resconstruction.
Definition: regrid_edge_values.F90:2
mom_remapping::remapping_ppm_ih4
integer, parameter remapping_ppm_ih4
O(h^3) remapping scheme.
Definition: MOM_remapping.F90:50
mom_remapping::hneglect_dflt
real, parameter hneglect_dflt
A thickness [H ~> m or kg m-2] that can be added to thicknesses in a denominator without changing the...
Definition: MOM_remapping.F90:76
mom_string_functions::uppercase
character(len=len(input_string)) function, public uppercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:42
mom_remapping
Provides column-wise vertical remapping functions.
Definition: MOM_remapping.F90:2
mom_remapping::integration_plm
integer, parameter integration_plm
Piecewise Linear Method.
Definition: MOM_remapping.F90:55
mom_remapping::remappingdefaultscheme
character(len=3), public remappingdefaultscheme
Default remapping method.
Definition: MOM_remapping.F90:69
mom_remapping::remappingschemesdoc
character(len=256), public remappingschemesdoc
Documentation for external callers.
Definition: MOM_remapping.F90:62
regrid_edge_values::edge_values_implicit_h4
subroutine, public edge_values_implicit_h4(N, h, u, edge_val, h_neglect, answers_2018)
Compute ih4 edge values (implicit fourth order accurate) in the same units as h.
Definition: regrid_edge_values.F90:492
mom_remapping::test_answer
logical function test_answer(verbose, n, u, u_true, label)
Returns true if any cell of u and u_true are not identical. Returns false otherwise.
Definition: MOM_remapping.F90:1880
mom_remapping::check_reconstructions_1d
subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, ppoly_r_coefs, ppoly_r_E, ppoly_r_S)
Checks that edge values and reconstructions satisfy bounds.
Definition: MOM_remapping.F90:442
mom_remapping::integraterecononinterval
subroutine integraterecononinterval(n0, h0, u0, ppoly0_E, ppoly0_coefs, method, xL, xR, hC, uAve, jStart, xStart, h_neglect)
Integrate the reconstructed column profile over a single cell.
Definition: MOM_remapping.F90:1223
pcm_functions::pcm_reconstruction
subroutine, public pcm_reconstruction(N, u, ppoly_E, ppoly_coef)
Reconstruction by constant polynomials within each cell. There is nothing to do but this routine is p...
Definition: PCM_functions.F90:19
pqm_functions
Piecewise quartic reconstruction functions.
Definition: PQM_functions.F90:2
mom_remapping::remapbydeltaz
subroutine remapbydeltaz(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, method, u1, h1, h_neglect)
Remaps column of values u0 on grid h0 to implied grid h1 where the interfaces of h1 differ from those...
Definition: MOM_remapping.F90:1144
mom_remapping::remapping_plm
integer, parameter remapping_plm
O(h^2) remapping scheme.
Definition: MOM_remapping.F90:48
pqm_functions::pqm_boundary_extrapolation_v1
subroutine, public pqm_boundary_extrapolation_v1(N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
Reconstruction by parabolas within boundary cells.
Definition: PQM_functions.F90:508
mom_remapping::average_value_ppoly
real function, public average_value_ppoly(n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb)
Returns the average value of a reconstruction within a single source cell, i0, between the non-dimens...
Definition: MOM_remapping.F90:926
plm_functions::plm_reconstruction
subroutine, public plm_reconstruction(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Reconstruction by linear polynomials within each cell.
Definition: PLM_functions.F90:19
mom_remapping::remapping_pcm
integer, parameter remapping_pcm
O(h^1) remapping scheme.
Definition: MOM_remapping.F90:47
mom_remapping::end_remapping
subroutine, public end_remapping(CS)
Destrcutor for remapping control structure.
Definition: MOM_remapping.F90:1603
mom_remapping::dumpgrid
subroutine dumpgrid(n, h, x, u)
Convenience function for printing grid to screen.
Definition: MOM_remapping.F90:1907
mom_remapping::integration_pqm
integer, parameter integration_pqm
Piecewise Quartic Method.
Definition: MOM_remapping.F90:57
mom_remapping::old_algorithm
logical, parameter old_algorithm
Use the old "broken" algorithm. This is a temporary measure to assist debugging until we delete the o...
Definition: MOM_remapping.F90:81
mom_remapping::remap_via_sub_cells
subroutine remap_via_sub_cells(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, force_bounds_in_subcell, u1, uh_err, ah_sub, aisub_src, aiss, aise)
Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating the n0+n1+1 sub-inte...
Definition: MOM_remapping.F90:516
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
regrid_edge_values::edge_values_implicit_h6
subroutine, public edge_values_implicit_h6(N, h, u, edge_val, h_neglect, answers_2018)
Compute ih6 edge values (implicit sixth order accurate) in the same units as h.
Definition: regrid_edge_values.F90:661
mom_remapping::measure_output_bounds
subroutine measure_output_bounds(n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max)
Measure totals and bounds on destination grid.
Definition: MOM_remapping.F90:1061
mom_remapping::remapping_core_w
subroutine, public remapping_core_w(CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_edge)
Remaps column of values u0 on grid h0 to implied grid h1 where the interfaces of h1 differ from those...
Definition: MOM_remapping.F90:266
mom_remapping::remapping_set_param
subroutine, public remapping_set_param(CS, remapping_scheme, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
Set parameters within remapping object.
Definition: MOM_remapping.F90:90
regrid_edge_slopes
Routines that estimate edge slopes to be used in high-order reconstruction schemes.
Definition: regrid_edge_slopes.F90:3
mom_remapping::initialize_remapping
subroutine, public initialize_remapping(CS, remapping_scheme, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
Constructor for remapping control structure.
Definition: MOM_remapping.F90:1547
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
plm_functions
Piecewise linear reconstruction functions.
Definition: PLM_functions.F90:2
pqm_functions::pqm_reconstruction
subroutine, public pqm_reconstruction(N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
Reconstruction by quartic polynomials within each cell.
Definition: PQM_functions.F90:21
mom_remapping::dzfromh1h2
subroutine, public dzfromh1h2(n1, h1, n2, h2, dx)
Calculates the change in interface positions based on h1 and h2.
Definition: MOM_remapping.F90:1522
regrid_edge_values::edge_values_explicit_h4
subroutine, public edge_values_explicit_h4(N, h, u, edge_val, h_neglect, answers_2018)
Compute h4 edge values (explicit fourth order accurate) in the same units as h.
Definition: regrid_edge_values.F90:289
mom_remapping::integration_ppm
integer, parameter integration_ppm
Piecewise Parabolic Method.
Definition: MOM_remapping.F90:56
ppm_functions::ppm_reconstruction
subroutine, public ppm_reconstruction(N, h, u, ppoly_E, ppoly_coef, h_neglect)
Builds quadratic polynomials coefficients from cell mean and edge values.
Definition: PPM_functions.F90:29