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)
42  use mod_global_parameters
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)
74  use mod_global_parameters
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)
101  use mod_global_parameters
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)
136  use mod_global_parameters
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)
166  use mod_global_parameters
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)
257  use mod_global_parameters
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)
384  use mod_global_parameters
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)
501  use mod_global_parameters
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)
632  use mod_global_parameters
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)
717  use mod_global_parameters
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)
801  use mod_global_parameters
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)
885  use mod_global_parameters
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)
969  use mod_global_parameters
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)
1063  use mod_global_parameters
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 
1273  use mod_global_parameters
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 
1286  use mod_global_parameters
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
mod_global_parameters
This module contains definitions of global parameters and variables and some generic functions/subrou...
Definition: mod_global_parameters.t:4
mod_global_parameters::block
type(state), pointer block
Block pointer for using one block and its previous state.
Definition: mod_global_parameters.t:755
mod_global_parameters::kr
integer, dimension(3, 3) kr
Kronecker delta tensor.
Definition: mod_global_parameters.t:449
mod_global_parameters::ndim
integer, parameter ndim
Number of spatial dimensions for grid variables.
Definition: mod_global_parameters.t:62
mod_global_parameters::d
integer, dimension(:), allocatable, parameter d
Definition: mod_global_parameters.t:88
mod_global_parameters::l
double precision l
Definition: mod_global_parameters.t:54
mod_physics
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
mod_physics::phys_check_w
procedure(sub_check_w), pointer phys_check_w
Definition: mod_physics.t:76
mod_weno
Definition: mod_weno.t:1
mod_weno::fix_onelimiter
subroutine, public fix_onelimiter(ixIL, iLL, wCin, wCout)
Definition: mod_weno.t:42
mod_weno::weno5limiterl
subroutine, public weno5limiterl(ixIL, iLL, idims, w, wLC, var)
Definition: mod_weno.t:632
mod_weno::fix_onelimiter1
subroutine, public fix_onelimiter1(ixIL, iLL, wCin, wCout)
Definition: mod_weno.t:74
mod_weno::weno5nmlimiter
subroutine, public weno5nmlimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:501
mod_weno::teno5adlimiter
subroutine, public teno5adlimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC)
Definition: mod_weno.t:384
mod_weno::weno5limiterr
subroutine, public weno5limiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:717
mod_weno::weno5nmlimiterl
subroutine, public weno5nmlimiterl(ixIL, iLL, idims, w, wLC, var)
Definition: mod_weno.t:801
mod_weno::fix_limiter
subroutine, public fix_limiter(ixIL, iLL, wLCin, wRCin, wLCout, wRCout)
Definition: mod_weno.t:101
mod_weno::weno5cu6limiter
subroutine, public weno5cu6limiter(ixIL, iLL, idims, w, wLC, wRC)
Definition: mod_weno.t:969
mod_weno::weno7limiter
subroutine, public weno7limiter(ixIL, iLL, idims, w, wLC, wRC, var)
Definition: mod_weno.t:1063
mod_weno::fix_limiter1
subroutine, public fix_limiter1(ixIL, iLL, wLCin, wRCin, wLCout, wRCout)
Definition: mod_weno.t:136
mod_weno::weno3limiter
subroutine, public weno3limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:166
mod_weno::minmod
subroutine minmod(ixIL, ixOL, a, b, minm)
Definition: mod_weno.t:1272
mod_weno::weno5limiter
subroutine, public weno5limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:257
mod_weno::weno5nmlimiterr
subroutine, public weno5nmlimiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:885
mod_weno::median
subroutine median(ixIL, ixOL, a, b, c, med)
Definition: mod_weno.t:1285
Generated by doxygen 1.9.1