subroutine abst_2d(x,y,dis)
  ! 2 次元ベクトルの絶対値を計算するルーチン
  implicit none
  real, intent(in) :: x(:,:)  ! x 方向のベクトル成分
  real, intent(in) :: y(size(x,1),size(x,2))  ! y 方向のベクトル成分
  real, intent(inout) :: dis(size(x,1),size(x,2))  ! 各点での絶対値ベクトル
  integer :: i, j, nx, ny

  nx=size(x,1)
  ny=size(x,2)

  do j=1,ny
     do i=1,nx
        dis(i,j)=sqrt(x(i,j)**2+y(i,j)**2)
     end do
  end do

end subroutine abst_2d

subroutine abst_3d(x,y,z,dis)  ! 3 次元ベクトルの絶対値を計算するルーチン
  ! 配列数を調整することにより, 2 次元での計算も可能.
  implicit none
  real, intent(in) :: x(:,:,:)  ! x 方向のベクトル成分
  real, intent(in) :: y(size(x,1),size(x,2),size(x,3))  ! y 方向のベクトル成分
  real, intent(in) :: z(size(x,1),size(x,2),size(x,3))  ! z 方向のベクトル成分
  real, intent(inout) :: dis(size(x,1),size(x,2),size(x,3))  ! 各点での絶対値ベクトル
  integer :: i, j, k, nx, ny, nz

  nx=size(x,1)
  ny=size(x,2)
  nz=size(x,3)

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j,k)
  do k=1,nz
     do j=1,ny
        do i=1,nx
           dis(i,j,k)=sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
        end do
     end do
  end do
!$omp end do
!$omp end parallel

end subroutine abst_3d

subroutine calc_radius(xp,yp,zp,x,y,z,rad)
  ! ある位置からの距離を計算するルーチン
  ! 配列数を調整することにより, 2 次元での計算も可能.
  implicit none
  real, intent(in) :: xp  ! 中心位置座標の x 成分
  real, intent(in) :: yp  ! 中心位置座標の y 成分
  real, intent(in) :: zp  ! 中心位置座標の z 成分
  real, intent(in) :: x(:)  ! x 方向の位置座標
  real, intent(in) :: y(:)  ! y 方向の位置座標
  real, intent(in) :: z(:)  ! z 方向の位置座標
  real, intent(inout) :: rad(size(x),size(y),size(z))  ! 距離
  integer :: i, j, k, nx, ny, nz

  nx=size(x)
  ny=size(y)
  nz=size(z)

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j,k)
  do k=1,nz
     do j=1,ny
        do i=1,nx
           rad(i,j,k)=sqrt((x(i)-xp)**2+(y(j)-yp)**2+(z(k)-zp)**2)
        end do
     end do
  end do
!$omp end do
!$omp end parallel


end subroutine calc_radius

subroutine dot_prod_2d(x,y,u,v,dot,undeff)
  ! 2ベクトルの内積計算ルーチン
  ! dot_prod の 2 次元版
  implicit none
  real, intent(in) :: x(:,:)  ! x 方向のベクトル成分
  real, intent(in) :: y(size(x,1),size(x,2))  ! y 方向のベクトル成分
  real, intent(in) :: u(size(x,1),size(x,2))  ! x 方向のベクトル成分
  real, intent(in) :: v(size(x,1),size(x,2))  ! y 方向のベクトル成分
  real, intent(inout) :: dot(size(x,1),size(x,2))  ! 内積
  real, intent(in), optional :: undeff        ! 未定義値
  integer :: i, j, nx, ny

  nx=size(x,1)
  ny=size(x,2)

  if(present(undeff))then

     do j=1,ny
        do i=1,nx
           if(x(i,j)==undeff.or.u(i,j)==undeff.or. y(i,j)==undeff.or.v(i,j)==undeff)then
              dot(i,j)=undeff
           else
              dot(i,j)=x(i,j)*u(i,j)+y(i,j)*v(i,j)
           end if
        end do
     end do

  else

     do j=1,ny
        do i=1,nx
           dot(i,j)=x(i,j)*u(i,j)+y(i,j)*v(i,j)
        end do
     end do

  end if

