mod_fix_conserve.t

Go to the documentation of this file.

!> Module for flux conservation near refinement boundaries
module mod_fix_conserve
  implicit none
  private
 
  type fluxalloc
     double precision, dimension(:^D&,:), pointer:: flux => null()
     double precision, dimension(:^D&,:), pointer:: edge => null()
  end type fluxalloc
  !> store flux to fix conservation
  type(fluxalloc), dimension(:,:,:), allocatable, public :: pflux
 
  double precision, allocatable, save  :: recvbuffer(:), sendbuffer(:)
  ! buffer for corner coarse
  double precision, allocatable, save  :: recvbuffer_cc(:), sendbuffer_cc(:)
  integer, dimension(:), allocatable   :: fc_recvreq, fc_sendreq
  integer, dimension(:,:), allocatable :: fc_recvstat, fc_sendstat
  integer, dimension(^ND), save        :: isize
  integer, save                        :: nrecv, nsend
  integer                              :: ibuf, ibuf_send
  ! ct for corner total
  integer, save                        :: nrecv_ct, nsend_ct
  integer, dimension(:), allocatable   :: cc_recvreq, cc_sendreq
  integer, dimension(:,:), allocatable :: cc_recvstat, cc_sendstat
  integer, dimension(^ND), save        :: isize_stg
  integer                              :: ibuf_cc, ibuf_cc_send
  integer                              :: itag, itag_cc, isend, isend_cc, irecv, irecv_cc
 
  public :: init_comm_fix_conserve
  public :: allocatebflux
  public :: deallocatebflux
  public :: sendflux
  public :: recvflux
  public :: store_flux
  public :: store_edge
  public :: fix_conserve
  public :: fix_edges
 
 contains
 
   subroutine init_comm_fix_conserve(idim^LIM,nwfluxin)
     use mod_global_parameters
 
     integer, intent(in) :: idim^lim,nwfluxin
 
     integer :: iigrid, igrid, idims, iside, i^d, nxco^d
     integer :: ic^d, inc^d, ipe_neighbor
     integer :: recvsize, sendsize
     integer :: recvsize_cc, sendsize_cc
 
     nsend    = 0
     nrecv    = 0
     recvsize = 0
     sendsize = 0
     if(stagger_grid) then
       ! Special communication for diagonal 'coarse corners'
       ! nrecv/send_cc (for 'coarse corners' is a dim=ndim-1 array which
       ! stores the faces that must be communicated in each direction.
       ! nrecv/send_ct (for 'corners total' is the total number of
       ! necessary communications. These special cases have their own
       ! send and receive buffers (send/recvbuffer_cc), their tags, etc.
       nsend_ct=0
       nrecv_ct=0
       recvsize_cc=0
       sendsize_cc=0
     end if
 
     do idims= idim^lim
       select case (idims)
         {case (^d)
         nrecv=nrecv+nrecv_fc(^d)
         nsend=nsend+nsend_fc(^d)
         nxco^d=1^d%nxCo^dd=ixghi^dd/2-nghostcells;
         isize(^d)={nxco^dd*}*(nwfluxin)
         recvsize=recvsize+nrecv_fc(^d)*isize(^d)
         sendsize=sendsize+nsend_fc(^d)*isize(^d)
         if(stagger_grid) then
           ! This does not consider the 'coarse corner' case
           nxco^d=1^d%nxCo^dd=ixghi^dd/2-nghostcells+1;
           isize_stg(^d)={nxco^dd*}*(^nd-1)
           ! the whole size is used (cell centered and staggered)
           isize(^d)=isize(^d)+isize_stg(^d)      
           recvsize=recvsize+nrecv_fc(^d)*isize_stg(^d)
           sendsize=sendsize+nsend_fc(^d)*isize_stg(^d)
           ! Coarse corner case
           nrecv_ct=nrecv_ct+nrecv_cc(^d)
           nsend_ct=nsend_ct+nsend_cc(^d)
           recvsize_cc=recvsize_cc+nrecv_cc(^d)*isize_stg(^d)
           sendsize_cc=sendsize_cc+nsend_cc(^d)*isize_stg(^d)
         end if
         \}
       end select
     end do
 
     ! Reallocate buffers when size differs
     if (allocated(recvbuffer)) then
       if (recvsize /= size(recvbuffer)) then
         deallocate(recvbuffer)
         allocate(recvbuffer(recvsize))
       end if
     else
       allocate(recvbuffer(recvsize))
     end if
 
     if (allocated(fc_recvreq)) then
       if (nrecv /= size(fc_recvreq)) then
         deallocate(fc_recvreq, fc_recvstat)
         allocate(fc_recvstat(mpi_status_size,nrecv), fc_recvreq(nrecv))
       end if
     else
       allocate(fc_recvstat(mpi_status_size,nrecv), fc_recvreq(nrecv))
     end if
 
     if (allocated(fc_sendreq)) then
       if (nsend /= size(fc_sendreq)) then
         deallocate(fc_sendreq, fc_sendstat)
         allocate(fc_sendstat(mpi_status_size,nsend), fc_sendreq(nsend))
       end if
     else
       allocate(fc_sendstat(mpi_status_size,nsend), fc_sendreq(nsend))
     end if
 
     if(stagger_grid) then
 
       if (allocated(sendbuffer)) then
         if (sendsize /= size(sendbuffer)) then
           deallocate(sendbuffer)
           allocate(sendbuffer(sendsize))
         end if
       else
         allocate(sendbuffer(sendsize))
       end if
 
       if (allocated(recvbuffer_cc)) then
         if (recvsize_cc /= size(recvbuffer_cc)) then
           deallocate(recvbuffer_cc)
           allocate(recvbuffer_cc(recvsize_cc))
         end if
       else
         allocate(recvbuffer_cc(recvsize_cc))
       end if
 
       if (allocated(cc_recvreq)) then
         if (nrecv_ct /= size(cc_recvreq)) then
           deallocate(cc_recvreq, cc_recvstat)
           allocate(cc_recvstat(mpi_status_size,nrecv_ct), cc_recvreq(nrecv_ct))
         end if
       else
         allocate(cc_recvstat(mpi_status_size,nrecv_ct), cc_recvreq(nrecv_ct))
       end if
 
       if (allocated(sendbuffer_cc)) then
         if (sendsize_cc /= size(sendbuffer_cc)) then
           deallocate(sendbuffer_cc)
           allocate(sendbuffer_cc(sendsize_cc))
         end if
       else
         allocate(sendbuffer_cc(sendsize_cc))
       end if
 
       if (allocated(cc_sendreq)) then
         if (nsend_ct /= size(cc_sendreq)) then
           deallocate(cc_sendreq, cc_sendstat)
           allocate(cc_sendstat(mpi_status_size,nsend_ct), cc_sendreq(nsend_ct))
         end if
       else
         allocate(cc_sendstat(mpi_status_size,nsend_ct), cc_sendreq(nsend_ct))
       end if
     end if
 
   end subroutine init_comm_fix_conserve
 
   subroutine recvflux(idim^LIM)
     use mod_global_parameters
 
     integer, intent(in) :: idim^lim
 
     integer :: iigrid, igrid, idims, iside, i^d, nxco^d
     integer :: ic^d, inc^d, ipe_neighbor
     integer :: pi^d,mi^d,ph^d,mh^d,idir
 
     if (nrecv>0) then
       fc_recvreq=mpi_request_null
       ibuf=1
       irecv=0
 
       do iigrid=1,igridstail; igrid=igrids(iigrid);
         do idims= idim^lim
           do iside=1,2
             i^d=kr(^d,idims)*(2*iside-3);
 
             if (neighbor_pole(i^d,igrid)/=0) cycle
 
             if (neighbor_type(i^d,igrid)/=4) cycle
             {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
               inc^db=2*i^db+ic^db\}
               ipe_neighbor=neighbor_child(2,inc^d,igrid)
               if (ipe_neighbor/=mype) then
                 irecv=irecv+1
                 itag=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
                 call mpi_irecv(recvbuffer(ibuf),isize(idims), &
                      mpi_double_precision,ipe_neighbor,itag, &
                      icomm,fc_recvreq(irecv),ierrmpi)
                 ibuf=ibuf+isize(idims)
               end if
             {end do\}
           end do
         end do
       end do
     end if
 
     if(stagger_grid) then
     ! receive corners
       if (nrecv_ct>0) then
         cc_recvreq=mpi_request_null
         ibuf_cc=1
         irecv_cc=0
 
         do iigrid=1,igridstail; igrid=igrids(iigrid);
           do idims= idim^lim
             do iside=1,2
               i^d=kr(^d,idims)*(2*iside-3);
               ! Check if there are special corners
               ! (Coarse block diagonal to a fine block)
               ! If there are, receive.
               ! Tags are calculated in the same way as for
               ! normal fluxes, but should not overlap because
               ! inc^D are different
               if (neighbor_type(i^d,igrid)==3) then
                 do idir=idims+1,ndim
                   pi^d=i^d+kr(idir,^d);
                   mi^d=i^d-kr(idir,^d);
                   ph^d=pi^d-kr(idims,^d)*(2*iside-3);
                   mh^d=mi^d-kr(idims,^d)*(2*iside-3);
 
                   if (neighbor_type(pi^d,igrid)==4.and.&
                       neighbor_type(ph^d,igrid)==3.and.&
                       neighbor_pole(pi^d,igrid)==0) then
                      ! Loop on children (several in 3D)
                    {do ic^db=1+int((1-pi^db)/2),2-int((1+pi^db)/2)
                       inc^db=2*pi^db+ic^db\}
                       ipe_neighbor=neighbor_child(2,inc^d,igrid)
                       if (mype/=ipe_neighbor) then
                         irecv_cc=irecv_cc+1
                         itag_cc=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
                         call mpi_irecv(recvbuffer_cc(ibuf_cc),isize_stg(idims),&
                                        mpi_double_precision,ipe_neighbor,itag_cc,&
                                        icomm,cc_recvreq(irecv_cc),ierrmpi)
                         ibuf_cc=ibuf_cc+isize_stg(idims)
                       end if
                    {end do\}
                   end if
 
                   if (neighbor_type(mi^d,igrid)==4.and.&
                       neighbor_type(mh^d,igrid)==3.and.&
                       neighbor_pole(mi^d,igrid)==0) then
                      ! Loop on children (several in 3D)
                    {do ic^db=1+int((1-mi^db)/2),2-int((1+mi^db)/2)
                        inc^db=2*mi^db+ic^db\}
                       ipe_neighbor=neighbor_child(2,inc^d,igrid)
                       if (mype/=ipe_neighbor) then
                         irecv_cc=irecv_cc+1
                         itag_cc=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
                         call mpi_irecv(recvbuffer_cc(ibuf_cc),isize_stg(idims),&
                                        mpi_double_precision,ipe_neighbor,itag_cc,&
                                        icomm,cc_recvreq(irecv_cc),ierrmpi)
                         ibuf_cc=ibuf_cc+isize_stg(idims)
                       end if
                    {end do\}
                   end if
                 end do
               end if
             end do
           end do
         end do
       end if
     end if ! end if stagger grid
 
   end subroutine recvflux
 
   subroutine sendflux(idim^LIM)
     use mod_global_parameters
 
     integer, intent(in) :: idim^lim
 
     integer :: idims, iside, i^d, ic^d, inc^d, ix^d, ixco^d, nxco^d, iw
     integer :: ineighbor, ipe_neighbor, igrid, iigrid, ibuf_send_next
     integer :: idir, ibuf_cc_send_next, pi^d, ph^d, mi^d, mh^d
 
     fc_sendreq = mpi_request_null
     isend      = 0
     if(stagger_grid) then
       ibuf_send  = 1
       cc_sendreq=mpi_request_null
       isend_cc=0
       ibuf_cc_send=1
     end if
 
     do iigrid=1,igridstail; igrid=igrids(iigrid);
       do idims = idim^lim
         select case (idims)
        {case (^d)
           do iside=1,2
             i^dd=kr(^dd,^d)*(2*iside-3);
 
             if (neighbor_pole(i^dd,igrid)/=0) cycle
 
             if (neighbor_type(i^dd,igrid)==neighbor_coarse) then
               ! send flux to coarser neighbor
               ineighbor=neighbor(1,i^dd,igrid)
               ipe_neighbor=neighbor(2,i^dd,igrid)
               if (ipe_neighbor/=mype) then
                 ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
                 inc^d=-2*i^d+ic^d^d%inc^dd=ic^dd;
                 itag=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
                 isend=isend+1
 
                 if(stagger_grid) then
                   ibuf_send_next=ibuf_send+isize(^d)
                   sendbuffer(ibuf_send:ibuf_send_next-isize_stg(^d)-1)=&
                   reshape(pflux(iside,^d,igrid)%flux,(/isize(^d)-isize_stg(^d)/))
 
                   sendbuffer(ibuf_send_next-isize_stg(^d):ibuf_send_next-1)=&
                   reshape(pflux(iside,^d,igrid)%edge,(/isize_stg(^d)/))
                   call mpi_isend(sendbuffer(ibuf_send),isize(^d), &
                        mpi_double_precision,ipe_neighbor,itag, &
                        icomm,fc_sendreq(isend),ierrmpi)
                   ibuf_send=ibuf_send_next
                 else
                   call mpi_isend(pflux(iside,^d,igrid)%flux,isize(^d), &
                        mpi_double_precision,ipe_neighbor,itag, &
                        icomm,fc_sendreq(isend),ierrmpi)
                 end if
               end if
 
               if(stagger_grid) then
                 ! If we are in a fine block surrounded by coarse blocks
                 do idir=idims+1,ndim
                   pi^dd=i^dd+kr(idir,^dd);
                   mi^dd=i^dd-kr(idir,^dd);
                   ph^dd=pi^dd-kr(idims,^dd)*(2*iside-3);
                   mh^dd=mi^dd-kr(idims,^dd)*(2*iside-3);
 
                   if (neighbor_type(pi^dd,igrid)==2.and.&
                       neighbor_type(ph^dd,igrid)==2.and.&
                       mype/=neighbor(2,pi^dd,igrid).and.&
                       neighbor_pole(pi^dd,igrid)==0) then
                     ! Get relative position in the grid for tags
                     ineighbor=neighbor(1,pi^dd,igrid)
                     ipe_neighbor=neighbor(2,pi^dd,igrid)
                     ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
                     inc^dd=-2*pi^dd+ic^dd;
                     itag_cc=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
                     ! Reshape to buffer and send
                     isend_cc=isend_cc+1
                     ibuf_cc_send_next=ibuf_cc_send+isize_stg(^d)
                     sendbuffer_cc(ibuf_cc_send:ibuf_cc_send_next-1)=&
                     reshape(pflux(iside,^d,igrid)%edge,shape=(/isize_stg(^d)/))
                     call mpi_isend(sendbuffer_cc(ibuf_cc_send),isize_stg(^d),&
                                    mpi_double_precision,ipe_neighbor,itag_cc,&
                                    icomm,cc_sendreq(isend_cc),ierrmpi)
                     ibuf_cc_send=ibuf_cc_send_next
                   end if
       
                   if (neighbor_type(mi^dd,igrid)==2.and.&
                       neighbor_type(mh^dd,igrid)==2.and.&
                       mype/=neighbor(2,mi^dd,igrid).and.&
                       neighbor_pole(mi^dd,igrid)==0) then
                     ! Get relative position in the grid for tags
                     ineighbor=neighbor(1,mi^dd,igrid)
                     ipe_neighbor=neighbor(2,mi^dd,igrid)
                     ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
                     inc^dd=-2*pi^dd+ic^dd;
                     inc^dd=-2*mi^dd+ic^dd;
                     itag_cc=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
                     ! Reshape to buffer and send
                     isend_cc=isend_cc+1
                     ibuf_cc_send_next=ibuf_cc_send+isize_stg(^d)
                     sendbuffer_cc(ibuf_cc_send:ibuf_cc_send_next-1)=&
                     reshape(pflux(iside,^d,igrid)%edge,shape=(/isize_stg(^d)/))
                     call mpi_isend(sendbuffer_cc(ibuf_cc_send),isize_stg(^d),&
                                    mpi_double_precision,ipe_neighbor,itag_cc,&
                                    icomm,cc_sendreq(isend_cc),ierrmpi)
                     ibuf_cc_send=ibuf_cc_send_next
                   end if
                 end do
               end if ! end if stagger grid
 
             end if
           end do\}
         end select
       end do
     end do
   end subroutine sendflux
 
   subroutine allocatebflux
     use mod_global_parameters
 
     integer :: iigrid, igrid, iside, i^d, nx^d, nxco^d
     integer :: idir,idim,pi^d, mi^d, ph^d, mh^d ! To detect corners
 
     nx^d=ixmhi^d-ixmlo^d+1;
     nxco^d=nx^d/2;
 
     do iigrid=1,igridstail; igrid=igrids(iigrid);
       ! For every grid,
       ! arrays for the fluxes are allocated for every face direction(^D)
       ! and every side (1=left, 2=right)
       {do iside=1,2
         i^dd=kr(^dd,^d)*(2*iside-3);
 
         if (neighbor_pole(i^dd,igrid)/=0) cycle
 
         select case (neighbor_type(i^dd,igrid))
         case(neighbor_fine)
           allocate(pflux(iside,^d,igrid)%flux(1^d%1:nx^dd,1:nwflux))
           if(stagger_grid) allocate(pflux(iside,^d,igrid)%edge(1^d%0:nx^dd,1:ndim-1))
         case(neighbor_coarse)
           allocate(pflux(iside,^d,igrid)%flux(1^d%1:nxco^dd,1:nwflux))
           if(stagger_grid) allocate(pflux(iside,^d,igrid)%edge(1^d%0:nxco^dd,1:ndim-1))
         case(neighbor_sibling)
           if(stagger_grid) then
             idim=^d
             do idir=idim+1,ndim
             !do idir=min(idim+1,ndim),ndim
               pi^dd=i^dd+kr(idir,^dd);
               mi^dd=i^dd-kr(idir,^dd);
               ph^dd=pi^dd-kr(^d,^dd)*(2*iside-3);
               mh^dd=mi^dd-kr(^d,^dd)*(2*iside-3);
               if ((neighbor_type(pi^dd,igrid)==4&
                   .and.neighbor_type(ph^dd,igrid)==3)&
                   .or.(neighbor_type(mi^dd,igrid)==4&
                   .and.neighbor_type(mh^dd,igrid)==3)) then
                 allocate(pflux(iside,^d,igrid)%edge(1^d%0:nx^dd,1:ndim-1))
                 exit
               end if
             end do
           end if
         end select
       end do\}
     end do
 
   end subroutine allocatebflux
 
   subroutine deallocatebflux
     use mod_global_parameters
 
     integer :: iigrid, igrid, iside
 
     do iigrid=1,igridstail; igrid=igrids(iigrid);
       {do iside=1,2
         if (associated(pflux(iside,^d,igrid)%flux)) then
           deallocate(pflux(iside,^d,igrid)%flux)
           nullify(pflux(iside,^d,igrid)%flux)
         end if
         if (associated(pflux(iside,^d,igrid)%edge)) then
           deallocate(pflux(iside,^d,igrid)%edge)
           nullify(pflux(iside,^d,igrid)%edge)
         end if
       end do\}
     end do
 
   end subroutine deallocatebflux
 
   subroutine fix_conserve(psb,idim^LIM,nw0,nwfluxin)
     use mod_global_parameters
 
     integer, intent(in) :: idim^lim, nw0, nwfluxin
     type(state) :: psb(max_blocks)
 
     integer :: iigrid, igrid, idims, iside, iotherside, i^d, ic^d, inc^d, ix^l
     integer :: nxco^d, iw, ix, ipe_neighbor, ineighbor, nbuf, ibufnext, nw1
     double precision :: cofiratio
 
     nw1=nw0-1+nwfluxin
     if (slab_uniform) then
       ! The flux is divided by volume of fine cell. We need, however,
       ! to divide by volume of coarse cell => muliply by volume ratio
       cofiratio=one/dble(2**ndim)
     end if
 
     if (nrecv>0) then
       call mpi_waitall(nrecv,fc_recvreq,fc_recvstat,ierrmpi)
       ibuf=1
     end if
 
     nxco^d=(ixmhi^d-ixmlo^d+1)/2;
 
     ! for all grids: perform flux update at Coarse-Fine interfaces
     do iigrid=1,igridstail; igrid=igrids(iigrid);
       do idims= idim^lim
         select case (idims)
           {case (^d)
           do iside=1,2
             i^dd=kr(^dd,^d)*(2*iside-3);
 
             if (neighbor_pole(i^dd,igrid)/=0) cycle
 
             if (neighbor_type(i^dd,igrid)/=4) cycle
 
             ! opedit: skip over active/passive interface since flux for passive ones is 
             ! not computed, keep the buffer counter up to date:
             if (.not.neighbor_active(i^dd,igrid).or.&
                  .not.neighbor_active(0^dd&,igrid) ) then
               {do ic^ddb=1+int((1-i^ddb)/2),2-int((1+i^ddb)/2)
               inc^ddb=2*i^ddb+ic^ddb\}
               ipe_neighbor=neighbor_child(2,inc^dd,igrid)
               if (ipe_neighbor/=mype) then
                 ibufnext=ibuf+isize(^d)
                 ibuf=ibufnext
               end if
               {end do\}
               cycle
             end if
             !
 
             select case (iside)
             case (1)
               ix=ixmlo^d
             case (2)
               ix=ixmhi^d
             end select
 
             ! remove coarse flux
             if (slab_uniform) then
               psb(igrid)%w(ix^d%ixM^t,nw0:nw1) &
                    = psb(igrid)%w(ix^d%ixM^t,nw0:nw1) &
                    -pflux(iside,^d,igrid)%flux(1^d%:^dd&,1:nwfluxin)
             else
               do iw=nw0,nw1
                 psb(igrid)%w(ix^d%ixM^t,iw)=psb(igrid)%w(ix^d%ixM^t,iw)&
                      -pflux(iside,^d,igrid)%flux(1^d%:^dd&,iw-nw0+1) &
                      /ps(igrid)%dvolume(ix^d%ixM^t)
               end do
             end if
 
 
             ! add fine flux
            {do ic^ddb=1+int((1-i^ddb)/2),2-int((1+i^ddb)/2)
               inc^ddb=2*i^ddb+ic^ddb\}
               ineighbor=neighbor_child(1,inc^dd,igrid)
               ipe_neighbor=neighbor_child(2,inc^dd,igrid)
               ixmin^d=ix^d%ixmin^dd=ixmlo^dd+(ic^dd-1)*nxco^dd;
               ixmax^d=ix^d%ixmax^dd=ixmin^dd-1+nxco^dd;
               if (ipe_neighbor==mype) then
                 iotherside=3-iside
                 if (slab_uniform) then
                   psb(igrid)%w(ix^s,nw0:nw1) &
                        = psb(igrid)%w(ix^s,nw0:nw1) &
                        + pflux(iotherside,^d,ineighbor)%flux(:^dd&,1:nwfluxin)&
                        * cofiratio
                 else
                   do iw=nw0,nw1
                     psb(igrid)%w(ix^s,iw)=psb(igrid)%w(ix^s,iw) &
                          +pflux(iotherside,^d,ineighbor)%flux(:^dd&,iw-nw0+1) &
                          /ps(igrid)%dvolume(ix^s)
                   end do
                 end if
               else
                 if (slab_uniform) then
                   ibufnext=ibuf+isize(^d)
                   if(stagger_grid) ibufnext=ibufnext-isize_stg(^d)
                   psb(igrid)%w(ix^s,nw0:nw1) &
                        = psb(igrid)%w(ix^s,nw0:nw1)+cofiratio &
                        *reshape(source=recvbuffer(ibuf:ibufnext-1), &
                        shape=shape(psb(igrid)%w(ix^s,nw0:nw1)))
                   ibuf=ibuf+isize(^d)
                 else
                   ibufnext=ibuf+isize(^d)
                   if(stagger_grid) then
                     nbuf=(isize(^d)-isize_stg(^d))/nwfluxin
                   else
                     nbuf=isize(^d)/nwfluxin
                   end if
                   do iw=nw0,nw1
                     psb(igrid)%w(ix^s,iw)=psb(igrid)%w(ix^s,iw) &
                          +reshape(source=recvbuffer(ibuf:ibufnext-1), &
                          shape=shape(psb(igrid)%w(ix^s,iw))) &
                          /ps(igrid)%dvolume(ix^s)
                     ibuf=ibuf+nbuf
                   end do
                   ibuf=ibufnext
                 end if
               end if
            {end do\}
           end do\}
         end select
       end do
     end do
 
     if (nsend>0) then
       call mpi_waitall(nsend,fc_sendreq,fc_sendstat,ierrmpi)
     end if
 
   end subroutine fix_conserve
 
   subroutine store_flux(igrid,fC,idim^LIM,nwfluxin)
     use mod_global_parameters
 
     integer, intent(in)          :: igrid, idim^lim, nwfluxin
     double precision, intent(in) :: fc(ixg^t,1:nwfluxin,1:ndim)
 
     integer :: idims, iside, i^d, ic^d, inc^d, ix^d, ixco^d, nxco^d, iw
 
     do idims = idim^lim
       select case (idims)
         {case (^d)
         do iside=1,2
           i^dd=kr(^dd,^d)*(2*iside-3);
 
           if (neighbor_pole(i^dd,igrid)/=0) cycle
 
           select case (neighbor_type(i^dd,igrid))
           case (neighbor_fine)
             select case (iside)
             case (1)
               pflux(iside,^d,igrid)%flux(1^d%:^dd&,1:nwfluxin) = &
                    -fc(nghostcells^d%ixM^t,1:nwfluxin,^d)
             case (2)
               pflux(iside,^d,igrid)%flux(1^d%:^dd&,1:nwfluxin) = &
                    fc(ixmhi^d^d%ixM^t,1:nwfluxin,^d)
             end select
           case (neighbor_coarse)
             nxco^d=1^d%nxCo^dd=ixghi^dd/2-nghostcells;
             select case (iside)
             case (1)
               do iw=1,nwfluxin
                {do ixco^ddb=1,nxco^ddb\}
                   ix^d=nghostcells^d%ix^dd=ixmlo^dd+2*(ixco^dd-1);
                   pflux(iside,^d,igrid)%flux(ixco^dd,iw) &
                        = {^noonedsum}(fc(ix^d^d%ix^dd:ix^dd+1,iw,^d))
                {end do\}
               end do
             case (2)
               do iw=1,nwfluxin
                {do ixco^ddb=1,nxco^ddb\}
                   ix^d=ixmhi^d^d%ix^dd=ixmlo^dd+2*(ixco^dd-1);
                   pflux(iside,^d,igrid)%flux(ixco^dd,iw) &
                        =-{^noonedsum}(fc(ix^d^d%ix^dd:ix^dd+1,iw,^d))
                {end do\}
               end do
             end select
           end select
         end do\}
       end select
     end do
 
   end subroutine store_flux
 
   subroutine store_edge(igrid,ixI^L,fE,idim^LIM)
     use mod_global_parameters
     
     integer, intent(in)          :: igrid, ixi^l, idim^lim
     double precision, intent(in) :: fe(ixi^s,sdim:3)
     
     integer :: idims, idir, iside, i^d
     integer :: pi^d, mi^d, ph^d, mh^d ! To detect corners
     integer :: ixmc^l
 
     do idims = idim^lim  !loop over face directions
       !! Loop over block faces
       do iside=1,2 
         i^d=kr(^d,idims)*(2*iside-3);
         if (neighbor_pole(i^d,igrid)/=0) cycle
         select case (neighbor_type(i^d,igrid))
         case (neighbor_fine)
           ! The neighbour is finer
           ! Face direction, side (left or right), restrict ==ired?, fE
           call flux_to_edge(igrid,ixi^l,idims,iside,.false.,fe)
         case(neighbor_coarse)
           ! The neighbour is coarser
           call flux_to_edge(igrid,ixi^l,idims,iside,.true.,fe)
         case(neighbor_sibling)
           ! If the neighbour is at the same level,
           ! check if there are corners
           ! If there is any corner, store the fluxes from that side
           do idir=idims+1,ndim
             pi^d=i^d+kr(idir,^d);
             mi^d=i^d-kr(idir,^d);
             ph^d=pi^d-kr(idims,^d)*(2*iside-3);
             mh^d=mi^d-kr(idims,^d)*(2*iside-3);
             if (neighbor_type(pi^d,igrid)==4&
               .and.neighbor_type(ph^d,igrid)==3) then
               call flux_to_edge(igrid,ixi^l,idims,iside,.false.,fe)
             end if
             if (neighbor_type(mi^d,igrid)==4&
               .and.neighbor_type(mh^d,igrid)==3) then
               call flux_to_edge(igrid,ixi^l,idims,iside,.false.,fe)
             end if
           end do
         end select
       end do
     end do
 
   end subroutine store_edge
 
   subroutine flux_to_edge(igrid,ixI^L,idims,iside,restrict,fE)
     use mod_global_parameters
 
     integer                      :: igrid,ixi^l,idims,iside
     logical                      :: restrict
     double precision, intent(in) :: fe(ixi^s,sdim:3)
 
     integer                      :: idir1,idir2
     integer                      :: ixe^l,ixf^l{^ifthreed, jxf^l,}, nx^d,nxco^d
 
     nx^d=ixmhi^d-ixmlo^d+1;
     nxco^d=nx^d/2;
     ! ixE are the indices on the 'edge' array.
     ! ixF are the indices on the 'fE' array
     ! jxF are indices advanced to perform the flux restriction (sum) in 3D
     ! A line integral of the electric field on the coarse side
     ! lies over two edges on the fine side. So, in 3D we restrict by summing
     ! over two cells on the fine side.
 
     do idir1=1,ndim-1
      {^ifthreed ! 3D: rotate indices among 1 and 2 to save space 
       idir2=mod(idir1+idims-1,3)+1}
      {^iftwod ! Assign the only field component (3) to the only index (1)
       idir2=3}
 
       if (restrict) then
         ! Set up indices for restriction
         ixfmin^d=ixmlo^d-1+kr(^d,idir2);
         ixfmax^d=ixmhi^d-kr(^d,idir2);
         {^ifthreed
         jxf^l=ixf^l+kr(^d,idir2);}
 
         ixemin^d=0+kr(^d,idir2);
         ixemax^d=nxco^d;
         select case(idims)
        {case(^d)
           ixemin^d=1;ixemax^d=1;
           select case(iside)
           case(1)
             ixfmax^d=ixfmin^d
             {^ifthreedjxfmax^d=ixfmin^d}
           case(2)
             ixfmin^d=ixfmax^d
             {^ifthreedjxfmin^d=ixfmax^d}
           end select
        \}
         end select
 
       pflux(iside,idims,igrid)%edge(ixe^s,idir1)=&
         fe(ixfmin^d:ixfmax^d:2,idir2){^ifthreed +&
         fe(jxfmin^d:jxfmax^d:2,idir2)};
 
       else
         ! Set up indices for copying 
         ixfmin^d=ixmlo^d-1+kr(^d,idir2);
         ixfmax^d=ixmhi^d;
         ixemin^d=0+kr(^d,idir2);
         ixemax^d=nx^d;
 
         select case(idims)
        {case(^d)
           ixemin^d=1;ixemax^d=1;
           select case(iside)
           case(1)
             ixfmax^d=ixfmin^d
           case(2)
             ixfmin^d=ixfmax^d
           end select
        \}
         end select
 
         pflux(iside,idims,igrid)%edge(ixe^s,idir1)=fe(ixf^s,idir2)
 
       end if
 
     end do
 
   end subroutine flux_to_edge
 
   subroutine fix_edges(psuse,idim^LIM)
     use mod_global_parameters
 
     type(state) :: psuse(max_blocks)
     integer, intent(in) :: idim^lim
 
     integer :: iigrid, igrid, idims, iside, iotherside, i^d, ic^d, inc^d, ixmc^l
     integer :: nbuf, ibufnext
     integer :: ibufnext_cc
     integer :: pi^d, mi^d, ph^d, mh^d ! To detect corners
     integer :: ixe^l(1:3), ixte^l, ixf^l(1:ndim), ixfe^l(1:3)
     integer :: nx^d, idir, ix, ipe_neighbor, ineighbor
     logical :: pcorner(1:ndim),mcorner(1:ndim)
 
     if (nrecv_ct>0) then
        call mpi_waitall(nrecv_ct,cc_recvreq,cc_recvstat,ierrmpi)
     end if
 
     ! Initialise buffer counter again
     ibuf=1
     ibuf_cc=1
     do iigrid=1,igridstail; igrid=igrids(iigrid);
       do idims= idim^lim
         do iside=1,2
           i^d=kr(^d,idims)*(2*iside-3);
           if (neighbor_pole(i^d,igrid)/=0) cycle
           select case(neighbor_type(i^d,igrid))
           case(neighbor_fine)
             ! The first neighbour is finer
             if (.not.neighbor_active(i^d,igrid).or.&
                 .not.neighbor_active(0^d&,igrid) ) then
               {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
                  inc^db=2*i^db+ic^db\}
                  ipe_neighbor=neighbor_child(2,inc^d,igrid)
                  !! When the neighbour is in a different process
                  if (ipe_neighbor/=mype) then
                     ibufnext=ibuf+isize(idims)
                     ibuf=ibufnext
                     end if
               {end do\}
                cycle
             end if
 
             ! Check if there are corners
             pcorner=.false.
             mcorner=.false.
             do idir=1,ndim
               pi^d=i^d+kr(idir,^d);
               mi^d=i^d-kr(idir,^d);
               ph^d=pi^d-kr(idims,^d)*(2*iside-3);
               mh^d=mi^d-kr(idims,^d)*(2*iside-3);
               if (neighbor_type(ph^d,igrid)==neighbor_fine) pcorner(idir)=.true.
               if (neighbor_type(mh^d,igrid)==neighbor_fine) mcorner(idir)=.true.
             end do
             ! Calculate indices range
             call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.false.,.false.,0^d&,pcorner,mcorner)
             ! Remove coarse part of circulation
             call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iside,idims,igrid)%edge,idims,iside,.false.,psuse(igrid))
             ! Add fine part of the circulation
            {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
               inc^db=2*i^db+ic^db\}
               ineighbor=neighbor_child(1,inc^d,igrid)
               ipe_neighbor=neighbor_child(2,inc^d,igrid)
               iotherside=3-iside
               nx^d=(ixmhi^d-ixmlo^d+1)/2;
               call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.true.,.false.,inc^d,pcorner,mcorner)
               if (ipe_neighbor==mype) then
                 call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iotherside,idims,ineighbor)%edge,idims,iside,.true.,psuse(igrid))
               else
                 ibufnext=ibuf+isize(idims)
                 call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,&
                   reshape(source=recvbuffer(ibufnext-isize_stg(idims):ibufnext-1),&
                   shape=(/ ixtemax^d-ixtemin^d+1 ,^nd-1 /)),&
                   idims,iside,.true.,psuse(igrid))
                 ibuf=ibufnext
               end if
            {end do\}
 
           case(neighbor_sibling)
             ! The first neighbour is at the same level
             ! Check if there are corners
             do idir=idims+1,ndim
               pcorner=.false.
               mcorner=.false.
               pi^d=i^d+kr(idir,^d);
               mi^d=i^d-kr(idir,^d);
               ph^d=pi^d-kr(idims,^d)*(2*iside-3);
               mh^d=mi^d-kr(idims,^d)*(2*iside-3);
               if (neighbor_type(pi^d,igrid)==neighbor_fine&
                 .and.neighbor_type(ph^d,igrid)==neighbor_sibling&
                 .and.neighbor_pole(pi^d,igrid)==0) then
                 pcorner(idir)=.true.
                 ! Remove coarse part
                 ! Set indices
                 call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.false.,.true.,0^d&,pcorner,mcorner)
                 ! Remove
                 call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iside,idims,igrid)%edge,idims,iside,.false.,psuse(igrid))
                 ! Add fine part
                 ! Find relative position of finer grid
      {^ifthreed do ix=1,2}
                 inc^d=kr(idims,^d)*3*(iside-1)+3*kr(idir,^d){^ifthreed+kr(6-idir-idims,^d)*ix};
                 ineighbor=neighbor_child(1,inc^d,igrid)
                 ipe_neighbor=neighbor_child(2,inc^d,igrid)
                 iotherside=3-iside
                 ! Set indices
                 call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.true.,.true.,inc^d,pcorner,mcorner)
                 ! add
                 if (ipe_neighbor==mype) then
                   call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iotherside,idims,ineighbor)%edge,idims,iside,.true.,psuse(igrid))
                 else
                   ibufnext_cc=ibuf_cc+isize_stg(idims)
                   call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,&
                     reshape(source=recvbuffer_cc(ibuf_cc:ibufnext_cc-1),&
                     shape=(/ ixtemax^d-ixtemin^d+1 ,^nd-1 /)),&
                     idims,iside,.true.,psuse(igrid))
                   ibuf_cc=ibufnext_cc
                 end if
      {^ifthreed end do}
               ! Set CoCorner to false again for next step
                 pcorner(idir)=.false.
               end if
 
               if (neighbor_type(mi^d,igrid)==neighbor_fine&
                 .and.neighbor_type(mh^d,igrid)==neighbor_sibling&
                 .and.neighbor_pole(mi^d,igrid)==0) then
                   mcorner(idir)=.true.
                   ! Remove coarse part
                   ! Set indices
                   call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.false.,.true.,0^d&,pcorner,mcorner)
                   ! Remove
                   call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iside,idims,igrid)%edge,idims,iside,.false.,psuse(igrid))
                   ! Add fine part
                   ! Find relative position of finer grid
        {^ifthreed do ix=1,2}
                   inc^d=kr(idims,^d)*3*(iside-1){^ifthreed+kr(6-idir-idims,^d)*ix};
                   ineighbor=neighbor_child(1,inc^d,igrid)
                   ipe_neighbor=neighbor_child(2,inc^d,igrid)
                   iotherside=3-iside
                   ! Set indices
                   call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.true.,.true.,inc^d,pcorner,mcorner)
                   ! add
                   if (ipe_neighbor==mype) then
                     call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iotherside,idims,ineighbor)%edge,idims,iside,.true.,psuse(igrid))
                   else
                     ibufnext_cc=ibuf_cc+isize_stg(idims)
                     call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,&
                       reshape(source=recvbuffer_cc(ibuf_cc:ibufnext_cc-1),&
                       shape=(/ ixtemax^d-ixtemin^d+1 ,^nd-1 /)),&
                       idims,iside,.true.,psuse(igrid))
                     ibuf_cc=ibufnext_cc
                   end if
        {^ifthreed end do}
                 ! Set CoCorner to false again for next step
                  mcorner(idir)=.false.
               end if
             end do
           end select
         end do
       end do
     end do
 
     if (nsend_ct>0) call mpi_waitall(nsend_ct,cc_sendreq,cc_sendstat,ierrmpi)
 
   end subroutine fix_edges
 
   !> This routine sets the indexes for the correction
   !> of the circulation according to several different
   !> cases, as grids located in different cpus,
   !> presence of corners, and different relative locations
   !> of the fine grid respect to the coarse one
   subroutine set_ix_circ(ixF^L,ixtE^L,ixE^L,ixfE^L,igrid,idims,iside,add,CoCorner,inc^D,pcorner,mcorner)
     use mod_global_parameters
     
     integer,intent(in)    :: igrid,idims,iside,inc^d
     logical,intent(in)    :: add,cocorner
     logical,intent(inout) :: pcorner(1:ndim),mcorner(1:ndim)
     integer,intent(out)   :: ixf^l(1:ndim),ixte^l,ixe^l(1:3),ixfe^l(1:3) ! Indices for faces and edges
     integer               :: icor^d,idim1,idir,nx^d,middle^d
     integer               :: ixtfe^l
 
     ! ixF -> Indices for the _F_aces, and
     ! depends on the field component
     ! ixtE -> are the _t_otal range of the 'edge' array
     ! ixE -> are the ranges of the edge array,
     ! depending on the component
     ! ixfE -> are the ranges of the fE array (3D),
     ! and also depend on the component
 
     ! ... General ...
     ! Assign indices for the size of the E field array
 
     ixtfemin^d=ixmlo^d-1;
     ixtfemax^d=ixmhi^d;
 
     if(add) then
       nx^d=(ixmhi^d-ixmlo^d+1)/2;
     else
       nx^d=ixmhi^d-ixmlo^d+1;
     end if
 
     do idim1=1,ndim
       ixtemin^d=0;
       ixtemax^d=nx^d;
       select case(idims)
       {case(^d)
         ixtemin^d=1;ixtemax^d=1;
         if (iside==1) ixtfemax^d=ixtfemin^d;
         if (iside==2) ixtfemin^d=ixtfemax^d;
       \}
       end select
     end do
 
     ! Assign indices, considering only the face
     ! (idims and iside)
     do idim1=1,ndim
       ixfmin^d(idim1)=ixmlo^d-kr(idim1,^d);
       ixfmax^d(idim1)=ixmhi^d;
       select case(idims)
       {case(^d)
          select case(iside)
          case(1)
          ixfmax^d(idim1)=ixfmin^d(idim1)
          case(2)
          ixfmin^d(idim1)=ixfmax^d(idim1)
          end select
       \}
       end select
     end do
     ! ... Relative position ...
     ! Restrict range using relative position
     if(add) then
       middle^d=(ixmhi^d+ixmlo^d)/2;
       {
       if(inc^d==1) then
         ixfmax^d(:)=middle^d
         ixtfemax^d=middle^d
       end if
       if(inc^d==2) then
         ixfmin^d(:)=middle^d+1
         ixtfemin^d=middle^d
       end if
       \}
     end if
     ! ... Adjust ranges of edges according to direction ...
     do idim1=1,3
       ixfemax^d(idim1)=ixtfemax^d;
       ixemax^d(idim1)=ixtemax^d;
       ixfemin^d(idim1)=ixtfemin^d+kr(idim1,^d);
       ixemin^d(idim1)=ixtemin^d+kr(idim1,^d);
     end do
     ! ... Corners ...
     ! 'Coarse' corners
     if (cocorner) then
       do idim1=idims+1,ndim
         if (pcorner(idim1)) then
           do idir=1,3!Index arrays have size ndim
             if (idir==6-idim1-idims) then
              !!! Something here has to change
              !!! Array ixfE must have size 3, while
              !!! ixE must have size ndim
              {if (^d==idim1) then
                 ixfemin^d(idir)=ixfemax^d(idir)
                 if (add) then
                   ixemax^d(idir) =ixemin^d(idir)
                 else
                   ixemin^d(idir) =ixemax^d(idir)
                 end if
               end if\}
             else
               ixemin^d(idir)=1;
               ixemax^d(idir)=0;
               ixfemin^d(idir)=1;
               ixfemax^d(idir)=0;
             end if
           end do
         end if
         if (mcorner(idim1)) then
           do idir=1,3
             if (idir==6-idim1-idims) then
              {if (^d==idim1) then
                 ixfemax^d(idir)=ixfemin^d(idir)
                 if (add) then
                   ixemin^d(idir) =ixemax^d(idir)
                 else
                   ixemax^d(idir) =ixemin^d(idir)
                 end if
               end if\}
             else
               ixemin^d(idir)=1;
               ixemax^d(idir)=0;
               ixfemin^d(idir)=1;
               ixfemax^d(idir)=0;
             end if
           end do
         end if
       end do
     else
     ! Other kinds of corners
     ! Crop ranges to account for corners
     ! When the fine fluxes are added, we consider 
     ! whether they come from the same cpu or from
     ! a different one, in order to minimise the 
     ! amount of communication
     ! Case for different processors still not implemented!!!
      {if((idims.gt.^d).and.pcorner(^d)) then
         if((.not.add).or.(inc^d==2)) then
           !ixFmax^DD(:)=ixFmax^DD(:)-kr(^D,^DD);
           do idir=1,3
             if ((idir==idims).or.(idir==^d)) cycle
               ixfemax^d(idir)=ixfemax^d(idir)-1
               ixemax^d(idir)=ixemax^d(idir)-1
           end do
         end if
       end if\}
      {if((idims>^d).and.mcorner(^d)) then
         if((.not.add).or.(inc^d==1)) then
           !ixFmin^DD(:)=ixFmin^DD(:)+kr(^D,^DD);
           do idir=1,3
             if ((idir==idims).or.(idir==^d)) cycle
               ixfemin^d(idir)=ixfemin^d(idir)+1
               ixemin^d(idir)=ixemin^d(idir)+1
           end do
         end if
       end if\}
     end if
 
   end subroutine set_ix_circ
 
   subroutine add_sub_circ(ixF^L,ixtE^L,ixE^L,ixfE^L,edge,idims,iside,add,s)
     use mod_global_parameters
 
     type(state)        :: s
     integer,intent(in) :: idims,iside
     integer            :: ixf^l(1:ndim),ixte^l,ixe^l(1:3),ixfe^l(1:3)
     double precision   :: edge(ixte^s,1:ndim-1)
     logical,intent(in) :: add
 
     double precision   :: fe(ixg^t,sdim:3)
     double precision   :: circ(ixg^t,1:ndim)
     integer            :: idim1,idim2,idir,middle^d
     integer            :: ixfec^l,ixec^l
     integer            :: ix^l,hx^l,ixc^l,hxc^l ! 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(:^d&,:)=zero
     do idim1=1,ndim-1
       {^ifthreed ! 3D: rotate indices (see routine flux_to_edge)
       idir=mod(idim1+idims-1,3)+1}
       {^iftwod   ! 2D: move E back to directon 3
       idir=3}
       ixfec^l=ixfe^l(idir);
       ixec^l=ixe^l(idir);
       fe(ixfec^s,idir)=edge(ixec^s,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
             ixc^l=ixf^l(idim1);
             hxc^l=ixc^l-kr(idim2,^d);
             if(idim1==idims) then
               circ(ixc^s,idim1)=circ(ixc^s,idim1)+lvc(idim1,idim2,idir)&
                                 *(fe(ixc^s,idir)-fe(hxc^s,idir))
             else
               select case(iside)
               case(2)
                 circ(ixc^s,idim1)=circ(ixc^s,idim1)+lvc(idim1,idim2,idir)*fe(ixc^s,idir)
               case(1)
                 circ(ixc^s,idim1)=circ(ixc^s,idim1)-lvc(idim1,idim2,idir)*fe(hxc^s,idir)
               end select
             end if
           end do
        end do
     end do
 
     ! Divide circulation by surface and add
     do idim1=1,ndim
        ixc^l=ixf^l(idim1);
        where(s%surfaceC(ixc^s,idim1)>1.0d-9*s%dvolume(ixc^s))
          circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
        elsewhere
          circ(ixc^s,idim1)=zero
        end where
        ! Add/subtract to field at face
        if (add) then
          s%ws(ixc^s,idim1)=s%ws(ixc^s,idim1)-circ(ixc^s,idim1)
        else
          s%ws(ixc^s,idim1)=s%ws(ixc^s,idim1)+circ(ixc^s,idim1)
        end if
     end do
 
   end subroutine add_sub_circ
 
end module mod_fix_conserve