9 implicit none ;
private
22 #undef __DO_SAFETY_CHECKS__
48 integer,
intent(in) :: n
49 real,
dimension(:),
intent(in) :: h
50 real,
dimension(:),
intent(in) :: u
51 real,
dimension(:,:),
intent(inout) :: edge_val
52 real,
optional,
intent(in) :: h_neglect
59 real :: sigma_l, sigma_c, sigma_r
64 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
78 elseif ( k == n )
then
100 sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
101 sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
102 sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
104 if ( (sigma_l * sigma_r) > 0.0 )
then
105 slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
114 slope = slope * h_c * 0.5
116 if ( (u_l-u0_l)*(u0_l-u_c) < 0.0 )
then
117 u0_l = u_c - sign( min( abs(slope), abs(u0_l-u_c) ), slope )
120 if ( (u_r-u0_r)*(u0_r-u_c) < 0.0 )
then
121 u0_r = u_c + sign( min( abs(slope), abs(u0_r-u_c) ), slope )
125 u0_l = max( min( u0_l, max(u_l, u_c) ), min(u_l, u_c) )
126 u0_r = max( min( u0_r, max(u_r, u_c) ), min(u_r, u_c) )
141 integer,
intent(in) :: n
142 real,
dimension(:,:),
intent(inout) :: edge_val
154 u0_minus = edge_val(k,2)
157 u0_plus = edge_val(k+1,1)
159 if ( u0_minus /= u0_plus )
then
160 u0_avg = 0.5 * ( u0_minus + u0_plus )
161 edge_val(k,2) = u0_avg
162 edge_val(k+1,1) = u0_avg
174 integer,
intent(in) :: n
175 real,
dimension(:),
intent(in) :: u
176 real,
dimension(:,:),
intent(inout) :: edge_val
189 u0_minus = edge_val(k,2)
192 u0_plus = edge_val(k+1,1)
200 if ( (u0_plus - u0_minus)*(um_plus - um_minus) < 0.0 )
then
201 u0_avg = 0.5 * ( u0_minus + u0_plus )
202 u0_avg = max( min( u0_avg, max(um_minus, um_plus) ), min(um_minus, um_plus) )
203 edge_val(k,2) = u0_avg
204 edge_val(k+1,1) = u0_avg
226 integer,
intent(in) :: n
227 real,
dimension(:),
intent(in) :: h
228 real,
dimension(:),
intent(in) :: u
229 real,
dimension(:,:),
intent(inout) :: edge_val
230 real,
optional,
intent(in) :: h_neglect
255 edge_val(k,1) = ( u0*h1 + u1*h0 ) / ( h0 + h1 )
259 edge_val(k-1,2) = edge_val(k,1)
289 integer,
intent(in) :: n
290 real,
dimension(:),
intent(in) :: h
291 real,
dimension(:),
intent(in) :: u
292 real,
dimension(:,:),
intent(inout) :: edge_val
293 real,
optional,
intent(in) :: h_neglect
294 logical,
optional,
intent(in) :: answers_2018
298 real :: u0, u1, u2, u3
299 real :: h0, h1, h2, h3
302 real,
dimension(5) :: x
303 real,
parameter :: c1_12 = 1.0 / 12.0
305 real,
dimension(4,4) :: a
306 real,
dimension(4) :: b, c
308 logical :: use_2018_answers
310 use_2018_answers = .true. ;
if (
present(answers_2018)) use_2018_answers = answers_2018
322 if (h0+h1==0. .or. h1+h2==0. .or. h2+h3==0.)
then
323 f1 = max( hneglect, h0+h1+h2+h3 )
335 f1 = (h0+h1) * (h2+h3) / (h1+h2)
336 f2 = u1 * h2 + u2 * h1
337 f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3)
341 f1 = h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) )
342 f2 = u1*(h0+2.0*h1) - u0*h1
346 f1 = h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) )
347 f2 = u2*(2.0*h2+h3) - u3*h2
351 e = e / ( h0 + h1 + h2 + h3)
356 #ifdef __DO_SAFETY_CHECKS__
358 write(0,*)
'NaN in explicit_edge_h4 at k=',i
359 write(0,*)
'u0-u3=',u0,u1,u2,u3
360 write(0,*)
'h0-h3=',h0,h1,h2,h3
361 write(0,*)
'f1-f3=',f1,f2,f3
362 stop
'Nan during edge_values_explicit_h4'
369 f1 = max( hneglect,
hminfrac*sum(h(1:4)) )
372 x(i) = x(i-1) + max(f1, h(i-1))
377 if (use_2018_answers)
then
378 do j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ;
enddo
380 xavg = 0.5 * (x(i+1) + x(i))
383 a(i,3) = dx * (xavg**2 + c1_12*dx**2)
384 a(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
398 edge_val(2,1) = edge_val(1,2)
400 #ifdef __DO_SAFETY_CHECKS__
401 if (edge_val(1,1) /= edge_val(1,1) .or. edge_val(1,2) /= edge_val(1,2))
then
402 write(0,*)
'NaN in explicit_edge_h4 at k=',1
406 write(0,*)
'h(1:4)=',h(1:4)
408 stop
'Nan during edge_values_explicit_h4'
413 f1 = max( hneglect,
hminfrac*sum(h(n-3:n)) )
416 x(i) = x(i-1) + max(f1, h(n-5+i))
420 dx = max(f1, h(n-4+i) )
421 if (use_2018_answers)
then
422 do j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ;
enddo
424 xavg = 0.5 * (x(i+1) + x(i))
427 a(i,3) = dx * (xavg**2 + c1_12*dx**2)
428 a(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
442 edge_val(n-1,2) = edge_val(n,1)
444 #ifdef __DO_SAFETY_CHECKS__
445 if (edge_val(n,1) /= edge_val(n,1) .or. edge_val(n,2) /= edge_val(n,2))
then
446 write(0,*)
'NaN in explicit_edge_h4 at k=',n
450 a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j)
453 b(i) = u(n-4+i) * ( h(n-4+i) )
457 write(0,*)
'h(:N)=',h(n-3:n)
459 stop
'Nan during edge_values_explicit_h4'
492 integer,
intent(in) :: n
493 real,
dimension(:),
intent(in) :: h
494 real,
dimension(:),
intent(in) :: u
495 real,
dimension(:,:),
intent(inout) :: edge_val
496 real,
optional,
intent(in) :: h_neglect
497 logical,
optional,
intent(in) :: answers_2018
502 real :: h0_2, h1_2, h0h1
506 real,
dimension(5) :: x
507 real,
parameter :: c1_12 = 1.0 / 12.0
509 real,
dimension(4,4) :: asys
510 real,
dimension(4) :: bsys, csys
511 real,
dimension(N+1) :: tri_l, &
517 logical :: use_2018_answers
519 use_2018_answers = .true. ;
if (
present(answers_2018)) use_2018_answers = answers_2018
545 a = 2.0 * h1_2 * ( h1_2 + 2.0 * h0_2 + 3.0 * h0h1 ) / d4
546 b = 2.0 * h0_2 * ( h0_2 + 2.0 * h1_2 + 3.0 * h0h1 ) / d4
552 tri_b(i+1) = a * u(i) + b * u(i+1)
557 h0 = max( hneglect,
hminfrac*sum(h(1:4)) )
560 x(i) = x(i-1) + max( h0, h(i-1) )
565 if (use_2018_answers)
then
566 do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;
enddo
568 xavg = 0.5 * (x(i+1) + x(i))
570 asys(i,2) = dx * xavg
571 asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
572 asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
586 h0 = max( hneglect,
hminfrac*sum(h(n-3:n)) )
589 x(i) = x(i-1) + max( h0, h(n-5+i) )
593 dx = max(h0, h(n-4+i) )
594 if (use_2018_answers)
then
595 do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;
enddo
597 xavg = 0.5 * (x(i+1) + x(i))
599 asys(i,2) = dx * xavg
600 asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
601 asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
603 bsys(i) = u(n-4+i) * dx
617 edge_val(i,1) = tri_x(i)
618 edge_val(i-1,2) = tri_x(i)
620 edge_val(1,1) = tri_x(1)
621 edge_val(n,2) = tri_x(n+1)
661 integer,
intent(in) :: n
662 real,
dimension(:),
intent(in) :: h
663 real,
dimension(:),
intent(in) :: u
664 real,
dimension(:,:),
intent(inout) :: edge_val
665 real,
optional,
intent(in) :: h_neglect
666 logical,
optional,
intent(in) :: answers_2018
670 real :: h0, h1, h2, h3
672 real :: g_4, g_5, g_6
673 real :: d2, d3, d4, d5, d6
674 real :: n2, n3, n4, n5, n6
680 real :: h0ph1, h0ph1_2
681 real :: h0ph1_3, h0ph1_4
682 real :: h2ph3, h2ph3_2
683 real :: h2ph3_3, h2ph3_4
684 real :: h0ph1_5, h2ph3_5
687 real,
dimension(7) :: x
688 real,
parameter :: c1_12 = 1.0 / 12.0
689 real,
parameter :: c5_6 = 5.0 / 6.0
691 real,
dimension(6,6) :: asys
692 real,
dimension(6) :: bsys, csys
693 real,
dimension(N+1) :: tri_l, &
699 logical :: use_2018_answers
701 use_2018_answers = .true. ;
if (
present(answers_2018)) use_2018_answers = answers_2018
715 g = max( hneglect, h0+h1+h2+h3 )
742 d2 = ( h1_2 - g_2 ) / h0
743 d3 = ( h1_3 - g_3 ) / h0
744 d4 = ( h1_4 - g_4 ) / h0
745 d5 = ( h1_5 - g_5 ) / h0
746 d6 = ( h1_6 - g_6 ) / h0
755 n2 = ( g_2 - h2_2 ) / h3
756 n3 = ( g_3 - h2_3 ) / h3
757 n4 = ( g_4 - h2_4 ) / h3
758 n5 = ( g_5 - h2_5 ) / h3
759 n6 = ( g_6 - h2_6 ) / h3
771 asys(2,3) = -0.5 * d2
773 asys(2,5) = -0.5 * h2
774 asys(2,6) = -0.5 * n2
776 asys(3,1) = 0.5 * h1_2
777 asys(3,2) = 0.5 * h2_2
779 asys(3,4) = - h1_2 / 6.0
780 asys(3,5) = - h2_2 / 6.0
781 asys(3,6) = - n3 / 6.0
783 asys(4,1) = - h1_3 / 6.0
784 asys(4,2) = h2_3 / 6.0
785 asys(4,3) = - d4 / 24.0
786 asys(4,4) = h1_3 / 24.0
787 asys(4,5) = - h2_3 / 24.0
788 asys(4,6) = - n4 / 24.0
790 asys(5,1) = h1_4 / 24.0
791 asys(5,2) = h2_4 / 24.0
792 asys(5,3) = d5 / 120.0
793 asys(5,4) = - h1_4 / 120.0
794 asys(5,5) = - h2_4 / 120.0
795 asys(5,6) = - n5 / 120.0
797 asys(6,1) = - h1_5 / 120.0
798 asys(6,2) = h2_5 / 120.0
799 asys(6,3) = - d6 / 720.0
800 asys(6,4) = h1_5 / 720.0
801 asys(6,5) = - h2_5 / 720.0
802 asys(6,6) = - n6 / 720.0
804 bsys(:) = (/ -1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
818 tri_b(k+1) = a * u(k-1) + b * u(k) + c * u(k+1) + d * u(k+2)
832 g = max( hneglect, h0+h1+h2+h3 )
860 h0ph1_2 = h0ph1 * h0ph1
861 h0ph1_3 = h0ph1_2 * h0ph1
862 h0ph1_4 = h0ph1_2 * h0ph1_2
863 h0ph1_5 = h0ph1_3 * h0ph1_2
865 d2 = ( h1_2 - g_2 ) / h0
866 d3 = ( h1_3 - g_3 ) / h0
867 d4 = ( h1_4 - g_4 ) / h0
868 d5 = ( h1_5 - g_5 ) / h0
869 d6 = ( h1_6 - g_6 ) / h0
878 n2 = ( g_2 - h2_2 ) / h3
879 n3 = ( g_3 - h2_3 ) / h3
880 n4 = ( g_4 - h2_4 ) / h3
881 n5 = ( g_5 - h2_5 ) / h3
882 n6 = ( g_6 - h2_6 ) / h3
894 asys(2,3) = -0.5 * d2
896 asys(2,5) = -0.5 * h2
897 asys(2,6) = -0.5 * n2
899 asys(3,1) = 0.5 * h0ph1_2
902 asys(3,4) = - h1_2 / 6.0
903 asys(3,5) = - h2_2 / 6.0
904 asys(3,6) = - n3 / 6.0
906 asys(4,1) = - h0ph1_3 / 6.0
908 asys(4,3) = - d4 / 24.0
909 asys(4,4) = h1_3 / 24.0
910 asys(4,5) = - h2_3 / 24.0
911 asys(4,6) = - n4 / 24.0
913 asys(5,1) = h0ph1_4 / 24.0
915 asys(5,3) = d5 / 120.0
916 asys(5,4) = - h1_4 / 120.0
917 asys(5,5) = - h2_4 / 120.0
918 asys(5,6) = - n5 / 120.0
920 asys(6,1) = - h0ph1_5 / 120.0
922 asys(6,3) = - d6 / 720.0
923 asys(6,4) = h1_5 / 720.0
924 asys(6,5) = - h2_5 / 720.0
925 asys(6,6) = - n6 / 720.0
927 bsys(:) = (/ -1.0, h1, -0.5*h1_2, h1_3/6.0, -h1_4/24.0, h1_5/120.0 /)
941 tri_b(2) = a * u(1) + b * u(2) + c * u(3) + d * u(4)
944 g = max( hneglect,
hminfrac*sum(h(1:6)) )
947 x(i) = x(i-1) + max( g, h(i-1) )
952 if (use_2018_answers)
then
953 do j = 1,6 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;
enddo
955 xavg = 0.5 * (x(i+1) + x(i))
957 asys(i,2) = dx * xavg
958 asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
959 asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
960 asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4)
961 asys(i,6) = dx * xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4)
984 g = max( hneglect, h0+h1+h2+h3 )
1012 h2ph3_2 = h2ph3 * h2ph3
1013 h2ph3_3 = h2ph3_2 * h2ph3
1014 h2ph3_4 = h2ph3_2 * h2ph3_2
1015 h2ph3_5 = h2ph3_3 * h2ph3_2
1017 d2 = ( h1_2 - g_2 ) / h0
1018 d3 = ( h1_3 - g_3 ) / h0
1019 d4 = ( h1_4 - g_4 ) / h0
1020 d5 = ( h1_5 - g_5 ) / h0
1021 d6 = ( h1_6 - g_6 ) / h0
1030 n2 = ( g_2 - h2_2 ) / h3
1031 n3 = ( g_3 - h2_3 ) / h3
1032 n4 = ( g_4 - h2_4 ) / h3
1033 n5 = ( g_5 - h2_5 ) / h3
1034 n6 = ( g_6 - h2_6 ) / h3
1046 asys(2,3) = -0.5 * d2
1047 asys(2,4) = 0.5 * h1
1048 asys(2,5) = -0.5 * h2
1049 asys(2,6) = -0.5 * n2
1052 asys(3,2) = 0.5 * h2ph3_2
1053 asys(3,3) = d3 / 6.0
1054 asys(3,4) = - h1_2 / 6.0
1055 asys(3,5) = - h2_2 / 6.0
1056 asys(3,6) = - n3 / 6.0
1059 asys(4,2) = h2ph3_3 / 6.0
1060 asys(4,3) = - d4 / 24.0
1061 asys(4,4) = h1_3 / 24.0
1062 asys(4,5) = - h2_3 / 24.0
1063 asys(4,6) = - n4 / 24.0
1066 asys(5,2) = h2ph3_4 / 24.0
1067 asys(5,3) = d5 / 120.0
1068 asys(5,4) = - h1_4 / 120.0
1069 asys(5,5) = - h2_4 / 120.0
1070 asys(5,6) = - n5 / 120.0
1073 asys(6,2) = h2ph3_5 / 120.0
1074 asys(6,3) = - d6 / 720.0
1075 asys(6,4) = h1_5 / 720.0
1076 asys(6,5) = - h2_5 / 720.0
1077 asys(6,6) = - n6 / 720.0
1079 bsys(:) = (/ -1.0, -h2, -0.5*h2_2, -h2_3/6.0, -h2_4/24.0, -h2_5/120.0 /)
1093 tri_b(n) = a * u(n-3) + b * u(n-2) + c * u(n-1) + d * u(n)
1096 g = max( hneglect,
hminfrac*sum(h(n-5:n)) )
1099 x(i) = x(i-1) + max( g, h(n-7+i) )
1103 dx = max( g, h(n-6+i) )
1104 if (use_2018_answers)
then
1105 do j = 1,6 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;
enddo
1107 xavg = 0.5 * (x(i+1) + x(i))
1109 asys(i,2) = dx * xavg
1110 asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
1111 asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
1112 asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4)
1113 asys(i,6) = dx * xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4)
1115 bsys(i) = u(n-6+i) * dx
1130 edge_val(i,1) = tri_x(i)
1131 edge_val(i-1,2) = tri_x(i)
1133 edge_val(1,1) = tri_x(1)
1134 edge_val(n,2) = tri_x(n+1)