divvectorS Subroutine

  subroutine divvectorS(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
     ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,divq)
    use mod_global_parameters
    use mod_limiter

    integer, intent(in)                :: ixImin1,ixImin2,ixImin3,ixImax1,&
       ixImax2,ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
    double precision, intent(in)       :: qvec(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3,1:ndir)
    double precision, intent(inout)    :: divq(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3)
    double precision, dimension(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3) :: qL,qR,dqC,ldq,rdq

    double precision :: invdx(1:ndim)
    integer          :: hxOmin1,hxOmin2,hxOmin3,hxOmax1,hxOmax2,hxOmax3,&
       ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,jxCmin1,jxCmin2,jxCmin3,&
       jxCmax1,jxCmax2,jxCmax3,idims,ixmin1,ixmin2,ixmin3,ixmax1,ixmax2,ixmax3,&
       gxCmin1,gxCmin2,gxCmin3,gxCmax1,gxCmax2,gxCmax3,hxCmin1,hxCmin2,hxCmin3,&
       hxCmax1,hxCmax2,hxCmax3

    ixmin1=ixOmin1-2;ixmin2=ixOmin2-2;ixmin3=ixOmin3-2;ixmax1=ixOmax1+2
    ixmax2=ixOmax2+2;ixmax3=ixOmax3+2;

    if (ixImin1>ixmin1.or.ixImax1<ixmax1.or.ixImin2>ixmin2.or.ixImax2<ixmax2.or.&
       ixImin3>ixmin3.or.ixImax3<ixmax3) call &
       mpistop("Error in divvectorS: Non-conforming input limits")

    invdx=1.d0/dxlevel
    divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)=zero
    do idims=1,ndim
      hxOmin1=ixOmin1-kr(idims,1);hxOmin2=ixOmin2-kr(idims,2)
      hxOmin3=ixOmin3-kr(idims,3);hxOmax1=ixOmax1-kr(idims,1)
      hxOmax2=ixOmax2-kr(idims,2);hxOmax3=ixOmax3-kr(idims,3);
      ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3;ixCmax1=ixOmax1
      ixCmax2=ixOmax2;ixCmax3=ixOmax3;
      jxCmin1=ixCmin1+kr(idims,1);jxCmin2=ixCmin2+kr(idims,2)
      jxCmin3=ixCmin3+kr(idims,3);jxCmax1=ixCmax1+kr(idims,1)
      jxCmax2=ixCmax2+kr(idims,2);jxCmax3=ixCmax3+kr(idims,3);
      gxCmin1=ixCmin1-kr(idims,1);gxCmin2=ixCmin2-kr(idims,2)
      gxCmin3=ixCmin3-kr(idims,3);gxCmax1=jxCmax1;gxCmax2=jxCmax2
      gxCmax3=jxCmax3;
      hxCmin1=gxCmin1+kr(idims,1);hxCmin2=gxCmin2+kr(idims,2)
      hxCmin3=gxCmin3+kr(idims,3);hxCmax1=gxCmax1+kr(idims,1)
      hxCmax2=gxCmax2+kr(idims,2);hxCmax3=gxCmax3+kr(idims,3);

      qR(gxCmin1:gxCmax1,gxCmin2:gxCmax2,&
         gxCmin3:gxCmax3) = qvec(hxCmin1:hxCmax1,hxCmin2:hxCmax2,&
         hxCmin3:hxCmax3,idims)
      qL(gxCmin1:gxCmax1,gxCmin2:gxCmax2,&
         gxCmin3:gxCmax3) = qvec(gxCmin1:gxCmax1,gxCmin2:gxCmax2,&
         gxCmin3:gxCmax3,idims)
      dqC(gxCmin1:gxCmax1,gxCmin2:gxCmax2,gxCmin3:gxCmax3)= qR(gxCmin1:gxCmax1,&
         gxCmin2:gxCmax2,gxCmin3:gxCmax3)-qL(gxCmin1:gxCmax1,gxCmin2:gxCmax2,&
         gxCmin3:gxCmax3)
      call dwlimiter2(dqC,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
         gxCmin1,gxCmin2,gxCmin3,gxCmax1,gxCmax2,gxCmax3,idims,&
         type_gradient_limiter(block%level),ldw=ldq,rdw=rdq)
      qL(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3) = qL(ixCmin1:ixCmax1,&
         ixCmin2:ixCmax2,ixCmin3:ixCmax3) + half*ldq(ixCmin1:ixCmax1,&
         ixCmin2:ixCmax2,ixCmin3:ixCmax3)
      qR(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3) = qR(ixCmin1:ixCmax1,&
         ixCmin2:ixCmax2,ixCmin3:ixCmax3) - half*rdq(jxCmin1:jxCmax1,&
         jxCmin2:jxCmax2,jxCmin3:jxCmax3)

      if (slab_uniform) then
        divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)+half*(qR(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)-qL(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
           hxOmin3:hxOmax3))*invdx(idims)
      else
        qR(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
           ixCmin3:ixCmax3)=block%surfaceC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
           ixCmin3:ixCmax3,idims)*qR(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
           ixCmin3:ixCmax3)
        qL(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
           ixCmin3:ixCmax3)=block%surfaceC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
           ixCmin3:ixCmax3,idims)*qL(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
           ixCmin3:ixCmax3)
        divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)+qR(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)-qL(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
           hxOmin3:hxOmax3)
      end if
    end do
    if(.not.slab_uniform) divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
       ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
       ixOmin3:ixOmax3)/block%dvolume(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
       ixOmin3:ixOmax3)

  end subroutine divvectorS