24 integer,
intent(in) :: ixi^
l, ix^
l, idims
25 double precision,
intent(in) :: q(ixi^s),qct(ixi^s)
27 double precision,
intent(inout) :: qrc(ixi^s),qlc(ixi^s)
29 double precision,
dimension(ixI^S) :: dqc,d2qc,ldq
30 double precision,
dimension(ixI^S) :: qmin,qmax,tmp
32 integer :: lxc^
l,lxr^
l,lxl^
l
33 integer :: ixl^
l,ixo^
l,ixr^
l,ixol^
l,ixor^
l
34 integer :: hxl^
l,hxc^
l,hxr^
l
35 integer :: kxl^
l,kxc^
l,kxr^
l
37 ixomin^
d=ixmin^
d-
kr(idims,^
d);ixomax^
d=ixmax^
d+
kr(idims,^
d);
38 ixolmin^
d=ixomin^
d-
kr(idims,^
d);ixolmax^
d=ixomax^
d;
39 ixor^
l=ixol^
l+
kr(idims,^
d);
40 ixl^
l=ixo^
l-
kr(idims,^
d);
41 ixr^
l=ixo^
l+
kr(idims,^
d);
43 hxcmin^
d=ixomin^
d;hxcmax^
d=ixmax^
d;
44 hxl^
l=hxc^
l-
kr(idims,^
d);
45 hxr^
l=hxc^
l+
kr(idims,^
d);
47 kxcmin^
d=ixlmin^
d-1; kxcmax^
d=ixrmax^
d;
48 kxr^
l=kxc^
l+
kr(idims,^
d);
50 lxcmin^
d=ixomin^
d-
kr(idims,^
d);lxcmax^
d=ixomax^
d+
kr(idims,^
d);
51 lxl^
l=lxc^
l-
kr(idims,^
d);
52 lxr^
l=lxc^
l+
kr(idims,^
d);
54 dqc(kxc^s)=q(kxr^s)-q(kxc^s)
56 d2qc(kxc^s)=half*(q(lxr^s)-q(lxl^s))
57 where(dqc(lxc^s)*dqc(lxl^s)>zero)
59 qmin(kxc^s)= sign(one,d2qc(lxc^s))
61 ldq(lxc^s)= qmin(lxc^s)*min(dabs(d2qc(lxc^s)),&
62 2.0d0*dabs(dqc(lxl^s)),&
63 2.0d0*dabs(dqc(lxc^s)))
69 qlc(ixol^s)=qlc(ixol^s)+half*dqc(ixol^s)&
70 +(ldq(ixol^s)-ldq(ixor^s))/6.0d0
74 call extremaq(ixi^
l,ixo^
l,qct,1,qmax,qmin)
77 qrc(ixl^s)=max(qmin(ixo^s),min(qmax(ixo^s),qrc(ixl^s)))
78 qlc(ixo^s)=max(qmin(ixo^s),min(qmax(ixo^s),qlc(ixo^s)))
81 where((qrc(ixl^s)-qct(ixo^s))*(qct(ixo^s)-qlc(ixo^s))<=zero)
86 qmax(ixo^s)=(qlc(ixo^s)-qrc(ixl^s))&
87 *(qct(ixo^s)-half*(qlc(ixo^s)+qrc(ixl^s)))
88 qmin(ixo^s)=(qlc(ixo^s)-qrc(ixl^s))**2/6.0d0
92 where(qmax(hxr^s)>qmin(hxr^s))
93 qrc(hxc^s)= 3.0d0*qct(hxr^s)-2.0d0*qlc(hxr^s)
96 where(qmax(hxc^s)<-qmin(hxc^s))
97 qlc(hxc^s)= 3.0d0*qct(hxc^s)-2.0d0*tmp(hxl^s)
114 integer,
intent(in) :: ixi^
l, ix^
l, idims
115 double precision,
intent(in) :: w(ixi^s,1:nw),wct(ixi^s,1:nw)
117 double precision,
intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
119 double precision,
dimension(ixI^S,1:nw_recon) :: dwc,d2wc,ldw
120 double precision,
dimension(ixI^S,1:nw_recon) :: wmin,wmax,tmp
121 double precision,
dimension(ixI^S) :: aa, ab, ac, ki
123 double precision,
parameter :: betamin=0.75d0, betamax=0.85d0,&
124 zmin=0.25d0, zmax=0.75d0,&
125 eta1=20.0d0,eta2=0.05d0,eps=0.01d0,kappa=0.1d0
126 integer :: lxc^
l,lxl^
l,lxr^
l
127 integer :: ixll^
l,ixl^
l,ixo^
l,ixr^
l,ixrr^
l,ixol^
l,ixor^
l
128 integer :: hxl^
l,hxc^
l,hxr^
l
129 integer :: kxll^
l,kxl^
l,kxc^
l,kxr^
l,kxrr^
l
130 integer :: iw, idimss
132 ixomin^
d=ixmin^
d-
kr(idims,^
d);ixomax^
d=ixmax^
d+
kr(idims,^
d);
133 ixolmin^
d=ixomin^
d-
kr(idims,^
d);ixolmax^
d=ixomax^
d;
134 ixor^
l=ixol^
l+
kr(idims,^
d);
135 ixl^
l=ixo^
l-
kr(idims,^
d);
136 ixr^
l=ixo^
l+
kr(idims,^
d);
138 hxcmin^
d=ixomin^
d;hxcmax^
d=ixmax^
d;
139 hxl^
l=hxc^
l-
kr(idims,^
d);
140 hxr^
l=hxc^
l+
kr(idims,^
d);
142 kxcmin^
d=ixlmin^
d-1; kxcmax^
d=ixrmax^
d;
143 kxr^
l=kxc^
l+
kr(idims,^
d);
145 lxcmin^
d=ixomin^
d-
kr(idims,^
d);lxcmax^
d=ixomax^
d+
kr(idims,^
d);
146 lxl^
l=lxc^
l-
kr(idims,^
d);
147 lxr^
l=lxc^
l+
kr(idims,^
d);
149 dwc(kxc^s,1:nw_recon)=w(kxr^s,1:nw_recon)-w(kxc^s,1:nw_recon)
151 d2wc(lxc^s,1:nw_recon)=half*(w(lxr^s,1:nw_recon)-w(lxl^s,1:nw_recon))
152 where(dwc(lxc^s,1:nw_recon)*dwc(lxl^s,1:nw_recon)>zero)
154 wmin(lxc^s,1:nw_recon)= sign(one,d2wc(lxc^s,1:nw_recon))
156 ldw(lxc^s,1:nw_recon)= wmin(lxc^s,1:nw_recon)*min(dabs(d2wc(lxc^s,1:nw_recon)),&
157 2.0d0*dabs(dwc(lxl^s,1:nw_recon)),&
158 2.0d0*dabs(dwc(lxc^s,1:nw_recon)))
160 ldw(lxc^s,1:nw_recon)=zero
164 wlc(ixol^s,1:nw_recon)=wlc(ixol^s,1:nw_recon)+half*dwc(ixol^s,1:nw_recon)&
165 +(ldw(ixol^s,1:nw_recon)-ldw(ixor^s,1:nw_recon))/6.0d0
166 wrc(ixl^s,1:nw_recon)=wlc(ixl^s,1:nw_recon)
169 call extremaw(ixi^
l,ixo^
l,w,1,wmax,wmin)
172 wrc(ixl^s,1:nw_recon)=max(wmin(ixo^s,1:nw_recon)&
173 ,min(wmax(ixo^s,1:nw_recon),wrc(ixl^s,1:nw_recon)))
174 wlc(ixo^s,1:nw_recon)=max(wmin(ixo^s,1:nw_recon)&
175 ,min(wmax(ixo^s,1:nw_recon),wlc(ixo^s,1:nw_recon)))
178 where((wrc(ixl^s,1:nw_recon)-wct(ixo^s,1:nw_recon))&
179 *(wct(ixo^s,1:nw_recon)-wlc(ixo^s,1:nw_recon))<=zero)
180 wrc(ixl^s,1:nw_recon)=wct(ixo^s,1:nw_recon)
181 wlc(ixo^s,1:nw_recon)=wct(ixo^s,1:nw_recon)
184 wmax(ixo^s,1:nw_recon)=(wlc(ixo^s,1:nw_recon)-wrc(ixl^s,1:nw_recon))*&
185 (wct(ixo^s,1:nw_recon)-half*(wlc(ixo^s,1:nw_recon)+wrc(ixl^s,1:nw_recon)))
186 wmin(ixo^s,1:nw_recon)=(wlc(ixo^s,1:nw_recon)-wrc(ixl^s,1:nw_recon))**2/6.0d0
187 tmp(hxl^s,1:nw_recon)=wrc(hxl^s,1:nw_recon)
189 where(wmax(hxr^s,1:nw_recon)>wmin(hxr^s,1:nw_recon))
190 wrc(hxc^s,1:nw_recon)= 3.0d0*wct(hxr^s,1:nw_recon)&
191 -2.0d0*wlc(hxr^s,1:nw_recon)
194 where(wmax(hxc^s,1:nw_recon)<-wmin(hxc^s,1:nw_recon))
195 wlc(hxc^s,1:nw_recon)= 3.0d0*wct(hxc^s,1:nw_recon)&
196 -2.0d0*tmp(hxl^s,1:nw_recon)
201 ixrr^
l=ixr^
l+
kr(idims,^
d);
202 kxl^
l=kxc^
l-
kr(idims,^
d);
203 call ppm_flatcd(ixi^
l,kxc^
l,kxl^
l,kxr^
l,wct,d2wc,aa,ab)
204 if(any(kappa*aa(kxc^s)>=ab(kxc^s)))
then
206 where(kappa*aa(kxc^s)>=ab(kxc^s).and. dabs(dwc(kxc^s,iw))>smalldouble)
207 wmax(kxc^s,iw) = wct(kxr^s,iw)-2.0d0*wct(kxc^s,iw)+wct(kxl^s,iw)
210 where(wmax(ixr^s,iw)*wmax(ixl^s,iw)<zero .and.&
211 dabs(wct(ixr^s,iw)-wct(ixl^s,iw))&
212 -eps*min(dabs(wct(ixr^s,iw)),dabs(wct(ixl^s,iw)))>zero &
213 .and. kappa*aa(ixo^s)>=ab(ixo^s)&
214 .and. dabs(dwc(ixo^s,iw))>smalldouble)
216 ac(ixo^s)=(wct(ixll^s,iw)-wct(ixrr^s,iw)+4.0d0*dwc(ixo^s,iw))&
217 /(12.0d0*dwc(ixo^s,iw))
218 wmin(ixo^s,iw)=max(zero,min(eta1*(ac(ixo^s)-eta2),one))
223 where(wmin(hxc^s,iw)>zero)
224 wlc(hxc^s,iw) = wlc(hxc^s,iw)*(one-wmin(hxc^s,iw))&
225 +(wct(hxc^s,iw)+half*ldw(hxc^s,iw))*wmin(hxc^s,iw)
227 where(wmin(hxr^s,iw)>zero)
228 wrc(hxc^s,iw) = wrc(hxc^s,iw)*(one-wmin(hxr^s,iw))&
229 +(wct(hxr^s,iw)-half*ldw(hxr^s,iw))*wmin(hxr^s,iw)
238 kxcmin^
d=ixmin^
d-2; kxcmax^
d=ixmax^
d+2;
241 kxl^
l=kxc^
l-
kr(idimss,^
d);
242 kxr^
l=kxc^
l+
kr(idimss,^
d);
243 kxll^
l=kxl^
l-
kr(idimss,^
d);
244 kxrr^
l=kxr^
l+
kr(idimss,^
d);
246 call ppm_flatsh(ixi^
l,kxc^
l,kxll^
l,kxl^
l,kxr^
l,kxrr^
l,idimss,wct,aa,ab)
249 ac(kxc^s) = max(zero,min(one,(betamax-aa(kxc^s))/(betamax-betamin)))
252 where(wct(kxr^s,iw_mom(idimss))<wct(kxl^s,iw_mom(idimss)))
253 aa(kxc^s) = max(ac(kxc^s), min(one,(zmax-ab(kxc^s))/(zmax-zmin)))
257 ixl^
l=ixo^
l-
kr(idimss,^
d);
258 ixr^
l=ixo^
l+
kr(idimss,^
d);
259 ki(ixo^s)=min(ki(ixo^s),aa(ixl^s),aa(ixo^s),aa(ixr^s))
263 where(dabs(ab(ixo^s)-one)>smalldouble)
264 wmax(ixo^s,iw) = (one-ab(ixo^s))*wct(ixo^s,iw)
267 where(dabs(ab(hxc^s)-one)>smalldouble)
268 wlc(hxc^s,iw) = ab(hxc^s)*wlc(hxc^s,iw)+wmax(hxc^s,iw)
271 where(dabs(ab(hxr^s)-one)>smalldouble)
272 wrc(hxc^s,iw) = ab(hxr^s)*wrc(hxc^s,iw)+wmax(hxr^s,iw)