MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_weno.t
Go to the documentation of this file.
1 module mod_weno
2  ! All kinds of (W)ENO schemes
3  !
4  ! 2019.9.19 WENO(-JS)5 transplant from the BHAC code;
5  ! 2019.9.20 WENO3;
6  ! 2019.9.21 WENO-Z5;
7  ! 2019.9.22 WENO-Z+5 transplant from the BHAC code;
8  ! 2019.10.30 WENO(-JS)7;
9  ! 2019.10.31 MPWENO7;
10  ! 2019.11.1 exENO7;
11  ! 2019.11.7 clean up the code, comment out the interpolation variation;
12  ! 2019.12.9 WENO-YC3;
13  ! 2020.1.15 new WENO5 limiter WENO5NM for stretched grid.
14  ! 2020.4.15 WENO5-CU6: hybrid sixth-order linear & WENO5
15 
16  implicit none
17  private
18 
19  public :: weno3limiter
20  public :: weno5limiter
21  public :: weno5nmlimiter
22  public :: weno5limiterl
23  public :: weno5nmlimiterl
24  public :: weno5limiterr
25  public :: weno5nmlimiterr
26  public :: teno5adlimiter
27  public :: weno5cu6limiter
28  public :: weno7limiter
29  public :: exeno7limiter
30 
31 contains
32 
33  subroutine weno3limiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
35 
36  integer, intent(in) :: ixI^L, iL^L, idims, var
37  double precision, intent(in) :: dxdim
38  double precision, intent(in) :: w(ixi^s,1:nw)
39  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
40  !> local
41  integer :: iLm^L, iLp^L, iLpp^L
42  double precision :: f_array(ixi^s,1:nw,2), d_array(2)
43  double precision :: beta(ixi^s,1:nw,2),tau(ixi^s,1:nw)
44  double precision :: u1_coeff(2), u2_coeff(2)
45  double precision :: alpha_array(ixi^s,1:nw,2), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
46  integer :: i, iw
47  double precision, parameter :: weno_eps_machine = 1.0d-18
48 
49  ! iL^L holds the indices of interfaces to reconstruct to. Convention is that a center index holds the _right-side_ interface.
50  ilm^l=il^l-kr(idims,^d);
51  ilp^l=il^l+kr(idims,^d);
52  ilpp^l=ilp^l+kr(idims,^d);
53  d_array(1:2) = (/ 1.0d0/3.0d0, 2.0d0/3.0d0 /)
54  u1_coeff(1:2) = (/ -1.d0/2.d0, 3.d0/2.d0 /)
55  u2_coeff(1:2) = (/ 1.d0/2.d0, 1.d0/2.d0 /)
56 
57  !> left side
58  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilm^s,1:nwflux) + u1_coeff(2) * w(il^s,1:nwflux)
59  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(il^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux)
60 
61  beta(il^s,1:nwflux,1) = (w(il^s,1:nwflux) - w(ilm^s,1:nwflux))**2
62  beta(il^s,1:nwflux,2) = (w(ilp^s,1:nwflux) - w(il^s,1:nwflux))**2
63 
64  alpha_sum(il^s,1:nwflux) = 0.0d0
65  select case(var)
66  case(1)
67  do i = 1,2
68  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
69  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
70  end do
71  case(2)
72  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,2) - beta(il^s,1:nwflux,1))
73  do i = 1,2
74  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
75  (beta(il^s,1:nwflux,i) + dxdim**2)))
76  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
77  end do
78  end select
79 
80  flux(il^s,1:nwflux) = 0.0d0
81  do i = 1,2
82  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
83  end do
84 
85  !> left value at right interface
86  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
87 
88  !> right side
89  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilpp^s,1:nwflux) + u1_coeff(2) * w(ilp^s,1:nwflux)
90  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilp^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux)
91 
92  beta(il^s,1:nwflux,1) = (w(ilpp^s,1:nwflux) - w(ilp^s,1:nwflux))**2
93  beta(il^s,1:nwflux,2) = (w(ilp^s,1:nwflux) - w(il^s,1:nwflux))**2
94 
95  alpha_sum(il^s,1:nwflux) = 0.0d0
96  select case(var)
97  case(1)
98  do i = 1,2
99  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
100  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
101  end do
102  case(2)
103  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,2) - beta(il^s,1:nwflux,1))
104  do i = 1,2
105  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
106  (beta(il^s,1:nwflux,i) + dxdim**2)))
107  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
108  end do
109  end select
110 
111  flux(il^s,1:nwflux) = 0.0d0
112  do i = 1,2
113  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
114  end do
115 
116  !> right value at right interface
117  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
118 
119  end subroutine weno3limiter
120 
121  subroutine weno5limiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
123 
124  integer, intent(in) :: ixI^L, iL^L, idims, var
125  double precision, intent(in) :: dxdim
126  double precision, intent(in) :: w(ixi^s,1:nw)
127  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
128  !> local
129  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
130  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
131  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
132  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
133  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
134  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
135  integer :: i
136  double precision, parameter :: weno_eps_machine = 1.0d-18
137  double precision :: lambda
138  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
139 
140  ilm^l=il^l-kr(idims,^d);
141  ilmm^l=ilm^l-kr(idims,^d);
142  ilp^l=il^l+kr(idims,^d);
143  ilpp^l=ilp^l+kr(idims,^d);
144  ilppp^l=ilpp^l+kr(idims,^d);
145  lambda = dxdim**weno_dx_exp
146  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
147 ! reconstruction variation
148  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
149  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
150  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
151  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
152 ! interpolation variation
153 ! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
154 ! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
155 ! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
156 ! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
157 
158  !> left side
159  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilmm^s,1:nwflux) + u1_coeff(2) * w(ilm^s,1:nwflux) + u1_coeff(3) * w(il^s,1:nwflux)
160  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux)
161  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(ilpp^s,1:nwflux)
162 
163  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilmm^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0*w(ilm^s,1:nwflux))**2 &
164  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
165  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilm^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
166  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
167  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(il^s,1:nwflux) + w(ilpp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
168  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**2
169 
170  alpha_sum(il^s,1:nwflux) = 0.0d0
171  select case(var)
172  ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
173  case(1)
174  do i = 1,3
175  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
176  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
177  end do
178  case(2)
179  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
180  do i = 1,3
181  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
182  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
183  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
184  end do
185  case(3)
186  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
187  do i = 1,3
188  tmp(il^s,1:nwflux) = (tau(il^s,1:nwflux) + weno_eps_machine) / (beta(il^s,1:nwflux,i) + weno_eps_machine)
189  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nwflux)**2 + lambda/tmp(il^s,1:nwflux))
190  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
191  end do
192  end select
193 
194  flux(il^s,1:nwflux) = 0.0d0
195  do i = 1,3
196  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
197  end do
198 
199  !> left value at right interface
200  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
201 
202  !> right side
203  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilppp^s,1:nwflux) + u1_coeff(2) * w(ilpp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
204  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilpp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux)
205  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilp^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilm^s,1:nwflux)
206 
207  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilppp^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0*w(ilpp^s,1:nwflux))**2 &
208  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
209  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilpp^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
210  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
211  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(ilp^s,1:nwflux) + w(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
212  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
213 
214  alpha_sum(il^s,1:nwflux) = 0.0d0
215  select case(var)
216  case(1)
217  do i = 1,3
218  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
219  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
220  end do
221  case(2)
222  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
223  do i = 1,3
224  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
225  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
226  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
227  end do
228  case(3)
229  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
230  do i = 1,3
231  tmp(il^s,1:nwflux) = (tau(il^s,1:nwflux) + weno_eps_machine) / (beta(il^s,1:nwflux,i) + weno_eps_machine)
232  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nwflux)**2 + lambda/tmp(il^s,1:nwflux))
233  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
234  end do
235  end select
236  flux(il^s,1:nwflux) = 0.0d0
237  do i = 1,3
238  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
239  end do
240 
241  !> right value at right interface
242  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
243 
244  end subroutine weno5limiter
245 
246 
247  subroutine teno5adlimiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC)
249  integer, intent(in) :: ixI^L, iL^L, idims
250  double precision, intent(in) :: dxdim
251  double precision, intent(in) :: w(ixi^s,1:nw)
252  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
253  !> local
254  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
255  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
256  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
257  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
258  double precision :: gam_sum(ixi^s,1:nw),tau(ixi^s,1:nw),delta_sum(ixi^s,1:nw)
259  double precision :: gam(ixi^s,1:nw,3), kai(ixi^s,1:nw,3), delta(ixi^s,1:nw,3)
260  double precision :: flux(ixi^s,1:nw), kai1(ixi^s,1:nw,3), theta(ixi^s,1:nw)
261  integer :: marray(ixi^s,1:nw)
262  integer :: i
263  double precision, parameter :: weno_eps_machine = 1.0d-18
264 
265  ilm^l=il^l-kr(idims,^d);
266  ilmm^l=ilm^l-kr(idims,^d);
267  ilp^l=il^l+kr(idims,^d);
268  ilpp^l=ilp^l+kr(idims,^d);
269  ilppp^l=ilpp^l+kr(idims,^d);
270  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
271  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
272  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
273  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
274  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
275 
276  !> left side
277  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilmm^s,1:nwflux) + u1_coeff(2) * w(ilm^s,1:nwflux) + u1_coeff(3) * w(il^s,1:nwflux)
278  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux)
279  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(ilpp^s,1:nwflux)
280 
281  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilmm^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0*w(ilm^s,1:nwflux))**2 &
282  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
283  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilm^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
284  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
285  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(il^s,1:nwflux) + w(ilpp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
286  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**2
287 
288  gam_sum(il^s,1:nwflux) = 0.0d0
289  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
290  do i = 1,3
291  kai1(il^s,1:nwflux,i) = (tau(il^s,1:nwflux) / (beta(il^s,1:nwflux,i) + weno_eps_machine))
292  gam(il^s,1:nwflux,i) = (1.d0 + kai1(il^s,1:nwflux,i))**2
293  gam_sum(il^s,1:nwflux) = gam_sum(il^s,1:nwflux) + gam(il^s,1:nwflux,i)
294  end do
295  theta(il^s,1:nwflux) = one / (one + maxval(kai1(il^s,1:nwflux,1:3)/10.d0, dim=ndim+2))
296  marray(il^s,1:nwflux)=-floor(4.d0 + theta(il^s,1:nwflux)*6.d0)
297  do i = 1,3
298  kai(il^s,1:nwflux,i) = gam(il^s,1:nwflux,i) / gam_sum(il^s,1:nwflux)
299  where(kai(il^s,1:nwflux,i) .lt. 10**marray(il^s,1:nwflux))
300  delta(il^s,1:nwflux,i)=zero
301  else where
302  delta(il^s,1:nwflux,i)=one
303  end where
304  end do
305  delta_sum=zero
306  do i = 1,3
307  delta_sum(il^s,1:nwflux)=delta_sum(il^s,1:nwflux)+delta(il^s,1:nwflux,i)*d_array(i)
308  end do
309  flux(il^s,1:nwflux)=0.0d0
310  do i = 1,3
311  flux(il^s,1:nwflux)=flux(il^s,1:nwflux)+f_array(il^s,1:nwflux,i)*delta(il^s,1:nwflux,i)*d_array(i)/(delta_sum(il^s,1:nwflux))
312  end do
313 
314  !> left value at right interface
315  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
316 
317  !> right side
318  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilppp^s,1:nwflux) + u1_coeff(2) * w(ilpp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
319  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilpp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux)
320  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilp^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilm^s,1:nwflux)
321 
322  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilppp^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0*w(ilpp^s,1:nwflux))**2 &
323  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
324  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilpp^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
325  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
326  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(ilp^s,1:nwflux) + w(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
327  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
328 
329 
330  gam_sum(il^s,1:nwflux)=0.0d0
331  tau(il^s,1:nwflux)=abs(beta(il^s,1:nwflux,1)-beta(il^s,1:nwflux,3))
332  do i=1,3
333  kai1(il^s,1:nwflux,i)=(tau(il^s,1:nwflux)/(beta(il^s,1:nwflux,i)+weno_eps_machine))
334  gam(il^s,1:nwflux,i)=(1.d0+kai1(il^s,1:nwflux,i))**6
335  gam_sum(il^s,1:nwflux)=gam_sum(il^s,1:nwflux)+gam(il^s,1:nwflux,i)
336  end do
337  theta(il^s,1:nwflux)=one/(one+maxval(kai1(il^s,1:nwflux,1:3)/10.d0,dim=ndim+2))
338  marray(il^s,1:nwflux)=-floor(4.d0+theta(il^s,1:nwflux)*6.d0)
339  do i=1,3
340  kai(il^s,1:nwflux,i) = gam(il^s,1:nwflux,i)/gam_sum(il^s,1:nwflux)
341  where(kai(il^s,1:nwflux,i) .lt. 10**marray(il^s,1:nwflux))
342  delta(il^s,1:nwflux,i)=zero
343  else where
344  delta(il^s,1:nwflux,i)=one
345  end where
346  end do
347  delta_sum=zero
348  do i = 1,3
349  delta_sum(il^s,1:nwflux)=delta_sum(il^s,1:nwflux)+delta(il^s,1:nwflux,i)*d_array(i)
350  end do
351  flux(il^s,1:nwflux)=0.0d0
352  do i = 1,3
353  flux(il^s,1:nwflux)=flux(il^s,1:nwflux)+f_array(il^s,1:nwflux,i)*delta(il^s,1:nwflux,i)*d_array(i)/(delta_sum(il^s,1:nwflux))
354  end do
355 
356  !> right value at right interface
357  wrc(il^s,1:nwflux)=flux(il^s,1:nwflux)
358  end subroutine teno5adlimiter
359 
360  subroutine weno5nmlimiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
362 
363  integer, intent(in) :: ixI^L,iL^L,idims,var
364  double precision, intent(in) :: dxdim
365  double precision, intent(in) :: w(ixi^s,1:nw)
366  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
367  !> local
368  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
369  integer :: iM^L, iMp^L
370  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
371  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
372  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
373  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
374  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
375  double precision :: wc(ixi^s,1:nw), wd(ixi^s,1:nw), we(ixi^s,1:nw)
376  integer :: i,j
377  double precision, parameter :: weno_eps_machine = 1.0d-18
378  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
379  double precision :: lambda(ixi^s)
380 
381  ilm^l=il^l-kr(idims,^d);
382  ilmm^l=ilm^l-kr(idims,^d);
383  ilp^l=il^l+kr(idims,^d);
384  ilpp^l=ilp^l+kr(idims,^d);
385  ilppp^l=ilpp^l+kr(idims,^d);
386  immin^d=ilmmin^d;
387  immax^d=ilpmax^d;
388  imp^l=im^l+kr(idims,^d);
389  lambda(il^s)=block%dx(il^s,idims)**weno_dx_exp
390  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
391  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
392  u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
393  u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
394  u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
395  do i = 1, nwflux
396  wc(im^s,i) = (block%dx(imp^s,idims) * w(im^s,i) + block%dx(im^s,idims) * w(imp^s,i)) / &
397  (block%dx(imp^s,idims) + block%dx(im^s,idims))
398  wd(il^s,i) = ((2.d0 * block%dx(ilm^s,idims) + block%dx(ilmm^s,idims)) * w(ilm^s,i) - block%dx(ilm^s,idims) * w(ilmm^s,i)) / &
399  (block%dx(ilmm^s,idims) + block%dx(ilm^s,idims))
400  we(il^s,i) = ((2.d0 * block%dx(ilpp^s,idims) + block%dx(ilppp^s,idims)) * w(ilpp^s,i) - block%dx(ilpp^s,idims) * w(ilppp^s,i)) / &
401  (block%dx(ilppp^s,idims) + block%dx(ilpp^s,idims))
402  enddo
403  !> left side
404  f_array(il^s,1:nwflux,1) = u1_coeff(1) * wd(il^s,1:nwflux) + u1_coeff(2) * wc(ilm^s,1:nwflux)+ u1_coeff(3) * w(il^s,1:nwflux)
405  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
406  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * wc(ilp^s,1:nwflux)
407 
408  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2 &
409  + beta_coeff(2) * (2.d0 * w(il^s,1:nwflux) - wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2
410  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilm^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
411  + beta_coeff(2) * (wc(ilm^s,1:nwflux) - wc(il^s,1:nwflux))**2
412  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
413  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + wc(ilp^s,1:nwflux))**2
414  alpha_sum(il^s,1:nwflux) = 0.0d0
415  select case(var)
416  ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
417  case(1)
418  do i = 1,3
419  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
420  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
421  end do
422  case(2)
423  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
424  do i = 1,3
425  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
426  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
427  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
428  end do
429  case(3)
430  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
431  do i=1,3
432  do j=1,nwflux
433  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
434  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
435  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
436  end do
437  end do
438  end select
439  flux(il^s,1:nwflux) = 0.0d0
440  do i = 1,3
441  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
442  end do
443  !> left value at right interface
444  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
445  !> right side
446  f_array(il^s,1:nwflux,1) = u1_coeff(1) * we(il^s,1:nwflux) + u1_coeff(2) * wc(ilp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
447  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
448  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * wc(ilm^s,1:nwflux)
449  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2 &
450  + beta_coeff(2) * (2.d0 * w(ilp^s,1:nwflux) - wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2
451  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilp^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
452  + beta_coeff(2) * (wc(ilp^s,1:nwflux) - wc(il^s,1:nwflux))**2
453  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
454  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + wc(ilm^s,1:nwflux))**2
455  alpha_sum(il^s,1:nwflux) = 0.0d0
456  select case(var)
457  case(1)
458  do i = 1,3
459  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
460  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
461  end do
462  case(2)
463  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
464  do i = 1,3
465  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
466  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
467  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
468  end do
469  case(3)
470  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
471  do i = 1,3
472  do j = 1,nwflux
473  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
474  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
475  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
476  end do
477  end do
478  end select
479  flux(il^s,1:nwflux) = 0.0d0
480  do i = 1,3
481  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
482  end do
483  !> right value at right interface
484  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
485 
486  end subroutine weno5nmlimiter
487 
488  subroutine weno5limiterl(ixI^L,iL^L,idims,w,wLC,var)
490 
491  integer, intent(in) :: ixI^L, iL^L, idims
492  integer, intent(in) :: var
493  double precision, intent(in) :: w(ixi^s,1:nw)
494  double precision, intent(inout) :: wLC(ixi^s,1:nw)
495  !> local
496  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
497  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
498  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
499  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
500  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
501  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
502  integer :: i,j
503  double precision :: lambda(ixi^s)
504  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
505  double precision, parameter :: weno_eps_machine = 1.0d-18
506 
507  ilm^l=il^l-kr(idims,^d);
508  ilmm^l=ilm^l-kr(idims,^d);
509  ilp^l=il^l+kr(idims,^d);
510  ilpp^l=ilp^l+kr(idims,^d);
511  lambda=block%dx(il^s,idims)**weno_dx_exp
512  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
513 ! reconstruction variation
514  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
515  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
516  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
517  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
518 ! interpolation variation
519 ! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
520 ! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
521 ! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
522 ! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
523 
524  !> left side
525  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilmm^s,1:nwflux) + u1_coeff(2) * w(ilm^s,1:nwflux) + u1_coeff(3) * w(il^s,1:nwflux)
526  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux)
527  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(ilpp^s,1:nwflux)
528 
529  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilmm^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0*w(ilm^s,1:nwflux))**2 &
530  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
531  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilm^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
532  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
533  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(il^s,1:nwflux) + w(ilpp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
534  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**2
535 
536  alpha_sum(il^s,1:nwflux) = 0.0d0
537  select case(var)
538  ! case1 for wenojs, case2 for wenoz
539  case(1)
540  do i = 1,3
541  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
542  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
543  end do
544  case(2)
545  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
546  do i = 1,3
547  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
548  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
549  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
550  end do
551  case(3)
552  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
553  do i=1,3
554  do j=1,nwflux
555  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
556  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
557  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
558  end do
559  end do
560  end select
561  flux(il^s,1:nwflux) = 0.0d0
562  do i = 1,3
563  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
564  end do
565 
566  !> left value at right interface
567  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
568 
569  end subroutine weno5limiterl
570 
571  subroutine weno5limiterr(ixI^L,iL^L,idims,w,wRC,var)
573 
574  integer, intent(in) :: ixI^L, iL^L, idims
575  integer, intent(in) :: var
576  double precision, intent(in) :: w(ixi^s,1:nw)
577  double precision, intent(inout) :: wRC(ixi^s,1:nw)
578  !> local
579  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
580  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
581  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
582  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
583  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
584  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
585  integer :: i,j
586  double precision :: lambda(ixi^s)
587  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
588  double precision, parameter :: weno_eps_machine = 1.0d-18
589 
590  ilm^l=il^l-kr(idims,^d);
591  ilp^l=il^l+kr(idims,^d);
592  ilpp^l=ilp^l+kr(idims,^d);
593  ilppp^l=ilpp^l+kr(idims,^d);
594  lambda=block%dx(il^s,idims)**weno_dx_exp
595  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
596 ! reconstruction variation
597  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
598  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
599  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
600  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
601 ! interpolation variation
602 ! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
603 ! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
604 ! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
605 ! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
606 
607  !> right side
608  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilppp^s,1:nwflux) + u1_coeff(2) * w(ilpp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
609  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilpp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux)
610  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilp^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilm^s,1:nwflux)
611 
612  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilppp^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0*w(ilpp^s,1:nwflux))**2 &
613  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
614  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilpp^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
615  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
616  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(ilp^s,1:nwflux) + w(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
617  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
618 
619  alpha_sum(il^s,1:nwflux) = 0.0d0
620  select case(var)
621  case(1)
622  do i = 1,3
623  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
624  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
625  end do
626  case(2)
627  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
628  do i = 1,3
629  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
630  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
631  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
632  end do
633  case(3)
634  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
635  do i = 1,3
636  do j = 1,nwflux
637  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
638  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
639  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
640  end do
641  end do
642  end select
643  flux(il^s,1:nwflux) = 0.0d0
644  do i = 1,3
645  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
646  end do
647 
648  !> right value at right interface
649  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
650 
651  end subroutine weno5limiterr
652 
653  subroutine weno5nmlimiterl(ixI^L,iL^L,idims,w,wLC,var)
655 
656  integer, intent(in) :: ixI^L,iL^L,idims,var
657  double precision, intent(in) :: w(ixi^s,1:nw)
658  double precision, intent(inout) :: wLC(ixi^s,1:nw)
659  !> local
660  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L
661  integer :: iM^L, iMp^L
662  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
663  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
664  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
665  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
666  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
667  double precision :: wc(ixi^s,1:nw), wd(ixi^s,1:nw)
668  integer :: i,j
669  double precision, parameter :: weno_eps_machine = 1.0d-18
670  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
671  double precision :: lambda(ixi^s)
672 
673  ilm^l=il^l-kr(idims,^d);
674  ilmm^l=ilm^l-kr(idims,^d);
675  ilp^l=il^l+kr(idims,^d);
676  ilpp^l=ilp^l+kr(idims,^d);
677  immin^d=ilmmin^d;
678  immax^d=ilpmax^d;
679  imp^l=im^l+kr(idims,^d);
680  lambda(il^s)=block%dx(il^s,idims)**weno_dx_exp
681  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
682  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
683  u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
684  u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
685  u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
686  do i = 1, nwflux
687  wc(im^s,i) = (block%dx(imp^s,idims) * w(im^s,i) + block%dx(im^s,idims) * w(imp^s,i)) / &
688  (block%dx(imp^s,idims) + block%dx(im^s,idims))
689  wd(il^s,i) = ((2.d0 * block%dx(ilm^s,idims) + block%dx(ilmm^s,idims)) * w(ilm^s,i) - block%dx(ilm^s,idims) * w(ilmm^s,i)) / &
690  (block%dx(ilmm^s,idims) + block%dx(ilm^s,idims))
691  enddo
692  f_array(il^s,1:nwflux,1) = u1_coeff(1) * wd(il^s,1:nwflux) + u1_coeff(2) * wc(ilm^s,1:nwflux)+ u1_coeff(3) * w(il^s,1:nwflux)
693  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
694  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * wc(ilp^s,1:nwflux)
695  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2 &
696  + beta_coeff(2) * (2.d0 * w(il^s,1:nwflux) - wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2
697  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilm^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
698  + beta_coeff(2) * (wc(ilm^s,1:nwflux) - wc(il^s,1:nwflux))**2
699  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
700  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + wc(ilp^s,1:nwflux))**2
701  alpha_sum(il^s,1:nwflux) = 0.0d0
702  select case(var)
703  ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
704  case(1)
705  do i = 1,3
706  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
707  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
708  end do
709  case(2)
710  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
711  do i = 1,3
712  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
713  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
714  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
715  end do
716  case(3)
717  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
718  do i=1,3
719  do j=1,nwflux
720  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
721  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
722  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
723  end do
724  end do
725  end select
726  flux(il^s,1:nwflux) = 0.0d0
727  do i = 1,3
728  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
729  end do
730  !> left value at right interface
731  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
732 
733  end subroutine weno5nmlimiterl
734 
735  subroutine weno5nmlimiterr(ixI^L,iL^L,idims,w,wRC,var)
737 
738  integer, intent(in) :: ixI^L,iL^L,idims,var
739  double precision, intent(in) :: w(ixi^s,1:nw)
740  double precision, intent(inout) :: wRC(ixi^s,1:nw)
741  !> local
742  integer :: iLm^L,iLp^L,iLpp^L,iLppp^L
743  integer :: iM^L, iMp^L
744  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
745  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
746  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
747  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
748  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
749  double precision :: wc(ixi^s,1:nw), we(ixi^s,1:nw)
750  integer :: i,j
751  double precision, parameter :: weno_eps_machine = 1.0d-18
752  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
753  double precision :: lambda(ixi^s)
754 
755  ilm^l=il^l-kr(idims,^d);
756  ilp^l=il^l+kr(idims,^d);
757  ilpp^l=ilp^l+kr(idims,^d);
758  ilppp^l=ilpp^l+kr(idims,^d);
759  immin^d=ilmmin^d;
760  immax^d=ilpmax^d;
761  imp^l=im^l+kr(idims,^d);
762  lambda(il^s)=block%dx(il^s,idims)**weno_dx_exp
763  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
764  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
765  u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
766  u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
767  u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
768  do i = 1, nwflux
769  wc(im^s,i) = (block%dx(imp^s,idims) * w(im^s,i) + block%dx(im^s,idims) * w(imp^s,i)) / &
770  (block%dx(imp^s,idims) + block%dx(im^s,idims))
771  we(il^s,i) = ((2.d0 * block%dx(ilpp^s,idims) + block%dx(ilppp^s,idims)) * w(ilpp^s,i) - block%dx(ilpp^s,idims) * w(ilppp^s,i)) / &
772  (block%dx(ilppp^s,idims) + block%dx(ilpp^s,idims))
773  enddo
774  !> right side
775  f_array(il^s,1:nwflux,1) = u1_coeff(1) * we(il^s,1:nwflux) + u1_coeff(2) * wc(ilp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
776  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
777  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * wc(ilm^s,1:nwflux)
778  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2 &
779  + beta_coeff(2) * (2.d0 * w(ilp^s,1:nwflux) - wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2
780  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilp^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
781  + beta_coeff(2) * (wc(ilp^s,1:nwflux) - wc(il^s,1:nwflux))**2
782  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
783  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + wc(ilm^s,1:nwflux))**2
784  alpha_sum(il^s,1:nwflux) = 0.0d0
785  select case(var)
786  case(1)
787  do i = 1,3
788  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
789  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
790  end do
791  case(2)
792  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
793  do i = 1,3
794  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
795  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
796  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
797  end do
798  case(3)
799  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
800  do i = 1,3
801  do j = 1,nwflux
802  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
803  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
804  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
805  end do
806  end do
807  end select
808  flux(il^s,1:nwflux) = 0.0d0
809  do i = 1,3
810  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
811  end do
812  !> right value at right interface
813  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
814 
815  end subroutine weno5nmlimiterr
816 
817 
818 
819 
820  subroutine weno5cu6limiter(ixI^L,iL^L,idims,w,wLC,wRC)
822 
823  integer, intent(in) :: ixI^L, iL^L, idims
824  double precision, intent(in) :: w(ixi^s,1:nw)
825  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
826  !> local
827  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
828  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
829  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
830  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
831  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw)
832  double precision :: theta2(ixi^s,1:nw)
833  integer :: i
834  double precision, parameter :: weno_eps_machine=1.0d-18, theta_limit=0.7d0
835 
836  ilm^l=il^l-kr(idims,^d);
837  ilmm^l=ilm^l-kr(idims,^d);
838  ilp^l=il^l+kr(idims,^d);
839  ilpp^l=ilp^l+kr(idims,^d);
840  ilppp^l=ilpp^l+kr(idims,^d);
841 
842  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
843  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
844  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
845  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
846  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
847 
848  !> left side
849  beta(il^s,1:nwflux,1)=beta_coeff(1)*(w(ilmm^s,1:nwflux)+w(il^s,1:nwflux)-2.0d0*w(ilm^s,1:nwflux))**2&
850  +beta_coeff(2)*(w(ilmm^s,1:nwflux)-4.0d0*w(ilm^s,1:nwflux)+3.0d0*w(il^s,1:nwflux))**2
851  beta(il^s,1:nwflux,2)=beta_coeff(1)*(w(ilm^s,1:nwflux)+w(ilp^s,1:nwflux)-2.0d0*w(il^s,1:nwflux))**2&
852  +beta_coeff(2)*(w(ilm^s,1:nwflux)-w(ilp^s,1:nwflux))**2
853  beta(il^s,1:nwflux,3)=beta_coeff(1)*(w(il^s,1:nwflux)+w(ilpp^s,1:nwflux)-2.0d0*w(ilp^s,1:nwflux))**2&
854  +beta_coeff(2)*(3.0d0*w(il^s,1:nwflux)-4.0d0*w(ilp^s,1:nwflux)+w(ilpp^s,1:nwflux))**2
855  alpha_sum(il^s,1:nwflux)=zero
856  do i=1,3
857  alpha_array(il^s,1:nwflux,i)=d_array(i)/(beta(il^s,1:nwflux,i)+weno_eps_machine)**2
858  alpha_sum(il^s,1:nwflux)=alpha_sum(il^s,1:nwflux)+alpha_array(il^s,1:nwflux,i)
859  end do
860  do i=1,3
861  alpha_array(il^s,1:nwflux,i)=alpha_array(il^s,1:nwflux,i)/alpha_sum(il^s,1:nwflux)
862  end do
863  theta2(il^s,1:nwflux)=((alpha_array(il^s,1:nwflux,1)/d_array(1)-1.d0)**2&
864  +(alpha_array(il^s,1:nwflux,2)/d_array(2)-1.d0)**2&
865  +(alpha_array(il^s,1:nwflux,3)/d_array(3)-1.d0)**2)/83.d0
866  where(theta2(il^s,1:nwflux) .gt. theta_limit)
867  f_array(il^s,1:nwflux,1)=u1_coeff(1)*w(ilmm^s,1:nwflux)+u1_coeff(2)*w(ilm^s,1:nwflux)+u1_coeff(3)*w(il^s,1:nwflux)
868  f_array(il^s,1:nwflux,2)=u2_coeff(1)*w(ilm^s,1:nwflux)+u2_coeff(2)*w(il^s,1:nwflux)+u2_coeff(3)*w(ilp^s,1:nwflux)
869  f_array(il^s,1:nwflux,3)=u3_coeff(1)*w(il^s,1:nwflux)+u3_coeff(2)*w(ilp^s,1:nwflux)+u3_coeff(3)*w(ilpp^s,1:nwflux)
870  wlc(il^s,1:nwflux)=f_array(il^s,1:nwflux,1)*alpha_array(il^s,1:nwflux,1)&
871  +f_array(il^s,1:nwflux,2)*alpha_array(il^s,1:nwflux,2)&
872  +f_array(il^s,1:nwflux,3)*alpha_array(il^s,1:nwflux,3)
873  else where
874  wlc(il^s,1:nwflux)=1.d0/60.d0*(w(ilmm^s,1:nwflux)-8.d0*w(ilm^s,1:nwflux)+37.d0*w(il^s,1:nwflux)&
875  +37.d0*w(ilp^s,1:nwflux)-8.d0*w(ilpp^s,1:nwflux)+w(ilppp^s,1:nwflux))
876  end where
877 
878  !> right side
879  beta(il^s,1:nwflux,1)=beta_coeff(1)*(w(ilppp^s,1:nwflux)+w(ilp^s,1:nwflux)-2.0d0*w(ilpp^s,1:nwflux))**2&
880  +beta_coeff(2)*(w(ilppp^s,1:nwflux)-4.0d0*w(ilpp^s,1:nwflux)+3.0d0*w(ilp^s,1:nwflux))**2
881  beta(il^s,1:nwflux,2)=beta_coeff(1)*(w(ilpp^s,1:nwflux)+w(il^s,1:nwflux)-2.0d0*w(ilp^s,1:nwflux))**2&
882  +beta_coeff(2)*(w(ilpp^s,1:nwflux)-w(il^s,1:nwflux))**2
883  beta(il^s,1:nwflux,3)=beta_coeff(1)*(w(ilp^s,1:nwflux)+w(ilm^s,1:nwflux)-2.0d0*w(il^s,1:nwflux))**2&
884  +beta_coeff(2)*(3.0d0*w(ilp^s,1:nwflux)-4.0d0*w(il^s,1:nwflux)+w(ilm^s,1:nwflux))**2
885  alpha_sum(il^s,1:nwflux)=zero
886  do i=1,3
887  alpha_array(il^s,1:nwflux,i)=d_array(i)/(beta(il^s,1:nwflux,i)+weno_eps_machine)**2
888  alpha_sum(il^s,1:nwflux)=alpha_sum(il^s,1:nwflux)+alpha_array(il^s,1:nwflux,i)
889  end do
890  do i=1,3
891  alpha_array(il^s,1:nwflux,i)=alpha_array(il^s,1:nwflux,i)/alpha_sum(il^s,1:nwflux)
892  end do
893  theta2(il^s,1:nwflux)=((alpha_array(il^s,1:nwflux,1)/d_array(1)-1.d0)**2&
894  +(alpha_array(il^s,1:nwflux,2)/d_array(2)-1.d0)**2&
895  +(alpha_array(il^s,1:nwflux,3)/d_array(3)-1.d0)**2)/83.d0
896  where(theta2(il^s,1:nwflux) .gt. theta_limit)
897  f_array(il^s,1:nwflux,1)=u1_coeff(1)*w(ilppp^s,1:nwflux)+u1_coeff(2)*w(ilpp^s,1:nwflux)+u1_coeff(3)*w(ilp^s,1:nwflux)
898  f_array(il^s,1:nwflux,2)=u2_coeff(1)*w(ilpp^s,1:nwflux)+u2_coeff(2)*w(ilp^s,1:nwflux)+u2_coeff(3)*w(il^s,1:nwflux)
899  f_array(il^s,1:nwflux,3)=u3_coeff(1)*w(ilp^s,1:nwflux)+u3_coeff(2)*w(il^s,1:nwflux)+u3_coeff(3)*w(ilm^s,1:nwflux)
900  wrc(il^s,1:nwflux)=f_array(il^s,1:nwflux,1)*alpha_array(il^s,1:nwflux,1)&
901  +f_array(il^s,1:nwflux,2)*alpha_array(il^s,1:nwflux,2)&
902  +f_array(il^s,1:nwflux,3)*alpha_array(il^s,1:nwflux,3)
903  else where
904  wrc(il^s,1:nwflux)=1.d0/60.d0*(w(ilppp^s,1:nwflux)-8.d0*w(ilpp^s,1:nwflux)+37.d0*w(ilp^s,1:nwflux)&
905  +37.d0*w(il^s,1:nwflux)-8.d0*w(ilm^s,1:nwflux)+w(ilmm^s,1:nwflux))
906  end where
907  end subroutine weno5cu6limiter
908 
909  subroutine weno7limiter(ixI^L,iL^L,idims,w,wLC,wRC,var)
911 
912  integer, intent(in) :: ixI^L, iL^L, idims, var
913  double precision, intent(in) :: w(ixi^s,1:nw)
914  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
915  !> local
916  integer :: iLm^L, iLmm^L, iLmmm^L
917  integer :: iLp^L, iLpp^L, iLppp^L, iLpppp^L
918  integer :: id^L, idp^L, idpp^L, idm^L, ie^L, iem^L, iep^L, iepp^L
919 
920  double precision, dimension(4) :: d_array, u1_coeff, u2_coeff, u3_coeff, u4_coeff
921  double precision, dimension(ixI^S,1:nw,4) :: f_array, beta, alpha_array
922  double precision, dimension(ixI^S) :: a, b, c, tmp, tmp2, tmp3
923  double precision, dimension(ixI^S,1:nw) :: alpha_sum, d, dm4
924  double precision, dimension(ixI^S,1:nw) :: flux, flux_min, flux_max, flux_ul, flux_md, flux_lc
925  integer :: i,iw
926  double precision, parameter :: mpalpha = 2.d0, mpbeta = 4.d0
927  double precision, parameter :: weno_eps_machine = 1.0d-18
928 
929  ilm^l=il^l-kr(idims,^d);
930  ilmm^l=ilm^l-kr(idims,^d);
931  ilmmm^l=ilmm^l-kr(idims,^d);
932  ilp^l=il^l+kr(idims,^d);
933  ilpp^l=ilp^l+kr(idims,^d);
934  ilppp^l=ilpp^l+kr(idims,^d);
935  ilpppp^l=ilppp^l+kr(idims,^d);
936 
937  d_array(1:4) = (/ 1.d0/35.d0, 12.d0/35.d0, 18.d0/35.d0, 4.d0/35.d0 /)
938  u1_coeff(1:4) = (/ -1.d0/4.d0, 13.d0/12.d0, -23.d0/12.d0, 25.d0/12.d0 /)
939  u2_coeff(1:4) = (/ 1.d0/12.d0, -5.d0/12.d0, 13.d0/12.d0, 1.d0/4.d0 /)
940  u3_coeff(1:4) = (/ -1.d0/12.d0, 7.d0/12.d0, 7.d0/12.d0, -1.d0/12.d0 /)
941  u4_coeff(1:4) = (/ 1.d0/4.d0, 13.d0/12.d0, -5.d0/12.d0, 1.d0/12.d0 /)
942 
943  !> left side
944  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilmmm^s,1:nwflux) + u1_coeff(2) * w(ilmm^s,1:nwflux) + u1_coeff(3) * w(ilm^s,1:nwflux) &
945  + u1_coeff(4) * w(il^s,1:nwflux)
946  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilmm^s,1:nwflux) + u2_coeff(2) * w(ilm^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux) &
947  + u2_coeff(4) * w(ilp^s,1:nwflux)
948  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilm^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilp^s,1:nwflux) &
949  + u3_coeff(4) * w(ilpp^s,1:nwflux)
950  f_array(il^s,1:nwflux,4) = u4_coeff(1) * w(il^s,1:nwflux) + u4_coeff(2) * w(ilp^s,1:nwflux) + u4_coeff(3) * w(ilpp^s,1:nwflux) &
951  + u4_coeff(4) * w(ilppp^s,1:nwflux)
952 
953  beta(il^s,1:nwflux,1) = w(ilmmm^s,1:nwflux) * (547.d0 * w(ilmmm^s,1:nwflux) - 3882.d0 * w(ilmm^s,1:nwflux) + 4642.d0 * w(ilm^s,1:nwflux) &
954  - 1854.d0 * w(il^s,1:nwflux)) &
955  + w(ilmm^s,1:nwflux) * (7043.d0 * w(ilmm^s,1:nwflux) - 17246.d0 * w(ilm^s,1:nwflux) + 7042.d0 * w(il^s,1:nwflux)) &
956  + w(ilm^s,1:nwflux) * (11003.d0 * w(ilm^s,1:nwflux) - 9402.d0 * w(il^s,1:nwflux)) + 2107.d0 * w(il^s,1:nwflux)**2
957  beta(il^s,1:nwflux,2) = w(ilmm^s,1:nwflux) * (267.d0 * w(ilmm^s,1:nwflux) - 1642.d0 * w(ilm^s,1:nwflux) + 1602.d0 * w(il^s,1:nwflux) &
958  - 494.d0 * w(ilp^s,1:nwflux)) &
959  + w(ilm^s,1:nwflux) * (2843.d0 * w(ilm^s,1:nwflux) - 5966.d0 * w(il^s,1:nwflux) + 1922.d0 * w(ilp^s,1:nwflux)) &
960  + w(il^s,1:nwflux) * (3443.d0 * w(il^s,1:nwflux) - 2522.d0 * w(ilp^s,1:nwflux)) + 547.d0 * w(ilp^s,1:nwflux) ** 2
961  beta(il^s,1:nwflux,3) = w(ilm^s,1:nwflux) * (547.d0 * w(ilm^s,1:nwflux) - 2522.d0 * w(il^s,1:nwflux) + 1922.d0 * w(ilp^s,1:nwflux) &
962  - 494.d0 * w(ilpp^s,1:nwflux)) &
963  + w(il^s,1:nwflux) * (3443.d0 * w(il^s,1:nwflux) - 5966.d0 * w(ilp^s,1:nwflux) + 1602.d0 * w(ilpp^s,1:nwflux)) &
964  + w(ilp^s,1:nwflux) * (2843.d0 * w(ilp^s,1:nwflux) - 1642.d0 * w(ilpp^s,1:nwflux)) + 267.d0 * w(ilpp^s,1:nwflux) ** 2
965  beta(il^s,1:nwflux,4) = w(il^s,1:nwflux) * (2107.d0 * w(il^s,1:nwflux) - 9402.d0 * w(ilp^s,1:nwflux) + 7042.d0 * w(ilpp^s,1:nwflux) &
966  - 1854.d0 * w(ilppp^s,1:nwflux)) &
967  + w(ilp^s,1:nwflux) * (11003.d0 * w(ilp^s,1:nwflux) - 17246.d0 * w(ilpp^s,1:nwflux) + 4642.d0 * w(ilppp^s,1:nwflux)) &
968  + w(ilpp^s,1:nwflux) * (7043.d0 * w(ilpp^s,1:nwflux) - 3882.d0 * w(ilppp^s,1:nwflux)) &
969  + 547.d0 * w(ilppp^s,1:nwflux) ** 2
970 
971  alpha_sum(il^s,1:nwflux) = 0.0d0
972  do i = 1,4
973  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
974  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
975  end do
976  flux(il^s,1:nwflux) = 0.0d0
977  do i = 1,4
978  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
979  end do
980 
981  select case(var)
982  ! case1 for wenojs, case2 for mpweno
983  case(1)
984  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
985  case(2)
986  idmax^d=ilmax^d; idmin^d=ilmin^d-kr(idims,^d);
987  idm^l=id^l-kr(idims,^d);
988  idp^l=id^l+kr(idims,^d);
989 
990  iemax^d=idmax^d+kr(idims,^d); iemin^d=idmin^d;
991  iem^l=ie^l-kr(idims,^d);
992  iep^l=ie^l+kr(idims,^d);
993 
994  d(ie^s,1:nwflux) = w(iep^s,1:nwflux)-2.0d0*w(ie^s,1:nwflux)+w(iem^s,1:nwflux)
995 
996  do iw=1,nwflux
997  a(id^s) = 4.0d0*d(id^s,iw)-d(idp^s,iw)
998  b(id^s) = 4.0d0*d(idp^s,iw)-d(id^s,iw)
999  call minmod(ixi^l,id^l,a,b,tmp)
1000  a(id^s) = d(id^s,iw)
1001  b(id^s) = d(idp^s,iw)
1002  call minmod(ixi^l,id^l,a,b,tmp2)
1003  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
1004  dm4(id^s,iw) = tmp3(id^s)
1005  end do
1006 
1007  flux_ul(il^s,1:nwflux) = w(il^s,1:nwflux) + mpalpha * (w(il^s,1:nwflux) - w(ilm^s,1:nwflux))
1008  flux_md(il^s,1:nwflux) = half * (w(il^s,1:nwflux) + w(ilp^s,1:nwflux) - dm4(il^s,1:nwflux))
1009  flux_lc(il^s,1:nwflux) = half * (3.d0 * w(il^s,1:nwflux) - w(ilm^s,1:nwflux)) + mpbeta / 3.d0 * dm4(ilm^s,1:nwflux)
1010 
1011  flux_min(il^s,1:nwflux) = max(min(w(il^s,1:nwflux), w(ilp^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1012  min(w(il^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1013 
1014  flux_max(il^s,1:nwflux) = min(max(w(il^s,1:nwflux), w(ilp^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1015  max(w(il^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1016  do iw=1,nwflux
1017  a(il^s) = flux(il^s,iw)
1018  b(il^s) = flux_min(il^s,iw)
1019  c(il^s) = flux_max(il^s,iw)
1020  call median(ixi^l, il^l, a, b, c, tmp)
1021  wlc(il^s,iw) = tmp(il^s)
1022  end do
1023  end select
1024 
1025  !> right side
1026  !>> mmm -> pppp
1027  !>> mm -> ppp
1028  !>> m -> pp
1029  !>> 0 -> p
1030  !>> p -> 0
1031  !>> pp -> m
1032  !>> ppp -> mm
1033  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilpppp^s,1:nwflux) + u1_coeff(2) * w(ilppp^s,1:nwflux) + u1_coeff(3) * w(ilpp^s,1:nwflux) &
1034  + u1_coeff(4) * w(ilp^s,1:nwflux)
1035  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilppp^s,1:nwflux) + u2_coeff(2) * w(ilpp^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux) &
1036  + u2_coeff(4) * w(il^s,1:nwflux)
1037  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilpp^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(il^s,1:nwflux) &
1038  + u3_coeff(4) * w(ilm^s,1:nwflux)
1039  f_array(il^s,1:nwflux,4) = u4_coeff(1) * w(ilp^s,1:nwflux) + u4_coeff(2) * w(il^s,1:nwflux) + u4_coeff(3) * w(ilm^s,1:nwflux) &
1040  + u4_coeff(4) * w(ilmm^s,1:nwflux)
1041 
1042  beta(il^s,1:nwflux,1) = w(ilpppp^s,1:nwflux) * (547.d0 * w(ilpppp^s,1:nwflux) - 3882.d0 * w(ilppp^s,1:nwflux) + 4642.d0 * w(ilpp^s,1:nwflux) &
1043  - 1854.d0 * w(ilp^s,1:nwflux)) &
1044  + w(ilppp^s,1:nwflux) * (7043.d0 * w(ilppp^s,1:nwflux) - 17246.d0 * w(ilpp^s,1:nwflux) + 7042.d0 * w(ilp^s,1:nwflux)) &
1045  + w(ilpp^s,1:nwflux) * (11003.d0 * w(ilpp^s,1:nwflux) - 9402.d0 * w(ilp^s,1:nwflux)) + 2107.d0 * w(ilp^s,1:nwflux)**2
1046  beta(il^s,1:nwflux,2) = w(ilppp^s,1:nwflux) * (267.d0 * w(ilppp^s,1:nwflux) - 1642.d0 * w(ilpp^s,1:nwflux) + 1602.d0 * w(ilp^s,1:nwflux) &
1047  - 494.d0 * w(il^s,1:nwflux)) &
1048  + w(ilpp^s,1:nwflux) * (2843.d0 * w(ilpp^s,1:nwflux) - 5966.d0 * w(ilp^s,1:nwflux) + 1922.d0 * w(il^s,1:nwflux)) &
1049  + w(ilp^s,1:nwflux) * (3443.d0 * w(ilp^s,1:nwflux) - 2522.d0 * w(il^s,1:nwflux)) + 547.d0 * w(il^s,1:nwflux) ** 2
1050  beta(il^s,1:nwflux,3) = w(ilpp^s,1:nwflux) * (547.d0 * w(ilpp^s,1:nwflux) - 2522.d0 * w(ilp^s,1:nwflux) + 1922.d0 * w(il^s,1:nwflux) &
1051  - 494.d0 * w(ilm^s,1:nwflux)) &
1052  + w(ilp^s,1:nwflux) * (3443.d0 * w(ilp^s,1:nwflux) - 5966.d0 * w(il^s,1:nwflux) + 1602.d0 * w(ilm^s,1:nwflux)) &
1053  + w(il^s,1:nwflux) * (2843.d0 * w(il^s,1:nwflux) - 1642.d0 * w(ilm^s,1:nwflux)) + 267.d0 * w(ilm^s,1:nwflux) ** 2
1054  beta(il^s,1:nwflux,4) = w(ilp^s,1:nwflux) * (2107.d0 * w(ilp^s,1:nwflux) - 9402.d0 * w(il^s,1:nwflux) + 7042.d0 * w(ilm^s,1:nwflux) &
1055  - 1854.d0 * w(ilmm^s,1:nwflux)) &
1056  + w(il^s,1:nwflux) * (11003.d0 * w(il^s,1:nwflux) - 17246.d0 * w(ilm^s,1:nwflux) + 4642.d0 * w(ilmm^s,1:nwflux)) &
1057  + w(ilm^s,1:nwflux) * (7043.d0 * w(ilm^s,1:nwflux) - 3882.d0 * w(ilmm^s,1:nwflux)) + 547.d0 * w(ilmm^s,1:nwflux) ** 2
1058 
1059  alpha_sum(il^s,1:nwflux) = 0.0d0
1060  do i = 1,4
1061  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
1062  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
1063  end do
1064  flux(il^s,1:nwflux) = 0.0d0
1065  do i = 1,4
1066  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
1067  end do
1068 
1069  select case(var)
1070  case(1)
1071  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
1072  case(2)
1073  idmax^d=ilmax^d+kr(idims,^d); idmin^d=ilmin^d;
1074  idm^l=id^l-kr(idims,^d);
1075  idp^l=id^l+kr(idims,^d);
1076 
1077  iemax^d=idmax^d; iemin^d=idmin^d-kr(idims,^d);
1078  iem^l=ie^l-kr(idims,^d);
1079  iep^l=ie^l+kr(idims,^d);
1080  iepp^l=iep^l+kr(idims,^d);
1081 
1082  d(ie^s,1:nwflux) = w(ie^s,1:nwflux)-2.0d0*w(iep^s,1:nwflux)+w(iepp^s,1:nwflux)
1083 
1084  do iw = 1,nwflux
1085  a(id^s) = 4.0d0*d(id^s,iw)-d(idm^s,iw)
1086  b(id^s) = 4.0d0*d(idm^s,iw)-d(id^s,iw)
1087  call minmod(ixi^l,id^l,a,b,tmp)
1088  a(id^s) = d(id^s,iw)
1089  b(id^s) = d(idm^s,iw)
1090  call minmod(ixi^l,id^l,a,b,tmp2)
1091  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
1092  dm4(id^s,iw) = tmp3(id^s)
1093  end do
1094 
1095  flux_ul(il^s,1:nwflux) = w(ilp^s,1:nwflux) + mpalpha * (w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux))
1096  flux_md(il^s,1:nwflux) = half * (w(il^s,1:nwflux) + w(ilp^s,1:nwflux) - dm4(il^s,1:nwflux))
1097  flux_lc(il^s,1:nwflux) = half * (3.d0 * w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux)) + mpbeta / 3.d0 * dm4(ilp^s,1:nwflux)
1098 
1099  flux_min(il^s,1:nwflux) = max(min(w(ilp^s,1:nwflux), w(il^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1100  min(w(ilp^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1101 
1102  flux_max(il^s,1:nwflux) = min(max(w(ilp^s,1:nwflux), w(il^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1103  max(w(ilp^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1104  do iw=1,nwflux
1105  a(il^s) = flux(il^s,iw)
1106  b(il^s) = flux_min(il^s,iw)
1107  c(il^s) = flux_max(il^s,iw)
1108  call median(ixi^l, il^l, a, b, c, tmp)
1109  wrc(il^s,iw) = tmp(il^s)
1110  end do
1111  end select
1112  end subroutine weno7limiter
1113 
1114  subroutine exeno7limiter(ixI^L,iL^L,idims,w,wLC,wRC)
1116 
1117  integer, intent(in) :: ixI^L, iL^L, idims
1118  double precision, intent(in) :: w(ixi^s,1:nw)
1119  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
1120  !> local
1121  integer :: i, iw
1122  integer :: iLm^L, iLmm^L, iLmmm^L
1123  integer :: iLp^L, iLpp^L, iLppp^L, iLpppp^L
1124  integer :: iM^L, iMm^L, iMmm^L
1125  integer :: iMp^L, iMpp^L, iMppp^L
1126  integer :: id^L, idp^L, idpp^L, idm^L, ie^L, iem^L, iep^L, iepp^L
1127  integer, dimension(ixI^S,1:nw) :: delta_sum
1128  integer, dimension(ixI^S,1:nw,3):: delta
1129  double precision, dimension(2) :: beta_coeff
1130  double precision, dimension(3) :: d_array
1131  double precision, dimension(ixI^S,1:nw) :: gamma_sum, flux
1132  double precision, dimension(ixI^S,1:nw,3) :: beta, gamma_array, kai_array
1133  double precision, parameter :: exeno_ct = 1.d-1
1134  double precision, parameter :: weno_eps_machine = 1.d-18
1135 
1136  ilm^l=il^l-kr(idims,^d);
1137  ilmm^l=ilm^l-kr(idims,^d);
1138  ilmmm^l=ilmm^l-kr(idims,^d);
1139  ilp^l=il^l+kr(idims,^d);
1140  ilpp^l=ilp^l+kr(idims,^d);
1141  ilppp^l=ilpp^l+kr(idims,^d);
1142  ilpppp^l=ilppp^l+kr(idims,^d);
1143 
1144  immin^d=ilmin^d-kr(idims,^d);
1145  immax^d=ilmax^d+kr(idims,^d);
1146  imm^l=im^l-kr(idims,^d);
1147  immm^l=imm^l-kr(idims,^d);
1148  imp^l=im^l+kr(idims,^d);
1149  impp^l=imp^l+kr(idims,^d);
1150  imppp^l=impp^l+kr(idims,^d);
1151 
1152  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
1153  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
1154 
1155  !>> left side
1156  beta(im^s,1:nwflux,1) = beta_coeff(1) * (w(immm^s,1:nwflux) + w(im^s,1:nwflux) - 2.0d0 * w(imm^s,1:nwflux))**2 &
1157  + beta_coeff(2) * (w(immm^s,1:nwflux) - 4.0d0 * w(imm^s,1:nwflux) + 3.0d0 * w(im^s,1:nwflux))**2
1158  beta(im^s,1:nwflux,2) = beta_coeff(1) * (w(imm^s,1:nwflux) + w(imp^s,1:nwflux) - 2.0d0 * w(im^s,1:nwflux))**2 &
1159  + beta_coeff(2) * (w(imm^s,1:nwflux) - w(imp^s,1:nwflux))**2
1160  beta(im^s,1:nwflux,3) = beta_coeff(1) * (w(im^s,1:nwflux) + w(impp^s,1:nwflux) - 2.0d0 * w(imp^s,1:nwflux))**2 &
1161  + beta_coeff(2) * (3.0d0 * w(im^s, 1:nwflux) - 4.0d0 * w(imp^s,1:nwflux) + w(impp^s,1:nwflux))**2
1162 
1163  gamma_sum(im^s,1:nwflux) = 0.0d0
1164  do i = 1,3
1165  gamma_array(im^s,1:nwflux,i) = d_array(i) / (beta(im^s,1:nwflux,i) + weno_eps_machine)**2
1166  gamma_sum(im^s,1:nwflux) = gamma_sum(im^s,1:nwflux) + gamma_array(im^s,1:nwflux,i)
1167  end do
1168  do i = 1,3
1169  kai_array(im^s,1:nwflux,i) = gamma_array(im^s,1:nwflux,i) / gamma_sum(im^s,1:nwflux)
1170  where(kai_array(im^s,1:nwflux,i) .lt. exeno_ct)
1171  delta(im^s,1:nwflux,i) = 0
1172  elsewhere
1173  delta(im^s,1:nwflux,i) = 1
1174  endwhere
1175  end do
1176 
1177  delta_sum(il^s,1:nwflux) = delta(ilm^s,1:nwflux,1) * 1 + delta(ilp^s,1:nwflux,3) * 2 + delta(il^s,1:nwflux,1) * 4 &
1178  + delta(il^s,1:nwflux,2) * 8 + delta(il^s,1:nwflux,3) * 16
1179 
1180  !> f3
1181  where(delta_sum(il^s,1:nwflux) .eq. 31)
1182  flux(il^s,1:nwflux) = (- 3.d0 * w(ilmmm^s,1:nwflux) + 25.d0 * w(ilmm^s,1:nwflux) - 101.d0 * w(ilm^s,1:nwflux) + 319.d0 * w(il^s,1:nwflux) &
1183  + 214.d0 * w(ilp^s,1:nwflux) - 38.d0 * w(ilpp^s,1:nwflux) + 4.d0 * w(ilppp^s,1:nwflux)) / 420.d0
1184  !> f4
1185  elsewhere(delta_sum(il^s,1:nwflux) .eq. 30)
1186  flux(il^s,1:nwflux) = (w(ilmm^s,1:nwflux) - 8.d0 * w(ilm^s,1:nwflux) + 37.d0 * w(il^s,1:nwflux) + 37.d0 * w(ilp^s,1:nwflux) &
1187  - 8.d0 * w(ilpp^s,1:nwflux) + w(ilppp^s,1:nwflux)) / 60.d0
1188  !> f5
1189  elsewhere(delta_sum(il^s,1:nwflux) .eq. 29)
1190  flux(il^s,1:nwflux) = (- w(ilmmm^s,1:nwflux) + 7.d0 * w(ilmm^s,1:nwflux) - 23.d0 * w(ilm^s,1:nwflux) + 57.d0 * w(il^s,1:nwflux) &
1191  + 22.d0 * w(ilp^s,1:nwflux) - 2.d0 * w(ilpp^s,1:nwflux)) / 60.d0
1192  !> f6
1193  elsewhere(delta_sum(il^s,1:nwflux) .eq. 28)
1194  flux(il^s,1:nwflux) = (2.d0 * w(ilmm^s,1:nwflux) - 13.d0 * w(ilm^s,1:nwflux) + 47.d0 * w(il^s,1:nwflux) + 27.d0 * w(ilp^s,1:nwflux) &
1195  - 3.d0 * w(ilpp^s,1:nwflux)) / 60.d0
1196  !> f7
1197  elsewhere(delta_sum(il^s,1:nwflux) .ge. 24)
1198  flux(il^s,1:nwflux) = (- w(ilm^s,1:nwflux) + 7.d0 * w(il^s,1:nwflux) + 7.d0 * w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux)) / 12.d0
1199  !> f9
1200  elsewhere(delta_sum(il^s,1:nwflux) .ge. 16)
1201  flux(il^s,1:nwflux) = (2.d0 * w(il^s,1:nwflux) + 5.d0 * w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux)) / 6.d0
1202  !> f8
1203  elsewhere(delta_sum(il^s,1:nwflux) .ge. 12)
1204  flux(il^s,1:nwflux) = (w(ilmm^s,1:nwflux) - 5.d0 * w(ilm^s,1:nwflux) + 13.d0 * w(il^s,1:nwflux) + 3.d0 * w(ilp^s,1:nwflux)) / 12.d0
1205  !> f10
1206  elsewhere(delta_sum(il^s,1:nwflux) .ge. 8)
1207  flux(il^s,1:nwflux) = (- w(ilm^s,1:nwflux) + 5.d0 * w(il^s,1:nwflux) + 2.d0 * w(ilp^s,1:nwflux)) / 6.d0
1208  !> f11
1209  elsewhere
1210  flux(il^s,1:nwflux) = (2.d0 * w(ilmm^s,1:nwflux) - 7.d0 * w(ilm^s,1:nwflux) + 11.d0 * w(il^s,1:nwflux)) / 6.d0
1211  endwhere
1212 
1213  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
1214 
1215  !> right side
1216  beta(im^s,1:nwflux,1) = beta_coeff(1) * (w(imppp^s,1:nwflux) + w(imp^s,1:nwflux) - 2.0d0 * w(impp^s,1:nwflux))**2 &
1217  + beta_coeff(2) * (w(imppp^s,1:nwflux) - 4.0d0 * w(impp^s,1:nwflux) + 3.0d0 * w(imp^s,1:nwflux))**2
1218  beta(im^s,1:nwflux,2) = beta_coeff(1) * (w(impp^s,1:nwflux) + w(im^s,1:nwflux) - 2.0d0 * w(imp^s,1:nwflux))**2 &
1219  + beta_coeff(2) * (w(impp^s,1:nwflux) - w(im^s,1:nwflux))**2
1220  beta(im^s,1:nwflux,3) = beta_coeff(1) * (w(imp^s,1:nwflux) + w(imm^s,1:nwflux) - 2.0d0 * w(im^s,1:nwflux))**2 &
1221  + beta_coeff(2) * (3.0d0 * w(imp^s, 1:nwflux) - 4.0d0 * w(im^s,1:nwflux) + w(imm^s,1:nwflux))**2
1222 
1223  gamma_sum(im^s,1:nwflux) = 0.0d0
1224  do i = 1,3
1225  gamma_array(im^s,1:nwflux,i) = d_array(i) / (beta(im^s,1:nwflux,i) + weno_eps_machine)**2
1226  gamma_sum(im^s,1:nwflux) = gamma_sum(im^s,1:nwflux) + gamma_array(im^s,1:nwflux,i)
1227  end do
1228  do i = 1,3
1229  kai_array(im^s,1:nwflux,i) = gamma_array(im^s,1:nwflux,i) / gamma_sum(im^s,1:nwflux)
1230  where(kai_array(im^s,1:nwflux,i) .lt. exeno_ct)
1231  delta(im^s,1:nwflux,i) = 0
1232  elsewhere
1233  delta(im^s,1:nwflux,i) = 1
1234  endwhere
1235  end do
1236 
1237  delta_sum(il^s,1:nwflux) = delta(ilp^s,1:nwflux,1) * 1 + delta(ilm^s,1:nwflux,3) * 2 + delta(il^s,1:nwflux,1) * 4 &
1238  + delta(il^s,1:nwflux,2) * 8 + delta(il^s,1:nwflux,3) * 16
1239 
1240  where(delta_sum(il^s,1:nwflux) .eq. 31)
1241  flux(il^s,1:nwflux) = (- 3.d0 * w(ilpppp^s,1:nwflux) + 25.d0 * w(ilppp^s,1:nwflux) - 101.d0 * w(ilpp^s,1:nwflux) &
1242  + 319.d0 * w(ilp^s,1:nwflux) + 214.d0 * w(il^s,1:nwflux) - 38.d0 * w(ilm^s,1:nwflux) &
1243  + 4.d0 * w(ilmm^s,1:nwflux)) / 420.d0
1244  elsewhere(delta_sum(il^s,1:nwflux) .eq. 30)
1245  flux(il^s,1:nwflux) = (w(ilppp^s,1:nwflux) - 8.d0 * w(ilpp^s,1:nwflux) + 37.d0 * w(ilp^s,1:nwflux) + 37.d0 * w(il^s,1:nwflux) &
1246  - 8.d0 * w(ilm^s,1:nwflux) + w(ilmm^s,1:nwflux)) / 60.d0
1247  elsewhere(delta_sum(il^s,1:nwflux) .eq. 29)
1248  flux(il^s,1:nwflux) = (- w(ilpppp^s,1:nwflux) + 7.d0 * w(ilppp^s,1:nwflux) - 23.d0 * w(ilpp^s,1:nwflux) + 57.d0 * w(ilp^s,1:nwflux) &
1249  + 22.d0 * w(il^s,1:nwflux) - 2.d0 * w(ilm^s,1:nwflux)) / 60.d0
1250  elsewhere(delta_sum(il^s,1:nwflux) .eq. 28)
1251  flux(il^s,1:nwflux) = (2.d0 * w(ilppp^s,1:nwflux) - 13.d0 * w(ilpp^s,1:nwflux) + 47.d0 * w(ilp^s,1:nwflux) + 27.d0 * w(il^s,1:nwflux) &
1252  - 3.d0 * w(ilm^s,1:nwflux)) / 60.d0
1253  elsewhere(delta_sum(il^s,1:nwflux) .ge. 24)
1254  flux(il^s,1:nwflux) = (- w(ilpp^s,1:nwflux) + 7.d0 * w(ilp^s,1:nwflux) + 7.d0 * w(il^s,1:nwflux) - w(ilm^s,1:nwflux)) / 12.d0
1255  elsewhere(delta_sum(il^s,1:nwflux) .ge. 16)
1256  flux(il^s,1:nwflux) = (2.d0 * w(ilp^s,1:nwflux) + 5.d0 * w(il^s,1:nwflux) - w(ilm^s,1:nwflux)) / 6.d0
1257  elsewhere(delta_sum(il^s,1:nwflux) .ge. 12)
1258  flux(il^s,1:nwflux) = (w(ilppp^s,1:nwflux) - 5.d0 * w(ilpp^s,1:nwflux) + 13.d0 * w(ilp^s,1:nwflux) + 3.d0 * w(il^s,1:nwflux)) / 12.d0
1259  elsewhere(delta_sum(il^s,1:nwflux) .ge. 8)
1260  flux(il^s,1:nwflux) = (- w(ilpp^s,1:nwflux) + 5.d0 * w(ilp^s,1:nwflux) + 2.d0 * w(il^s,1:nwflux)) / 6.d0
1261  elsewhere
1262  flux(il^s,1:nwflux) = (2.d0 * w(ilppp^s,1:nwflux) - 7.d0 * w(ilpp^s,1:nwflux) + 11.d0 * w(ilp^s,1:nwflux)) / 6.d0
1263  endwhere
1264 
1265  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
1266 
1267  end subroutine exeno7limiter
1268 
1269  subroutine minmod(ixI^L,ixO^L,a,b,minm)
1272 
1273  integer, intent(in) :: ixI^L, ixO^L
1274  double precision, intent(in) :: a(ixi^s), b(ixi^s)
1275  double precision, intent(out):: minm(ixi^s)
1276 
1277  minm(ixo^s) = (sign(one,a(ixo^s))+sign(one,b(ixo^s)))/2.0d0 * &
1278  min(abs(a(ixo^s)),abs(b(ixo^s)))
1279 
1280  end subroutine minmod
1281 
1282  subroutine median(ixI^L,ixO^L,a,b,c,med)
1285 
1286  integer, intent(in) :: ixI^L, ixO^L
1287  double precision, intent(in) :: a(ixi^s), b(ixi^s), c(ixi^s)
1288  double precision, intent(out):: med(ixi^s)
1289 
1290  double precision :: tmp1(ixi^s),tmp2(ixi^s)
1291 
1292  tmp1(ixo^s) = b(ixo^s) - a(ixo^s); tmp2(ixo^s) = c(ixo^s) - a(ixo^s)
1293 
1294  med(ixo^s) = a(ixo^s) + (sign(one,tmp1(ixo^s))+sign(one,tmp2(ixo^s)))/2.0d0 * &
1295  min(abs(tmp1(ixo^s)),abs(tmp2(ixo^s)))
1296 
1297  end subroutine median
1298 end module mod_weno
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine, public weno3limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:34
subroutine, public weno5nmlimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:361
subroutine, public teno5adlimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC)
Definition: mod_weno.t:248
subroutine, public exeno7limiter(ixIL, iLL, idims, w, wLC, wRC)
Definition: mod_weno.t:1115
subroutine, public weno5limiterl(ixIL, iLL, idims, w, wLC, var)
Definition: mod_weno.t:489
subroutine median(ixIL, ixOL, a, b, c, med)
Definition: mod_weno.t:1283
subroutine minmod(ixIL, ixOL, a, b, minm)
Definition: mod_weno.t:1270
subroutine, public weno5limiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:572
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
subroutine, public weno5cu6limiter(ixIL, iLL, idims, w, wLC, wRC)
Definition: mod_weno.t:821
subroutine, public weno7limiter(ixIL, iLL, idims, w, wLC, wRC, var)
Definition: mod_weno.t:910
subroutine, public weno5nmlimiterl(ixIL, iLL, idims, w, wLC, var)
Definition: mod_weno.t:654
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine, public weno5limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:122
subroutine, public weno5nmlimiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:736