prolong_2nd_stg Subroutine

  subroutine prolong_2nd_stg(sCo,sFi,ixComin1in,ixComin2in,ixComin3in,&
     ixComax1in,ixComax2in,ixComax3in,ixFimin1in,ixFimin2in,ixFimin3in,&
     ixFimax1in,ixFimax2in,ixFimax3in,dxCo1,dxCo2,dxCo3,xComin1,xComin2,&
     xComin3,dxFi1,dxFi2,dxFi3,xFimin1,xFimin2,xFimin3,ghost,fine_min1in,&
     fine_min2in,fine_min3in,fine_max1in,fine_max2in,fine_max3in)
    use mod_global_parameters
    use mod_physics

    logical, intent(in)          :: ghost
    integer, intent(in)          :: ixComin1in,ixComin2in,ixComin3in,&
       ixComax1in,ixComax2in,ixComax3in, ixFimin1in,ixFimin2in,ixFimin3in,&
       ixFimax1in,ixFimax2in,ixFimax3in
    double precision, intent(in) :: dxCo1,dxCo2,dxCo3, xComin1,xComin2,xComin3,&
        dxFi1,dxFi2,dxFi3, xFimin1,xFimin2,xFimin3
    type(state), intent(in)      :: sCo
    type(state), intent(inout)   :: sFi

    logical, optional :: fine_min1in,fine_min2in,fine_min3in,fine_max1in,&
       fine_max2in,fine_max3in
    logical           :: fine_min1,fine_min2,fine_min3,fine_max1,fine_max2,&
       fine_max3

    double precision :: eta1,eta2,eta3, invdxCo1,invdxCo2,invdxCo3
    integer :: ixComin1,ixComin2,ixComin3,ixComax1,ixComax2,ixComax3,ixFimin1,&
       ixFimin2,ixFimin3,ixFimax1,ixFimax2,ixFimax3
    integer :: idim1,idim2,ix2,ix3,idim3,ixFismin1,ixFismin2,ixFismin3,&
       ixFismax1,ixFismax2,ixFismax3,ixGsmin1,ixGsmin2,ixGsmin3,ixGsmax1,&
       ixGsmax2,ixGsmax3,ixCosmin1,ixCosmin2,ixCosmin3,ixCosmax1,ixCosmax2,&
       ixCosmax3,ixFisCmin1,ixFisCmin2,ixFisCmin3,ixFisCmax1,ixFisCmax2,&
       ixFisCmax3
    integer :: ixCosVmin1(1:ndim),ixCosVmin2(1:ndim),ixCosVmin3(1:ndim),&
       ixCosVmax1(1:ndim),ixCosVmax2(1:ndim),ixCosVmax3(1:ndim),&
       ixFisVmin1(1:ndim),ixFisVmin2(1:ndim),ixFisVmin3(1:ndim),&
       ixFisVmax1(1:ndim),ixFisVmax2(1:ndim),ixFisVmax3(1:ndim)
    integer :: hxCosmin1,hxCosmin2,hxCosmin3,hxCosmax1,hxCosmax2,hxCosmax3,&
       jxCosmin1,jxCosmin2,jxCosmin3,jxCosmax1,jxCosmax2,jxCosmax3,ixCosEmin1,&
       ixCosEmin2,ixCosEmin3,ixCosEmax1,ixCosEmax2,ixCosEmax3,ixFisEmin1,&
       ixFisEmin2,ixFisEmin3,ixFisEmax1,ixFisEmax2,ixFisEmax3,hxFisCmin1,&
       hxFisCmin2,hxFisCmin3,hxFisCmax1,hxFisCmax2,hxFisCmax3,jxFisCmin1,&
       jxFisCmin2,jxFisCmin3,jxFisCmax1,jxFisCmax2,jxFisCmax3,ipxFisCmin1,&
       ipxFisCmin2,ipxFisCmin3,ipxFisCmax1,ipxFisCmax2,ipxFisCmax3,ixCosCmin1,&
       ixCosCmin2,ixCosCmin3,ixCosCmax1,ixCosCmax2,ixCosCmax3,imxFisCmin1,&
       imxFisCmin2,imxFisCmin3,imxFisCmax1,imxFisCmax2,imxFisCmax3,jpxFisCmin1,&
       jpxFisCmin2,jpxFisCmin3,jpxFisCmax1,jpxFisCmax2,jpxFisCmax3,jmxFisCmin1,&
       jmxFisCmin2,jmxFisCmin3,jmxFisCmax1,jmxFisCmax2,jmxFisCmax3,hpxFisCmin1,&
       hpxFisCmin2,hpxFisCmin3,hpxFisCmax1,hpxFisCmax2,hpxFisCmax3
    integer :: hxFimin1,hxFimin2,hxFimin3,hxFimax1,hxFimax2,hxFimax3,jxFimin1,&
       jxFimin2,jxFimin3,jxFimax1,jxFimax2,jxFimax3,hijxFimin1,hijxFimin2,&
       hijxFimin3,hijxFimax1,hijxFimax2,hijxFimax3,hjixFimin1,hjixFimin2,&
       hjixFimin3,hjixFimax1,hjixFimax2,hjixFimax3,hjjxFimin1,hjjxFimin2,&
       hjjxFimin3,hjjxFimax1,hjjxFimax2,hjjxFimax3
    integer :: iihxFimin1,iihxFimin2,iihxFimin3,iihxFimax1,iihxFimax2,&
       iihxFimax3,iijxFimin1,iijxFimin2,iijxFimin3,iijxFimax1,iijxFimax2,&
       iijxFimax3,ijhxFimin1,ijhxFimin2,ijhxFimin3,ijhxFimax1,ijhxFimax2,&
       ijhxFimax3,ijjxFimin1,ijjxFimin2,ijjxFimin3,ijjxFimax1,ijjxFimax2,&
       ijjxFimax3,ihixFimin1,ihixFimin2,ihixFimin3,ihixFimax1,ihixFimax2,&
       ihixFimax3,ijixFimin1,ijixFimin2,ijixFimin3,ijixFimax1,ijixFimax2,&
       ijixFimax3,ihjxFimin1,ihjxFimin2,ihjxFimin3,ihjxFimax1,ihjxFimax2,&
       ihjxFimax3
    integer :: jihxFimin1,jihxFimin2,jihxFimin3,jihxFimax1,jihxFimax2,&
       jihxFimax3,jijxFimin1,jijxFimin2,jijxFimin3,jijxFimax1,jijxFimax2,&
       jijxFimax3,jjhxFimin1,jjhxFimin2,jjhxFimin3,jjhxFimax1,jjhxFimax2,&
       jjhxFimax3,jjjxFimin1,jjjxFimin2,jjjxFimin3,jjjxFimax1,jjjxFimax2,&
       jjjxFimax3,jhixFimin1,jhixFimin2,jhixFimin3,jhixFimax1,jhixFimax2,&
       jhixFimax3,jjixFimin1,jjixFimin2,jjixFimin3,jjixFimax1,jjixFimax2,&
       jjixFimax3,jhjxFimin1,jhjxFimin2,jhjxFimin3,jhjxFimax1,jhjxFimax2,&
       jhjxFimax3
    double precision :: bfluxCo(sCo%ixGsmin1:sCo%ixGsmax1,&
       sCo%ixGsmin2:sCo%ixGsmax2,sCo%ixGsmin3:sCo%ixGsmax3,nws),&
       bfluxFi(sFi%ixGsmin1:sFi%ixGsmax1,sFi%ixGsmin2:sFi%ixGsmax2,&
       sFi%ixGsmin3:sFi%ixGsmax3,nws)
    double precision :: slopes(sCo%ixGsmin1:sCo%ixGsmax1,&
       sCo%ixGsmin2:sCo%ixGsmax2,sCo%ixGsmin3:sCo%ixGsmax3,ndim),&
       B_energy_change(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3)

    ! Directional bias, see pdf
 !These arrays should have ranges ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3
 !Since they are still not known, we give ranges ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,
    ! even though it is wasteful
    double precision :: sigmau(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3),&
       sigmad(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3),sigma(ixGlo1:ixGhi1,&
       ixGlo2:ixGhi2,ixGlo3:ixGhi3,1:ndim), alpha(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
       ixGlo3:ixGhi3,1:ndim)
    ! Auxiliary arrays for magnetic fluxes
    double precision :: F1(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3),&
       F2(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3),F3(ixGlo1:ixGhi1,&
       ixGlo2:ixGhi2,ixGlo3:ixGhi3),F4(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
       ixGlo3:ixGhi3)





    ! Note on the indices:
    ! ixCo  Cells where
    !       divergence-preserving prolongation will be applied.
    ! ixFi  Fine cells which correspond to that extent.
    ! ixCoE For 'expanded', faces that need to be used to calculate
    !       the slopes for the prolongation
    ! ixCos
    ! ixFis
    ! ixCosV And
 !ixFisV For 'variable', since the ranges depend on the direction of the faces
    ! ixCosC For 'copy',
 !ixFisC For 'copy', to fill the information from the fine grid needed in the internal faces prolongation step.
    ! At the end, ixFisC is used to copy the assign the magnetic
    ! fields components to the state structure.

    fine_min1=.false.;fine_min2=.false.;fine_min3=.false.;fine_max1=.false.
    fine_max2=.false.;fine_max3=.false.;

    if(present(fine_min1in)) fine_min1=fine_min1in
    if(present(fine_min2in)) fine_min2=fine_min2in
    if(present(fine_min3in)) fine_min3=fine_min3in;
    if(present(fine_max1in)) fine_max1=fine_max1in
    if(present(fine_max2in)) fine_max2=fine_max2in
    if(present(fine_max3in)) fine_max3=fine_max3in;

 !When filling ghost cells, ixFimin1,ixFimin2,ixFimin3,ixFimax1,ixFimax2,ixFimax3 are given.
 !When refining the block, ixComin1,ixComin2,ixComin3,ixComax1,ixComax2,ixComax3 are given.
    if (ghost) then
      ixFimin1=ixFimin1in;ixFimin2=ixFimin2in;ixFimin3=ixFimin3in
      ixFimax1=ixFimax1in;ixFimax2=ixFimax2in;ixFimax3=ixFimax3in;
      ixComin1=int((ixFimin1+nghostcells+1)/2)
      ixComin2=int((ixFimin2+nghostcells+1)/2)
      ixComin3=int((ixFimin3+nghostcells+1)/2);
      ixComax1=int((ixFimax1+nghostcells+1)/2)
      ixComax2=int((ixFimax2+nghostcells+1)/2)
      ixComax3=int((ixFimax3+nghostcells+1)/2);
    else
      ixComin1=ixComin1in;ixComin2=ixComin2in;ixComin3=ixComin3in
      ixComax1=ixComax1in;ixComax2=ixComax2in;ixComax3=ixComax3in;
      ixFimin1=ixMlo1;ixFimin2=ixMlo2;ixFimin3=ixMlo3;ixFimax1=ixMhi1
      ixFimax2=ixMhi2;ixFimax3=ixMhi3;
    end if

    ! Expanded range for staggered variables
    ixCosmin1=ixComin1-1;ixCosmin2=ixComin2-1;ixCosmin3=ixComin3-1;
    ixCosmax1=ixComax1;ixCosmax2=ixComax2;ixCosmax3=ixComax3;

    ixFismin1=ixFimin1-1;ixFismin2=ixFimin2-1;ixFismin3=ixFimin3-1;
    ixFismax1=ixFimax1;ixFismax2=ixFimax2;ixFismax3=ixFimax3;


    ! Assemble general indices
    ixGsmin1=sFi%ixGsmin1;ixGsmin2=sFi%ixGsmin2;ixGsmin3=sFi%ixGsmin3;
    ixGsmax1=sFi%ixGsmax1;ixGsmax2=sFi%ixGsmax2;ixGsmax3=sFi%ixGsmax3;

    do idim1=1,ndim
      ixCosVmin1(idim1)=ixComin1-kr(1,idim1)
      ixCosVmin2(idim1)=ixComin2-kr(2,idim1)
      ixCosVmin3(idim1)=ixComin3-kr(3,idim1);
      ixCosVmax1(idim1)=ixComax1;ixCosVmax2(idim1)=ixComax2
      ixCosVmax3(idim1)=ixComax3;
      ixFisVmin1(idim1)=ixFimin1-kr(1,idim1)
      ixFisVmin2(idim1)=ixFimin2-kr(2,idim1)
      ixFisVmin3(idim1)=ixFimin3-kr(3,idim1);
      ixFisVmax1(idim1)=ixFimax1;ixFisVmax2(idim1)=ixFimax2
      ixFisVmax3(idim1)=ixFimax3;
    end do

    ! Initialize auxiliary arrays at zero
    bfluxCo = zero
    bfluxFi = zero
    slopes  = zero


    invdxCo1=1.d0/dxCo1;invdxCo2=1.d0/dxCo2;invdxCo3=1.d0/dxCo3;

    ! Fill coarse magnetic flux array
    do idim1=1,ndim
      ! Fill information from parts already at the fine level
      ! First set up indices
      if (ghost) then
        ixFisEmin1=max(1-kr(1,idim1),ixFisVmin1(idim1)-2*(1-kr(1,idim1)))
        ixFisEmin2=max(1-kr(2,idim1),ixFisVmin2(idim1)-2*(1-kr(2,idim1)))
        ixFisEmin3=max(1-kr(3,idim1),ixFisVmin3(idim1)-2*(1-kr(3,idim1)));
        ixFisEmax1=min(ixGsmax1,ixFisVmax1(idim1)+2*(1-kr(1,idim1)))
        ixFisEmax2=min(ixGsmax2,ixFisVmax2(idim1)+2*(1-kr(2,idim1)))
        ixFisEmax3=min(ixGsmax3,ixFisVmax3(idim1)+2*(1-kr(3,idim1)));
        ixCosEmin1=int((ixFisEmin1+nghostcells+1)/2)
        ixCosEmin2=int((ixFisEmin2+nghostcells+1)/2)
        ixCosEmin3=int((ixFisEmin3+nghostcells+1)/2)
        ixCosEmax1=int((ixFisEmax1+nghostcells+1)/2)
        ixCosEmax2=int((ixFisEmax2+nghostcells+1)/2)
        ixCosEmax3=int((ixFisEmax3+nghostcells+1)/2);
      else
        ixCosEmin1=ixCosVmin1(idim1)-(1-kr(idim1,1))
        ixCosEmin2=ixCosVmin2(idim1)-(1-kr(idim1,2))
        ixCosEmin3=ixCosVmin3(idim1)-(1-kr(idim1,3))
        ixCosEmax1=ixCosVmax1(idim1)+(1-kr(idim1,1))
        ixCosEmax2=ixCosVmax2(idim1)+(1-kr(idim1,2))
        ixCosEmax3=ixCosVmax3(idim1)+(1-kr(idim1,3));
        ixFisEmin1=ixFisVmin1(idim1)-2*(1-kr(idim1,1))
        ixFisEmin2=ixFisVmin2(idim1)-2*(1-kr(idim1,2))
        ixFisEmin3=ixFisVmin3(idim1)-2*(1-kr(idim1,3))
        ixFisEmax1=ixFisVmax1(idim1)+2*(1-kr(idim1,1))
        ixFisEmax2=ixFisVmax2(idim1)+2*(1-kr(idim1,2))
        ixFisEmax3=ixFisVmax3(idim1)+2*(1-kr(idim1,3));
      end if
      ! Convert fine fields to fluxes
      bfluxFi(ixFisEmin1:ixFisEmax1,ixFisEmin2:ixFisEmax2,&
         ixFisEmin3:ixFisEmax3,idim1)=sFi%ws(ixFisEmin1:ixFisEmax1,&
         ixFisEmin2:ixFisEmax2,ixFisEmin3:ixFisEmax3,&
         idim1)*sFi%surfaceC(ixFisEmin1:ixFisEmax1,ixFisEmin2:ixFisEmax2,&
         ixFisEmin3:ixFisEmax3,idim1)


      idim2=1+mod(idim1,3)
      idim3=1+mod(idim1+1,3)

      bfluxCo(ixCosEmin1:ixCosEmax1,ixCosEmin2:ixCosEmax2,&
         ixCosEmin3:ixCosEmax3,idim1) = zero
      ! Add fine fluxes sharing the same fine face
     do ix2=0,1
    do ix3=0,1
        ixFisCmin1=ixFisEmin1+ix2*kr(idim2,1)
        ixFisCmin2=ixFisEmin2+ix2*kr(idim2,2)
        ixFisCmin3=ixFisEmin3+ix2*kr(idim2,3)
        ixFisCmax1=ixFisEmax1+ix2*kr(idim2,1)
        ixFisCmax2=ixFisEmax2+ix2*kr(idim2,2)
        ixFisCmax3=ixFisEmax3+ix2*kr(idim2,3);

        ixFisCmin1=ixFisCmin1+ix3*kr(idim3,1)
        ixFisCmin2=ixFisCmin2+ix3*kr(idim3,2)
        ixFisCmin3=ixFisCmin3+ix3*kr(idim3,3)
        ixFisCmax1=ixFisCmax1+ix3*kr(idim3,1)
        ixFisCmax2=ixFisCmax2+ix3*kr(idim3,2)
        ixFisCmax3=ixFisCmax3+ix3*kr(idim3,3);

        bfluxCo(ixCosEmin1:ixCosEmax1,ixCosEmin2:ixCosEmax2,&
           ixCosEmin3:ixCosEmax3,idim1)=bfluxCo(ixCosEmin1:ixCosEmax1,&
           ixCosEmin2:ixCosEmax2,ixCosEmin3:ixCosEmax3,&
           idim1)+bfluxFi(ixFisCmin1:ixFisCmax1:2,ixFisCmin2:ixFisCmax2:2,&
           ixFisCmin3:ixFisCmax3:2,idim1)
     end do
    end do
    end do

    ! Omit indices for already refined face, if any

    if (fine_min1) then
      ixCosVmin1(1)=ixCosVmin1(1)+1
      ixFisVmin1(1)=ixFisVmin1(1)+2
    end if
    if (fine_max1) then
      ixCosVmax1(1)=ixCosVmax1(1)-1
      ixFisVmax1(1)=ixFisVmax1(1)-2
    end if


    if (fine_min2) then
      ixCosVmin2(2)=ixCosVmin2(2)+1
      ixFisVmin2(2)=ixFisVmin2(2)+2
    end if
    if (fine_max2) then
      ixCosVmax2(2)=ixCosVmax2(2)-1
      ixFisVmax2(2)=ixFisVmax2(2)-2
    end if


    if (fine_min3) then
      ixCosVmin3(3)=ixCosVmin3(3)+1
      ixFisVmin3(3)=ixFisVmin3(3)+2
    end if
    if (fine_max3) then
      ixCosVmax3(3)=ixCosVmax3(3)-1
      ixFisVmax3(3)=ixFisVmax3(3)-2
    end if


    do idim1=1,ndim
      ixCosEmin1=ixCosVmin1(idim1);ixCosEmin2=ixCosVmin2(idim1)
      ixCosEmin3=ixCosVmin3(idim1);ixCosEmax1=ixCosVmax1(idim1)
      ixCosEmax2=ixCosVmax2(idim1);ixCosEmax3=ixCosVmax3(idim1);
      ! Omit part already refined

      if (1/=idim1) then
        if ((.not.fine_min1).or.(.not.ghost)) then
         ixCosEmin1=ixCosVmin1(idim1)-1
        end if
        if ((.not.fine_max1).or.(.not.ghost)) then
         ixCosEmax1=ixCosVmax1(idim1)+1
        end if
      end if


      if (2/=idim1) then
        if ((.not.fine_min2).or.(.not.ghost)) then
         ixCosEmin2=ixCosVmin2(idim1)-1
        end if
        if ((.not.fine_max2).or.(.not.ghost)) then
         ixCosEmax2=ixCosVmax2(idim1)+1
        end if
      end if


      if (3/=idim1) then
        if ((.not.fine_min3).or.(.not.ghost)) then
         ixCosEmin3=ixCosVmin3(idim1)-1
        end if
        if ((.not.fine_max3).or.(.not.ghost)) then
         ixCosEmax3=ixCosVmax3(idim1)+1
        end if
      end if

      ! Fill coarse flux array from coarse field
      bfluxCo(ixCosEmin1:ixCosEmax1,ixCosEmin2:ixCosEmax2,&
         ixCosEmin3:ixCosEmax3,idim1)=sCo%ws(ixCosEmin1:ixCosEmax1,&
         ixCosEmin2:ixCosEmax2,ixCosEmin3:ixCosEmax3,&
         idim1)*sCo%surfaceC(ixCosEmin1:ixCosEmax1,ixCosEmin2:ixCosEmax2,&
         ixCosEmin3:ixCosEmax3,idim1)
    end do
    ! Finished filling coarse flux array

    ! Distribute coarse fluxes among fine fluxes
    ! There are too many loops here, perhaps optimize later
    do idim1=1,ndim
       ixCosmin1=ixCosVmin1(idim1);ixCosmin2=ixCosVmin2(idim1)
       ixCosmin3=ixCosVmin3(idim1);ixCosmax1=ixCosVmax1(idim1)
       ixCosmax2=ixCosVmax2(idim1);ixCosmax3=ixCosVmax3(idim1);
       do idim2=1,ndim
         if(idim1==idim2) cycle
         ! Calculate slope in direction idim2
         ! Set up indices
         jxCosmin1=ixCosmin1+kr(idim2,1);jxCosmin2=ixCosmin2+kr(idim2,2)
         jxCosmin3=ixCosmin3+kr(idim2,3);jxCosmax1=ixCosmax1+kr(idim2,1)
         jxCosmax2=ixCosmax2+kr(idim2,2);jxCosmax3=ixCosmax3+kr(idim2,3);
         hxCosmin1=ixCosmin1-kr(idim2,1);hxCosmin2=ixCosmin2-kr(idim2,2)
         hxCosmin3=ixCosmin3-kr(idim2,3);hxCosmax1=ixCosmax1-kr(idim2,1)
         hxCosmax2=ixCosmax2-kr(idim2,2);hxCosmax3=ixCosmax3-kr(idim2,3);
         slopes(ixCosmin1:ixCosmax1,ixCosmin2:ixCosmax2,ixCosmin3:ixCosmax3,&
            idim2)=0.125d0*(bfluxCo(jxCosmin1:jxCosmax1,jxCosmin2:jxCosmax2,&
            jxCosmin3:jxCosmax3,idim1)-bfluxCo(hxCosmin1:hxCosmax1,&
            hxCosmin2:hxCosmax2,hxCosmin3:hxCosmax3,idim1))


         do idim3=1,ndim
           ! Calculate slope in direction idim3
           ! Set up indices
           jxCosmin1=ixCosmin1+kr(idim3,1);jxCosmin2=ixCosmin2+kr(idim3,2)
           jxCosmin3=ixCosmin3+kr(idim3,3);jxCosmax1=ixCosmax1+kr(idim3,1)
           jxCosmax2=ixCosmax2+kr(idim3,2);jxCosmax3=ixCosmax3+kr(idim3,3);
           hxCosmin1=ixCosmin1-kr(idim3,1);hxCosmin2=ixCosmin2-kr(idim3,2)
           hxCosmin3=ixCosmin3-kr(idim3,3);hxCosmax1=ixCosmax1-kr(idim3,1)
           hxCosmax2=ixCosmax2-kr(idim3,2);hxCosmax3=ixCosmax3-kr(idim3,3);
           slopes(ixCosmin1:ixCosmax1,ixCosmin2:ixCosmax2,ixCosmin3:ixCosmax3,&
              idim3)=0.125d0*(bfluxCo(jxCosmin1:jxCosmax1,jxCosmin2:jxCosmax2,&
              jxCosmin3:jxCosmax3,idim1)-bfluxCo(hxCosmin1:hxCosmax1,&
              hxCosmin2:hxCosmax2,hxCosmin3:hxCosmax3,idim1))
           if(lvc(idim1,idim2,idim3)<1) cycle
           do ix2=0,1
             do ix3=0,1
              ixFisCmin1=ixFisVmin1(idim1)+ix2*kr(1,idim2)+ix3*kr(1,idim3)
              ixFisCmin2=ixFisVmin2(idim1)+ix2*kr(2,idim2)+ix3*kr(2,idim3)
              ixFisCmin3=ixFisVmin3(idim1)+ix2*kr(3,idim2)+ix3*kr(3,idim3);
              ixFisCmax1=ixFisVmax1(idim1)+ix2*kr(1,idim2)+ix3*kr(1,idim3)
              ixFisCmax2=ixFisVmax2(idim1)+ix2*kr(2,idim2)+ix3*kr(2,idim3)
              ixFisCmax3=ixFisVmax3(idim1)+ix2*kr(3,idim2)+ix3*kr(3,idim3);
               bfluxFi(ixFisCmin1:ixFisCmax1:2,ixFisCmin2:ixFisCmax2:2,&
                  ixFisCmin3:ixFisCmax3:2,&
                  idim1)=quarter*bfluxCo(ixCosmin1:ixCosmax1,&
                  ixCosmin2:ixCosmax2,ixCosmin3:ixCosmax3,&
                  idim1)+quarter*(2*ix2-1)*slopes(ixCosmin1:ixCosmax1,&
                  ixCosmin2:ixCosmax2,ixCosmin3:ixCosmax3,&
                  idim2)+quarter*(2*ix3-1)*slopes(ixCosmin1:ixCosmax1,&
                  ixCosmin2:ixCosmax2,ixCosmin3:ixCosmax3,idim3)
             end do
           end do
         end do

       end do
    end do

    ! Calculate interior fine fluxes

    ! Directional bias for nonlinear prolongation
    do idim1=1,ndim
      do idim2=1,ndim
        do idim3=1,ndim
          if (lvc(idim1,idim2,idim3)<1) cycle
          ! Set up indices
          hxFimin1=ixFimin1-kr(idim1,1);hxFimin2=ixFimin2-kr(idim1,2)
          hxFimin3=ixFimin3-kr(idim1,3);hxFimax1=ixFimax1-kr(idim1,1)
          hxFimax2=ixFimax2-kr(idim1,2);hxFimax3=ixFimax3-kr(idim1,3);
          jxFimin1=ixFimin1+kr(idim1,1);jxFimin2=ixFimin2+kr(idim1,2)
          jxFimin3=ixFimin3+kr(idim1,3);jxFimax1=ixFimax1+kr(idim1,1)
          jxFimax2=ixFimax2+kr(idim1,2);jxFimax3=ixFimax3+kr(idim1,3);

          hijxFimin1=hxFimin1+kr(idim3,1);hijxFimin2=hxFimin2+kr(idim3,2)
          hijxFimin3=hxFimin3+kr(idim3,3);hijxFimax1=hxFimax1+kr(idim3,1)
          hijxFimax2=hxFimax2+kr(idim3,2);hijxFimax3=hxFimax3+kr(idim3,3);
          hjixFimin1=hxFimin1+kr(idim2,1);hjixFimin2=hxFimin2+kr(idim2,2)
          hjixFimin3=hxFimin3+kr(idim2,3);hjixFimax1=hxFimax1+kr(idim2,1)
          hjixFimax2=hxFimax2+kr(idim2,2);hjixFimax3=hxFimax3+kr(idim2,3);
          hjjxFimin1=hijxFimin1+kr(idim2,1);hjjxFimin2=hijxFimin2+kr(idim2,2)
          hjjxFimin3=hijxFimin3+kr(idim2,3);hjjxFimax1=hijxFimax1+kr(idim2,1)
          hjjxFimax2=hijxFimax2+kr(idim2,2);hjjxFimax3=hijxFimax3+kr(idim2,3);

          iihxFimin1=ixFimin1-kr(idim3,1);iihxFimin2=ixFimin2-kr(idim3,2)
          iihxFimin3=ixFimin3-kr(idim3,3);iihxFimax1=ixFimax1-kr(idim3,1)
          iihxFimax2=ixFimax2-kr(idim3,2);iihxFimax3=ixFimax3-kr(idim3,3);
          iijxFimin1=ixFimin1+kr(idim3,1);iijxFimin2=ixFimin2+kr(idim3,2)
          iijxFimin3=ixFimin3+kr(idim3,3);iijxFimax1=ixFimax1+kr(idim3,1)
          iijxFimax2=ixFimax2+kr(idim3,2);iijxFimax3=ixFimax3+kr(idim3,3);
          ihixFimin1=ixFimin1-kr(idim2,1);ihixFimin2=ixFimin2-kr(idim2,2)
          ihixFimin3=ixFimin3-kr(idim2,3);ihixFimax1=ixFimax1-kr(idim2,1)
          ihixFimax2=ixFimax2-kr(idim2,2);ihixFimax3=ixFimax3-kr(idim2,3);
          ihjxFimin1=ihixFimin1+kr(idim3,1);ihjxFimin2=ihixFimin2+kr(idim3,2)
          ihjxFimin3=ihixFimin3+kr(idim3,3);ihjxFimax1=ihixFimax1+kr(idim3,1)
          ihjxFimax2=ihixFimax2+kr(idim3,2);ihjxFimax3=ihixFimax3+kr(idim3,3);
          ijixFimin1=ixFimin1+kr(idim2,1);ijixFimin2=ixFimin2+kr(idim2,2)
          ijixFimin3=ixFimin3+kr(idim2,3);ijixFimax1=ixFimax1+kr(idim2,1)
          ijixFimax2=ixFimax2+kr(idim2,2);ijixFimax3=ixFimax3+kr(idim2,3);
          ijhxFimin1=ijixFimin1-kr(idim3,1);ijhxFimin2=ijixFimin2-kr(idim3,2)
          ijhxFimin3=ijixFimin3-kr(idim3,3);ijhxFimax1=ijixFimax1-kr(idim3,1)
          ijhxFimax2=ijixFimax2-kr(idim3,2);ijhxFimax3=ijixFimax3-kr(idim3,3);
          ijjxFimin1=ijixFimin1+kr(idim3,1);ijjxFimin2=ijixFimin2+kr(idim3,2)
          ijjxFimin3=ijixFimin3+kr(idim3,3);ijjxFimax1=ijixFimax1+kr(idim3,1)
          ijjxFimax2=ijixFimax2+kr(idim3,2);ijjxFimax3=ijixFimax3+kr(idim3,3);

          jihxFimin1=jxFimin1-kr(idim3,1);jihxFimin2=jxFimin2-kr(idim3,2)
          jihxFimin3=jxFimin3-kr(idim3,3);jihxFimax1=jxFimax1-kr(idim3,1)
          jihxFimax2=jxFimax2-kr(idim3,2);jihxFimax3=jxFimax3-kr(idim3,3);
          jijxFimin1=jxFimin1+kr(idim3,1);jijxFimin2=jxFimin2+kr(idim3,2)
          jijxFimin3=jxFimin3+kr(idim3,3);jijxFimax1=jxFimax1+kr(idim3,1)
          jijxFimax2=jxFimax2+kr(idim3,2);jijxFimax3=jxFimax3+kr(idim3,3);
          jhixFimin1=jxFimin1-kr(idim2,1);jhixFimin2=jxFimin2-kr(idim2,2)
          jhixFimin3=jxFimin3-kr(idim2,3);jhixFimax1=jxFimax1-kr(idim2,1)
          jhixFimax2=jxFimax2-kr(idim2,2);jhixFimax3=jxFimax3-kr(idim2,3);
          jhjxFimin1=jhixFimin1+kr(idim3,1);jhjxFimin2=jhixFimin2+kr(idim3,2)
          jhjxFimin3=jhixFimin3+kr(idim3,3);jhjxFimax1=jhixFimax1+kr(idim3,1)
          jhjxFimax2=jhixFimax2+kr(idim3,2);jhjxFimax3=jhixFimax3+kr(idim3,3);
          jjixFimin1=jxFimin1+kr(idim2,1);jjixFimin2=jxFimin2+kr(idim2,2)
          jjixFimin3=jxFimin3+kr(idim2,3);jjixFimax1=jxFimax1+kr(idim2,1)
          jjixFimax2=jxFimax2+kr(idim2,2);jjixFimax3=jxFimax3+kr(idim2,3);
          jjhxFimin1=jjixFimin1-kr(idim3,1);jjhxFimin2=jjixFimin2-kr(idim3,2)
          jjhxFimin3=jjixFimin3-kr(idim3,3);jjhxFimax1=jjixFimax1-kr(idim3,1)
          jjhxFimax2=jjixFimax2-kr(idim3,2);jjhxFimax3=jjixFimax3-kr(idim3,3);
          jjjxFimin1=jjixFimin1+kr(idim3,1);jjjxFimin2=jjixFimin2+kr(idim3,2)
          jjjxFimin3=jjixFimin3+kr(idim3,3);jjjxFimax1=jjixFimax1+kr(idim3,1)
          jjjxFimax2=jjixFimax2+kr(idim3,2);jjjxFimax3=jjixFimax3+kr(idim3,3);

          sigmau(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)=abs(bfluxFi(jhixFimin1:jhixFimax1:2,&
             jhixFimin2:jhixFimax2:2,jhixFimin3:jhixFimax3:2,&
             idim2))+ abs(bfluxFi(jhjxFimin1:jhjxFimax1:2,&
             jhjxFimin2:jhjxFimax2:2,jhjxFimin3:jhjxFimax3:2,&
             idim2))+ abs(bfluxFi(jjixFimin1:jjixFimax1:2,&
             jjixFimin2:jjixFimax2:2,jjixFimin3:jjixFimax3:2,&
             idim2))+ abs(bfluxFi(jjjxFimin1:jjjxFimax1:2,&
             jjjxFimin2:jjjxFimax2:2,jjjxFimin3:jjjxFimax3:2,&
             idim2))+ abs(bfluxFi(jihxFimin1:jihxFimax1:2,&
             jihxFimin2:jihxFimax2:2,jihxFimin3:jihxFimax3:2,&
             idim3))+ abs(bfluxFi(jijxFimin1:jijxFimax1:2,&
             jijxFimin2:jijxFimax2:2,jijxFimin3:jijxFimax3:2,&
             idim3))+ abs(bfluxFi(jjhxFimin1:jjhxFimax1:2,&
             jjhxFimin2:jjhxFimax2:2,jjhxFimin3:jjhxFimax3:2,&
             idim3))+ abs(bfluxFi(jjjxFimin1:jjjxFimax1:2,&
             jjjxFimin2:jjjxFimax2:2,jjjxFimin3:jjjxFimax3:2,idim3))

          sigmad(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)=abs(bfluxFi(ihixFimin1:ihixFimax1:2,&
             ihixFimin2:ihixFimax2:2,ihixFimin3:ihixFimax3:2,&
             idim2))+ abs(bfluxFi(ihjxFimin1:ihjxFimax1:2,&
             ihjxFimin2:ihjxFimax2:2,ihjxFimin3:ihjxFimax3:2,&
             idim2))+ abs(bfluxFi(ijixFimin1:ijixFimax1:2,&
             ijixFimin2:ijixFimax2:2,ijixFimin3:ijixFimax3:2,&
             idim2))+ abs(bfluxFi(ijjxFimin1:ijjxFimax1:2,&
             ijjxFimin2:ijjxFimax2:2,ijjxFimin3:ijjxFimax3:2,&
             idim2))+ abs(bfluxFi(iihxFimin1:iihxFimax1:2,&
             iihxFimin2:iihxFimax2:2,iihxFimin3:iihxFimax3:2,&
             idim3))+ abs(bfluxFi(iijxFimin1:iijxFimax1:2,&
             iijxFimin2:iijxFimax2:2,iijxFimin3:iijxFimax3:2,&
             idim3))+ abs(bfluxFi(ijhxFimin1:ijhxFimax1:2,&
             ijhxFimin2:ijhxFimax2:2,ijhxFimin3:ijhxFimax3:2,&
             idim3))+ abs(bfluxFi(ijjxFimin1:ijjxFimax1:2,&
             ijjxFimin2:ijjxFimax2:2,ijjxFimin3:ijjxFimax3:2,idim3))

          sigma(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
             idim1)=sigmau(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+sigmad(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)
          where(sigma(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
             idim1)/=zero)
            sigma(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
               idim1)=abs(sigmau(ixComin1:ixComax1,ixComin2:ixComax2,&
               ixComin3:ixComax3)-sigmad(ixComin1:ixComax1,ixComin2:ixComax2,&
               ixComin3:ixComax3))/sigma(ixComin1:ixComax1,ixComin2:ixComax2,&
               ixComin3:ixComax3,idim1)
          elsewhere
            sigma(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
               idim1)=zero
          end where

        end do
      end do
    end do

    ! Directional bias for Toth-Roe prolongation
    do idim1=1,ndim
      do idim2=1,ndim
        do idim3=1,ndim
          if (lvc(idim1,idim2,idim3)<1) cycle
    !       Nonlinear
 !alpha(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,idim1)=sigma(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,idim2)-sigma(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,idim3)
    !       Homogeneous
 !alpha(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,idim1)=0.d0
    !       Toth-Roe
            alpha(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
               idim1)=(sCo%dx(ixComin1:ixComax1,ixComin2:ixComax2,&
               ixComin3:ixComax3,idim2)-sCo%dx(ixComin1:ixComax1,&
               ixComin2:ixComax2,ixComin3:ixComax3,&
               idim3))/(sCo%dx(ixComin1:ixComax1,ixComin2:ixComax2,&
               ixComin3:ixComax3,idim2)+sCo%dx(ixComin1:ixComax1,&
               ixComin2:ixComax2,ixComin3:ixComax3,idim3))
        end do
      end do
    end do


    do idim1=1,ndim
      do idim2=1,ndim


        do idim3=1,ndim
          if (lvc(idim1,idim2,idim3)<1) cycle
          ! Set up indices
          hxFimin1=ixFimin1-kr(idim1,1);hxFimin2=ixFimin2-kr(idim1,2)
          hxFimin3=ixFimin3-kr(idim1,3);hxFimax1=ixFimax1-kr(idim1,1)
          hxFimax2=ixFimax2-kr(idim1,2);hxFimax3=ixFimax3-kr(idim1,3);
          jxFimin1=ixFimin1+kr(idim1,1);jxFimin2=ixFimin2+kr(idim1,2)
          jxFimin3=ixFimin3+kr(idim1,3);jxFimax1=ixFimax1+kr(idim1,1)
          jxFimax2=ixFimax2+kr(idim1,2);jxFimax3=ixFimax3+kr(idim1,3);

          hijxFimin1=hxFimin1+kr(idim3,1);hijxFimin2=hxFimin2+kr(idim3,2)
          hijxFimin3=hxFimin3+kr(idim3,3);hijxFimax1=hxFimax1+kr(idim3,1)
          hijxFimax2=hxFimax2+kr(idim3,2);hijxFimax3=hxFimax3+kr(idim3,3);
          hjixFimin1=hxFimin1+kr(idim2,1);hjixFimin2=hxFimin2+kr(idim2,2)
          hjixFimin3=hxFimin3+kr(idim2,3);hjixFimax1=hxFimax1+kr(idim2,1)
          hjixFimax2=hxFimax2+kr(idim2,2);hjixFimax3=hxFimax3+kr(idim2,3);
          hjjxFimin1=hijxFimin1+kr(idim2,1);hjjxFimin2=hijxFimin2+kr(idim2,2)
          hjjxFimin3=hijxFimin3+kr(idim2,3);hjjxFimax1=hijxFimax1+kr(idim2,1)
          hjjxFimax2=hijxFimax2+kr(idim2,2);hjjxFimax3=hijxFimax3+kr(idim2,3);

          iihxFimin1=ixFimin1-kr(idim3,1);iihxFimin2=ixFimin2-kr(idim3,2)
          iihxFimin3=ixFimin3-kr(idim3,3);iihxFimax1=ixFimax1-kr(idim3,1)
          iihxFimax2=ixFimax2-kr(idim3,2);iihxFimax3=ixFimax3-kr(idim3,3);
          iijxFimin1=ixFimin1+kr(idim3,1);iijxFimin2=ixFimin2+kr(idim3,2)
          iijxFimin3=ixFimin3+kr(idim3,3);iijxFimax1=ixFimax1+kr(idim3,1)
          iijxFimax2=ixFimax2+kr(idim3,2);iijxFimax3=ixFimax3+kr(idim3,3);
          ihixFimin1=ixFimin1-kr(idim2,1);ihixFimin2=ixFimin2-kr(idim2,2)
          ihixFimin3=ixFimin3-kr(idim2,3);ihixFimax1=ixFimax1-kr(idim2,1)
          ihixFimax2=ixFimax2-kr(idim2,2);ihixFimax3=ixFimax3-kr(idim2,3);
          ihjxFimin1=ihixFimin1+kr(idim3,1);ihjxFimin2=ihixFimin2+kr(idim3,2)
          ihjxFimin3=ihixFimin3+kr(idim3,3);ihjxFimax1=ihixFimax1+kr(idim3,1)
          ihjxFimax2=ihixFimax2+kr(idim3,2);ihjxFimax3=ihixFimax3+kr(idim3,3);
          ijixFimin1=ixFimin1+kr(idim2,1);ijixFimin2=ixFimin2+kr(idim2,2)
          ijixFimin3=ixFimin3+kr(idim2,3);ijixFimax1=ixFimax1+kr(idim2,1)
          ijixFimax2=ixFimax2+kr(idim2,2);ijixFimax3=ixFimax3+kr(idim2,3);
          ijhxFimin1=ijixFimin1-kr(idim3,1);ijhxFimin2=ijixFimin2-kr(idim3,2)
          ijhxFimin3=ijixFimin3-kr(idim3,3);ijhxFimax1=ijixFimax1-kr(idim3,1)
          ijhxFimax2=ijixFimax2-kr(idim3,2);ijhxFimax3=ijixFimax3-kr(idim3,3);
          ijjxFimin1=ijixFimin1+kr(idim3,1);ijjxFimin2=ijixFimin2+kr(idim3,2)
          ijjxFimin3=ijixFimin3+kr(idim3,3);ijjxFimax1=ijixFimax1+kr(idim3,1)
          ijjxFimax2=ijixFimax2+kr(idim3,2);ijjxFimax3=ijixFimax3+kr(idim3,3);

          jihxFimin1=jxFimin1-kr(idim3,1);jihxFimin2=jxFimin2-kr(idim3,2)
          jihxFimin3=jxFimin3-kr(idim3,3);jihxFimax1=jxFimax1-kr(idim3,1)
          jihxFimax2=jxFimax2-kr(idim3,2);jihxFimax3=jxFimax3-kr(idim3,3);
          jijxFimin1=jxFimin1+kr(idim3,1);jijxFimin2=jxFimin2+kr(idim3,2)
          jijxFimin3=jxFimin3+kr(idim3,3);jijxFimax1=jxFimax1+kr(idim3,1)
          jijxFimax2=jxFimax2+kr(idim3,2);jijxFimax3=jxFimax3+kr(idim3,3);
          jhixFimin1=jxFimin1-kr(idim2,1);jhixFimin2=jxFimin2-kr(idim2,2)
          jhixFimin3=jxFimin3-kr(idim2,3);jhixFimax1=jxFimax1-kr(idim2,1)
          jhixFimax2=jxFimax2-kr(idim2,2);jhixFimax3=jxFimax3-kr(idim2,3);
          jhjxFimin1=jhixFimin1+kr(idim3,1);jhjxFimin2=jhixFimin2+kr(idim3,2)
          jhjxFimin3=jhixFimin3+kr(idim3,3);jhjxFimax1=jhixFimax1+kr(idim3,1)
          jhjxFimax2=jhixFimax2+kr(idim3,2);jhjxFimax3=jhixFimax3+kr(idim3,3);
          jjixFimin1=jxFimin1+kr(idim2,1);jjixFimin2=jxFimin2+kr(idim2,2)
          jjixFimin3=jxFimin3+kr(idim2,3);jjixFimax1=jxFimax1+kr(idim2,1)
          jjixFimax2=jxFimax2+kr(idim2,2);jjixFimax3=jxFimax3+kr(idim2,3);
          jjhxFimin1=jjixFimin1-kr(idim3,1);jjhxFimin2=jjixFimin2-kr(idim3,2)
          jjhxFimin3=jjixFimin3-kr(idim3,3);jjhxFimax1=jjixFimax1-kr(idim3,1)
          jjhxFimax2=jjixFimax2-kr(idim3,2);jjhxFimax3=jjixFimax3-kr(idim3,3);
          jjjxFimin1=jjixFimin1+kr(idim3,1);jjjxFimin2=jjixFimin2+kr(idim3,2)
          jjjxFimin3=jjixFimin3+kr(idim3,3);jjjxFimax1=jjixFimax1+kr(idim3,1)
          jjjxFimax2=jjixFimax2+kr(idim3,2);jjjxFimax3=jjixFimax3+kr(idim3,3);

          ! Prolongation formulas
          F1(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)=bfluxFi(ihixFimin1:ihixFimax1:2,&
             ihixFimin2:ihixFimax2:2,ihixFimin3:ihixFimax3:2,&
             idim2)-bfluxFi(jhixFimin1:jhixFimax1:2,jhixFimin2:jhixFimax2:2,&
             jhixFimin3:jhixFimax3:2,idim2)-bfluxFi(ijixFimin1:ijixFimax1:2,&
             ijixFimin2:ijixFimax2:2,ijixFimin3:ijixFimax3:2,&
             idim2)+bfluxFi(jjixFimin1:jjixFimax1:2,jjixFimin2:jjixFimax2:2,&
             jjixFimin3:jjixFimax3:2,idim2)

          F2(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)=bfluxFi(ihjxFimin1:ihjxFimax1:2,&
             ihjxFimin2:ihjxFimax2:2,ihjxFimin3:ihjxFimax3:2,&
             idim2)-bfluxFi(jhjxFimin1:jhjxFimax1:2,jhjxFimin2:jhjxFimax2:2,&
             jhjxFimin3:jhjxFimax3:2,idim2)-bfluxFi(ijjxFimin1:ijjxFimax1:2,&
             ijjxFimin2:ijjxFimax2:2,ijjxFimin3:ijjxFimax3:2,&
             idim2)+bfluxFi(jjjxFimin1:jjjxFimax1:2,jjjxFimin2:jjjxFimax2:2,&
             jjjxFimin3:jjjxFimax3:2,idim2)

          F3(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)=bfluxFi(iihxFimin1:iihxFimax1:2,&
             iihxFimin2:iihxFimax2:2,iihxFimin3:iihxFimax3:2,&
             idim3)-bfluxFi(jihxFimin1:jihxFimax1:2,jihxFimin2:jihxFimax2:2,&
             jihxFimin3:jihxFimax3:2,idim3)-bfluxFi(iijxFimin1:iijxFimax1:2,&
             iijxFimin2:iijxFimax2:2,iijxFimin3:iijxFimax3:2,&
             idim3)+bfluxFi(jijxFimin1:jijxFimax1:2,jijxFimin2:jijxFimax2:2,&
             jijxFimin3:jijxFimax3:2,idim3)

          F4(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)=bfluxFi(ijhxFimin1:ijhxFimax1:2,&
             ijhxFimin2:ijhxFimax2:2,ijhxFimin3:ijhxFimax3:2,&
             idim3)-bfluxFi(jjhxFimin1:jjhxFimax1:2,jjhxFimin2:jjhxFimax2:2,&
             jjhxFimin3:jjhxFimax3:2,idim3)-bfluxFi(ijjxFimin1:ijjxFimax1:2,&
             ijjxFimin2:ijjxFimax2:2,ijjxFimin3:ijjxFimax3:2,&
             idim3)+bfluxFi(jjjxFimin1:jjjxFimax1:2,jjjxFimin2:jjjxFimax2:2,&
             jjjxFimin3:jjjxFimax3:2,idim3)

          bfluxFi(ixFimin1:ixFimax1:2,ixFimin2:ixFimax2:2,ixFimin3:ixFimax3:2,&
             idim1)=half*(bfluxFi(hxFimin1:hxFimax1:2,hxFimin2:hxFimax2:2,&
             hxFimin3:hxFimax3:2,idim1)+bfluxFi(jxFimin1:jxFimax1:2,&
             jxFimin2:jxFimax2:2,jxFimin3:jxFimax3:2,&
             idim1))+6.25d-2*((3.d0+alpha(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3,idim2))*F1(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(1.d0-alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim2))*F2(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3)+&
             (3.d0-alpha(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
             idim3))*F3(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(1.d0+alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim3))*F4(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3))

          bfluxFi(ijixFimin1:ijixFimax1:2,ijixFimin2:ijixFimax2:2,&
             ijixFimin3:ijixFimax3:2,idim1)=half*(bfluxFi(&
             hjixFimin1:hjixFimax1:2,hjixFimin2:hjixFimax2:2,&
             hjixFimin3:hjixFimax3:2,idim1)+bfluxFi(jjixFimin1:jjixFimax1:2,&
             jjixFimin2:jjixFimax2:2,jjixFimin3:jjixFimax3:2,&
             idim1))+6.25d-2*((3.d0+alpha(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3,idim2))*F1(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(1.d0-alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim2))*F2(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3)+(1.d0+&
             alpha(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
             idim3))*F3(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(3.d0-alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim3))*F4(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3))

          bfluxFi(iijxFimin1:iijxFimax1:2,iijxFimin2:iijxFimax2:2,&
             iijxFimin3:iijxFimax3:2,idim1)=half*(bfluxFi(&
             hijxFimin1:hijxFimax1:2,hijxFimin2:hijxFimax2:2,&
             hijxFimin3:hijxFimax3:2,idim1)+bfluxFi(jijxFimin1:jijxFimax1:2,&
             jijxFimin2:jijxFimax2:2,jijxFimin3:jijxFimax3:2,&
             idim1))+6.25d-2*((1.d0-alpha(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3,idim2))*F1(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(3.d0+alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim2))*F2(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3)+&
             (3.d0-alpha(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
             idim3))*F3(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(1.d0+alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim3))*F4(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3))

          bfluxFi(ijjxFimin1:ijjxFimax1:2,ijjxFimin2:ijjxFimax2:2,&
             ijjxFimin3:ijjxFimax3:2,idim1)=half*(bfluxFi(&
             hjjxFimin1:hjjxFimax1:2,hjjxFimin2:hjjxFimax2:2,&
             hjjxFimin3:hjjxFimax3:2,idim1)+bfluxFi(jjjxFimin1:jjjxFimax1:2,&
             jjjxFimin2:jjjxFimax2:2,jjjxFimin3:jjjxFimax3:2,&
             idim1))+6.25d-2*((1.d0-alpha(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3,idim2))*F1(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(3.d0+alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim2))*F2(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3)+(1.d0+&
             alpha(ixComin1:ixComax1,ixComin2:ixComax2,ixComin3:ixComax3,&
             idim3))*F3(ixComin1:ixComax1,ixComin2:ixComax2,&
             ixComin3:ixComax3)+(3.d0-alpha(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3,idim3))*F4(ixComin1:ixComax1,&
             ixComin2:ixComax2,ixComin3:ixComax3))
        end do

      end do
    end do

    ! Go back to magnetic fields
    do idim1=1,ndim
      ixFisCmax1=ixFimax1;ixFisCmax2=ixFimax2;ixFisCmax3=ixFimax3;
      ixFisCmin1=ixFimin1-kr(1,idim1);ixFisCmin2=ixFimin2-kr(2,idim1)
      ixFisCmin3=ixFimin3-kr(3,idim1);
      where(sFi%surfaceC(ixFisCmin1:ixFisCmax1,ixFisCmin2:ixFisCmax2,&
         ixFisCmin3:ixFisCmax3,idim1)/=zero)
        sFi%ws(ixFisCmin1:ixFisCmax1,ixFisCmin2:ixFisCmax2,ixFisCmin3:ixFisCmax3,&
           idim1)=bfluxFi(ixFisCmin1:ixFisCmax1,ixFisCmin2:ixFisCmax2,&
           ixFisCmin3:ixFisCmax3,idim1)/sFi%surfaceC(ixFisCmin1:ixFisCmax1,&
           ixFisCmin2:ixFisCmax2,ixFisCmin3:ixFisCmax3,idim1)
      elsewhere
        sFi%ws(ixFisCmin1:ixFisCmax1,ixFisCmin2:ixFisCmax2,ixFisCmin3:ixFisCmax3,&
           idim1)=zero
      end where
    end do

    if(phys_total_energy.and. .not.prolongprimitive) then
      B_energy_change(ixFimin1:ixFimax1,ixFimin2:ixFimax2,&
         ixFimin3:ixFimax3)=0.5d0*sum(sFi%w(ixFimin1:ixFimax1,ixFimin2:ixFimax2,&
         ixFimin3:ixFimax3,iw_mag(:))**2,dim=ndim+1)
    end if
    call phys_face_to_center(ixFimin1,ixFimin2,ixFimin3,ixFimax1,ixFimax2,&
       ixFimax3,sFi)
    if(phys_total_energy.and. .not.prolongprimitive) then
      B_energy_change(ixFimin1:ixFimax1,ixFimin2:ixFimax2,&
         ixFimin3:ixFimax3)=0.5d0*sum(sFi%w(ixFimin1:ixFimax1,ixFimin2:ixFimax2,&
         ixFimin3:ixFimax3,iw_mag(:))**2,&
         dim=ndim+1)-B_energy_change(ixFimin1:ixFimax1,ixFimin2:ixFimax2,&
         ixFimin3:ixFimax3)
      sFi%w(ixFimin1:ixFimax1,ixFimin2:ixFimax2,ixFimin3:ixFimax3,&
         iw_e)=sFi%w(ixFimin1:ixFimax1,ixFimin2:ixFimax2,ixFimin3:ixFimax3,&
         iw_e)+B_energy_change(ixFimin1:ixFimax1,ixFimin2:ixFimax2,&
         ixFimin3:ixFimax3)
    end if

    ! END NOONED
  end subroutine prolong_2nd_stg