17 implicit none ;
private
19 #include <MOM_memory.h>
25 integer :: remapping_scheme = -911
29 logical :: boundary_extrapolation = .true.
31 logical :: check_reconstruction = .false.
33 logical :: check_remapping = .false.
35 logical :: force_bounds_in_subcell = .false.
37 logical :: answers_2018 = .true.
59 character(len=40) ::
mdl =
"MOM_remapping"
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"
74 #define __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
89 check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
91 character(len=*),
optional,
intent(in) :: remapping_scheme
92 logical,
optional,
intent(in) :: boundary_extrapolation
93 logical,
optional,
intent(in) :: check_reconstruction
94 logical,
optional,
intent(in) :: check_remapping
95 logical,
optional,
intent(in) :: force_bounds_in_subcell
96 logical,
optional,
intent(in) :: answers_2018
98 if (
present(remapping_scheme))
then
101 if (
present(boundary_extrapolation))
then
102 cs%boundary_extrapolation = boundary_extrapolation
104 if (
present(check_reconstruction))
then
105 cs%check_reconstruction = check_reconstruction
107 if (
present(check_remapping))
then
108 cs%check_remapping = check_remapping
110 if (
present(force_bounds_in_subcell))
then
111 cs%force_bounds_in_subcell = force_bounds_in_subcell
113 if (
present(answers_2018))
then
114 cs%answers_2018 = answers_2018
119 check_remapping, force_bounds_in_subcell)
121 integer,
optional,
intent(out) :: remapping_scheme
122 integer,
optional,
intent(out) :: degree
123 logical,
optional,
intent(out) :: boundary_extrapolation
124 logical,
optional,
intent(out) :: check_reconstruction
125 logical,
optional,
intent(out) :: check_remapping
127 logical,
optional,
intent(out) :: force_bounds_in_subcell
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
140 integer,
intent(in) :: nz
141 real,
dimension(nz),
intent(in) :: h
142 real,
dimension(nz+1),
intent(inout) :: x
162 integer,
intent(in) :: n1
163 integer,
intent(in) :: n2
164 real,
intent(in) :: sum1
165 real,
intent(in) :: sum2
168 real :: sumerr, allowederr, eps
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)
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
187 subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
189 integer,
intent(in) :: n0
190 real,
dimension(n0),
intent(in) :: h0
191 real,
dimension(n0),
intent(in) :: u0
192 integer,
intent(in) :: n1
193 real,
dimension(n1),
intent(in) :: h1
194 real,
dimension(n1),
intent(out) :: u1
195 real,
optional,
intent(in) :: h_neglect
198 real,
optional,
intent(in) :: h_neglect_edge
203 real,
dimension(n0,2) :: ppoly_r_e
204 real,
dimension(n0,2) :: ppoly_r_s
205 real,
dimension(n0,CS%degree+1) :: ppoly_r_coefs
207 real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
208 real :: hneglect, hneglect_edge
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
214 hneglect, hneglect_edge )
217 cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
221 cs%force_bounds_in_subcell, u1, uh_err )
223 if (cs%check_remapping)
then
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'
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)
246 write(0,
'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
248 write(0,
'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
251 write(0,
'(a3,2a24)')
'k',
'u0',
'Polynomial coefficients'
253 write(0,
'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
255 call mom_error( fatal,
'MOM_remapping, remapping_core_h: '//&
256 'Remapping result is inconsistent!' )
265 subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_edge )
267 integer,
intent(in) :: n0
268 real,
dimension(n0),
intent(in) :: h0
269 real,
dimension(n0),
intent(in) :: u0
270 integer,
intent(in) :: n1
271 real,
dimension(n1+1),
intent(in) :: dx
272 real,
dimension(n1),
intent(out) :: u1
273 real,
optional,
intent(in) :: h_neglect
276 real,
optional,
intent(in) :: h_neglect_edge
281 real,
dimension(n0,2) :: ppoly_r_e
282 real,
dimension(n0,2) :: ppoly_r_s
283 real,
dimension(n0,CS%degree+1) :: ppoly_r_coefs
285 real :: eps, h0tot, h0err, h1tot, h1err
286 real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
287 real,
dimension(n1) :: h1
288 real :: hneglect, hneglect_edge
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
294 hneglect, hneglect_edge )
297 cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
302 h1(k) = max( 0., h0(k) + ( dx(k+1) - dx(k) ) )
304 h1(k) = max( 0., dx(k+1) - dx(k) )
308 cs%force_bounds_in_subcell,u1, uh_err )
312 if (cs%check_remapping)
then
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'
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)
335 write(0,
'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
337 write(0,
'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
340 write(0,
'(a3,2a24)')
'k',
'u0',
'Polynomial coefficients'
342 write(0,
'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
344 call mom_error( fatal,
'MOM_remapping, remapping_core_w: '//&
345 'Remapping result is inconsistent!' )
354 ppoly_r_E, ppoly_r_S, iMethod, h_neglect, &
357 integer,
intent(in) :: n0
358 real,
dimension(n0),
intent(in) :: h0
359 real,
dimension(n0),
intent(in) :: u0
360 real,
dimension(n0,CS%degree+1), &
361 intent(out) :: ppoly_r_coefs
362 real,
dimension(n0,2),
intent(out) :: ppoly_r_e
363 real,
dimension(n0,2),
intent(out) :: ppoly_r_s
364 integer,
intent(out) :: imethod
365 real,
optional,
intent(in) :: h_neglect
368 real,
optional,
intent(in) :: h_neglect_edge
372 integer :: local_remapping_scheme
373 integer :: remapping_scheme
374 logical :: boundary_extrapolation
379 ppoly_r_coefs(:,:) = 0.0
382 local_remapping_scheme = cs%remapping_scheme
386 local_remapping_scheme = min( local_remapping_scheme,
remapping_plm )
390 select case ( local_remapping_scheme )
396 if ( cs%boundary_extrapolation )
then
401 call edge_values_explicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
403 if ( cs%boundary_extrapolation )
then
410 if ( cs%boundary_extrapolation )
then
418 if ( cs%boundary_extrapolation )
then
420 ppoly_r_coefs, h_neglect )
427 if ( cs%boundary_extrapolation )
then
429 ppoly_r_coefs, h_neglect )
433 call mom_error( fatal,
'MOM_remapping, build_reconstructions_1d: '//&
434 'The selected remapping method is invalid' )
441 ppoly_r_coefs, ppoly_r_E, ppoly_r_S)
442 integer,
intent(in) :: n0
443 real,
dimension(n0),
intent(in) :: h0
444 real,
dimension(n0),
intent(in) :: u0
445 integer,
intent(in) :: deg
446 logical,
intent(in) :: boundary_extrapolation
447 real,
dimension(n0,deg+1),
intent(out) :: ppoly_r_coefs
448 real,
dimension(n0,2),
intent(out) :: ppoly_r_E
449 real,
dimension(n0,2),
intent(out) :: ppoly_r_S
452 real :: u_l, u_c, u_r
454 logical :: problem_detected
456 problem_detected = .false.
458 u_l = u0(max(1,i0-1))
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.
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.
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.
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.
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.
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'
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,:)
504 call mom_error(fatal,
'MOM_remapping, check_reconstructions_1d: '// &
505 'Edge values or polynomial coefficients were inconsistent!')
515 force_bounds_in_subcell, u1, uh_err, ah_sub, aisub_src, aiss, aise )
516 integer,
intent(in) :: n0
517 real,
intent(in) :: h0(n0)
518 real,
intent(in) :: u0(n0)
519 real,
intent(in) :: ppoly0_E(n0,2)
520 real,
intent(in) :: ppoly0_coefs(:,:)
521 integer,
intent(in) :: n1
522 real,
intent(in) :: h1(n1)
523 integer,
intent(in) :: method
524 logical,
intent(in) :: force_bounds_in_subcell
525 real,
intent(out) :: u1(n1)
526 real,
intent(out) :: uh_err
527 real,
optional,
intent(out) :: ah_sub(n0+n1+1)
528 integer,
optional,
intent(out) :: aisub_src(n0+n1+1)
529 integer,
optional,
intent(out) :: aiss(n0)
530 integer,
optional,
intent(out) :: aise(n0)
539 real,
dimension(n0+n1+1) :: h_sub
540 real,
dimension(n0+n1+1) :: uh_sub
541 real,
dimension(n0+n1+1) :: u_sub
542 integer,
dimension(n0+n1+1) :: isub_src
543 integer,
dimension(n0) :: isrc_start
544 integer,
dimension(n0) :: isrc_end
545 integer,
dimension(n0) :: isrc_max
546 real,
dimension(n0) :: h0_eff
547 real,
dimension(n0) :: u0_min
548 real,
dimension(n0) :: u0_max
549 integer,
dimension(n1) :: itgt_start
550 integer,
dimension(n1) :: itgt_end
552 real :: h0_supply, h1_supply
557 logical,
parameter :: force_bounds_in_target = .true.
558 logical,
parameter :: adjust_thickest_subcell = .true.
559 logical,
parameter :: debug_bounds = .false.
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
564 logical :: tgt_has_volume
568 i0_last_thick_cell = 0
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
578 src_has_volume = .true.
579 tgt_has_volume = .true.
581 i_start0 = 1 ; i_start1 = 1
594 do i_sub = 2, n0+n1+1
598 dh = min(h0_supply, h1_supply)
603 dh0_eff = dh0_eff + min(dh, h0_supply)
612 if (dh >= dh_max)
then
619 if (h0_supply <= h1_supply .and. src_has_volume)
then
622 h1_supply = h1_supply - dh
625 isrc_start(i0) = i_start0
636 if (i0 < i0_last_thick_cell)
then
640 do while (h0_supply==0. .and. i0<i0_last_thick_cell)
655 src_has_volume = .false.
658 elseif (h0_supply >= h1_supply .and. tgt_has_volume)
then
661 h0_supply = h0_supply - dh
664 itgt_start(i1) = i_start1
676 tgt_has_volume = .false.
679 elseif (src_has_volume)
then
681 h_sub(i_sub) = h0_supply
683 isrc_start(i0) = i_start0
698 src_has_volume = .false.
700 elseif (tgt_has_volume)
then
702 h_sub(i_sub) = h1_supply
704 itgt_start(i1) = i_start1
713 tgt_has_volume = .false.
716 stop
'remap_via_sub_cells: THIS SHOULD NEVER HAPPEN!'
725 u_sub(1) = ppoly0_e(1,1)
738 dh0_eff = dh0_eff + dh
739 if (h0_eff(i0)>0.)
then
740 xb = dh0_eff / h0_eff(i0)
745 u_sub(i_sub) = u0(i0)
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!' )
760 if (force_bounds_in_subcell)
then
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 )
769 uh_sub(i_sub) = dh * u_sub(i_sub)
771 if (isub_src(i_sub+1) /= i0)
then
780 u_sub(n0+n1+1) = ppoly0_e(n0,2)
781 uh_sub(n0+n1+1) = ppoly0_e(n0,2) * h_sub(n0+n1+1)
783 if (adjust_thickest_subcell)
then
787 do i0 = 1, i0_last_thick_cell
789 dh_max = h_sub(i_max)
790 if (dh_max > 0.)
then
793 do i_sub = isrc_start(i0), isrc_end(i0)
794 if (i_sub /= i_max) duh = duh + uh_sub(i_sub)
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) )
805 if (h1(i1) > 0.)
then
807 i_sub = itgt_start(i1)
808 if (force_bounds_in_target)
then
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))
817 dh = dh + h_sub(i_sub)
818 duh = duh + uh_sub(i_sub)
820 uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh)
824 uh_err = uh_err + abs(duh)*epsilon(duh)
825 if (force_bounds_in_target)
then
827 u1(i1) = max(u1min, min(u1max, u1(i1)))
829 uh_err = uh_err + dh*abs( u1(i1)-u_orig )
832 u1(i1) = u_sub(itgt_start(i1))
837 if (debug_bounds)
then
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'
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)
879 write(0,
'(i3,96x,1p2e24.16,2i3)') k,h1(k),u1(k),itgt_start(k),itgt_end(k)
881 write(0,
'(i3,1p4e24.16)') k,h0(k),ppoly0_e(k,1),u0(k),ppoly0_e(k,2)
884 write(0,
'(a3,2a24)')
'k',
'u0',
'Polynomial coefficients'
886 write(0,
'(i3,1p6e24.16)') k,u0(k),ppoly0_coefs(k,:)
888 write(0,
'(a3,3a24,a3,2a24)')
'k',
'Sub-cell h',
'Sub-cell u',
'Sub-cell hu',
'i0',
'xa',
'xb'
894 dh0_eff = dh0_eff + dh
895 xb = dh0_eff / h0_eff(i0)
897 write(0,
'(i3,1p3e24.16,i3,1p2e24.16)') k,h_sub(k),u_sub(k),uh_sub(k),i0,xa,xb
899 if (isub_src(k+1) /= i0)
then
900 dh0_eff = 0.; xa = 0.
906 call mom_error( fatal,
'MOM_remapping, remap_via_sub_cells: '//&
907 'Remapping result is inconsistent!' )
913 uh_err = uh_err + u02_err
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)
926 integer,
intent(in) :: n0
927 real,
intent(in) :: u0(:)
928 real,
intent(in) :: ppoly0_e(:,:)
929 real,
intent(in) :: ppoly0_coefs(:,:)
930 integer,
intent(in) :: method
931 integer,
intent(in) :: i0
932 real,
intent(in) :: xa
933 real,
intent(in) :: xb
935 real :: u_ave, xa_2, xb_2, xa2pxb2, xapxb
936 real,
parameter :: r_3 = 1.0/3.0
938 real :: mx, a_l, a_r, u_c, ya, yb, my, xa2b2ab, ya2b2ab, a_c
941 select case ( method )
947 + ppoly0_coefs(i0,2) * 0.5 * ( xb + xa ) )
949 mx = 0.5 * ( xa + xb )
950 a_l = ppoly0_e(i0, 1)
951 a_r = ppoly0_e(i0, 2)
953 a_c = 0.5 * ( ( u_c - a_l ) + ( u_c - a_r ) )
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 ) )
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 ) )
971 xa2pxb2 = xa_2 + xb_2
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 ) ) ) ) )
980 call mom_error( fatal,
'The selected integration method is invalid' )
983 select case ( method )
985 u_ave = ppoly0_coefs(i0,1)
989 a_l = ppoly0_e(i0, 1)
990 a_r = ppoly0_e(i0, 2)
993 u_ave = a_l + xa * ( a_r - a_l )
995 u_ave = a_r + ya * ( a_l - a_r )
1001 a_l = ppoly0_e(i0, 1)
1002 a_r = ppoly0_e(i0, 2)
1004 a_c = 3. * ( ( u_c - a_l ) + ( u_c - a_r ) )
1007 u_ave = a_l + xa * ( ( a_r - a_l ) + a_c * ya )
1009 u_ave = a_r + ya * ( ( a_l - a_r ) + a_c * xa )
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) ) ) )
1018 call mom_error( fatal,
'The selected integration method is invalid' )
1026 subroutine measure_input_bounds( n0, h0, u0, ppoly_E, h0tot, h0err, u0tot, u0err, u0min, u0max )
1027 integer,
intent(in) :: n0
1028 real,
dimension(n0),
intent(in) :: h0
1029 real,
dimension(n0),
intent(in) :: u0
1030 real,
dimension(n0,2),
intent(in) :: ppoly_E
1031 real,
intent(out) :: h0tot
1032 real,
intent(out) :: h0err
1033 real,
intent(out) :: u0tot
1034 real,
intent(out) :: u0err
1035 real,
intent(out) :: u0min
1036 real,
intent(out) :: u0max
1041 eps = epsilon(h0(1))
1044 u0tot = h0(1) * u0(1)
1046 u0min = min( ppoly_e(1,1), ppoly_e(1,2) )
1047 u0max = max( ppoly_e(1,1), ppoly_e(1,2) )
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) )
1061 integer,
intent(in) :: n1
1062 real,
dimension(n1),
intent(in) :: h1
1063 real,
dimension(n1),
intent(in) :: u1
1064 real,
intent(out) :: h1tot
1065 real,
intent(out) :: h1err
1066 real,
intent(out) :: u1tot
1067 real,
intent(out) :: u1err
1068 real,
intent(out) :: u1min
1069 real,
intent(out) :: u1max
1074 eps = epsilon(h1(1))
1077 u1tot = h1(1) * u1(1)
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))
1095 n1, h1, method, u1, h_neglect )
1096 integer,
intent(in) :: n0
1097 real,
intent(in) :: h0(:)
1098 real,
intent(in) :: u0(:)
1099 real,
intent(in) :: ppoly0_E(:,:)
1100 real,
intent(in) :: ppoly0_coefs(:,:)
1101 integer,
intent(in) :: n1
1102 real,
intent(in) :: h1(:)
1103 integer,
intent(in) :: method
1104 real,
intent(out) :: u1(:)
1105 real,
optional,
intent(in) :: h_neglect
1123 xr = xl + h1(itarget)
1126 xl, xr, h1(itarget), u1(itarget), jstart, xstart, h_neglect )
1142 subroutine remapbydeltaz( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, &
1143 method, u1, h1, h_neglect )
1144 integer,
intent(in) :: n0
1145 real,
dimension(:),
intent(in) :: h0
1146 real,
dimension(:),
intent(in) :: u0
1147 real,
dimension(:,:),
intent(in) :: ppoly0_E
1148 real,
dimension(:,:),
intent(in) :: ppoly0_coefs
1149 integer,
intent(in) :: n1
1150 real,
dimension(:),
intent(in) :: dx1
1151 integer,
intent(in) :: method
1152 real,
dimension(:),
intent(out) :: u1
1153 real,
dimension(:), &
1154 optional,
intent(out) :: h1
1155 real,
optional,
intent(in) :: h_neglect
1161 real :: xOld, hOld, uOld
1162 real :: xNew, hNew, h_err
1163 real :: uhNew, hFlux, uAve, fluxL, fluxR
1179 if (itarget == 0)
then
1183 elseif (itarget <= n0)
then
1184 xold = xold + h0(itarget)
1187 h_err = h_err + epsilon(hold) * max(hold, xold)
1192 xnew = xold + dx1(itarget+1)
1193 xl = min( xold, xnew )
1194 xr = max( xold, xnew )
1197 hflux = abs(dx1(itarget+1))
1199 xl, xr, hflux, uave, jstart, xstart )
1201 fluxr = dx1(itarget+1)*uave
1204 hnew = hold + ( dx1(itarget+1) - dx1(itarget) )
1205 hnew = max( 0., hnew )
1206 uhnew = ( uold * hold ) + ( fluxr - fluxl )
1208 u1(itarget) = uhnew / hnew
1212 if (
present(h1)) h1(itarget) = hnew
1222 xL, xR, hC, uAve, jStart, xStart, h_neglect )
1223 integer,
intent(in) :: n0
1224 real,
dimension(:),
intent(in) :: h0
1225 real,
dimension(:),
intent(in) :: u0
1226 real,
dimension(:,:),
intent(in) :: ppoly0_E
1227 real,
dimension(:,:),
intent(in) :: ppoly0_coefs
1228 integer,
intent(in) :: method
1229 real,
intent(in) :: xL
1230 real,
intent(in) :: xR
1231 real,
intent(in) :: hC
1232 real,
intent(out) :: uAve
1233 integer,
intent(inout) :: jStart
1235 real,
intent(inout) :: xStart
1237 real,
optional,
intent(in) :: h_neglect
1247 real :: x0jLl, x0jLr
1248 real :: x0jRl, x0jRr
1251 real :: x0_2, x1_2, x02px12, x0px1
1253 real,
parameter :: r_3 = 1.0/3.0
1255 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
1266 x0jlr = x0jll + h0(j)
1268 if ( ( xl >= x0jll ) .AND. ( xl <= x0jlr ) )
then
1285 'MOM_remapping, integrateReconOnInterval: '//&
1286 'The location of the left-most cell could not be found')
1298 if ( abs(xr - xl) == 0.0 )
then
1304 if ( h0(jl) == 0.0 )
then
1305 uave = 0.5 * ( ppoly0_e(jl,1) + ppoly0_e(jl,2) )
1308 xi0 = xl / ( h0(jl) + hneglect ) - x0jll / ( h0(jl) + hneglect )
1310 select case ( method )
1312 uave = ppoly0_coefs(jl,1)
1314 uave = ppoly0_coefs(jl,1) &
1315 + xi0 * ppoly0_coefs(jl,2)
1317 uave = ppoly0_coefs(jl,1) &
1318 + xi0 * ( ppoly0_coefs(jl,2) &
1319 + xi0 * ppoly0_coefs(jl,3) )
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) ) ) )
1327 call mom_error( fatal,
'The selected integration method is invalid' )
1340 x0jrr = x0jrl + h0(j)
1342 if ( ( xr >= x0jrl ) .AND. ( xr <= x0jrr ) )
then
1351 if (xr>x0jrr) jr = n0
1357 if ( jl == jr )
then
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 ) ) )
1370 xi0 = xl / h0(jl) - x0jll / ( h0(jl) + hneglect )
1371 xi1 = xr / h0(jl) - x0jll / ( h0(jl) + hneglect )
1374 hact = h0(jl) * ( xi1 - xi0 )
1379 select case ( method )
1381 q = ( xr - xl ) * ppoly0_coefs(jl,1)
1383 q = ( xr - xl ) * ( &
1384 ppoly0_coefs(jl,1) &
1385 + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
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 ) ) )
1394 x02px12 = x0_2 + x1_2
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 ) ) ) )
1403 call mom_error( fatal,
'The selected integration method is invalid' )
1422 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1423 xi0 = max( 0., min( 1., ( xl - x0jll ) / ( h0(jl) + hneglect ) ) )
1425 xi0 = (xl - x0jll) / ( h0(jl) + hneglect )
1429 hact = h0(jl) * ( xi1 - xi0 )
1431 select case ( method )
1433 q = q + ( x0jlr - xl ) * ppoly0_coefs(jl,1)
1435 q = q + ( x0jlr - xl ) * ( &
1436 ppoly0_coefs(jl,1) &
1437 + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
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 ) ) )
1446 x02px12 = x0_2 + x1_2
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 ) ) ) )
1455 call mom_error( fatal,
'The selected integration method is invalid' )
1459 if ( jr > (jl+1) )
then
1461 q = q + h0(k) * u0(k)
1468 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1469 xi1 = max( 0., min( 1., ( xr - x0jrl ) / ( h0(jr) + hneglect ) ) )
1471 xi1 = (xr - x0jrl) / ( h0(jr) + hneglect )
1474 hact = hact + h0(jr) * ( xi1 - xi0 )
1476 select case ( method )
1478 q = q + ( xr - x0jrl ) * ppoly0_coefs(jr,1)
1480 q = q + ( xr - x0jrl ) * ( &
1481 ppoly0_coefs(jr,1) &
1482 + ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) )
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 ) ) )
1491 x02px12 = x0_2 + x1_2
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 ) ) ) )
1500 call mom_error( fatal,
'The selected integration method is invalid' )
1506 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1508 uave = ppoly0_coefs(jl,1)
1522 integer,
intent(in) :: n1
1523 real,
dimension(:),
intent(in) :: h1
1524 integer,
intent(in) :: n2
1525 real,
dimension(:),
intent(in) :: h2
1526 real,
dimension(:),
intent(out) :: dx
1534 do k = 1, max(n1,n2)
1535 if (k <= n1) x1 = x1 + h1(k)
1546 check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
1549 character(len=*),
intent(in) :: remapping_scheme
1550 logical,
optional,
intent(in) :: boundary_extrapolation
1551 logical,
optional,
intent(in) :: check_reconstruction
1552 logical,
optional,
intent(in) :: check_remapping
1553 logical,
optional,
intent(in) :: force_bounds_in_subcell
1554 logical,
optional,
intent(in) :: answers_2018
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)
1568 character(len=*),
intent(in) :: string
1593 call mom_error(fatal,
"setReconstructionType: "//&
1594 "Unrecognized choice for REMAPPING_SCHEME ("//trim(string)//
").")
1613 logical,
intent(in) :: verbose
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./
1624 real,
allocatable,
dimension(:,:) :: ppoly0_e, ppoly0_s, ppoly0_coefs
1625 logical :: answers_2018
1627 real :: err, h_neglect, h_neglect_edge
1628 logical :: thistest, v
1632 h_neglect_edge = 1.0e-10
1633 answers_2018 = .true.
1635 write(*,*)
'==== MOM_remapping: remapping_unit_tests ================='
1641 err=x0(i)-0.75*real(i-1)
1642 if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
1644 if (thistest)
write(*,*)
'remapping_unit_tests: Failed buildGridFromH() 1'
1649 if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
1651 if (thistest)
write(*,*)
'remapping_unit_tests: Failed buildGridFromH() 2'
1656 if (verbose)
write(*,*)
'h0 (test data)'
1657 if (verbose)
call dumpgrid(n0,h0,x0,u0)
1660 call remapping_core_w( cs, n0, h0, u0, n1, dx1, u1, h_neglect, h_neglect_edge)
1662 err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1663 if (abs(err)>real(n1-1)*epsilon(err)) thistest = .true.
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()'
1671 allocate(ppoly0_e(n0,2))
1672 allocate(ppoly0_s(n0,2))
1673 allocate(ppoly0_coefs(n0,cs%degree+1))
1677 ppoly0_coefs(:,:) = 0.0
1679 call edge_values_explicit_h4( n0, h0, u0, ppoly0_e, h_neglect=1e-10, answers_2018=answers_2018 )
1686 err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1687 if (abs(err)>2.*epsilon(err)) thistest = .true.
1689 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remapByProjection()'
1695 n1, x1-x0(1:n1+1), &
1697 if (verbose)
write(*,*)
'h1 (by delta)'
1698 if (verbose)
call dumpgrid(n1,h1,x1,u1)
1701 err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1702 if (abs(err)>2.*epsilon(err)) thistest = .true.
1704 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remapByDeltaZ() 1'
1709 dx2(1:n0+1) = x2(1:n0+1) - x0
1710 dx2(n0+2:n2+1) = x2(n0+2:n2+1) - x0(n0+1)
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)
1720 err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
1721 if (abs(err)>2.*epsilon(err)) thistest = .true.
1723 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remapByDeltaZ() 2'
1726 if (verbose)
write(*,*)
'Via sub-cells'
1730 if (verbose)
call dumpgrid(n2,h2,x2,u2)
1733 err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
1734 if (abs(err)>2.*epsilon(err)) thistest = .true.
1736 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remap_via_sub_cells() 2'
1740 6, (/.125,.125,.125,.125,.125,.125/),
integration_ppm, .false., u2, err )
1741 if (verbose)
call dumpgrid(6,h2,x2,u2)
1745 if (verbose)
call dumpgrid(3,h2,x2,u2)
1749 write(*,*)
'===== MOM_remapping: new remapping_unit_tests =================='
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))
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')
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')
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')
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')
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')
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 )
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')
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 )
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')
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')
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,:), &
1870 test_answer(v, 2, u2, (/4.,2./),
'PLM: remapped h=0110->h=11')
1872 deallocate(ppoly0_e, ppoly0_s, ppoly0_coefs)
1879 logical function test_answer(verbose, n, u, u_true, label)
1880 logical,
intent(in) :: verbose
1881 integer,
intent(in) :: n
1882 real,
dimension(n),
intent(in) :: u
1883 real,
dimension(n),
intent(in) :: u_true
1884 character(len=*),
intent(in) :: label
1893 write(*,
'(a4,2a24,x,a)')
'k',
'Calculated value',
'Correct value',label
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'
1898 write(*,
'(i4,1p2e24.16)') k,u(k),u_true(k)
1907 integer,
intent(in) :: n
1908 real,
dimension(:),
intent(in) :: h
1909 real,
dimension(:),
intent(in) :: x
1910 real,
dimension(:),
intent(in) :: u
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)