end subroutine dot_prod_2d

subroutine dot_prod_3d(x,y,z,u,v,w,dot,undeff)
  ! 2ベクトルの内積計算ルーチン
  ! 配列を工夫すると, 2 次元での内積を計算することも可能
  implicit none
  real, intent(in) :: x(:,:,:)  ! x 方向のベクトル成分
  real, intent(in) :: y(size(x,1),size(x,2),size(x,3))  ! y 方向のベクトル成分
  real, intent(in) :: z(size(x,1),size(x,2),size(x,3))  ! z 方向のベクトル成分
  real, intent(in) :: u(size(x,1),size(x,2),size(x,3))  ! x 方向のベクトル成分
  real, intent(in) :: v(size(x,1),size(x,2),size(x,3))  ! y 方向のベクトル成分
  real, intent(in) :: w(size(x,1),size(x,2),size(x,3))  ! z 方向のベクトル成分
  real, intent(inout) :: dot(size(x,1),size(x,2),size(x,3))  ! 内積
  real, intent(in), optional :: undeff
  integer :: i, j, k, nx, ny, nz

  nx=size(x,1)
  ny=size(x,2)
  nz=size(x,3)

  if(present(undeff))then

!$omp parallel default(shared)
!$omp do schedule(dynamic) private(i,j,k)
     do k=1,nz
        do j=1,ny
           do i=1,nx
              if(x(i,j,k)==undeff.or.u(i,j,k)==undeff.or.y(i,j,k)==undeff.or. v(i,j,k)==undeff.or.z(i,j,k)==undeff.or.w(i,j,k)==undeff)then
                 dot(i,j,k)=undeff
              else
                 dot(i,j,k)=x(i,j,k)*u(i,j,k)+y(i,j,k)*v(i,j,k)+z(i,j,k)*w(i,j,k)
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j,k)
     do k=1,nz
        do j=1,ny
           do i=1,nx
              dot(i,j,k)=x(i,j,k)*u(i,j,k)+y(i,j,k)*v(i,j,k)+z(i,j,k)*w(i,j,k)
           end do
        end do
     end do
!$omp end do
!$omp end parallel
  end if

end subroutine dot_prod_3d

subroutine rectangle_int( x, y, bot, top, res, undeff )  ! 1 次元台形積分
  ! 不等間隔でも計算可能であるが, 精度は保証しない.
  implicit none
  real, intent(in) :: bot  ! 積分区間左端
  real, intent(in) :: top  ! 積分区間右端
  real, intent(in) :: x(:)  ! 積分変数
  real, intent(in) :: y(size(x))  ! 非積分関数
  real, intent(inout) :: res  ! 台形積分の積分値
  real, intent(in), optional :: undeff
  integer :: i, nx, i_bot, i_top
  real :: y_bot, y_top

  nx=size(x)

  res=0.0

!-- bot < top でなければエラー
  if(bot>top)then
     write(*,*) "#### ERROR (algebra:rectangle_int) ####"
     write(*,*) "integrated interval must be bot < top. STOP"
     stop
  end if

!-- 以下で積分域を決定
!-- 実際には, 積分境界に最近する点

  if(x(1)>=bot)then
     if(x(1)>bot)then
        write(*,*) "#### WARNING ####"
        write(*,*) "there is NOT the bot in the x(i)."  ! このときは, i_bot=1としておいても, x(1)=bot なので, 余分計算はゼロとなり影響はない.
     end if
     i_bot=1
  end if
  if(x(nx)<=top)then
     if(x(nx)<top)then
        write(*,*) "#### WARNING ####"
        write(*,*) "there is NOT the top in the x(i)."  ! このとき, i_top=nx としても, x(nx)=top なので, 余分計算はゼロとなる.
     end if
     i_top=nx
  end if

  do i=2,nx-1
     if(x(i)>=bot)then  ! i_bot は bot - top 内の最も bot に近づく格子
        i_bot=i
        exit
     end if
  end do

  do i=nx-1,2,-1
     if(x(i)<=top)then  ! i_top は bot - top 内の最も top に近づく格子
        i_top=i
        exit
     end if
  end do

