add_sub_circ Subroutine

   subroutine add_sub_circ(ixFmin1,ixFmin2,ixFmin3,ixFmax1,ixFmax2,ixFmax3,&
      ixtEmin1,ixtEmin2,ixtEmin3,ixtEmax1,ixtEmax2,ixtEmax3,ixEmin1,ixEmin2,&
      ixEmin3,ixEmax1,ixEmax2,ixEmax3,ixfEmin1,ixfEmin2,ixfEmin3,ixfEmax1,&
      ixfEmax2,ixfEmax3,edge,idims,iside,add,s)
     use mod_global_parameters

     type(state)        :: s
     integer,intent(in) :: idims,iside
     integer            :: ixFmin1(1:ndim),ixFmin2(1:ndim),ixFmin3(1:ndim),&
        ixFmax1(1:ndim),ixFmax2(1:ndim),ixFmax3(1:ndim),ixtEmin1,ixtEmin2,&
        ixtEmin3,ixtEmax1,ixtEmax2,ixtEmax3,ixEmin1(1:3),ixEmin2(1:3),&
        ixEmin3(1:3),ixEmax1(1:3),ixEmax2(1:3),ixEmax3(1:3),ixfEmin1(1:3),&
        ixfEmin2(1:3),ixfEmin3(1:3),ixfEmax1(1:3),ixfEmax2(1:3),ixfEmax3(1:3)
     double precision   :: edge(ixtEmin1:ixtEmax1,ixtEmin2:ixtEmax2,&
        ixtEmin3:ixtEmax3,1:ndim-1)
     logical,intent(in) :: add

     integer            :: idim1,idim2,idir,middle1,middle2,middle3
     integer            :: ixfECmin1,ixfECmin2,ixfECmin3,ixfECmax1,ixfECmax2,&
        ixfECmax3,ixECmin1,ixECmin2,ixECmin3,ixECmax1,ixECmax2,ixECmax3
     double precision   :: fE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
        sdim:3)
     double precision   :: circ(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
        1:ndim)
     integer            :: ixmin1,ixmin2,ixmin3,ixmax1,ixmax2,ixmax3,hxmin1,&
        hxmin2,hxmin3,hxmax1,hxmax2,hxmax3,ixCmin1,ixCmin2,ixCmin3,ixCmax1,&
        ixCmax2,ixCmax3,hxCmin1,hxCmin2,hxCmin3,hxCmax1,hxCmax2,hxCmax3 !Indices for edges

     ! ixF -> Indices for the faces, depends on the field component
     ! ixE -> Total range for the edges
     ! ixfE -> Edges in fE (3D) array
     ! ix,hx,ixC,hxC -> Auxiliary indices
     ! Assign quantities stored ad edges to make it as similar as
     ! possible to the routine updatefaces.
     fE(:,:,:,:)=zero
     do idim1=1,ndim-1
        ! 3D: rotate indices (see routine flux_to_edge)
       idir=mod(idim1+idims-1,3)+1

       ixfECmin1=ixfEmin1(idir);ixfECmin2=ixfEmin2(idir)
       ixfECmin3=ixfEmin3(idir);ixfECmax1=ixfEmax1(idir)
       ixfECmax2=ixfEmax2(idir);ixfECmax3=ixfEmax3(idir);
       ixECmin1=ixEmin1(idir);ixECmin2=ixEmin2(idir);ixECmin3=ixEmin3(idir)
       ixECmax1=ixEmax1(idir);ixECmax2=ixEmax2(idir);ixECmax3=ixEmax3(idir);
       fE(ixfECmin1:ixfECmax1,ixfECmin2:ixfECmax2,ixfECmin3:ixfECmax3,&
          idir)=edge(ixECmin1:ixECmax1,ixECmin2:ixECmax2,ixECmin3:ixECmax3,&
          idim1)
     end do

     ! Calculate part of circulation needed
     circ=zero
     do idim1=1,ndim
        do idim2=1,ndim
           do idir=sdim,3
             if (lvc(idim1,idim2,idir)==0) cycle
             ! Assemble indices
             ixCmin1=ixFmin1(idim1);ixCmin2=ixFmin2(idim1)
             ixCmin3=ixFmin3(idim1);ixCmax1=ixFmax1(idim1)
             ixCmax2=ixFmax2(idim1);ixCmax3=ixFmax3(idim1);
             hxCmin1=ixCmin1-kr(idim2,1);hxCmin2=ixCmin2-kr(idim2,2)
             hxCmin3=ixCmin3-kr(idim2,3);hxCmax1=ixCmax1-kr(idim2,1)
             hxCmax2=ixCmax2-kr(idim2,2);hxCmax3=ixCmax3-kr(idim2,3);
             if(idim1==idims) then
               circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                  idim1)=circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                  idim1)+lvc(idim1,idim2,idir)*(fE(ixCmin1:ixCmax1,&
                  ixCmin2:ixCmax2,ixCmin3:ixCmax3,idir)-fE(hxCmin1:hxCmax1,&
                  hxCmin2:hxCmax2,hxCmin3:hxCmax3,idir))
             else
               select case(iside)
               case(2)
                 circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                    idim1)=circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                    ixCmin3:ixCmax3,idim1)+lvc(idim1,idim2,&
                    idir)*fE(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                    idir)
               case(1)
                 circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                    idim1)=circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                    ixCmin3:ixCmax3,idim1)-lvc(idim1,idim2,&
                    idir)*fE(hxCmin1:hxCmax1,hxCmin2:hxCmax2,hxCmin3:hxCmax3,&
                    idir)
               end select
             end if
           end do
        end do
     end do

     ! Divide circulation by surface and add
     do idim1=1,ndim
        ixCmin1=ixFmin1(idim1);ixCmin2=ixFmin2(idim1);ixCmin3=ixFmin3(idim1)
        ixCmax1=ixFmax1(idim1);ixCmax2=ixFmax2(idim1);ixCmax3=ixFmax3(idim1);
        where(s%surfaceC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
           idim1)>1.0d-9*s%dvolume(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
           ixCmin3:ixCmax3))
          circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)=circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)/s%surfaceC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)
        elsewhere
          circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,idim1)=zero
        end where
        ! Add/subtract to field at face
        if (add) then
          s%ws(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)=s%ws(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)-circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)
        else
          s%ws(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)=s%ws(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)+circ(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
             idim1)
        end if
     end do

   end subroutine add_sub_circ