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