!-- 以下で格子に当てはまらない部分を内挿で補完
  if(present(undeff))then
     if(i_bot/=1)then
        if(y(i_bot)/=undeff.and.y(i_bot-1)/=undeff.and.x(i_bot)/=x(i_bot-1))then
           y_bot=y(i_bot-1) +((y(i_bot)-y(i_bot-1))/(x(i_bot)-x(i_bot-1)))*(bot-x(i_bot-1))
        else
           y_bot=-y(i_bot)  ! 最後の積分でゼロにするための措置
        end if
     else
        y_bot=-y(i_bot)  ! 最後の積分でゼロにするための措置
     end if
     if(i_top/=nx)then
        if(y(i_top)/=undeff.and.y(i_top+1)/=undeff.and.x(i_top+1)/=x(i_top))then
           y_top=y(i_top) +((y(i_top+1)-y(i_top))/(x(i_top+1)-x(i_top)))*(x(i_top+1)-top)
        else
           y_top=-y(i_top)
        end if
     else
        y_top=-y(i_top)  ! 最後の積分でゼロにするための措置
     end if
  else
     if(i_bot/=1)then
        y_bot=y(i_bot-1) +((y(i_bot)-y(i_bot-1))/(x(i_bot)-x(i_bot-1)))*(bot-x(i_bot-1))
     else
        y_bot=-y(i_bot)  ! 最後の積分でゼロにするための措置
     end if
     if(i_top/=nx)then
        y_top=y(i_top) +((y(i_top+1)-y(i_top))/(x(i_top+1)-x(i_top)))*(x(i_top+1)-top)
     else
        y_top=-y(i_top)  ! 最後の積分でゼロにするための措置
     end if
  end if

  if(i_bot<i_top)then  ! 積分区間内に格子が 2 つ以上あるとき

     if(present(undeff))then
        do i=i_bot,i_top-1
           if(y(i+1)/=undeff.and.y(i)/=undeff)then
              if(i==i_bot)then
                 res=res+0.5*(x(i)-bot)*(y(i)+y_bot) +0.5*(x(i+1)-x(i))*(y(i+1)+y(i))
                     ! 下端の余りと通常の短冊分
              else
                 res=res+0.5*(x(i+1)-x(i))*(y(i+1)+y(i))  ! 通常の短冊
              end if
           end if
        end do
        res=res+0.5*(top-x(i_top))*(y_top+y(i_top))  ! 上端の余りのみ
     else
        do i=i_bot,i_top-1
           if(i==i_bot)then
              res=res+0.5*(x(i)-bot)*(y(i)+y_bot) +0.5*(x(i+1)-x(i))*(y(i+1)+y(i))
                  ! 下端の余りと通常の短冊分
           else
              res=res+0.5*(x(i+1)-x(i))*(y(i+1)+y(i))  ! 通常の短冊
           end if
        end do
        res=res+0.5*(top-x(i_top))*(y_top+y(i_top))  ! 上端の余りのみ
     end if
  else
     if(present(undeff))then
        if(y(i_bot)/=undeff)then
           res=res+0.5*(x(i)-bot)*(y(i)+y_bot) +0.5*(top-x(i))*(y_top+y(i))
        else
           res=0.5*(top-bot)*(y_top+y_bot)
        end if
     else
        res=res+0.5*(x(i_bot)-bot)*(y(i_bot)+y_bot) +0.5*(top-x(i_bot))*(y_top+y(i_bot))
     end if
  end if

end subroutine rectangle_int

