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:nwflux)  :: dwc,d2wc,ldw
 
  120    double precision,
dimension(ixI^S,1:nwflux)  :: 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-
kr(idims,^
d); 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:nwflux)=w(kxr^s,1:nwflux)-w(kxc^s,1:nwflux)
 
  151    d2wc(lxc^s,1:nwflux)=half*(w(lxr^s,1:nwflux)-w(lxl^s,1:nwflux))
 
  152    where(dwc(lxc^s,1:nwflux)*dwc(lxl^s,1:nwflux)>zero)
 
  154       wmin(lxc^s,1:nwflux)= sign(one,d2wc(lxc^s,1:nwflux))
 
  156       ldw(lxc^s,1:nwflux)= wmin(lxc^s,1:nwflux)*min(dabs(d2wc(lxc^s,1:nwflux)),&
 
  157            2.0d0*dabs(dwc(lxl^s,1:nwflux)),&
 
  158            2.0d0*dabs(dwc(lxc^s,1:nwflux)))
 
  160       ldw(lxc^s,1:nwflux)=zero
 
  164    wlc(ixol^s,1:nwflux)=wlc(ixol^s,1:nwflux)+half*dwc(ixol^s,1:nwflux)&
 
  165         +(ldw(ixol^s,1:nwflux)-ldw(ixor^s,1:nwflux))/6.0d0
 
  166    wrc(ixl^s,1:nwflux)=wlc(ixl^s,1:nwflux)
 
  169    call extremaw(ixi^
l,ixo^
l,w,1,wmax,wmin)
 
  172    wrc(ixl^s,1:nwflux)=max(wmin(ixo^s,1:nwflux)&
 
  173         ,min(wmax(ixo^s,1:nwflux),wrc(ixl^s,1:nwflux))) 
 
  174    wlc(ixo^s,1:nwflux)=max(wmin(ixo^s,1:nwflux)&
 
  175         ,min(wmax(ixo^s,1:nwflux),wlc(ixo^s,1:nwflux)))
 
  178    where((wrc(ixl^s,1:nwflux)-wct(ixo^s,1:nwflux))&
 
  179         *(wct(ixo^s,1:nwflux)-wlc(ixo^s,1:nwflux))<=zero)
 
  180       wrc(ixl^s,1:nwflux)=wct(ixo^s,1:nwflux)
 
  181       wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
 
  184    wmax(ixo^s,1:nwflux)=(wlc(ixo^s,1:nwflux)-wrc(ixl^s,1:nwflux))*&
 
  185         (wct(ixo^s,1:nwflux)-half*(wlc(ixo^s,1:nwflux)+wrc(ixl^s,1:nwflux)))
 
  186    wmin(ixo^s,1:nwflux)=(wlc(ixo^s,1:nwflux)-wrc(ixl^s,1:nwflux))**2/6.0d0
 
  187    tmp(hxl^s,1:nwflux)=wrc(hxl^s,1:nwflux)
 
  189    where(wmax(hxr^s,1:nwflux)>wmin(hxr^s,1:nwflux))
 
  190       wrc(hxc^s,1:nwflux)= 3.0d0*wct(hxr^s,1:nwflux)&
 
  191            -2.0d0*wlc(hxr^s,1:nwflux)
 
  194    where(wmax(hxc^s,1:nwflux)<-wmin(hxc^s,1:nwflux))
 
  195       wlc(hxc^s,1:nwflux)= 3.0d0*wct(hxc^s,1:nwflux)&
 
  196            -2.0d0*tmp(hxl^s,1:nwflux)
 
  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)