mod_ppm.t

Go to the documentation of this file.

module mod_ppm
 
  implicit none
  private
 
  public :: ppmlimiter
  public :: ppmlimitervar
 
contains
 
  subroutine ppmlimitervar(ixI^L,ix^L,idims,q,qCT,qLC,qRC)
 
    ! references:
    ! Mignone et al 2005, ApJS 160, 199,
    ! Miller and Colella 2002, JCP 183, 26
    ! Fryxell et al. 2000 ApJ, 131, 273 (Flash)
    ! baciotti Phd (http://www.aei.mpg.de/~baiotti/Baiotti_PhD.pdf)
    ! old version : april 2009
    ! author: zakaria.meliani@wis.kuleuven.be
    ! current version : March 2023 by Chun Xia
 
    use mod_global_parameters
 
    integer, intent(in)             :: ixi^l, ix^l, idims
    double precision, intent(in)    :: q(ixi^s),qct(ixi^s)
 
    double precision, intent(inout) :: qrc(ixi^s),qlc(ixi^s)
 
    double precision,dimension(ixI^S)  :: dqc,d2qc,ldq
    double precision,dimension(ixI^S)  :: qmin,qmax,tmp
 
    integer   :: lxc^l,lxr^l,lxl^l
    integer   :: ixl^l,ixo^l,ixr^l,ixol^l,ixor^l
    integer   :: hxl^l,hxc^l,hxr^l
    integer   :: kxl^l,kxc^l,kxr^l
 
    ixomin^d=ixmin^d-kr(idims,^d);ixomax^d=ixmax^d+kr(idims,^d);!ixO[ixMmin-1,ixMmax+1]
    ixolmin^d=ixomin^d-kr(idims,^d);ixolmax^d=ixomax^d;         !ixOL[ixMmin-2,ixMmax+1]
    ixor^l=ixol^l+kr(idims,^d);                                 !ixOR=[iMmin-1,ixMmax+2]
    ixl^l=ixo^l-kr(idims,^d);                                   !ixL[ixMmin-2,ixMmax]
    ixr^l=ixo^l+kr(idims,^d);                                   !ixR=[iMmin,ixMmax+2]
 
    hxcmin^d=ixomin^d;hxcmax^d=ixmax^d; ! hxC = [ixMmin-1,ixMmax]
    hxl^l=hxc^l-kr(idims,^d);           ! hxL = [ixMmin-2,ixMmax-1]
    hxr^l=hxc^l+kr(idims,^d);           ! hxR = [ixMmin,ixMmax+1]
 
    kxcmin^d=ixlmin^d-1; kxcmax^d=ixrmax^d; ! kxC=[iMmin-3,ixMmax+2]
    kxr^l=kxc^l+kr(idims,^d);              ! kxR=[iMmin-2,ixMmax+3]
 
    lxcmin^d=ixomin^d-kr(idims,^d);lxcmax^d=ixomax^d+kr(idims,^d);! lxC=[iMmin-2,ixMmax+2]
    lxl^l=lxc^l-kr(idims,^d);                          ! lxL=[iMmin-3,ixMmax+1]
    lxr^l=lxc^l+kr(idims,^d);                          ! lxR=[iMmin-1,ixMmax+3]
 
    dqc(kxc^s)=q(kxr^s)-q(kxc^s)
    ! Eq. 64,  Miller and Colella 2002, JCP 183, 26 
    d2qc(kxc^s)=half*(q(lxr^s)-q(lxl^s))
    where(dqc(lxc^s)*dqc(lxl^s)>zero)
       ! Store the sign of dwC in wMin
       qmin(kxc^s)= sign(one,d2qc(lxc^s))
       ! Eq. 65,  Miller and Colella 2002, JCP 183, 26 
       ldq(lxc^s)= qmin(lxc^s)*min(dabs(d2qc(lxc^s)),&
            2.0d0*dabs(dqc(lxl^s)),&
            2.0d0*dabs(dqc(lxc^s)))
    else where
       ldq(lxc^s)=zero
    end where
 
    ! Eq. 66, Miller and Colella 2002, JCP 183, 26
    qlc(ixol^s)=qlc(ixol^s)+half*dqc(ixol^s)&
         +(ldq(ixol^s)-ldq(ixor^s))/6.0d0
    qrc(ixl^s)=qlc(ixl^s)
 
    ! make sure that min wCT(i)<wLC(i)<wCT(i+1) same for wRC(i)
    call extremaq(ixi^l,ixo^l,qct,1,qmax,qmin)
 
    ! Eq. B8, page 217, Mignone et al 2005, ApJS
    qrc(ixl^s)=max(qmin(ixo^s),min(qmax(ixo^s),qrc(ixl^s)))
    qlc(ixo^s)=max(qmin(ixo^s),min(qmax(ixo^s),qlc(ixo^s)))
 
    ! Eq. B9, page 217, Mignone et al 2005, ApJS
    where((qrc(ixl^s)-qct(ixo^s))*(qct(ixo^s)-qlc(ixo^s))<=zero)
       qrc(ixl^s)=qct(ixo^s)
       qlc(ixo^s)=qct(ixo^s)
    end where
 
    qmax(ixo^s)=(qlc(ixo^s)-qrc(ixl^s))&
         *(qct(ixo^s)-half*(qlc(ixo^s)+qrc(ixl^s)))
    qmin(ixo^s)=(qlc(ixo^s)-qrc(ixl^s))**2/6.0d0
 
    tmp(hxl^s)=qrc(hxl^s)
    ! Eq. B10, page 218, Mignone et al 2005, ApJS
    where(qmax(hxr^s)>qmin(hxr^s))
       qrc(hxc^s)= 3.0d0*qct(hxr^s)-2.0d0*qlc(hxr^s)
    end where
    ! Eq. B11, page 218, Mignone et al 2005, ApJS
    where(qmax(hxc^s)<-qmin(hxc^s))
       qlc(hxc^s)= 3.0d0*qct(hxc^s)-2.0d0*tmp(hxl^s)
    end where
 
  end subroutine ppmlimitervar
 
  subroutine ppmlimiter(ixI^L,ix^L,idims,w,wCT,wLC,wRC)
 
    ! references:
    ! Mignone et al 2005, ApJS 160, 199, 
    ! Miller and Colella 2002, JCP 183, 26 
    ! Fryxell et al. 2000 ApJ, 131, 273 (Flash)
    ! baciotti Phd (http://www.aei.mpg.de/~baiotti/Baiotti_PhD.pdf)
    ! old version : april 2009
    ! author: zakaria.meliani@wis.kuleuven.be
    ! current version : March 2023 by Chun Xia
    use mod_global_parameters
 
    integer, intent(in)             :: ixi^l, ix^l, idims
    double precision, intent(in)    :: w(ixi^s,1:nw),wct(ixi^s,1:nw)
 
    double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw) 
 
    double precision,dimension(ixI^S,1:nwflux)  :: dwc,d2wc,ldw
    double precision,dimension(ixI^S,1:nwflux)  :: wmin,wmax,tmp
    double precision,dimension(ixI^S) :: aa, ab, ac, ki
 
    integer   :: lxc^l,lxl^l,lxr^l
    integer   :: ixll^l,ixl^l,ixo^l,ixr^l,ixrr^l,ixol^l,ixor^l
    integer   :: hxl^l,hxc^l,hxr^l
    integer   :: kxll^l,kxl^l,kxc^l,kxr^l,kxrr^l
    integer   :: iw, idimss
 
    double precision, parameter :: betamin=0.75d0, betamax=0.85d0,&
         zmin=0.25d0, zmax=0.75d0,&
         eta1=20.0d0,eta2=0.05d0,eps=0.01d0,kappa=0.1d0
 
    ixomin^d=ixmin^d-kr(idims,^d);ixomax^d=ixmax^d+kr(idims,^d);!ixO[ixMmin-1,ixMmax+1]
    ixolmin^d=ixomin^d-kr(idims,^d);ixolmax^d=ixomax^d;         !ixOL[ixMmin-2,ixMmax+1]
    ixor^l=ixol^l+kr(idims,^d);                                 !ixOR=[iMmin-1,ixMmax+2]
    ixl^l=ixo^l-kr(idims,^d);                                   !ixL[ixMmin-2,ixMmax]
    ixr^l=ixo^l+kr(idims,^d);                                   !ixR=[iMmin,ixMmax+2]
 
    hxcmin^d=ixomin^d;hxcmax^d=ixmax^d; ! hxC = [ixMmin-1,ixMmax]
    hxl^l=hxc^l-kr(idims,^d);           ! hxL = [ixMmin-2,ixMmax-1]
    hxr^l=hxc^l+kr(idims,^d);           ! hxR = [ixMmin,ixMmax+1]
 
    kxcmin^d=ixlmin^d-1; kxcmax^d=ixrmax^d; ! kxC=[iMmin-3,ixMmax+2]
    kxr^l=kxc^l+kr(idims,^d);              ! kxR=[iMmin-2,ixMmax+3]
 
    lxcmin^d=ixomin^d-kr(idims,^d);lxcmax^d=ixomax^d+kr(idims,^d);! lxC=[iMmin-2,ixMmax+2]
    lxl^l=lxc^l-kr(idims,^d);                          ! lxL=[iMmin-3,ixMmax+1]
    lxr^l=lxc^l+kr(idims,^d);                          ! lxR=[iMmin-1,ixMmax+3]
 
    dwc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)-w(kxc^s,1:nwflux)
    ! Eq. 64,  Miller and Colella 2002, JCP 183, 26 
    d2wc(lxc^s,1:nwflux)=half*(w(lxr^s,1:nwflux)-w(lxl^s,1:nwflux))
    where(dwc(lxc^s,1:nwflux)*dwc(lxl^s,1:nwflux)>zero)
       ! Store the sign of dwC in wMin
       wmin(lxc^s,1:nwflux)= sign(one,d2wc(lxc^s,1:nwflux))
       ! Eq. 65,  Miller and Colella 2002, JCP 183, 26 
       ldw(lxc^s,1:nwflux)= wmin(lxc^s,1:nwflux)*min(dabs(d2wc(lxc^s,1:nwflux)),&
            2.0d0*dabs(dwc(lxl^s,1:nwflux)),&
            2.0d0*dabs(dwc(lxc^s,1:nwflux)))
    else where
       ldw(lxc^s,1:nwflux)=zero
    end where
 
    ! Eq. 66,  Miller and Colella 2002, JCP 183, 26 
    wlc(ixol^s,1:nwflux)=wlc(ixol^s,1:nwflux)+half*dwc(ixol^s,1:nwflux)&
         +(ldw(ixol^s,1:nwflux)-ldw(ixor^s,1:nwflux))/6.0d0
    wrc(ixl^s,1:nwflux)=wlc(ixl^s,1:nwflux)
 
    ! make sure that min w(i)<wLC(i)<w(i+1) same for wRC(i)
    call extremaw(ixi^l,ixo^l,w,1,wmax,wmin)
 
    ! Eq. B8, page 217, Mignone et al 2005, ApJS
    wrc(ixl^s,1:nwflux)=max(wmin(ixo^s,1:nwflux)&
         ,min(wmax(ixo^s,1:nwflux),wrc(ixl^s,1:nwflux))) 
    wlc(ixo^s,1:nwflux)=max(wmin(ixo^s,1:nwflux)&
         ,min(wmax(ixo^s,1:nwflux),wlc(ixo^s,1:nwflux)))
 
    ! Eq. B9, page 217, Mignone et al 2005, ApJS
    where((wrc(ixl^s,1:nwflux)-wct(ixo^s,1:nwflux))&
         *(wct(ixo^s,1:nwflux)-wlc(ixo^s,1:nwflux))<=zero)
       wrc(ixl^s,1:nwflux)=wct(ixo^s,1:nwflux)
       wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
    end where
 
    wmax(ixo^s,1:nwflux)=(wlc(ixo^s,1:nwflux)-wrc(ixl^s,1:nwflux))*&
         (wct(ixo^s,1:nwflux)-half*(wlc(ixo^s,1:nwflux)+wrc(ixl^s,1:nwflux)))
    wmin(ixo^s,1:nwflux)=(wlc(ixo^s,1:nwflux)-wrc(ixl^s,1:nwflux))**2/6.0d0
    tmp(hxl^s,1:nwflux)=wrc(hxl^s,1:nwflux)
    ! Eq. B10, page 218, Mignone et al 2005, ApJS
    where(wmax(hxr^s,1:nwflux)>wmin(hxr^s,1:nwflux))
       wrc(hxc^s,1:nwflux)= 3.0d0*wct(hxr^s,1:nwflux)&
            -2.0d0*wlc(hxr^s,1:nwflux)
    end where
    ! Eq. B11, page 218, Mignone et al 2005, ApJS
    where(wmax(hxc^s,1:nwflux)<-wmin(hxc^s,1:nwflux))
       wlc(hxc^s,1:nwflux)= 3.0d0*wct(hxc^s,1:nwflux)&
            -2.0d0*tmp(hxl^s,1:nwflux)
    end where
 
    ! flattening at the contact discontinuities
    if(flatcd)then
      ixrr^l=ixr^l+kr(idims,^d);               !ixRR=[iMmin+1,ixMmax+3]
      kxl^l=kxc^l-kr(idims,^d);                ! kxL=[iMmin-4,ixMmax+1]
      call ppm_flatcd(ixi^l,kxc^l,kxl^l,kxr^l,wct,d2wc,aa,ab)
      if(any(kappa*aa(kxc^s)>=ab(kxc^s)))then
        do iw=1,nwflux
          where(kappa*aa(kxc^s)>=ab(kxc^s).and. dabs(dwc(kxc^s,iw))>smalldouble)
            wmax(kxc^s,iw) = wct(kxr^s,iw)-2.0d0*wct(kxc^s,iw)+wct(kxl^s,iw)
          end where
 
          where(wmax(ixr^s,iw)*wmax(ixl^s,iw)<zero .and.&
               dabs(wct(ixr^s,iw)-wct(ixl^s,iw))&
               -eps*min(dabs(wct(ixr^s,iw)),dabs(wct(ixl^s,iw)))>zero &
               .and. kappa*aa(ixo^s)>=ab(ixo^s)&
               .and. dabs(dwc(ixo^s,iw))>smalldouble)
 
            ac(ixo^s)=(wct(ixll^s,iw)-wct(ixrr^s,iw)+4.0d0*dwc(ixo^s,iw))&
                 /(12.0d0*dwc(ixo^s,iw))
            wmin(ixo^s,iw)=max(zero,min(eta1*(ac(ixo^s)-eta2),one))
          else where
            wmin(ixo^s,iw)=zero
          end where
 
          where(wmin(hxc^s,iw)>zero)
            wlc(hxc^s,iw) = wlc(hxc^s,iw)*(one-wmin(hxc^s,iw))&
                 +(wct(hxc^s,iw)+half*ldw(hxc^s,iw))*wmin(hxc^s,iw)
          end where
          where(wmin(hxr^s,iw)>zero)
            wrc(hxc^s,iw) = wrc(hxc^s,iw)*(one-wmin(hxr^s,iw))&
                 +(wct(hxr^s,iw)-half*ldw(hxr^s,iw))*wmin(hxr^s,iw)
          end where
        end do
      end if
    end if
 
    ! flattening at the shocks
    if(flatsh)then
       ! following MILLER and COLELLA 2002 JCP 183, 26
       kxcmin^d=ixmin^d-2; kxcmax^d=ixmax^d+2; ! kxC=[ixMmin-2,ixMmax+2]
       ki=bigdouble
       do idimss=1,ndim
          kxl^l=kxc^l-kr(idimss,^d); ! kxL=[ixMmin-3,ixMmax+1]
          kxr^l=kxc^l+kr(idimss,^d); ! kxR=[ixMmin-1,ixMmax+3]
          kxll^l=kxl^l-kr(idimss,^d);! kxLL=[ixMmin-4,ixMmax]
          kxrr^l=kxr^l+kr(idimss,^d);! kxRR=[ixMmin,ixMmax+4]
 
          call ppm_flatsh(ixi^l,kxc^l,kxll^l,kxl^l,kxr^l,kxrr^l,idimss,wct,aa,ab)
 
          ! eq. B17, page 218, Mignone et al 2005, ApJS (Xi1 min)
          ac(kxc^s) = max(zero,min(one,(betamax-aa(kxc^s))/(betamax-betamin)))
          ! eq. B18, page 218, Mignone et al 2005, ApJS (Xi1)
          ! recycling aa(ixL^S)
          where(wct(kxr^s,iw_mom(idimss))<wct(kxl^s,iw_mom(idimss)))
             aa(kxc^s) = max(ac(kxc^s), min(one,(zmax-ab(kxc^s))/(zmax-zmin)))
          else where
             aa(kxc^s) = one
          end where
          ixl^l=ixo^l-kr(idimss,^d); ! ixL=[ixMmin-2,ixMmax]
          ixr^l=ixo^l+kr(idimss,^d); ! ixR=[ixMmin,ixMmax+2]
          ki(ixo^s)=min(ki(ixo^s),aa(ixl^s),aa(ixo^s),aa(ixr^s))
       end do
       ! recycling wMax
       do iw=1,nwflux
          where(dabs(ab(ixo^s)-one)>smalldouble)
             wmax(ixo^s,iw) = (one-ab(ixo^s))*wct(ixo^s,iw)
          end where
 
          where(dabs(ab(hxc^s)-one)>smalldouble)
             wlc(hxc^s,iw) = ab(hxc^s)*wlc(hxc^s,iw)+wmax(hxc^s,iw)
          end where
 
          where(dabs(ab(hxr^s)-one)>smalldouble)
             wrc(hxc^s,iw) = ab(hxr^s)*wrc(hxc^s,iw)+wmax(hxr^s,iw)
          end where
       end do
    end if
 
  end subroutine ppmlimiter
 
  subroutine ppm_flatcd(ixI^L,ixO^L,ixL^L,ixR^L,w,d2w,drho,dp)
    use mod_global_parameters
    use mod_physics, only: phys_gamma
 
    integer, intent(in)             :: ixI^L, ixO^L, ixL^L, ixR^L
    double precision, intent(in)    :: w(ixI^S, nw), d2w(ixG^T, 1:nwflux)
    double precision, intent(inout) :: drho(ixG^T), dp(ixG^T)
 
    drho(ixo^s) = phys_gamma*abs(d2w(ixo^s, iw_rho))&
         /min(w(ixl^s, iw_rho), w(ixr^s, iw_rho))
    dp(ixo^s) = abs(d2w(ixo^s, iw_e))/min(w(ixl^s, iw_e), w(ixr^s, iw_e))
 
  end subroutine ppm_flatcd
 
  subroutine ppm_flatsh(ixI^L,ixO^L,ixLL^L,ixL^L,ixR^L,ixRR^L,idims,w,drho,dp)
    use mod_global_parameters
    use mod_physics, only: phys_gamma
 
    integer, intent(in)             :: ixI^L, ixO^L, ixLL^L, ixL^L, ixR^L, ixRR^L
    integer, intent(in)             :: idims
    double precision, intent(in)    :: w(ixI^S, nw)
    double precision, intent(inout) :: drho(ixI^S), dp(ixI^S)
 
    ! eq. B15, page 218, Mignone and Bodo 2005, ApJS (beta1)
    where (abs(w(ixrr^s, iw_e)-w(ixll^s, iw_e))>smalldouble)
       drho(ixo^s) = abs((w(ixr^s, iw_e)-w(ixl^s, iw_e))&
            /(w(ixrr^s, iw_e)-w(ixll^s, iw_e)))
    else where
       drho(ixo^s) = zero
    end where
 
    ! eq. 76, page 48, Miller and Colella 2002, JCoPh, adjusted
    dp(ixo^s) = abs(w(ixr^s, iw_e)-w(ixl^s, iw_e))&
         /(phys_gamma*(w(ixr^s, iw_e)+w(ixl^s, iw_e)))
 
  end subroutine ppm_flatsh
 
  subroutine extremaq(ixI^L,ixO^L,q,nshift,qMax,qMin)
 
    use mod_global_parameters
 
    integer,intent(in)           :: ixI^L,ixO^L
    double precision, intent(in) :: q(ixI^S)
    integer,intent(in)           :: nshift
 
    double precision, intent(out) :: qMax(ixI^S),qMin(ixI^S)
 
    integer           :: ixs^L,ixsR^L,ixsL^L,idims,jdims,kdims,ishift,i,j
 
    do ishift=1,nshift
      idims=1
      ixsr^l=ixo^l+ishift*kr(idims,^d);
      ixsl^l=ixo^l-ishift*kr(idims,^d);
      if (ishift==1) then
        qmax(ixo^s)=max(q(ixo^s),q(ixsr^s),q(ixsl^s))
        qmin(ixo^s)=min(q(ixo^s),q(ixsr^s),q(ixsl^s))
      else
        qmax(ixo^s)=max(qmax(ixo^s),q(ixsr^s),q(ixsl^s))
        qmin(ixo^s)=min(qmin(ixo^s),q(ixsr^s),q(ixsl^s))
      end if
      {^nooned
      idims=1
      jdims=idims+1
      do i=-1,1
        ixs^l=ixo^l+i*ishift*kr(idims,^d);
        ixsr^l=ixs^l+ishift*kr(jdims,^d);
        ixsl^l=ixs^l-ishift*kr(jdims,^d);
        qmax(ixo^s)=max(qmax(ixo^s),q(ixsr^s),q(ixsl^s))
        qmin(ixo^s)=min(qmin(ixo^s),q(ixsr^s),q(ixsl^s))
      end do
      }
      {^ifthreed
      idims=1
      jdims=idims+1
      kdims=jdims+1
      do i=-1,1
        ixs^l=ixo^l+i*ishift*kr(idims,^d);
        do j=-1,1
          ixs^l=ixo^l+j*ishift*kr(jdims,^d);
          ixsr^l=ixs^l+ishift*kr(kdims,^d);
          ixsl^l=ixs^l-ishift*kr(kdims,^d);
          qmax(ixo^s)=max(qmax(ixo^s),q(ixsr^s),q(ixsl^s))
          qmin(ixo^s)=min(qmin(ixo^s),q(ixsr^s),q(ixsl^s))
        end do
      end do
      }
    enddo
 
  end subroutine  extremaq
 
  subroutine extremaw(ixI^L,ixO^L,w,nshift,wMax,wMin)
    use mod_global_parameters
 
    integer,intent(in)            :: ixI^L,ixO^L
    double precision, intent(in)  :: w(ixI^S,1:nw)
    integer,intent(in)            :: nshift
 
    double precision, intent(out) :: wMax(ixI^S,1:nwflux),wMin(ixI^S,1:nwflux)
 
    integer          :: ixs^L,ixsR^L,ixsL^L,idims,jdims,kdims,ishift,i,j
 
    do ishift=1,nshift
      idims=1
      ixsr^l=ixo^l+ishift*kr(idims,^d);
      ixsl^l=ixo^l-ishift*kr(idims,^d);
      if (ishift==1) then
        wmax(ixo^s,1:nwflux)= &
             max(w(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
        wmin(ixo^s,1:nwflux)= &
             min(w(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
      else
        wmax(ixo^s,1:nwflux)= &
             max(wmax(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
        wmin(ixo^s,1:nwflux)= &
             min(wmin(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
      end if
      {^nooned
      idims=1
      jdims=idims+1
      do i=-1,1
        ixs^l=ixo^l+i*ishift*kr(idims,^d);
        ixsr^l=ixs^l+ishift*kr(jdims,^d);
        ixsl^l=ixs^l-ishift*kr(jdims,^d);
        wmax(ixo^s,1:nwflux)= &
             max(wmax(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
        wmin(ixo^s,1:nwflux)= &
             min(wmin(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
      end do
      }
      {^ifthreed
      idims=1
      jdims=idims+1
      kdims=jdims+1
      do i=-1,1
        ixs^l=ixo^l+i*ishift*kr(idims,^d);
        do j=-1,1
          ixs^l=ixo^l+j*ishift*kr(jdims,^d);
          ixsr^l=ixs^l+ishift*kr(kdims,^d);
          ixsl^l=ixs^l-ishift*kr(kdims,^d);
          wmax(ixo^s,1:nwflux)= &
               max(wmax(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
          wmin(ixo^s,1:nwflux)= &
               min(wmin(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
        end do
      end do
      }
    enddo
 
  end subroutine  extremaw
 
end module mod_ppm