subroutine vec_prod_2d(x,y,u,v,vecz,undeff)
  ! 2ベクトルの外積計算ルーチン
  ! vec_prod の 2 次元版
  ! 外積なので, 2 次元で計算すると, 外積の 1 成分しか現れないので,
  ! intent(inout) の引数が 1 つしかないことに注意.
  implicit none
  real, intent(in) :: x(:,:)  ! x 方向のベクトル成分
  real, intent(in) :: y(size(x,1),size(x,2))  ! y 方向のベクトル成分
  real, intent(in) :: u(size(x,1),size(x,2))  ! x 方向のベクトル成分
  real, intent(in) :: v(size(x,1),size(x,2))  ! y 方向のベクトル成分
  real, intent(inout) :: vecz(size(x,1),size(x,2))  ! 外積の直交成分
  real, intent(in), optional :: undeff
  integer :: i, j, nx, ny

  nx=size(x,1)
  ny=size(x,2)

  if(present(undeff))then

        do j=1,ny
           do i=1,nx
              if(x(i,j)==undeff.or.u(i,j)==undeff.or. y(i,j)==undeff.or.v(i,j)==undeff)then
                 vecz(i,j)=undeff
              else
                 vecz(i,j)=x(i,j)*v(i,j)-y(i,j)*u(i,j)
              end if
           end do
        end do

  else

        do j=1,ny
           do i=1,nx
              vecz(i,j)=x(i,j)*v(i,j)-y(i,j)*u(i,j)
           end do
        end do

  end if

end subroutine vec_prod_2d

subroutine vec_prod_3d(x,y,z,u,v,w,vecx,vecy,vecz,undeff)
  ! 2ベクトルの外積計算ルーチン
  implicit none
  real, intent(in) :: x(:,:,:)  ! x 方向のベクトル成分
  real, intent(in) :: y(size(x,1),size(x,2),size(x,3))  ! y 方向のベクトル成分
  real, intent(in) :: z(size(x,1),size(x,2),size(x,3))  ! z 方向のベクトル成分
  real, intent(in) :: u(size(x,1),size(x,2),size(x,3))  ! x 方向のベクトル成分
  real, intent(in) :: v(size(x,1),size(x,2),size(x,3))  ! y 方向のベクトル成分
  real, intent(in) :: w(size(x,1),size(x,2),size(x,3))  ! z 方向のベクトル成分
  real, intent(inout) :: vecx(size(x,1),size(x,2),size(x,3))  ! 外積の x 成分
  real, intent(inout) :: vecy(size(x,1),size(x,2),size(x,3))  ! 外積の y 成分
  real, intent(inout) :: vecz(size(x,1),size(x,2),size(x,3))  ! 外積の z 成分
  real, intent(in), optional :: undeff
  integer :: i, j, k, nx, ny, nz

  nx=size(x,1)
  ny=size(x,2)
  nz=size(x,3)

  if(present(undeff))then

!$omp parallel default(shared)
!$omp do schedule(dynamic) private(i,j,k)
     do k=1,nz
        do j=1,ny
           do i=1,nx
              if(x(i,j,k)==undeff.or.u(i,j,k)==undeff.or.y(i,j,k)==undeff.or. v(i,j,k)==undeff.or.z(i,j,k)==undeff.or.w(i,j,k)==undeff)then
                 vecx(i,j,k)=undeff
                 vecy(i,j,k)=undeff
                 vecz(i,j,k)=undeff
              else
                 vecx(i,j,k)=y(i,j,k)*w(i,j,k)-z(i,j,k)*v(i,j,k)
                 vecy(i,j,k)=z(i,j,k)*u(i,j,k)-x(i,j,k)*w(i,j,k)
                 vecz(i,j,k)=x(i,j,k)*v(i,j,k)-y(i,j,k)*u(i,j,k)
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j,k)
     do k=1,nz
        do j=1,ny
           do i=1,nx
              vecx(i,j,k)=y(i,j,k)*w(i,j,k)-z(i,j,k)*v(i,j,k)
              vecy(i,j,k)=z(i,j,k)*u(i,j,k)-x(i,j,k)*w(i,j,k)
              vecz(i,j,k)=x(i,j,k)*v(i,j,k)-y(i,j,k)*u(i,j,k)
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  end if

end subroutine vec_prod_3d