13 use plm_functions,
only : plm_reconstruction, plm_boundary_extrapolation
14 use ppm_functions,
only : ppm_reconstruction, ppm_boundary_extrapolation
15 use pqm_functions,
only : pqm_reconstruction, pqm_boundary_extrapolation_v1
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.
39 public remapping_core_h, remapping_core_w
40 public initialize_remapping, end_remapping, remapping_set_param, extract_member_remapping_cs
41 public remapping_unit_tests, build_reconstructions_1d, average_value_ppoly
45 integer,
parameter :: remapping_pcm = 0
46 integer,
parameter :: remapping_plm = 1
47 integer,
parameter :: remapping_ppm_h4 = 2
48 integer,
parameter :: remapping_ppm_ih4 = 3
49 integer,
parameter :: remapping_pqm_ih4ih3 = 4
50 integer,
parameter :: remapping_pqm_ih6ih5 = 5
52 integer,
parameter :: integration_pcm = 0
53 integer,
parameter :: integration_plm = 1
54 integer,
parameter :: integration_ppm = 3
55 integer,
parameter :: integration_pqm = 5
57 character(len=40) :: mdl =
"MOM_remapping"
60 character(len=256),
public :: remappingschemesdoc = &
61 "PCM (1st-order accurate)\n"//&
62 "PLM (2nd-order accurate)\n"//&
63 "PPM_H4 (3rd-order accurate)\n"//&
64 "PPM_IH4 (3rd-order accurate)\n"//&
65 "PQM_IH4IH3 (4th-order accurate)\n"//&
66 "PQM_IH6IH5 (5th-order accurate)\n"
67 character(len=3),
public :: remappingdefaultscheme =
"PLM"
72 #define __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
74 real,
parameter :: hneglect_dflt = 1.e-30
79 logical,
parameter :: old_algorithm = .false.
86 subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, &
87 check_reconstruction, check_remapping, force_bounds_in_subcell)
89 character(len=*),
optional,
intent(in) :: remapping_scheme
90 logical,
optional,
intent(in) :: boundary_extrapolation
91 logical,
optional,
intent(in) :: check_reconstruction
92 logical,
optional,
intent(in) :: check_remapping
93 logical,
optional,
intent(in) :: force_bounds_in_subcell
95 if (
present(remapping_scheme))
then
96 call setreconstructiontype( remapping_scheme, cs )
98 if (
present(boundary_extrapolation))
then
99 cs%boundary_extrapolation = boundary_extrapolation
101 if (
present(check_reconstruction))
then
102 cs%check_reconstruction = check_reconstruction
104 if (
present(check_remapping))
then
105 cs%check_remapping = check_remapping
107 if (
present(force_bounds_in_subcell))
then
108 cs%force_bounds_in_subcell = force_bounds_in_subcell
110 end subroutine remapping_set_param
112 subroutine extract_member_remapping_cs(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, &
113 check_remapping, force_bounds_in_subcell)
115 integer,
optional,
intent(out) :: remapping_scheme
116 integer,
optional,
intent(out) :: degree
117 logical,
optional,
intent(out) :: boundary_extrapolation
118 logical,
optional,
intent(out) :: check_reconstruction
119 logical,
optional,
intent(out) :: check_remapping
121 logical,
optional,
intent(out) :: force_bounds_in_subcell
124 if (
present(remapping_scheme)) remapping_scheme = cs%remapping_scheme
125 if (
present(degree)) degree = cs%degree
126 if (
present(boundary_extrapolation)) boundary_extrapolation = cs%boundary_extrapolation
127 if (
present(check_reconstruction)) check_reconstruction = cs%check_reconstruction
128 if (
present(check_remapping)) check_remapping = cs%check_remapping
129 if (
present(force_bounds_in_subcell)) force_bounds_in_subcell = cs%force_bounds_in_subcell
131 end subroutine extract_member_remapping_cs
133 subroutine buildgridfromh(nz, h, x)
134 integer,
intent(in) :: nz
135 real,
dimension(nz),
intent(in) :: h
136 real,
dimension(nz+1),
intent(inout) :: x
145 end subroutine buildgridfromh
155 function ispossumerrsignificant(n1, sum1, n2, sum2)
156 integer,
intent(in) :: n1
157 integer,
intent(in) :: n2
158 real,
intent(in) :: sum1
159 real,
intent(in) :: sum2
160 logical :: ispossumerrsignificant
162 real :: sumerr, allowederr, eps
164 if (sum1<0.)
call mom_error(fatal,
'isPosSumErrSignificant: sum1<0 is not allowed!')
165 if (sum2<0.)
call mom_error(fatal,
'isPosSumErrSignificant: sum2<0 is not allowed!')
166 sumerr = abs(sum1-sum2)
168 allowederr = eps*0.5*(real(n1-1)*sum1+real(n2-1)*sum2)
169 if (sumerr>allowederr)
then
170 write(0,*)
'isPosSumErrSignificant: sum1,sum2=',sum1,sum2
171 write(0,*)
'isPosSumErrSignificant: eps=',eps
172 write(0,*)
'isPosSumErrSignificant: err,n*eps=',sumerr,allowederr
173 write(0,*)
'isPosSumErrSignificant: err/eps,n1,n2,n1+n2=',sumerr/eps,n1,n2,n1+n2
174 ispossumerrsignificant = .true.
176 ispossumerrsignificant = .false.
178 end function ispossumerrsignificant
181 subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
183 integer,
intent(in) :: n0
184 real,
dimension(n0),
intent(in) :: h0
185 real,
dimension(n0),
intent(in) :: u0
186 integer,
intent(in) :: n1
187 real,
dimension(n1),
intent(in) :: h1
188 real,
dimension(n1),
intent(out) :: u1
189 real,
optional,
intent(in) :: h_neglect
192 real,
optional,
intent(in) :: h_neglect_edge
197 real,
dimension(n0,2) :: ppoly_r_e
198 real,
dimension(n0,2) :: ppoly_r_s
199 real,
dimension(n0,CS%degree+1) :: ppoly_r_coefs
201 real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
202 real :: hneglect, hneglect_edge
204 hneglect = 1.0e-30 ;
if (
present(h_neglect)) hneglect = h_neglect
205 hneglect_edge = 1.0e-10 ;
if (
present(h_neglect_edge)) hneglect_edge = h_neglect_edge
207 call build_reconstructions_1d( cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod, &
208 hneglect, hneglect_edge )
210 if (cs%check_reconstruction)
call check_reconstructions_1d(n0, h0, u0, cs%degree, &
211 cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
214 call remap_via_sub_cells( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h1, imethod, &
215 cs%force_bounds_in_subcell, u1, uh_err )
217 if (cs%check_remapping)
then
219 call measure_input_bounds( n0, h0, u0, ppoly_r_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
220 call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
222 if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err .and. abs(h1tot-h0tot)<h0err+h1err) &
223 .or. (u1min<u0min .or. u1max>u0max) )
then
224 write(0,*)
'iMethod = ',imethod
225 write(0,*)
'H: h0tot=',h0tot,
'h1tot=',h1tot,
'dh=',h1tot-h0tot,
'h0err=',h0err,
'h1err=',h1err
226 if (abs(h1tot-h0tot)>h0err+h1err) &
227 write(0,*)
'H non-conservation difference=',h1tot-h0tot,
'allowed err=',h0err+h1err,
' <-----!'
228 write(0,*)
'UH: u0tot=',u0tot,
'u1tot=',u1tot,
'duh=',u1tot-u0tot,
'u0err=',u0err,
'u1err=',u1err,
'uh_err=',uh_err
229 if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
230 write(0,*)
'U non-conservation difference=',u1tot-u0tot,
'allowed err=',u0err+u1err+uh_err,
' <-----!'
231 write(0,*)
'U: u0min=',u0min,
'u1min=',u1min
232 if (u1min<u0min)
write(0,*)
'U minimum overshoot=',u1min-u0min,
' <-----!'
233 write(0,*)
'U: u0max=',u0max,
'u1max=',u1max
234 if (u1max>u0max)
write(0,*)
'U maximum overshoot=',u1max-u0max,
' <-----!'
235 write(0,
'(a3,6a24)')
'k',
'h0',
'left edge',
'u0',
'right edge',
'h1',
'u1'
237 if (k<=min(n0,n1))
then
238 write(0,
'(i3,1p6e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2),h1(k),u1(k)
240 write(0,
'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
242 write(0,
'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
245 write(0,
'(a3,2a24)')
'k',
'u0',
'Polynomial coefficients'
247 write(0,
'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
249 call mom_error( fatal,
'MOM_remapping, remapping_core_h: '//&
250 'Remapping result is inconsistent!' )
255 end subroutine remapping_core_h
259 subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_edge )
261 integer,
intent(in) :: n0
262 real,
dimension(n0),
intent(in) :: h0
263 real,
dimension(n0),
intent(in) :: u0
264 integer,
intent(in) :: n1
265 real,
dimension(n1+1),
intent(in) :: dx
266 real,
dimension(n1),
intent(out) :: u1
267 real,
optional,
intent(in) :: h_neglect
270 real,
optional,
intent(in) :: h_neglect_edge
275 real,
dimension(n0,2) :: ppoly_r_e
276 real,
dimension(n0,2) :: ppoly_r_s
277 real,
dimension(n0,CS%degree+1) :: ppoly_r_coefs
279 real :: eps, h0tot, h0err, h1tot, h1err
280 real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
281 real,
dimension(n1) :: h1
282 real :: hneglect, hneglect_edge
284 hneglect = 1.0e-30 ;
if (
present(h_neglect)) hneglect = h_neglect
285 hneglect_edge = 1.0e-10 ;
if (
present(h_neglect_edge)) hneglect_edge = h_neglect_edge
287 call build_reconstructions_1d( cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod,&
288 hneglect, hneglect_edge )
290 if (cs%check_reconstruction)
call check_reconstructions_1d(n0, h0, u0, cs%degree, &
291 cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
296 h1(k) = max( 0., h0(k) + ( dx(k+1) - dx(k) ) )
298 h1(k) = max( 0., dx(k+1) - dx(k) )
301 call remap_via_sub_cells( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h1, imethod, &
302 cs%force_bounds_in_subcell,u1, uh_err )
306 if (cs%check_remapping)
then
308 call measure_input_bounds( n0, h0, u0, ppoly_r_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
309 call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
311 if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err .and. abs(h1tot-h0tot)<h0err+h1err) &
312 .or. (u1min<u0min .or. u1max>u0max) )
then
313 write(0,*)
'iMethod = ',imethod
314 write(0,*)
'H: h0tot=',h0tot,
'h1tot=',h1tot,
'dh=',h1tot-h0tot,
'h0err=',h0err,
'h1err=',h1err
315 if (abs(h1tot-h0tot)>h0err+h1err) &
316 write(0,*)
'H non-conservation difference=',h1tot-h0tot,
'allowed err=',h0err+h1err,
' <-----!'
317 write(0,*)
'UH: u0tot=',u0tot,
'u1tot=',u1tot,
'duh=',u1tot-u0tot,
'u0err=',u0err,
'u1err=',u1err,
'uh_err=',uh_err
318 if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
319 write(0,*)
'U non-conservation difference=',u1tot-u0tot,
'allowed err=',u0err+u1err+uh_err,
' <-----!'
320 write(0,*)
'U: u0min=',u0min,
'u1min=',u1min
321 if (u1min<u0min)
write(0,*)
'U minimum overshoot=',u1min-u0min,
' <-----!'
322 write(0,*)
'U: u0max=',u0max,
'u1max=',u1max
323 if (u1max>u0max)
write(0,*)
'U maximum overshoot=',u1max-u0max,
' <-----!'
324 write(0,
'(a3,6a24)')
'k',
'h0',
'left edge',
'u0',
'right edge',
'h1',
'u1'
326 if (k<=min(n0,n1))
then
327 write(0,
'(i3,1p6e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2),h1(k),u1(k)
329 write(0,
'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
331 write(0,
'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
334 write(0,
'(a3,2a24)')
'k',
'u0',
'Polynomial coefficients'
336 write(0,
'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
338 call mom_error( fatal,
'MOM_remapping, remapping_core_w: '//&
339 'Remapping result is inconsistent!' )
344 end subroutine remapping_core_w
347 subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, &
348 ppoly_r_E, ppoly_r_S, iMethod, h_neglect, &
351 integer,
intent(in) :: n0
352 real,
dimension(n0),
intent(in) :: h0
353 real,
dimension(n0),
intent(in) :: u0
354 real,
dimension(n0,CS%degree+1), &
355 intent(out) :: ppoly_r_coefs
356 real,
dimension(n0,2),
intent(out) :: ppoly_r_e
357 real,
dimension(n0,2),
intent(out) :: ppoly_r_s
358 integer,
intent(out) :: imethod
359 real,
optional,
intent(in) :: h_neglect
362 real,
optional,
intent(in) :: h_neglect_edge
366 integer :: local_remapping_scheme
367 integer :: remapping_scheme
368 logical :: boundary_extrapolation
373 ppoly_r_coefs(:,:) = 0.0
376 local_remapping_scheme = cs%remapping_scheme
378 local_remapping_scheme = remapping_pcm
380 local_remapping_scheme = min( local_remapping_scheme, remapping_plm )
382 local_remapping_scheme = min( local_remapping_scheme, remapping_ppm_h4 )
384 select case ( local_remapping_scheme )
385 case ( remapping_pcm )
386 call pcm_reconstruction( n0, u0, ppoly_r_e, ppoly_r_coefs)
387 imethod = integration_pcm
388 case ( remapping_plm )
389 call plm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
390 if ( cs%boundary_extrapolation )
then
391 call plm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect)
393 imethod = integration_plm
394 case ( remapping_ppm_h4 )
395 call edge_values_explicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge )
396 call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
397 if ( cs%boundary_extrapolation )
then
398 call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
400 imethod = integration_ppm
401 case ( remapping_ppm_ih4 )
402 call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge )
403 call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
404 if ( cs%boundary_extrapolation )
then
405 call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
407 imethod = integration_ppm
408 case ( remapping_pqm_ih4ih3 )
409 call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge )
410 call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_s, h_neglect )
411 call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect )
412 if ( cs%boundary_extrapolation )
then
413 call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &
414 ppoly_r_coefs, h_neglect )
416 imethod = integration_pqm
417 case ( remapping_pqm_ih6ih5 )
418 call edge_values_implicit_h6( n0, h0, u0, ppoly_r_e, h_neglect_edge )
419 call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_s, h_neglect )
420 call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect )
421 if ( cs%boundary_extrapolation )
then
422 call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &
423 ppoly_r_coefs, h_neglect )
425 imethod = integration_pqm
427 call mom_error( fatal,
'MOM_remapping, build_reconstructions_1d: '//&
428 'The selected remapping method is invalid' )
431 end subroutine build_reconstructions_1d
434 subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, &
435 ppoly_r_coefs, ppoly_r_E, ppoly_r_S)
436 integer,
intent(in) :: n0
437 real,
dimension(n0),
intent(in) :: h0
438 real,
dimension(n0),
intent(in) :: u0
439 integer,
intent(in) :: deg
440 logical,
intent(in) :: boundary_extrapolation
441 real,
dimension(n0,deg+1),
intent(out) :: ppoly_r_coefs
442 real,
dimension(n0,2),
intent(out) :: ppoly_r_E
443 real,
dimension(n0,2),
intent(out) :: ppoly_r_S
446 real :: u_l, u_c, u_r
448 logical :: problem_detected
450 problem_detected = .false.
452 u_l = u0(max(1,i0-1))
454 u_r = u0(min(n0,i0+1))
455 if (i0 > 1 .or. .not. boundary_extrapolation)
then
456 u_min = min(u_l, u_c)
457 u_max = max(u_l, u_c)
458 if (ppoly_r_e(i0,1) < u_min)
then
459 write(0,
'(a,i4,5(x,a,1pe24.16))')
'Left edge undershoot at',i0,
'u(i0-1)=',u_l,
'u(i0)=',u_c, &
460 'edge=',ppoly_r_e(i0,1),
'err=',ppoly_r_e(i0,1)-u_min
461 problem_detected = .true.
463 if (ppoly_r_e(i0,1) > u_max)
then
464 write(0,
'(a,i4,5(x,a,1pe24.16))')
'Left edge overshoot at',i0,
'u(i0-1)=',u_l,
'u(i0)=',u_c, &
465 'edge=',ppoly_r_e(i0,1),
'err=',ppoly_r_e(i0,1)-u_max
466 problem_detected = .true.
469 if (i0 < n0 .or. .not. boundary_extrapolation)
then
470 u_min = min(u_c, u_r)
471 u_max = max(u_c, u_r)
472 if (ppoly_r_e(i0,2) < u_min)
then
473 write(0,
'(a,i4,5(x,a,1pe24.16))')
'Right edge undershoot at',i0,
'u(i0)=',u_c,
'u(i0+1)=',u_r, &
474 'edge=',ppoly_r_e(i0,2),
'err=',ppoly_r_e(i0,2)-u_min
475 problem_detected = .true.
477 if (ppoly_r_e(i0,2) > u_max)
then
478 write(0,
'(a,i4,5(x,a,1pe24.16))')
'Right edge overshoot at',i0,
'u(i0)=',u_c,
'u(i0+1)=',u_r, &
479 'edge=',ppoly_r_e(i0,2),
'err=',ppoly_r_e(i0,2)-u_max
480 problem_detected = .true.
484 if ( (u_c-u_l)*(ppoly_r_e(i0,1)-ppoly_r_e(i0-1,2)) < 0.)
then
485 write(0,
'(a,i4,5(x,a,1pe24.16))')
'Non-monotonic edges at',i0,
'u(i0-1)=',u_l,
'u(i0)=',u_c, &
486 'right edge=',ppoly_r_e(i0-1,2),
'left edge=',ppoly_r_e(i0,1)
487 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)
488 problem_detected = .true.
491 if (problem_detected)
then
492 write(0,
'(a,1p9e24.16)')
'Polynomial coeffs:',ppoly_r_coefs(i0,:)
493 write(0,
'(3(a,1pe24.16,x))')
'u_l=',u_l,
'u_c=',u_c,
'u_r=',u_r
494 write(0,
'(a4,10a24)')
'i0',
'h0(i0)',
'u0(i0)',
'left edge',
'right edge',
'Polynomial coefficients'
496 write(0,
'(i4,1p10e24.16)') n,h0(n),u0(n),ppoly_r_e(n,1),ppoly_r_e(n,2),ppoly_r_coefs(n,:)
498 call mom_error(fatal,
'MOM_remapping, check_reconstructions_1d: '// &
499 'Edge values or polynomial coefficients were inconsistent!')
503 end subroutine check_reconstructions_1d
508 subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, &
509 force_bounds_in_subcell, u1, uh_err, ah_sub, aisub_src, aiss, aise )
510 integer,
intent(in) :: n0
511 real,
intent(in) :: h0(n0)
512 real,
intent(in) :: u0(n0)
513 real,
intent(in) :: ppoly0_E(n0,2)
514 real,
intent(in) :: ppoly0_coefs(:,:)
515 integer,
intent(in) :: n1
516 real,
intent(in) :: h1(n1)
517 integer,
intent(in) :: method
518 logical,
intent(in) :: force_bounds_in_subcell
519 real,
intent(out) :: u1(n1)
520 real,
intent(out) :: uh_err
521 real,
optional,
intent(out) :: ah_sub(n0+n1+1)
522 integer,
optional,
intent(out) :: aisub_src(n0+n1+1)
523 integer,
optional,
intent(out) :: aiss(n0)
524 integer,
optional,
intent(out) :: aise(n0)
533 real,
dimension(n0+n1+1) :: h_sub
534 real,
dimension(n0+n1+1) :: uh_sub
535 real,
dimension(n0+n1+1) :: u_sub
536 integer,
dimension(n0+n1+1) :: isub_src
537 integer,
dimension(n0) :: isrc_start
538 integer,
dimension(n0) :: isrc_end
539 integer,
dimension(n0) :: isrc_max
540 real,
dimension(n0) :: h0_eff
541 real,
dimension(n0) :: u0_min
542 real,
dimension(n0) :: u0_max
543 integer,
dimension(n1) :: itgt_start
544 integer,
dimension(n1) :: itgt_end
546 real :: h0_supply, h1_supply
551 logical,
parameter :: force_bounds_in_target = .true.
552 logical,
parameter :: adjust_thickest_subcell = .true.
553 logical,
parameter :: debug_bounds = .false.
554 integer :: k, i0_last_thick_cell
555 real :: h0tot, h0err, h1tot, h1err, h2tot, h2err, u02_err
556 real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, u2tot, u2err, u2min, u2max, u_orig
557 logical :: src_has_volume
558 logical :: tgt_has_volume
560 if (old_algorithm) isrc_max(:)=1
562 i0_last_thick_cell = 0
564 u0_min(i0) = min(ppoly0_e(i0,1), ppoly0_e(i0,2))
565 u0_max(i0) = max(ppoly0_e(i0,1), ppoly0_e(i0,2))
566 if (h0(i0)>0.) i0_last_thick_cell = i0
572 src_has_volume = .true.
573 tgt_has_volume = .true.
575 i_start0 = 1 ; i_start1 = 1
588 do i_sub = 2, n0+n1+1
592 dh = min(h0_supply, h1_supply)
597 dh0_eff = dh0_eff + min(dh, h0_supply)
606 if (dh >= dh_max)
then
613 if (h0_supply <= h1_supply .and. src_has_volume)
then
616 h1_supply = h1_supply - dh
619 isrc_start(i0) = i_start0
629 if (old_algorithm)
then
630 if (i0 < i0_last_thick_cell)
then
634 do while (h0_supply==0. .and. i0<i0_last_thick_cell)
649 src_has_volume = .false.
652 elseif (h0_supply >= h1_supply .and. tgt_has_volume)
then
655 h0_supply = h0_supply - dh
658 itgt_start(i1) = i_start1
666 if (old_algorithm)
then
670 tgt_has_volume = .false.
673 elseif (src_has_volume)
then
675 h_sub(i_sub) = h0_supply
677 isrc_start(i0) = i_start0
692 src_has_volume = .false.
694 elseif (tgt_has_volume)
then
696 h_sub(i_sub) = h1_supply
698 itgt_start(i1) = i_start1
707 tgt_has_volume = .false.
710 stop
'remap_via_sub_cells: THIS SHOULD NEVER HAPPEN!'
719 u_sub(1) = ppoly0_e(1,1)
732 dh0_eff = dh0_eff + dh
733 if (h0_eff(i0)>0.)
then
734 xb = dh0_eff / h0_eff(i0)
736 u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_e, ppoly0_coefs, method, i0, xa, xb)
739 u_sub(i_sub) = u0(i0)
741 if (debug_bounds)
then
742 if (method<5 .and.(u_sub(i_sub)<u0_min(i0) .or. u_sub(i_sub)>u0_max(i0)))
then
743 write(0,*)
'Sub cell average is out of bounds',i_sub,
'method=',method
744 write(0,*)
'xa,xb: ',xa,xb
745 write(0,*)
'Edge values: ',ppoly0_e(i0,:),
'mean',u0(i0)
746 write(0,*)
'a_c: ',(u0(i0)-ppoly0_e(i0,1))+(u0(i0)-ppoly0_e(i0,2))
747 write(0,*)
'Polynomial coeffs: ',ppoly0_coefs(i0,:)
748 write(0,*)
'Bounds min=',u0_min(i0),
'max=',u0_max(i0)
749 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)
750 call mom_error( fatal,
'MOM_remapping, remap_via_sub_cells: '//&
751 'Sub-cell average is out of bounds!' )
754 if (force_bounds_in_subcell)
then
758 u_orig = u_sub(i_sub)
759 u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )
760 u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )
761 u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig )
763 uh_sub(i_sub) = dh * u_sub(i_sub)
765 if (isub_src(i_sub+1) /= i0)
then
774 u_sub(n0+n1+1) = ppoly0_e(n0,2)
775 uh_sub(n0+n1+1) = ppoly0_e(n0,2) * h_sub(n0+n1+1)
777 if (adjust_thickest_subcell)
then
781 do i0 = 1, i0_last_thick_cell
783 dh_max = h_sub(i_max)
784 if (dh_max > 0.)
then
787 do i_sub = isrc_start(i0), isrc_end(i0)
788 if (i_sub /= i_max) duh = duh + uh_sub(i_sub)
790 uh_sub(i_max) = u0(i0)*h0(i0) - duh
791 u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) )
799 if (h1(i1) > 0.)
then
801 i_sub = itgt_start(i1)
802 if (force_bounds_in_target)
then
806 do i_sub = itgt_start(i1), itgt_end(i1)
807 if (force_bounds_in_target)
then
808 u1min = min(u1min, u_sub(i_sub))
809 u1max = max(u1max, u_sub(i_sub))
811 dh = dh + h_sub(i_sub)
812 duh = duh + uh_sub(i_sub)
814 uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh)
818 uh_err = uh_err + abs(duh)*epsilon(duh)
819 if (force_bounds_in_target)
then
821 u1(i1) = max(u1min, min(u1max, u1(i1)))
823 uh_err = uh_err + dh*abs( u1(i1)-u_orig )
826 u1(i1) = u_sub(itgt_start(i1))
831 if (debug_bounds)
then
832 call measure_input_bounds( n0, h0, u0, ppoly0_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
833 call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
834 call measure_output_bounds( n0+n1+1, h_sub, u_sub, h2tot, h2err, u2tot, u2err, u2min, u2max )
836 if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err+u02_err .and. abs(h1tot-h0tot)<h0err+h1err) &
837 .or. (abs(u2tot-u0tot)>u0err+u2err+u02_err .and. abs(h2tot-h0tot)<h0err+h2err) &
838 .or. (u1min<u0min .or. u1max>u0max) )
then
839 write(0,*)
'method = ',method
840 write(0,*)
'Source to sub-cells:'
841 write(0,*)
'H: h0tot=',h0tot,
'h2tot=',h2tot,
'dh=',h2tot-h0tot,
'h0err=',h0err,
'h2err=',h2err
842 if (abs(h2tot-h0tot)>h0err+h2err) &
843 write(0,*)
'H non-conservation difference=',h2tot-h0tot,
'allowed err=',h0err+h2err,
' <-----!'
844 write(0,*)
'UH: u0tot=',u0tot,
'u2tot=',u2tot,
'duh=',u2tot-u0tot,
'u0err=',u0err,
'u2err=',u2err,&
845 'adjustment err=',u02_err
846 if (abs(u2tot-u0tot)>u0err+u2err) &
847 write(0,*)
'U non-conservation difference=',u2tot-u0tot,
'allowed err=',u0err+u2err,
' <-----!'
848 write(0,*)
'Sub-cells to target:'
849 write(0,*)
'H: h2tot=',h2tot,
'h1tot=',h1tot,
'dh=',h1tot-h2tot,
'h2err=',h2err,
'h1err=',h1err
850 if (abs(h1tot-h2tot)>h2err+h1err) &
851 write(0,*)
'H non-conservation difference=',h1tot-h2tot,
'allowed err=',h2err+h1err,
' <-----!'
852 write(0,*)
'UH: u2tot=',u2tot,
'u1tot=',u1tot,
'duh=',u1tot-u2tot,
'u2err=',u2err,
'u1err=',u1err,
'uh_err=',uh_err
853 if (abs(u1tot-u2tot)>u2err+u1err) &
854 write(0,*)
'U non-conservation difference=',u1tot-u2tot,
'allowed err=',u2err+u1err,
' <-----!'
855 write(0,*)
'Source to target:'
856 write(0,*)
'H: h0tot=',h0tot,
'h1tot=',h1tot,
'dh=',h1tot-h0tot,
'h0err=',h0err,
'h1err=',h1err
857 if (abs(h1tot-h0tot)>h0err+h1err) &
858 write(0,*)
'H non-conservation difference=',h1tot-h0tot,
'allowed err=',h0err+h1err,
' <-----!'
859 write(0,*)
'UH: u0tot=',u0tot,
'u1tot=',u1tot,
'duh=',u1tot-u0tot,
'u0err=',u0err,
'u1err=',u1err,
'uh_err=',uh_err
860 if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
861 write(0,*)
'U non-conservation difference=',u1tot-u0tot,
'allowed err=',u0err+u1err+uh_err,
' <-----!'
862 write(0,*)
'U: u0min=',u0min,
'u1min=',u1min,
'u2min=',u2min
863 if (u1min<u0min)
write(0,*)
'U minimum overshoot=',u1min-u0min,
' <-----!'
864 if (u2min<u0min)
write(0,*)
'U2 minimum overshoot=',u2min-u0min,
' <-----!'
865 write(0,*)
'U: u0max=',u0max,
'u1max=',u1max,
'u2max=',u2max
866 if (u1max>u0max)
write(0,*)
'U maximum overshoot=',u1max-u0max,
' <-----!'
867 if (u2max>u0max)
write(0,*)
'U2 maximum overshoot=',u2max-u0max,
' <-----!'
868 write(0,
'(a3,6a24,2a3)')
'k',
'h0',
'left edge',
'u0',
'right edge',
'h1',
'u1',
'is',
'ie'
870 if (k<=min(n0,n1))
then
871 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)
873 write(0,
'(i3,96x,1p2e24.16,2i3)') k,h1(k),u1(k),itgt_start(k),itgt_end(k)
875 write(0,
'(i3,1p4e24.16)') k,h0(k),ppoly0_e(k,1),u0(k),ppoly0_e(k,2)
878 write(0,
'(a3,2a24)')
'k',
'u0',
'Polynomial coefficients'
880 write(0,
'(i3,1p6e24.16)') k,u0(k),ppoly0_coefs(k,:)
882 write(0,
'(a3,3a24,a3,2a24)')
'k',
'Sub-cell h',
'Sub-cell u',
'Sub-cell hu',
'i0',
'xa',
'xb'
888 dh0_eff = dh0_eff + dh
889 xb = dh0_eff / h0_eff(i0)
891 write(0,
'(i3,1p3e24.16,i3,1p2e24.16)') k,h_sub(k),u_sub(k),uh_sub(k),i0,xa,xb
893 if (isub_src(k+1) /= i0)
then
894 dh0_eff = 0.; xa = 0.
900 call mom_error( fatal,
'MOM_remapping, remap_via_sub_cells: '//&
901 'Remapping result is inconsistent!' )
907 uh_err = uh_err + u02_err
909 if (
present(ah_sub)) ah_sub(1:n0+n1+1) = h_sub(1:n0+n1+1)
910 if (
present(aisub_src)) aisub_src(1:n0+n1+1) = isub_src(1:n0+n1+1)
911 if (
present(aiss)) aiss(1:n0) = isrc_start(1:n0)
912 if (
present(aise)) aise(1:n0) = isrc_end(1:n0)
914 end subroutine remap_via_sub_cells
919 real function average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb)
920 integer,
intent(in) :: n0
921 real,
intent(in) :: u0(:)
922 real,
intent(in) :: ppoly0_e(:,:)
923 real,
intent(in) :: ppoly0_coefs(:,:)
924 integer,
intent(in) :: method
925 integer,
intent(in) :: i0
926 real,
intent(in) :: xa
927 real,
intent(in) :: xb
929 real :: u_ave, xa_2, xb_2, xa2pxb2, xapxb
930 real,
parameter :: r_3 = 1.0/3.0
932 real :: mx, a_l, a_r, u_c, ya, yb, my, xa2b2ab, ya2b2ab, a_c
935 select case ( method )
936 case ( integration_pcm )
938 case ( integration_plm )
941 + ppoly0_coefs(i0,2) * 0.5 * ( xb + xa ) )
942 case ( integration_ppm )
943 mx = 0.5 * ( xa + xb )
944 a_l = ppoly0_e(i0, 1)
945 a_r = ppoly0_e(i0, 2)
947 a_c = 0.5 * ( ( u_c - a_l ) + ( u_c - a_r ) )
950 xa2b2ab = (xa*xa+xb*xb)+xa*xb
951 u_ave = a_l + ( ( a_r - a_l ) * mx &
952 + a_c * ( 3. * ( xb + xa ) - 2.*xa2b2ab ) )
957 my = 0.5 * ( ya + yb )
958 ya2b2ab = (ya*ya+yb*yb)+ya*yb
959 u_ave = a_r + ( ( a_l - a_r ) * my &
960 + a_c * ( 3. * ( yb + ya ) - 2.*ya2b2ab ) )
962 case ( integration_pqm )
965 xa2pxb2 = xa_2 + xb_2
969 + ( ppoly0_coefs(i0,2) * 0.5 * ( xapxb ) &
970 + ( ppoly0_coefs(i0,3) * r_3 * ( xa2pxb2 + xa*xb ) &
971 + ( ppoly0_coefs(i0,4) * 0.25* ( xa2pxb2 * xapxb ) &
972 + ppoly0_coefs(i0,5) * 0.2 * ( ( xb*xb_2 + xa*xa_2 ) * xapxb + xa_2*xb_2 ) ) ) ) )
974 call mom_error( fatal,
'The selected integration method is invalid' )
977 select case ( method )
978 case ( integration_pcm )
979 u_ave = ppoly0_coefs(i0,1)
980 case ( integration_plm )
983 a_l = ppoly0_e(i0, 1)
984 a_r = ppoly0_e(i0, 2)
987 u_ave = a_l + xa * ( a_r - a_l )
989 u_ave = a_r + ya * ( a_l - a_r )
991 case ( integration_ppm )
995 a_l = ppoly0_e(i0, 1)
996 a_r = ppoly0_e(i0, 2)
998 a_c = 3. * ( ( u_c - a_l ) + ( u_c - a_r ) )
1001 u_ave = a_l + xa * ( ( a_r - a_l ) + a_c * ya )
1003 u_ave = a_r + ya * ( ( a_l - a_r ) + a_c * xa )
1005 case ( integration_pqm )
1006 u_ave = ppoly0_coefs(i0,1) &
1007 + xa * ( ppoly0_coefs(i0,2) &
1008 + xa * ( ppoly0_coefs(i0,3) &
1009 + xa * ( ppoly0_coefs(i0,4) &
1010 + xa * ppoly0_coefs(i0,5) ) ) )
1012 call mom_error( fatal,
'The selected integration method is invalid' )
1015 average_value_ppoly = u_ave
1017 end function average_value_ppoly
1020 subroutine measure_input_bounds( n0, h0, u0, ppoly_E, h0tot, h0err, u0tot, u0err, u0min, u0max )
1021 integer,
intent(in) :: n0
1022 real,
dimension(n0),
intent(in) :: h0
1023 real,
dimension(n0),
intent(in) :: u0
1024 real,
dimension(n0,2),
intent(in) :: ppoly_E
1025 real,
intent(out) :: h0tot
1026 real,
intent(out) :: h0err
1027 real,
intent(out) :: u0tot
1028 real,
intent(out) :: u0err
1029 real,
intent(out) :: u0min
1030 real,
intent(out) :: u0max
1035 eps = epsilon(h0(1))
1038 u0tot = h0(1) * u0(1)
1040 u0min = min( ppoly_e(1,1), ppoly_e(1,2) )
1041 u0max = max( ppoly_e(1,1), ppoly_e(1,2) )
1043 h0tot = h0tot + h0(k)
1044 h0err = h0err + eps * max(h0tot, h0(k))
1045 u0tot = u0tot + h0(k) * u0(k)
1046 u0err = u0err + eps * max(abs(u0tot), abs(h0(k) * u0(k)))
1047 u0min = min( u0min, ppoly_e(k,1), ppoly_e(k,2) )
1048 u0max = max( u0max, ppoly_e(k,1), ppoly_e(k,2) )
1051 end subroutine measure_input_bounds
1054 subroutine measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
1055 integer,
intent(in) :: n1
1056 real,
dimension(n1),
intent(in) :: h1
1057 real,
dimension(n1),
intent(in) :: u1
1058 real,
intent(out) :: h1tot
1059 real,
intent(out) :: h1err
1060 real,
intent(out) :: u1tot
1061 real,
intent(out) :: u1err
1062 real,
intent(out) :: u1min
1063 real,
intent(out) :: u1max
1068 eps = epsilon(h1(1))
1071 u1tot = h1(1) * u1(1)
1076 h1tot = h1tot + h1(k)
1077 h1err = h1err + eps * max(h1tot, h1(k))
1078 u1tot = u1tot + h1(k) * u1(k)
1079 u1err = u1err + eps * max(abs(u1tot), abs(h1(k) * u1(k)))
1080 u1min = min(u1min, u1(k))
1081 u1max = max(u1max, u1(k))
1084 end subroutine measure_output_bounds
1088 subroutine remapbyprojection( n0, h0, u0, ppoly0_E, ppoly0_coefs, &
1089 n1, h1, method, u1, h_neglect )
1090 integer,
intent(in) :: n0
1091 real,
intent(in) :: h0(:)
1092 real,
intent(in) :: u0(:)
1093 real,
intent(in) :: ppoly0_E(:,:)
1094 real,
intent(in) :: ppoly0_coefs(:,:)
1095 integer,
intent(in) :: n1
1096 real,
intent(in) :: h1(:)
1097 integer,
intent(in) :: method
1098 real,
intent(out) :: u1(:)
1099 real,
optional,
intent(in) :: h_neglect
1117 xr = xl + h1(itarget)
1119 call integraterecononinterval( n0, h0, u0, ppoly0_e, ppoly0_coefs, method, &
1120 xl, xr, h1(itarget), u1(itarget), jstart, xstart, h_neglect )
1124 end subroutine remapbyprojection
1136 subroutine remapbydeltaz( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, &
1137 method, u1, h1, h_neglect )
1138 integer,
intent(in) :: n0
1139 real,
dimension(:),
intent(in) :: h0
1140 real,
dimension(:),
intent(in) :: u0
1141 real,
dimension(:,:),
intent(in) :: ppoly0_E
1142 real,
dimension(:,:),
intent(in) :: ppoly0_coefs
1143 integer,
intent(in) :: n1
1144 real,
dimension(:),
intent(in) :: dx1
1145 integer,
intent(in) :: method
1146 real,
dimension(:),
intent(out) :: u1
1147 real,
dimension(:), &
1148 optional,
intent(out) :: h1
1149 real,
optional,
intent(in) :: h_neglect
1155 real :: xOld, hOld, uOld
1156 real :: xNew, hNew, h_err
1157 real :: uhNew, hFlux, uAve, fluxL, fluxR
1173 if (itarget == 0)
then
1177 elseif (itarget <= n0)
then
1178 xold = xold + h0(itarget)
1181 h_err = h_err + epsilon(hold) * max(hold, xold)
1186 xnew = xold + dx1(itarget+1)
1187 xl = min( xold, xnew )
1188 xr = max( xold, xnew )
1191 hflux = abs(dx1(itarget+1))
1192 call integraterecononinterval( n0, h0, u0, ppoly0_e, ppoly0_coefs, method, &
1193 xl, xr, hflux, uave, jstart, xstart )
1195 fluxr = dx1(itarget+1)*uave
1198 hnew = hold + ( dx1(itarget+1) - dx1(itarget) )
1199 hnew = max( 0., hnew )
1200 uhnew = ( uold * hold ) + ( fluxr - fluxl )
1202 u1(itarget) = uhnew / hnew
1206 if (
present(h1)) h1(itarget) = hnew
1211 end subroutine remapbydeltaz
1215 subroutine integraterecononinterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, &
1216 xL, xR, hC, uAve, jStart, xStart, h_neglect )
1217 integer,
intent(in) :: n0
1218 real,
dimension(:),
intent(in) :: h0
1219 real,
dimension(:),
intent(in) :: u0
1220 real,
dimension(:,:),
intent(in) :: ppoly0_E
1221 real,
dimension(:,:),
intent(in) :: ppoly0_coefs
1222 integer,
intent(in) :: method
1223 real,
intent(in) :: xL
1224 real,
intent(in) :: xR
1225 real,
intent(in) :: hC
1226 real,
intent(out) :: uAve
1227 integer,
intent(inout) :: jStart
1229 real,
intent(inout) :: xStart
1231 real,
optional,
intent(in) :: h_neglect
1241 real :: x0jLl, x0jLr
1242 real :: x0jRl, x0jRr
1245 real :: x0_2, x1_2, x02px12, x0px1
1247 real,
parameter :: r_3 = 1.0/3.0
1249 hneglect = hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
1260 x0jlr = x0jll + h0(j)
1262 if ( ( xl >= x0jll ) .AND. ( xl <= x0jlr ) )
then
1278 if ( jl == -1 )
call mom_error(fatal, &
1279 'MOM_remapping, integrateReconOnInterval: '//&
1280 'The location of the left-most cell could not be found')
1292 if ( abs(xr - xl) == 0.0 )
then
1298 if ( h0(jl) == 0.0 )
then
1299 uave = 0.5 * ( ppoly0_e(jl,1) + ppoly0_e(jl,2) )
1302 xi0 = xl / ( h0(jl) + hneglect ) - x0jll / ( h0(jl) + hneglect )
1304 select case ( method )
1305 case ( integration_pcm )
1306 uave = ppoly0_coefs(jl,1)
1307 case ( integration_plm )
1308 uave = ppoly0_coefs(jl,1) &
1309 + xi0 * ppoly0_coefs(jl,2)
1310 case ( integration_ppm )
1311 uave = ppoly0_coefs(jl,1) &
1312 + xi0 * ( ppoly0_coefs(jl,2) &
1313 + xi0 * ppoly0_coefs(jl,3) )
1314 case ( integration_pqm )
1315 uave = ppoly0_coefs(jl,1) &
1316 + xi0 * ( ppoly0_coefs(jl,2) &
1317 + xi0 * ( ppoly0_coefs(jl,3) &
1318 + xi0 * ( ppoly0_coefs(jl,4) &
1319 + xi0 * ppoly0_coefs(jl,5) ) ) )
1321 call mom_error( fatal,
'The selected integration method is invalid' )
1334 x0jrr = x0jrl + h0(j)
1336 if ( ( xr >= x0jrl ) .AND. ( xr <= x0jrr ) )
then
1345 if (xr>x0jrr) jr = n0
1351 if ( jl == jr )
then
1360 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1361 xi0 = max( 0., min( 1., ( xl - x0jll ) / ( h0(jl) + hneglect ) ) )
1362 xi1 = max( 0., min( 1., ( xr - x0jll ) / ( h0(jl) + hneglect ) ) )
1364 xi0 = xl / h0(jl) - x0jll / ( h0(jl) + hneglect )
1365 xi1 = xr / h0(jl) - x0jll / ( h0(jl) + hneglect )
1368 hact = h0(jl) * ( xi1 - xi0 )
1373 select case ( method )
1374 case ( integration_pcm )
1375 q = ( xr - xl ) * ppoly0_coefs(jl,1)
1376 case ( integration_plm )
1377 q = ( xr - xl ) * ( &
1378 ppoly0_coefs(jl,1) &
1379 + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
1380 case ( integration_ppm )
1381 q = ( xr - xl ) * ( &
1382 ppoly0_coefs(jl,1) &
1383 + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1384 + ppoly0_coefs(jl,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
1385 case ( integration_pqm )
1388 x02px12 = x0_2 + x1_2
1390 q = ( xr - xl ) * ( &
1391 ppoly0_coefs(jl,1) &
1392 + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1393 + ( ppoly0_coefs(jl,3) * r_3 * ( x02px12 + xi0*xi1 ) &
1394 + ppoly0_coefs(jl,4) * 0.25* ( x02px12 * x0px1 ) &
1395 + ppoly0_coefs(jl,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
1397 call mom_error( fatal,
'The selected integration method is invalid' )
1416 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1417 xi0 = max( 0., min( 1., ( xl - x0jll ) / ( h0(jl) + hneglect ) ) )
1419 xi0 = (xl - x0jll) / ( h0(jl) + hneglect )
1423 hact = h0(jl) * ( xi1 - xi0 )
1425 select case ( method )
1426 case ( integration_pcm )
1427 q = q + ( x0jlr - xl ) * ppoly0_coefs(jl,1)
1428 case ( integration_plm )
1429 q = q + ( x0jlr - xl ) * ( &
1430 ppoly0_coefs(jl,1) &
1431 + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
1432 case ( integration_ppm )
1433 q = q + ( x0jlr - xl ) * ( &
1434 ppoly0_coefs(jl,1) &
1435 + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1436 + ppoly0_coefs(jl,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
1437 case ( integration_pqm )
1440 x02px12 = x0_2 + x1_2
1442 q = q + ( x0jlr - xl ) * ( &
1443 ppoly0_coefs(jl,1) &
1444 + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
1445 + ( ppoly0_coefs(jl,3) * r_3 * ( x02px12 + xi0*xi1 ) &
1446 + ppoly0_coefs(jl,4) * 0.25* ( x02px12 * x0px1 ) &
1447 + ppoly0_coefs(jl,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
1449 call mom_error( fatal,
'The selected integration method is invalid' )
1453 if ( jr > (jl+1) )
then
1455 q = q + h0(k) * u0(k)
1462 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1463 xi1 = max( 0., min( 1., ( xr - x0jrl ) / ( h0(jr) + hneglect ) ) )
1465 xi1 = (xr - x0jrl) / ( h0(jr) + hneglect )
1468 hact = hact + h0(jr) * ( xi1 - xi0 )
1470 select case ( method )
1471 case ( integration_pcm )
1472 q = q + ( xr - x0jrl ) * ppoly0_coefs(jr,1)
1473 case ( integration_plm )
1474 q = q + ( xr - x0jrl ) * ( &
1475 ppoly0_coefs(jr,1) &
1476 + ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) )
1477 case ( integration_ppm )
1478 q = q + ( xr - x0jrl ) * ( &
1479 ppoly0_coefs(jr,1) &
1480 + ( ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) &
1481 + ppoly0_coefs(jr,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
1482 case ( integration_pqm )
1485 x02px12 = x0_2 + x1_2
1487 q = q + ( xr - x0jrl ) * ( &
1488 ppoly0_coefs(jr,1) &
1489 + ( ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) &
1490 + ( ppoly0_coefs(jr,3) * r_3 * ( x02px12 + xi0*xi1 ) &
1491 + ppoly0_coefs(jr,4) * 0.25* ( x02px12 * x0px1 ) &
1492 + ppoly0_coefs(jr,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
1494 call mom_error( fatal,
'The selected integration method is invalid' )
1500 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
1502 uave = ppoly0_coefs(jl,1)
1512 end subroutine integraterecononinterval
1515 subroutine dzfromh1h2( n1, h1, n2, h2, dx )
1516 integer,
intent(in) :: n1
1517 real,
dimension(:),
intent(in) :: h1
1518 integer,
intent(in) :: n2
1519 real,
dimension(:),
intent(in) :: h2
1520 real,
dimension(:),
intent(out) :: dx
1528 do k = 1, max(n1,n2)
1529 if (k <= n1) x1 = x1 + h1(k)
1536 end subroutine dzfromh1h2
1539 subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, &
1540 check_reconstruction, check_remapping, force_bounds_in_subcell)
1543 character(len=*),
intent(in) :: remapping_scheme
1544 logical,
optional,
intent(in) :: boundary_extrapolation
1545 logical,
optional,
intent(in) :: check_reconstruction
1546 logical,
optional,
intent(in) :: check_remapping
1547 logical,
optional,
intent(in) :: force_bounds_in_subcell
1550 call remapping_set_param(cs, remapping_scheme=remapping_scheme, boundary_extrapolation=boundary_extrapolation, &
1551 check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
1552 force_bounds_in_subcell=force_bounds_in_subcell)
1554 end subroutine initialize_remapping
1560 subroutine setreconstructiontype(string,CS)
1561 character(len=*),
intent(in) :: string
1566 select case ( uppercase(trim(string)) )
1568 cs%remapping_scheme = remapping_pcm
1571 cs%remapping_scheme = remapping_plm
1574 cs%remapping_scheme = remapping_ppm_h4
1577 cs%remapping_scheme = remapping_ppm_ih4
1580 cs%remapping_scheme = remapping_pqm_ih4ih3
1583 cs%remapping_scheme = remapping_pqm_ih6ih5
1586 call mom_error(fatal,
"setReconstructionType: "//&
1587 "Unrecognized choice for REMAPPING_SCHEME ("//trim(string)//
").")
1592 end subroutine setreconstructiontype
1595 subroutine end_remapping(CS)
1600 end subroutine end_remapping
1605 logical function remapping_unit_tests(verbose)
1606 logical,
intent(in) :: verbose
1608 integer,
parameter :: n0 = 4, n1 = 3, n2 = 6
1609 real :: h0(n0), x0(n0+1), u0(n0)
1610 real :: h1(n1), x1(n1+1), u1(n1), hn1(n1), dx1(n1+1)
1611 real :: h2(n2), x2(n2+1), u2(n2), hn2(n2), dx2(n2+1)
1612 data u0 /9., 3., -3., -9./
1617 real,
allocatable,
dimension(:,:) :: ppoly0_e, ppoly0_s, ppoly0_coefs
1619 real :: err, h_neglect, h_neglect_edge
1620 logical :: thistest, v
1623 h_neglect = hneglect_dflt
1624 h_neglect_edge = 1.0e-10
1626 write(*,*)
'==== MOM_remapping: remapping_unit_tests ================='
1627 remapping_unit_tests = .false.
1630 call buildgridfromh(n0, h0, x0)
1632 err=x0(i)-0.75*real(i-1)
1633 if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
1635 if (thistest)
write(*,*)
'remapping_unit_tests: Failed buildGridFromH() 1'
1636 remapping_unit_tests = remapping_unit_tests .or. thistest
1637 call buildgridfromh(n1, h1, x1)
1640 if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
1642 if (thistest)
write(*,*)
'remapping_unit_tests: Failed buildGridFromH() 2'
1643 remapping_unit_tests = remapping_unit_tests .or. thistest
1646 call initialize_remapping(cs,
'PPM_H4')
1647 if (verbose)
write(*,*)
'h0 (test data)'
1648 if (verbose)
call dumpgrid(n0,h0,x0,u0)
1650 call dzfromh1h2( n0, h0, n1, h1, dx1 )
1651 call remapping_core_w( cs, n0, h0, u0, n1, dx1, u1, h_neglect, h_neglect_edge)
1653 err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1654 if (abs(err)>real(n1-1)*epsilon(err)) thistest = .true.
1656 if (verbose)
write(*,*)
'h1 (by projection)'
1657 if (verbose)
call dumpgrid(n1,h1,x1,u1)
1658 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remapping_core_w()'
1659 remapping_unit_tests = remapping_unit_tests .or. thistest
1662 allocate(ppoly0_e(n0,2))
1663 allocate(ppoly0_s(n0,2))
1664 allocate(ppoly0_coefs(n0,cs%degree+1))
1668 ppoly0_coefs(:,:) = 0.0
1670 call edge_values_explicit_h4( n0, h0, u0, ppoly0_e, h_neglect=1e-10 )
1671 call ppm_reconstruction( n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )
1672 call ppm_boundary_extrapolation( n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )
1674 call remapbyprojection( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1675 n1, h1, integration_ppm, u1, h_neglect )
1677 err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1678 if (abs(err)>2.*epsilon(err)) thistest = .true.
1680 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remapByProjection()'
1681 remapping_unit_tests = remapping_unit_tests .or. thistest
1685 call remapbydeltaz( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1686 n1, x1-x0(1:n1+1), &
1687 integration_ppm, u1, hn1, h_neglect )
1688 if (verbose)
write(*,*)
'h1 (by delta)'
1689 if (verbose)
call dumpgrid(n1,h1,x1,u1)
1692 err=u1(i)-8.*(0.5*real(1+n1)-real(i))
1693 if (abs(err)>2.*epsilon(err)) thistest = .true.
1695 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remapByDeltaZ() 1'
1696 remapping_unit_tests = remapping_unit_tests .or. thistest
1699 call buildgridfromh(n2, h2, x2)
1700 dx2(1:n0+1) = x2(1:n0+1) - x0
1701 dx2(n0+2:n2+1) = x2(n0+2:n2+1) - x0(n0+1)
1702 call remapbydeltaz( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1704 integration_ppm, u2, hn2, h_neglect )
1705 if (verbose)
write(*,*)
'h2'
1706 if (verbose)
call dumpgrid(n2,h2,x2,u2)
1707 if (verbose)
write(*,*)
'hn2'
1708 if (verbose)
call dumpgrid(n2,hn2,x2,u2)
1711 err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
1712 if (abs(err)>2.*epsilon(err)) thistest = .true.
1714 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remapByDeltaZ() 2'
1715 remapping_unit_tests = remapping_unit_tests .or. thistest
1717 if (verbose)
write(*,*)
'Via sub-cells'
1719 call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1720 n2, h2, integration_ppm, .false., u2, err )
1721 if (verbose)
call dumpgrid(n2,h2,x2,u2)
1724 err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
1725 if (abs(err)>2.*epsilon(err)) thistest = .true.
1727 if (thistest)
write(*,*)
'remapping_unit_tests: Failed remap_via_sub_cells() 2'
1728 remapping_unit_tests = remapping_unit_tests .or. thistest
1730 call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1731 6, (/.125,.125,.125,.125,.125,.125/), integration_ppm, .false., u2, err )
1732 if (verbose)
call dumpgrid(6,h2,x2,u2)
1734 call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
1735 3, (/2.25,1.5,1./), integration_ppm, .false., u2, err )
1736 if (verbose)
call dumpgrid(3,h2,x2,u2)
1738 if (.not. remapping_unit_tests)
write(*,*)
'Pass'
1740 write(*,*)
'===== MOM_remapping: new remapping_unit_tests =================='
1742 deallocate(ppoly0_e, ppoly0_s, ppoly0_coefs)
1743 allocate(ppoly0_coefs(5,6))
1744 allocate(ppoly0_e(5,2))
1745 allocate(ppoly0_s(5,2))
1747 call pcm_reconstruction(3, (/1.,2.,4./), ppoly0_e(1:3,:), &
1748 ppoly0_coefs(1:3,:) )
1749 remapping_unit_tests = remapping_unit_tests .or. &
1750 test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,4./),
'PCM: left edges')
1751 remapping_unit_tests = remapping_unit_tests .or. &
1752 test_answer(v, 3, ppoly0_e(:,2), (/1.,2.,4./),
'PCM: right edges')
1753 remapping_unit_tests = remapping_unit_tests .or. &
1754 test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,4./),
'PCM: P0')
1756 call plm_reconstruction(3, (/1.,1.,1./), (/1.,3.,5./), ppoly0_e(1:3,:), &
1757 ppoly0_coefs(1:3,:), h_neglect )
1758 remapping_unit_tests = remapping_unit_tests .or. &
1759 test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,5./),
'Unlim PLM: left edges')
1760 remapping_unit_tests = remapping_unit_tests .or. &
1761 test_answer(v, 3, ppoly0_e(:,2), (/1.,4.,5./),
'Unlim PLM: right edges')
1762 remapping_unit_tests = remapping_unit_tests .or. &
1763 test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,5./),
'Unlim PLM: P0')
1764 remapping_unit_tests = remapping_unit_tests .or. &
1765 test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./),
'Unlim PLM: P1')
1767 call plm_reconstruction(3, (/1.,1.,1./), (/1.,2.,7./), ppoly0_e(1:3,:), &
1768 ppoly0_coefs(1:3,:), h_neglect )
1769 remapping_unit_tests = remapping_unit_tests .or. &
1770 test_answer(v, 3, ppoly0_e(:,1), (/1.,1.,7./),
'Left lim PLM: left edges')
1771 remapping_unit_tests = remapping_unit_tests .or. &
1772 test_answer(v, 3, ppoly0_e(:,2), (/1.,3.,7./),
'Left lim PLM: right edges')
1773 remapping_unit_tests = remapping_unit_tests .or. &
1774 test_answer(v, 3, ppoly0_coefs(:,1), (/1.,1.,7./),
'Left lim PLM: P0')
1775 remapping_unit_tests = remapping_unit_tests .or. &
1776 test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./),
'Left lim PLM: P1')
1778 call plm_reconstruction(3, (/1.,1.,1./), (/1.,6.,7./), ppoly0_e(1:3,:), &
1779 ppoly0_coefs(1:3,:), h_neglect )
1780 remapping_unit_tests = remapping_unit_tests .or. &
1781 test_answer(v, 3, ppoly0_e(:,1), (/1.,5.,7./),
'Right lim PLM: left edges')
1782 remapping_unit_tests = remapping_unit_tests .or. &
1783 test_answer(v, 3, ppoly0_e(:,2), (/1.,7.,7./),
'Right lim PLM: right edges')
1784 remapping_unit_tests = remapping_unit_tests .or. &
1785 test_answer(v, 3, ppoly0_coefs(:,1), (/1.,5.,7./),
'Right lim PLM: P0')
1786 remapping_unit_tests = remapping_unit_tests .or. &
1787 test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./),
'Right lim PLM: P1')
1789 call plm_reconstruction(3, (/1.,2.,3./), (/1.,4.,9./), ppoly0_e(1:3,:), &
1790 ppoly0_coefs(1:3,:), h_neglect )
1791 remapping_unit_tests = remapping_unit_tests .or. &
1792 test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,9./),
'Non-uniform line PLM: left edges')
1793 remapping_unit_tests = remapping_unit_tests .or. &
1794 test_answer(v, 3, ppoly0_e(:,2), (/1.,6.,9./),
'Non-uniform line PLM: right edges')
1795 remapping_unit_tests = remapping_unit_tests .or. &
1796 test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,9./),
'Non-uniform line PLM: P0')
1797 remapping_unit_tests = remapping_unit_tests .or. &
1798 test_answer(v, 3, ppoly0_coefs(:,2), (/0.,4.,0./),
'Non-uniform line PLM: P1')
1800 call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_e, &
1803 thistest = test_answer(v, 5, ppoly0_e(:,1), (/0.,2.,4.,6.,8./),
'Line H4: left edges')
1804 thistest = test_answer(v, 5, ppoly0_e(:,2), (/2.,4.,6.,8.,10./),
'Line H4: right edges')
1805 ppoly0_e(:,1) = (/0.,2.,4.,6.,8./)
1806 ppoly0_e(:,2) = (/2.,4.,6.,8.,10./)
1807 call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_e(1:5,:), &
1808 ppoly0_coefs(1:5,:), h_neglect )
1809 remapping_unit_tests = remapping_unit_tests .or. &
1810 test_answer(v, 5, ppoly0_coefs(:,1), (/1.,2.,4.,6.,9./),
'Line PPM: P0')
1811 remapping_unit_tests = remapping_unit_tests .or. &
1812 test_answer(v, 5, ppoly0_coefs(:,2), (/0.,2.,2.,2.,0./),
'Line PPM: P1')
1813 remapping_unit_tests = remapping_unit_tests .or. &
1814 test_answer(v, 5, ppoly0_coefs(:,3), (/0.,0.,0.,0.,0./),
'Line PPM: P2')
1816 call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_e, &
1819 thistest = test_answer(v, 5, ppoly0_e(:,1), (/3.,0.,3.,12.,27./),
'Parabola H4: left edges')
1820 thistest = test_answer(v, 5, ppoly0_e(:,2), (/0.,3.,12.,27.,48./),
'Parabola H4: right edges')
1821 ppoly0_e(:,1) = (/0.,0.,3.,12.,27./)
1822 ppoly0_e(:,2) = (/0.,3.,12.,27.,48./)
1823 call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,1.,7.,19.,37./), ppoly0_e(1:5,:), &
1824 ppoly0_coefs(1:5,:), h_neglect )
1825 remapping_unit_tests = remapping_unit_tests .or. &
1826 test_answer(v, 5, ppoly0_e(:,1), (/0.,0.,3.,12.,37./),
'Parabola PPM: left edges')
1827 remapping_unit_tests = remapping_unit_tests .or. &
1828 test_answer(v, 5, ppoly0_e(:,2), (/0.,3.,12.,27.,37./),
'Parabola PPM: right edges')
1829 remapping_unit_tests = remapping_unit_tests .or. &
1830 test_answer(v, 5, ppoly0_coefs(:,1), (/0.,0.,3.,12.,37./),
'Parabola PPM: P0')
1831 remapping_unit_tests = remapping_unit_tests .or. &
1832 test_answer(v, 5, ppoly0_coefs(:,2), (/0.,0.,6.,12.,0./),
'Parabola PPM: P1')
1833 remapping_unit_tests = remapping_unit_tests .or. &
1834 test_answer(v, 5, ppoly0_coefs(:,3), (/0.,3.,3.,3.,0./),
'Parabola PPM: P2')
1836 ppoly0_e(:,1) = (/0.,0.,6.,10.,15./)
1837 ppoly0_e(:,2) = (/0.,6.,12.,17.,15./)
1838 call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,5.,7.,16.,15./), ppoly0_e(1:5,:), &
1839 ppoly0_coefs(1:5,:), h_neglect )
1840 remapping_unit_tests = remapping_unit_tests .or. &
1841 test_answer(v, 5, ppoly0_e(:,1), (/0.,3.,6.,16.,15./),
'Limits PPM: left edges')
1842 remapping_unit_tests = remapping_unit_tests .or. &
1843 test_answer(v, 5, ppoly0_e(:,2), (/0.,6.,9.,16.,15./),
'Limits PPM: right edges')
1844 remapping_unit_tests = remapping_unit_tests .or. &
1845 test_answer(v, 5, ppoly0_coefs(:,1), (/0.,3.,6.,16.,15./),
'Limits PPM: P0')
1846 remapping_unit_tests = remapping_unit_tests .or. &
1847 test_answer(v, 5, ppoly0_coefs(:,2), (/0.,6.,0.,0.,0./),
'Limits PPM: P1')
1848 remapping_unit_tests = remapping_unit_tests .or. &
1849 test_answer(v, 5, ppoly0_coefs(:,3), (/0.,-3.,3.,0.,0./),
'Limits PPM: P2')
1851 call plm_reconstruction(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_e(1:4,:), &
1852 ppoly0_coefs(1:4,:), h_neglect )
1853 remapping_unit_tests = remapping_unit_tests .or. &
1854 test_answer(v, 4, ppoly0_e(1:4,1), (/5.,5.,3.,1./),
'PPM: left edges h=0110')
1855 remapping_unit_tests = remapping_unit_tests .or. &
1856 test_answer(v, 4, ppoly0_e(1:4,2), (/5.,3.,1.,1./),
'PPM: right edges h=0110')
1857 call remap_via_sub_cells( 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_e(1:4,:), &
1858 ppoly0_coefs(1:4,:), &
1859 2, (/1.,1./), integration_plm, .false., u2, err )
1860 remapping_unit_tests = remapping_unit_tests .or. &
1861 test_answer(v, 2, u2, (/4.,2./),
'PLM: remapped h=0110->h=11')
1863 deallocate(ppoly0_e, ppoly0_s, ppoly0_coefs)
1865 if (.not. remapping_unit_tests)
write(*,*)
'Pass'
1867 end function remapping_unit_tests
1870 logical function test_answer(verbose, n, u, u_true, label)
1871 logical,
intent(in) :: verbose
1872 integer,
intent(in) :: n
1873 real,
dimension(n),
intent(in) :: u
1874 real,
dimension(n),
intent(in) :: u_true
1875 character(len=*),
intent(in) :: label
1879 test_answer = .false.
1881 if (u(k) /= u_true(k)) test_answer = .true.
1883 if (test_answer .or. verbose)
then
1884 write(*,
'(a4,2a24,x,a)')
'k',
'Calculated value',
'Correct value',label
1886 if (u(k) /= u_true(k))
then
1887 write(*,
'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),
' err=',u(k)-u_true(k),
' < wrong'
1889 write(*,
'(i4,1p2e24.16)') k,u(k),u_true(k)
1894 end function test_answer
1897 subroutine dumpgrid(n,h,x,u)
1898 integer,
intent(in) :: n
1899 real,
dimension(:),
intent(in) :: h
1900 real,
dimension(:),
intent(in) :: x
1901 real,
dimension(:),
intent(in) :: u
1903 write(*,
'("i=",20i10)') (i,i=1,n+1)
1904 write(*,
'("x=",20es10.2)') (x(i),i=1,n+1)
1905 write(*,
'("i=",5x,20i10)') (i,i=1,n)
1906 write(*,
'("h=",5x,20es10.2)') (h(i),i=1,n)
1907 write(*,
'("u=",5x,20es10.2)') (u(i),i=1,n)
1908 end subroutine dumpgrid