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