anchor.diff:          & -J_m(:,ix_m,iy_m) * 4.0*pi*(dt**2) - RR(:)*(c_light*dt)**2 )
       end do
     end do
  +!$omp end parallel do
  +
     ! Boundary condition
  -  Ac_new_m(:,:,NYvacB_m-1)=Ac_new_m(:,:,NYvacB_m)
  -  Ac_new_m(:,:,NYvacT_m+1)=0.0d0
  +!$omp parallel do default(none) &
  +!$    private(ix_m) &
  +!$    shared(Ac_new_m,NXvacL_m,NXvacR_m,NYvacB_m,NYvacT_m)
  +  do ix_m=NXvacL_m-1,NXvacR_m+1
  +    Ac_new_m(:,ix_m,NYvacB_m-1)=Ac_new_m(:,ix_m,NYvacB_m)
  +    Ac_new_m(:,ix_m,NYvacT_m+1)=0.0d0
  +  enddo
  +!$omp end parallel do
     return
   end subroutine dt_evolve_Ac_2dc
   !===========================================================
  
positive.code:   real(8) :: rc(3)  ! rot Ac
    ! calculate the magnetic field from the vector potential (2D case)
    ! (written by M.Uemoto on 2016-11-22)
    do iy_m=NYvacB_m, NYvacT_m
      do ix_m=NXvacL_m, NXvacR_m
        Rc(1) = + (Ac_m(3, ix_m, iy_m+1) - Ac_m(3, ix_m, iy_m-1)) / (2*HY_m)
  
negative.code:       Bmag(:,ix_m,iy_m) = Rc(:) * c_light
      end do
    end do
  !$omp end parallel do
  
    return
  end subroutine calc_bmag_field_2dc
  !===========================================================
  