Class Alge_Solv
In: alge_solv.f90

代数演算を用いて偏微分方程式を解くモジュール

Methods

Public Instance methods

Subroutine :
x(:) :real, intent(in)
: 領域の横座標
y(:) :real, intent(in)
: 領域の縦座標
rho(size(x),size(y)) :real, intent(in)
: ポアソン方程式の強制項 rho =0 でラプラス方程式も求積可能
eps :real, intent(in)
: 収束条件
boundary :character(4), intent(in)
: 境界条件 4 文字で各辺の境界条件を与える. 1 文字目 : x 下端, 2 文字目 : y 左端, 3 文字目 : x 上端, 4 文字目 : y 右端 boundary は 1 : 固定端境界, 2 : 自由端境界, 3 : 周期境界
psi(size(x),size(y)) :real, intent(inout)
: ポアソン方程式の解
bound_opt(size(x),size(y)) :real, intent(in), optional
: 境界での強制 ノイマン境界の場合 : フラックス値
a(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
b(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
c(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
d(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
e(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
undef :real, intent(in), optional
: 未定義値
inner_bound(size(x),size(y)) :integer, intent(in), optional
: 内部領域の境界. 値に応じてその格子点で境界値計算 1 = 固定端境界, 10 = 境界の内側. 2 = y 方向自由端境界 (フラックスは上向き) -2 = y 方向自由端境界 (フラックスは下向き) 4 = x 方向自由端境界 (フラックスは右向き) -4 = x 方向自由端境界 (フラックスは左向き) 3 = 周期境界, -3 = 隅領域で両方とも周期境界 8 = |_, ~| で両方とも自由境界条件 -8 = |~, _| で両方とも自由境界条件 この引数が与えられなければ全領域を計算する. 境界の内側格子点 (10) は反復計算を行わず, undef で設定された値もしくはゼロが入る. このときの境界値は bound_opt の値が用いられる.
init_flag :logical, intent(in), optional
: psi の値をゼロで初期化するか. .true. = 初期化する. .false. = 初期化しない. デフォルトでは初期化する.

ガウス=ザイデル法によるポアソン方程式の求積 各オプション配列は, ポアソン系の各微分項の値. デフォルトはゼロで設定される. $$a\dfrac{partial ^2\psi}{partial x^2} +b\dfrac{partial ^2\psi}{partial x\partial y} +c\dfrac{partial ^2\psi}{partial y^2} +d\dfrac{partial psi}{partial x} +e\dfrac{partial psi}{partial y} =rho $$ の各係数に対応している.

[Source]

subroutine Poisson_GauSei(x, y, rho, eps, boundary, psi, bound_opt, a, b, c, d, e, undef, inner_bound, init_flag )
! ガウス=ザイデル法によるポアソン方程式の求積
! 各オプション配列は, ポアソン系の各微分項の値. デフォルトはゼロで設定される.
! $$a\dfrac{\partial ^2\psi}{\partial x^2} +b\dfrac{\partial ^2\psi}{\partial x\partial y} +c\dfrac{\partial ^2\psi}{\partial y^2} +d\dfrac{\partial \psi}{\partial x} +e\dfrac{\partial \psi}{\partial y} =\rho $$
! の各係数に対応している.
  implicit none
  real, intent(in) :: x(:)  ! 領域の横座標
  real, intent(in) :: y(:)  ! 領域の縦座標
  real, intent(in) :: rho(size(x),size(y))  ! ポアソン方程式の強制項
                   ! rho =0 でラプラス方程式も求積可能
  real, intent(in) :: eps  ! 収束条件
  character(4), intent(in) :: boundary  ! 境界条件
                ! 4 文字で各辺の境界条件を与える.
                ! 1 文字目 : x 下端, 2 文字目 : y 左端, 3 文字目 : x 上端,
                ! 4 文字目 : y 右端
                ! boundary は 1 : 固定端境界, 2 : 自由端境界, 3 : 周期境界
  real, intent(in), optional :: bound_opt(size(x),size(y))  ! 境界での強制
                             ! ノイマン境界の場合 : フラックス値
  real, intent(in), optional :: a(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: b(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: c(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: d(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: e(size(x),size(y))  ! 各微分項の係数
  real, intent(inout) :: psi(size(x),size(y))  ! ポアソン方程式の解
  real, intent(in), optional :: undef  ! 未定義値
  integer, intent(in), optional :: inner_bound(size(x),size(y))
                             ! 内部領域の境界. 値に応じてその格子点で境界値計算
                             ! 1 = 固定端境界, 10 = 境界の内側.
                             ! 2 = y 方向自由端境界 (フラックスは上向き)
                             ! -2 = y 方向自由端境界 (フラックスは下向き)
                             ! 4 = x 方向自由端境界 (フラックスは右向き)
                             ! -4 = x 方向自由端境界 (フラックスは左向き)
                             ! 3 = 周期境界, -3 = 隅領域で両方とも周期境界
                             ! 8 = |_, ~| で両方とも自由境界条件
                             ! -8 = |~, _| で両方とも自由境界条件
                             ! この引数が与えられなければ全領域を計算する.
                             ! 境界の内側格子点 (10) は反復計算を行わず,
                             ! undef で設定された値もしくはゼロが入る.
                             ! このときの境界値は bound_opt の値が用いられる.
  logical, intent(in), optional :: init_flag  ! psi の値をゼロで初期化するか.
                             ! .true. = 初期化する. .false. = 初期化しない.
                             ! デフォルトでは初期化する.
  integer :: i, j
  integer :: nx  ! x 方向の配列要素
  integer :: ny  ! y 方向の配列要素
  integer :: signb  ! 各係数を計算するかどうか
  integer, dimension(size(x),size(y)) :: ib

  real :: defun
  real :: tmp, err, err_max
  real :: tmp_b
  real :: bnd(size(x),size(y))
  real :: dx(size(x)), dy(size(y)), dx2(size(x)), dy2(size(y))
  real, dimension(size(x),size(y)) :: dxdy
  real, dimension(size(x),size(y)) :: at, bt, ct, dt, et
  real, dimension(size(x),size(y)) :: adp, adm, cep, cem, ac, divi

  character(4) :: bound
  logical, dimension(size(x),size(y)) :: inner_flag

  bound(1:4)=boundary(1:4)

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

!-- 応答関数の初期化

  if(present(init_flag))then
     if(init_flag.eqv..true.)then
        psi = 0.0
     end if
  else
     psi = 0.0
  end if

!-- 内部境界の判別フラグの設定

  if(present(inner_bound))then
     call set_bound( bound, ib, inner_flag, inner_bound )
  else
     call set_bound( bound, ib, inner_flag )
  end if

!-- 領域・内部境界における境界値の設定

  if(present(bound_opt))then
     call setval_bound( ib, bnd, psi, bound_opt )
  else
     call setval_bound( ib, bnd, psi )
  end if

!-- 未定義値の設定

  if(present(undef))then
     defun=undef
  else
     defun=0.0
  end if

!-- 係数の代入
!-- a, c については, 値が入れられていなければ, 全配列に 1 を代入する.
  if(present(a))then
     call set_coe( at, ext=a )
  else
     call set_coe( at, def=1.0 )
  end if

  if(present(c))then
     call set_coe( ct, ext=c )
  else
     call set_coe( ct, def=1.0 )
  end if

  if(present(b))then
     call set_coe( bt, ext=b )
     signb=1
  else
     call set_coe( bt, def=0.0 )
     signb=0
  end if

  if(present(d))then
     call set_coe( dt, ext=d )
  else
     call set_coe( dt, def=0.0 )
  end if

  if(present(e))then
     call set_coe( et, ext=e )
  else
     call set_coe( et, def=0.0 )
  end if

!-- 以下で先に格子間隔等の 1 回計算でよいものを求めておく.
!-- これらは 1 方向のみで変化すればよい.
!-- 格子点間隔の計算
  do i=2,nx-1
     dx(i)=(x(i+1)-x(i-1))*0.5
     dx2(i)=dx(i)**2
  end do
  do j=2,ny-1
     dy(j)=(y(j+1)-y(j-1))*0.5
     dy2(j)=dy(j)**2
  end do

  dx(1)=(x(2)-x(1))
  dx(nx)=(x(nx)-x(nx-1))
  dy(1)=(y(2)-y(1))
  dy(ny)=(y(ny)-y(ny-1))

  do j=1,ny
     do i=1,nx
        dxdy(i,j)=dx(i)*dy(j)
     end do
  end do

!-- ポアソン係数の計算
!-- 実際にポアソン係数が使用されるのは, 境界を除いた 1 回り内側なので,
!-- 計算量削減のため, ループをそのようにしておく.

!-- 最高次数係数 ac の計算
!$omp parallel default(shared)
!$omp do private(i,j)

  do j=2,ny-1
     do i=2,nx-1
        ac(i,j)=2.0*(at(i,j)/dx2(i)+ct(i,j)/dy2(j))
        adp(i,j)=(at(i,j)/(dx2(i))+0.5*dt(i,j)/dx(i))/ac(i,j)
        adm(i,j)=(at(i,j)/(dx2(i))-0.5*dt(i,j)/dx(i))/ac(i,j)
        cep(i,j)=(ct(i,j)/(dy2(j))+0.5*et(i,j)/dy(j))/ac(i,j)
        cem(i,j)=(ct(i,j)/(dy2(j))-0.5*et(i,j)/dy(j))/ac(i,j)
        bt(i,j)=0.25*bt(i,j)/(dxdy(i,j))/ac(i,j)
     end do
  end do

!$omp end do
!$omp end parallel

  err_max=eps  ! while に入るための便宜的措置

!-- 実際のソルバ ---
  do while(err_max>=eps)
     err_max=0.0

     do j=2,ny-1
        do i=2,nx-1

           tmp=-rho(i,j)/ac(i,j)
           divi(i,j)=1.0

!-- 以降, 反復計算に必要な点の周囲 8 点についてそれぞれ
!-- inner_flag のチェックと同時に適切な値をそれぞれ逐次計算する,
!-- 下で, 各方向の格子計算を行っているが, ib=8,-8 は 4 隅格子のときにしか
!-- case select しない. なぜなら, 真横や上下に隅格子がくることはありえない.
           if(inner_flag(i,j).eqv..false.)then  ! .false. なら領域計算開始
           !-- 右横格子
              if(inner_flag(i+1,j).eqv..false.)then
                 tmp=tmp+adp(i,j)*psi(i+1,j)
              else
                 select case (ib(i+1,j))
                 case(1)
                    tmp=tmp+adp(i,j)*bnd(i+1,j)
                 case(-2)  ! 右横格子は 2 はありえない.
                    tmp=tmp+adp(i,j)*bnd(i+1,j)*dx(i+1)
                    divi(i,j)=divi(i,j)-adp(i,j)
                 case(3)
                    tmp=tmp+adp(i,j)*psi(2,j)
                 case(4)
                    tmp=tmp+adp(i,j)*(psi(i+1,j+1)-bnd(i+1,j)*dy(j))
                 case(-4)
                    tmp=tmp+adp(i,j)*(psi(i+1,j-1)+bnd(i+1,j)*dy(j))
                 end select
              end if

           !-- 左横格子
              if(inner_flag(i-1,j).eqv..false.)then
                 tmp=tmp+adm(i,j)*psi(i-1,j)
              else
                 select case (ib(i-1,j))
                 case(1)
                    tmp=tmp+adm(i,j)*bnd(i-1,j)
                 case(2)  ! 左横格子は -2 はありえない.
                    tmp=tmp-adm(i,j)*bnd(i-1,j)*dx(i-1)
                    divi(i,j)=divi(i,j)-adm(i,j)
                 case(3)
                    tmp=tmp+adm(i,j)*psi(nx-1,j)
                 case(4)
                    tmp=tmp+adm(i,j)*(psi(i-1,j+1)-bnd(i-1,j)*dy(j))
                 case(-4)
                    tmp=tmp+adm(i,j)*(psi(i-1,j-1)+bnd(i-1,j)*dy(j))
                 end select
              end if

           !-- 真上格子
              if(inner_flag(i,j+1).eqv..false.)then
                 tmp=tmp+cep(i,j)*psi(i,j+1)
              else
                 select case (ib(i,j+1))
                 case(1)
                    tmp=tmp+cep(i,j)*bnd(i,j+1)
                 case(2)
                    tmp=tmp+cep(i,j)*(psi(i+1,j+1)-bnd(i,j+1)*dx(i))
                 case(-2)
                    tmp=tmp+cep(i,j)*(psi(i-1,j+1)+bnd(i,j+1)*dx(i))
                 case(3)
                    tmp=tmp+cep(i,j)*psi(i,2)
                 case(-4)  ! 真上格子は 4 はありえない.
                    tmp=tmp+cep(i,j)*bnd(i,j+1)*dy(j+1)
                    divi(i,j)=divi(i,j)-cep(i,j)
                 end select
              end if

           !-- 真下格子
              if(inner_flag(i,j-1).eqv..false.)then
                 tmp=tmp+cem(i,j)*psi(i,j-1)
              else
                 select case (ib(i,j-1))
                 case(1)
                    tmp=tmp+cem(i,j)*bnd(i,j-1)
                 case(2)
                    tmp=tmp+cem(i,j)*(psi(i+1,j-1)-bnd(i,j-1)*dx(i))
                 case(-2)
                    tmp=tmp+cem(i,j)*(psi(i-1,j-1)+bnd(i,j-1)*dx(i))
                 case(3)
                    tmp=tmp+cem(i,j)*psi(i,ny-1)
                 case(4)  ! 真下格子は -4 はありえない.
                    tmp=tmp-cem(i,j)*bnd(i,j-1)*dy(j-1)
                    divi(i,j)=divi(i,j)-cem(i,j)
                 end select
              end if

           !-- 4 隅格子
              if(signb==1)then  ! そもそも bt = 0 なら計算しない.
                 tmp_b=0.0
                 select case (ib(i-1,j-1))  ! 左下格子
                 case(1)
                    tmp_b=tmp_b+bnd(i-1,j-1)
                 case(2)
                    tmp_b=tmp_b+psi(i,j-1)-bnd(i-1,j-1)*dx(i-1)
                 case(-2)
                    tmp_b=tmp_b+psi(i-2,j-1)+bnd(i-1,j-1)*dx(i-1)
                 case(4)
                    tmp_b=tmp_b+psi(i-1,j)-bnd(i-1,j-1)*dy(j-1)
                 case(-4)
                    tmp_b=tmp_b+psi(i-1,j-2)+bnd(i-1,j-1)*dy(j-1)
                 case(3)  ! 左下で周期境界なら, i==2 か j==2 しかありえない.
                    if(i==2)then
                       tmp_b=tmp_b+psi(nx-1,j-1)
                    else if(j==2)then
                       tmp_b=tmp_b+psi(i-1,ny-1)
                    end if
                 case(-3)  ! 左下で -3 なら, 領域の左下隅しかない.
                    tmp_b=tmp_b+psi(nx-1,ny-1)
                 case(8)  ! 左下では -8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b=tmp_b -0.5*(bnd(i,j-1)*dy(j-1)+bnd(i-1,j)*dx(i-1))

                 end select
                 
                 select case (ib(i+1,j+1))  ! 右上格子
                 case(1)
                    tmp_b=tmp_b+bnd(i+1,j+1)
                 case(2)
                    tmp_b=tmp_b+psi(i,j+2)-bnd(i+1,j+1)*dx(i+1)
                 case(-2)
                    tmp_b=tmp_b+psi(i,j+1)+bnd(i+1,j+1)*dx(i+1)
                 case(4)
                    tmp_b=tmp_b+psi(i+1,j+2)-bnd(i+1,j+1)*dy(j+1)
                 case(-4)
                    tmp_b=tmp_b+psi(i+1,j)+bnd(i+1,j+1)*dy(j+1)
                 case(3)  ! 右上の場合, 周期境界は i==nx-1 か j==ny-1 しかない.
                    if(i==nx-1)then
                       tmp_b=tmp_b+psi(2,j-1)
                    else if(j==ny-1)then
                       tmp_b=tmp_b+psi(i-1,2)
                    end if
                 case(-3)  ! 右上で -3 なら, 領域の右上隅しかない.
                    tmp_b=tmp_b+psi(2,2)
                 case(8)  ! 右上では -8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b=tmp_b +0.5*(bnd(i,j+1)*dy(j+1)+bnd(i+1,j)*dx(i+1))

                 end select

                 select case (ib(i-1,j+1))  ! 左上格子
                 case(1)
                    tmp_b=tmp_b+bnd(i-1,j+1)
                 case(2)
                    tmp_b=tmp_b+psi(i,j+1)-bnd(i-1,j+1)*dx(i-1)
                 case(-2)
                    tmp_b=tmp_b+psi(i-2,j+1)+bnd(i-1,j+1)*dx(i-1)
                 case(4)
                    tmp_b=tmp_b+psi(i-1,j+2)-bnd(i-1,j+1)*dy(j+1)
                 case(-4)
                    tmp_b=tmp_b+psi(i-1,j)+bnd(i-1,j+1)*dy(j+1)
                 case(3)  ! 左上で周期境界なら, i==2 か j==ny-1 しかありえない.
                    if(i==2)then
                       tmp_b=tmp_b+psi(nx-1,j-1)
                    else if(j==ny-1)then
                       tmp_b=tmp_b+psi(i-1,2)
                    end if
                 case(-3)  ! 左上で -3 なら, 領域の左上隅しかない.
                    tmp_b=tmp_b+psi(nx-1,2)
                 case(-8)  ! 左上では 8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b=tmp_b +0.5*(bnd(i,j+1)*dy(j+1)-bnd(i-1,j)*dx(i-1))

                 end select

                 select case (ib(i+1,j-1))  ! 右下格子
                 case(1)
                    tmp_b=tmp_b+bnd(i+1,j-1)
                 case(2)
                    tmp_b=tmp_b+psi(i+2,j-1)-bnd(i+1,j-1)*dx(i+1)
                 case(-2)
                    tmp_b=tmp_b+psi(i,j-1)+bnd(i+1,j-1)*dx(i+1)
                 case(4)
                    tmp_b=tmp_b+psi(i-1,j)-bnd(i+1,j-1)*dy(j-1)
                 case(-4)
                    tmp_b=tmp_b+psi(i-1,j-2)+bnd(i+1,j-1)*dy(j-1)
                 case(3)  ! 右下で周期境界なら, i==nx-1 か j==2 しかありえない.
                    if(i==nx-1)then
                       tmp_b=tmp_b+psi(2,j-1)
                    else if(j==2)then
                       tmp_b=tmp_b+psi(i+1,ny-1)
                    end if
                 case(-3)  ! 右下で -3 なら, 領域の右下隅しかない.
                    tmp_b=tmp_b+psi(2,ny-1)
                 case(-8)  ! 右下では 8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b=tmp_b +0.5*(-bnd(i,j-1)*dy(j-1)+bnd(i+1,j)*dx(i+1))

                 end select
              end if

              tmp=tmp+bt(i,j)*tmp_b
              tmp=tmp/divi(i,j)

           end if

           err=abs(tmp-psi(i,j))

!-- 最大誤差の更新
           if(err_max<=err)then
              err_max=err
           end if

           psi(i,j)=tmp

        end do
     end do

  end do

!-- 境界の設定

  call calculate_bound( ib, dx, dy, bnd, psi )

!-- 未定義領域には undef を代入する.

  do j=1,ny
     do i=1,nx
        if(ib(i,j)==10)then
           psi(i,j)=defun
        end if
     end do
  end do

end subroutine Poisson_GauSei
Subroutine :
x(:) :real, intent(in)
: 領域の横座標
y(:) :real, intent(in)
: 領域の縦座標
rho(size(x),size(y)) :real, intent(in)
: ポアソン方程式の強制項 rho =0 でラプラス方程式も求積可能
eps :real, intent(in)
: 収束条件
boundary :character(4), intent(in)
: 境界条件 4 文字で各辺の境界条件を与える. 1 文字目 : x 下端, 2 文字目 : y 左端, 3 文字目 : x 上端, 4 文字目 : y 右端 boundary は 1 : 固定端境界, 2 : 自由端境界, 3 : 周期境界
psi(size(x),size(y)) :real, intent(inout)
: ポアソン方程式の解
bound_opt(size(x),size(y)) :real, intent(in), optional
: 境界での強制 ノイマン境界の場合 : フラックス値
a(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
b(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
c(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
d(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
e(size(x),size(y)) :real, intent(in), optional
: 各微分項の係数
undef :real, intent(in), optional
: 未定義値
inner_bound(size(x),size(y)) :integer, intent(in), optional
: 内部領域の境界. 値に応じてその格子点で境界値計算 1 = 固定端境界, 10 = 境界の内側. 2 = y 方向自由端境界 (フラックスは上向き) -2 = y 方向自由端境界 (フラックスは下向き) 4 = x 方向自由端境界 (フラックスは右向き) -4 = x 方向自由端境界 (フラックスは左向き) 3 = 周期境界, -3 = 隅領域で両方とも周期境界 8 = |_, ~| で両方とも自由境界条件 -8 = |~, _| で両方とも自由境界条件 この引数が与えられなければ全領域を計算する. 境界の内側格子点 (10) は反復計算を行わず, undef で設定された値もしくはゼロが入る. このときの境界値は bound_opt の値が用いられる.
init_flag :logical, intent(in), optional
: psi の値をゼロで初期化するか. .true. = 初期化する. .false. = 初期化しない. デフォルトでは初期化する.

openmp によるスレッド並列が可能. ガウス・ザイデル法ではアルゴリズムの点から並列化が困難と思われたので, 並列計算によるポアソン方程式の求積が必要となるなら, ヤコビ法のものを使用されたい.

[Source]

subroutine Poisson_Jacobi(x, y, rho, eps, boundary, psi, bound_opt, a, b, c, d, e, undef, inner_bound, init_flag )
  ! openmp によるスレッド並列が可能.
  ! ガウス・ザイデル法ではアルゴリズムの点から並列化が困難と思われたので,
  ! 並列計算によるポアソン方程式の求積が必要となるなら,
  ! ヤコビ法のものを使用されたい.
  implicit none
  real, intent(in) :: x(:)  ! 領域の横座標
  real, intent(in) :: y(:)  ! 領域の縦座標
  real, intent(in) :: rho(size(x),size(y))  ! ポアソン方程式の強制項
                   ! rho =0 でラプラス方程式も求積可能
  real, intent(in) :: eps  ! 収束条件
  character(4), intent(in) :: boundary  ! 境界条件
                ! 4 文字で各辺の境界条件を与える.
                ! 1 文字目 : x 下端, 2 文字目 : y 左端, 3 文字目 : x 上端,
                ! 4 文字目 : y 右端
                ! boundary は 1 : 固定端境界, 2 : 自由端境界, 3 : 周期境界
  real, intent(in), optional :: bound_opt(size(x),size(y))  ! 境界での強制
                             ! ノイマン境界の場合 : フラックス値
  real, intent(in), optional :: a(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: b(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: c(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: d(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: e(size(x),size(y))  ! 各微分項の係数
  real, intent(in), optional :: undef  ! 未定義値
  integer, intent(in), optional :: inner_bound(size(x),size(y))
                             ! 内部領域の境界. 値に応じてその格子点で境界値計算
                             ! 1 = 固定端境界, 10 = 境界の内側.
                             ! 2 = y 方向自由端境界 (フラックスは上向き)
                             ! -2 = y 方向自由端境界 (フラックスは下向き)
                             ! 4 = x 方向自由端境界 (フラックスは右向き)
                             ! -4 = x 方向自由端境界 (フラックスは左向き)
                             ! 3 = 周期境界, -3 = 隅領域で両方とも周期境界
                             ! 8 = |_, ~| で両方とも自由境界条件
                             ! -8 = |~, _| で両方とも自由境界条件
                             ! この引数が与えられなければ全領域を計算する.
                             ! 境界の内側格子点 (10) は反復計算を行わず,
                             ! undef で設定された値もしくはゼロが入る.
                             ! このときの境界値は bound_opt の値が用いられる.
  real, intent(inout) :: psi(size(x),size(y))  ! ポアソン方程式の解
  logical, intent(in), optional :: init_flag  ! psi の値をゼロで初期化するか.
                             ! .true. = 初期化する. .false. = 初期化しない.
                             ! デフォルトでは初期化する.

  integer :: i, j
  integer :: nx  ! x 方向の配列要素
  integer :: ny  ! y 方向の配列要素
  integer :: signb  ! クロスターム b が存在するか
  integer, dimension(size(x),size(y)) :: ib

  real :: defun
  real :: err, err_max
  real :: bnd(size(x),size(y))
  real :: dx(size(x)), dy(size(y)), dx2(size(x)), dy2(size(y))
  real, dimension(size(x),size(y)) :: dxdy
  real, dimension(size(x),size(y)) :: at, bt, ct, dt, et
  real, dimension(size(x),size(y)) :: adp, adm, cep, cem, ac
  real, dimension(size(x),size(y)) :: tmp, tmp_b, divi

  character(4) :: bound
  logical, dimension(size(x),size(y)) :: inner_flag

  bound(1:4)=boundary(1:4)

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

!-- 応答関数の初期化

  if(present(init_flag))then
     if(init_flag.eqv..true.)then
        psi = 0.0
     end if
  else
     psi = 0.0
  end if

!-- 内部境界の判別フラグの設定

  if(present(inner_bound))then
     call set_bound( bound, ib, inner_flag, inner_bound )
  else
     call set_bound( bound, ib, inner_flag )
  end if

!-- 領域・内部境界における境界値の設定

  if(present(bound_opt))then
     call setval_bound( ib, bnd, psi, bound_opt )
  else
     call setval_bound( ib, bnd, psi )
  end if

!-- 未定義値の設定

  if(present(undef))then
     defun=undef
  else
     defun=0.0
  end if

!-- 係数の代入
!-- a, c については, 値が入れられていなければ, 全配列に 1 を代入する.

  if(present(a))then
     call set_coe( at, ext=a )
  else
     call set_coe( at, def=1.0 )
  end if

  if(present(c))then
     call set_coe( ct, ext=c )
  else
     call set_coe( ct, def=1.0 )
  end if

  if(present(b))then
     call set_coe( bt, ext=b )
     signb=1
  else
     call set_coe( bt, def=0.0 )
     signb=0
  end if

  if(present(d))then
     call set_coe( dt, ext=d )
  else
     call set_coe( dt, def=0.0 )
  end if

  if(present(e))then
     call set_coe( et, ext=e )
  else
     call set_coe( et, def=0.0 )
  end if

!-- 以下で先に格子間隔等の 1 回計算でよいものを求めておく.
!-- これらは 1 方向のみで変化すればよい.
!-- 格子点間隔の計算
  do i=2,nx-1
     dx(i)=(x(i+1)-x(i-1))*0.5
     dx2(i)=dx(i)**2
  end do
  do j=2,ny-1
     dy(j)=(y(j+1)-y(j-1))*0.5
     dy2(j)=dy(j)**2
  end do

  dx(1)=(x(2)-x(1))
  dx(nx)=(x(nx)-x(nx-1))
  dy(1)=(y(2)-y(1))
  dy(ny)=(y(ny)-y(ny-1))

  do j=1,ny
     do i=1,nx
        dxdy(i,j)=dx(i)*dy(j)
     end do
  end do

!-- ポアソン係数の計算
!-- 実際にポアソン係数が使用されるのは, 境界を除いた 1 回り内側なので,
!-- 計算量削減のため, ループをそのようにしておく.

!-- 最高次数係数 ac の計算
!$omp parallel default(shared)
!$omp do private(i,j)

  do j=2,ny-1
     do i=2,nx-1
        ac(i,j)=2.0*(at(i,j)/dx2(i)+ct(i,j)/dy2(j))
        adp(i,j)=(at(i,j)/(dx2(i))+0.5*dt(i,j)/dx(i))/ac(i,j)
        adm(i,j)=(at(i,j)/(dx2(i))-0.5*dt(i,j)/dx(i))/ac(i,j)
        cep(i,j)=(ct(i,j)/(dy2(j))+0.5*et(i,j)/dy(j))/ac(i,j)
        cem(i,j)=(ct(i,j)/(dy2(j))-0.5*et(i,j)/dy(j))/ac(i,j)
        bt(i,j)=0.25*bt(i,j)/(dxdy(i,j))/ac(i,j)
     end do
  end do

!$omp end do
!$omp end parallel

  err_max=eps  ! while に入るための便宜的措置

!-- 実際のソルバ ---
  do while(err_max>=eps)
     err_max=0.0
!$omp parallel default(shared)
!$omp do schedule(dynamic) private(i,j)
     do j=2,ny-1
        do i=2,nx-1

           tmp(i,j)=-rho(i,j)/ac(i,j)
           divi(i,j)=1.0

!-- 以降, 反復計算に必要な点の周囲 8 点についてそれぞれ
!-- inner_flag のチェックと同時に適切な値をそれぞれ逐次計算する,
!-- 下で, 各方向の格子計算を行っているが, ib=8,-8 は 4 隅格子のときにしか
!-- case select しない. なぜなら, 真横や上下に隅格子がくることはありえない.
           if(inner_flag(i,j).eqv..false.)then  ! .false. なら領域計算開始
           !-- 右横格子
              if(inner_flag(i+1,j).eqv..false.)then
                 tmp(i,j)=tmp(i,j)+adp(i,j)*psi(i+1,j)
              else
                 select case (ib(i+1,j))
                 case(1)
                    tmp(i,j)=tmp(i,j)+adp(i,j)*bnd(i+1,j)
                 case(-2)  ! 右横格子は 2 はありえない.
                    tmp(i,j)=tmp(i,j)+adp(i,j)*bnd(i+1,j)*dx(i+1)
                    divi(i,j)=divi(i,j)-adp(i,j)
                 case(3)
                    tmp(i,j)=tmp(i,j)+adp(i,j)*psi(2,j)
                 case(4)
                    tmp(i,j)=tmp(i,j)+adp(i,j)*(psi(i+1,j+1)-bnd(i+1,j)*dy(j))
                 case(-4)
                    tmp(i,j)=tmp(i,j)+adp(i,j)*(psi(i+1,j-1)+bnd(i+1,j)*dy(j))
                 end select
              end if

           !-- 左横格子
              if(inner_flag(i-1,j).eqv..false.)then
                 tmp(i,j)=tmp(i,j)+adm(i,j)*psi(i-1,j)
              else
                 select case (ib(i-1,j))
                 case(1)
                    tmp(i,j)=tmp(i,j)+adm(i,j)*bnd(i-1,j)
                 case(2)  ! 左横格子は -2 はありえない.
                    tmp(i,j)=tmp(i,j)-adm(i,j)*bnd(i-1,j)*dx(i-1)
                    divi(i,j)=divi(i,j)-adm(i,j)
                 case(3)
                    tmp(i,j)=tmp(i,j)+adm(i,j)*psi(nx-1,j)
                 case(4)
                    tmp(i,j)=tmp(i,j)+adm(i,j)*(psi(i-1,j+1)-bnd(i-1,j)*dy(j))
                 case(-4)
                    tmp(i,j)=tmp(i,j)+adm(i,j)*(psi(i-1,j-1)+bnd(i-1,j)*dy(j))
                 end select
              end if

           !-- 真上格子
              if(inner_flag(i,j+1).eqv..false.)then
                 tmp(i,j)=tmp(i,j)+cep(i,j)*psi(i,j+1)
              else
                 select case (ib(i,j+1))
                 case(1)
                    tmp(i,j)=tmp(i,j)+cep(i,j)*bnd(i,j+1)
                 case(2)
                    tmp(i,j)=tmp(i,j)+cep(i,j)*(psi(i+1,j+1)-bnd(i,j+1)*dx(i))
                 case(-2)
                    tmp(i,j)=tmp(i,j)+cep(i,j)*(psi(i-1,j+1)+bnd(i,j+1)*dx(i))
                 case(3)
                    tmp(i,j)=tmp(i,j)+cep(i,j)*psi(i,2)
                 case(-4)  ! 真上格子は 4 はありえない.
                    tmp(i,j)=tmp(i,j)+cep(i,j)*bnd(i,j+1)*dy(j+1)
                    divi(i,j)=divi(i,j)-cep(i,j)
                 end select
              end if

           !-- 真下格子
              if(inner_flag(i,j-1).eqv..false.)then
                 tmp(i,j)=tmp(i,j)+cem(i,j)*psi(i,j-1)
              else
                 select case (ib(i,j-1))
                 case(1)
                    tmp(i,j)=tmp(i,j)+cem(i,j)*bnd(i,j-1)
                 case(2)
                    tmp(i,j)=tmp(i,j)+cem(i,j)*(psi(i+1,j-1)-bnd(i,j-1)*dx(i))
                 case(-2)
                    tmp(i,j)=tmp(i,j)+cem(i,j)*(psi(i-1,j-1)+bnd(i,j-1)*dx(i))
                 case(3)
                    tmp(i,j)=tmp(i,j)+cem(i,j)*psi(i,ny-1)
                 case(4)  ! 真下格子は -4 はありえない.
                    tmp(i,j)=tmp(i,j)-cem(i,j)*bnd(i,j-1)*dy(j-1)
                    divi(i,j)=divi(i,j)-cem(i,j)
                 end select
              end if

           !-- 4 隅格子
              if(signb==1)then  ! そもそも bt = 0 なら計算しない.
                 tmp_b(i,j)=0.0
                 select case (ib(i-1,j-1))  ! 左下格子
                 case(1)
                    tmp_b(i,j)=tmp_b(i,j)+bnd(i-1,j-1)
                 case(2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i,j-1)-bnd(i-1,j-1)*dx(i-1)
                 case(-2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-2,j-1)+bnd(i-1,j-1)*dx(i-1)
                 case(4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j)-bnd(i-1,j-1)*dy(j-1)
                 case(-4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j-2)+bnd(i-1,j-1)*dy(j-1)
                 case(3)  ! 左下で周期境界なら, i==2 か j==2 しかありえない.
                    if(i==2)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,j-1)
                    else if(j==2)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(i-1,ny-1)
                    end if
                 case(-3)  ! 左下で -3 なら, 領域の左下隅しかない.
                    tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,ny-1)
                 case(8)  ! 左下では -8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b(i,j)=tmp_b(i,j) -0.5*(bnd(i,j-1)*dy(j-1)+bnd(i-1,j)*dx(i-1))

                 end select
                 
                 select case (ib(i+1,j+1))  ! 右上格子
                 case(1)
                    tmp_b(i,j)=tmp_b(i,j)+bnd(i+1,j+1)
                 case(2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i,j+2)-bnd(i+1,j+1)*dx(i+1)
                 case(-2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i,j+1)+bnd(i+1,j+1)*dx(i+1)
                 case(4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i+1,j+2)-bnd(i+1,j+1)*dy(j+1)
                 case(-4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i+1,j)+bnd(i+1,j+1)*dy(j+1)
                 case(3)  ! 右上の場合, 周期境界は i==nx-1 か j==ny-1 しかない.
                    if(i==nx-1)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(2,j-1)
                    else if(j==ny-1)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(i-1,2)
                    end if
                 case(-3)  ! 右上で -3 なら, 領域の右上隅しかない.
                    tmp_b(i,j)=tmp_b(i,j)+psi(2,2)
                 case(8)  ! 右上では -8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b(i,j)=tmp_b(i,j) +0.5*(bnd(i,j+1)*dy(j+1)+bnd(i+1,j)*dx(i+1))

                 end select

                 select case (ib(i-1,j+1))  ! 左上格子
                 case(1)
                    tmp_b(i,j)=tmp_b(i,j)+bnd(i-1,j+1)
                 case(2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i,j+1)-bnd(i-1,j+1)*dx(i-1)
                 case(-2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-2,j+1)+bnd(i-1,j+1)*dx(i-1)
                 case(4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j+2)-bnd(i-1,j+1)*dy(j+1)
                 case(-4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j)+bnd(i-1,j+1)*dy(j+1)
                 case(3)  ! 左上で周期境界なら, i==2 か j==ny-1 しかありえない.
                    if(i==2)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,j-1)
                    else if(j==ny-1)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(i-1,2)
                    end if
                 case(-3)  ! 左上で -3 なら, 領域の左上隅しかない.
                    tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,2)
                 case(-8)  ! 左上では 8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b(i,j)=tmp_b(i,j) +0.5*(bnd(i,j+1)*dy(j+1)-bnd(i-1,j)*dx(i-1))

                 end select

                 select case (ib(i+1,j-1))  ! 右下格子
                 case(1)
                    tmp_b(i,j)=tmp_b(i,j)+bnd(i+1,j-1)
                 case(2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i+2,j-1)-bnd(i+1,j-1)*dx(i+1)
                 case(-2)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i,j-1)+bnd(i+1,j-1)*dx(i+1)
                 case(4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j)-bnd(i+1,j-1)*dy(j-1)
                 case(-4)
                    tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j-2)+bnd(i+1,j-1)*dy(j-1)
                 case(3)  ! 右下で周期境界なら, i==nx-1 か j==2 しかありえない.
                    if(i==nx-1)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(2,j-1)
                    else if(j==2)then
                       tmp_b(i,j)=tmp_b(i,j)+psi(i+1,ny-1)
                    end if
                 case(-3)  ! 右下で -3 なら, 領域の右下隅しかない.
                    tmp_b(i,j)=tmp_b(i,j)+psi(2,ny-1)
                 case(-8)  ! 右下では 8 は存在しない
                    divi(i,j)=divi(i,j)-bt(i,j)
                    tmp_b(i,j)=tmp_b(i,j) +0.5*(-bnd(i,j-1)*dy(j-1)+bnd(i+1,j)*dx(i+1))

                 end select
              end if

              tmp(i,j)=tmp(i,j)+bt(i,j)*tmp_b(i,j)
              tmp(i,j)=tmp(i,j)/divi(i,j)

           end if
        end do
     end do
!$omp end do
!$omp end parallel

!!-- ここまでは if 文のインデントを調節しない.
!        end if
!
!        tmp=tmp/ac(i,j)
!-- 誤差の計算 ---
     do j=2,ny-1
        do i=2,nx-1
           err=abs(tmp(i,j)-psi(i,j))

!-- 最大誤差の更新
           if(err_max<=err)then
              err_max=err
           end if
        end do
     end do

!-- 一斉更新

     do j=2,ny-1
        do i=2,nx-1
           psi(i,j)=tmp(i,j)
        end do
     end do

  end do
  
!-- 境界の設定

  call calculate_bound( ib, dx, dy, bnd, psi )

!-- 未定義領域には undef を代入する.

  do j=1,ny
     do i=1,nx
        if(ib(i,j)==10)then
           psi(i,j)=defun
        end if
     end do
  end do

end subroutine Poisson_Jacobi
Subroutine :
ib(:,:) :integer, intent(in)
: 境界整数判別
dx(size(ib,1)) :real, intent(in)
: x 方向の格子解像度
dy(size(ib,2)) :real, intent(in)
: y 方向の格子解像度
bnd(size(ib,1),size(ib,2)) :real, intent(in)
: 境界値(ノイマン型のみ使用)
psi(size(ib,1),size(ib,2)) :real, intent(inout)
: 応答関数

ib 判別を元に, ノイマン型, 周期境界型について境界値を計算する.

[Source]

subroutine calculate_bound( ib, dx, dy, bnd, psi )
! ib 判別を元に, ノイマン型, 周期境界型について境界値を計算する.
  integer, intent(in) :: ib(:,:)  ! 境界整数判別
  real, intent(in) :: dx(size(ib,1))  ! x 方向の格子解像度
  real, intent(in) :: dy(size(ib,2))  ! y 方向の格子解像度
  real, intent(in) :: bnd(size(ib,1),size(ib,2))  ! 境界値(ノイマン型のみ使用)
  real, intent(inout) :: psi(size(ib,1),size(ib,2))  ! 応答関数
  integer :: i, j, nx, ny

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

  do j=1,ny
     do i=1,nx
        select case (ib(i,j))
        case (1)
           psi(i,j)=bnd(i,j)
        case (2)  ! x 方向にフラックス一定, 上側が参照値
           psi(i,j)=psi(i+1,j)-bnd(i,j)*dx(i)
        case (-2)  ! x 方向にフラックス一定, 下側が参照値
           psi(i,j)=psi(i-1,j)+bnd(i,j)*dx(i)
        case (4)  ! y 方向にフラックス一定, 右側が参照値
           psi(i,j)=psi(i,j+1)-bnd(i,j)*dy(j)
        case (-4)  ! y 方向にフラックス一定, 左側が参照値
           psi(i,j)=psi(i,j-1)+bnd(i,j)*dy(j)
        case (3)  ! 周期境界
           if(j==1)then
              psi(i,j)=psi(i,ny-1)
           else if(j==ny)then
              psi(i,j)=psi(i,2)
           else if(i==1)then
              psi(i,j)=psi(nx-1,j)
           else if(i==nx)then
              psi(i,j)=psi(2,j)
           end if
        case (-3)  ! 4 隅で同時周期境界
           if(i==1.and.j==1)then
              psi(i,j)=psi(nx-1,ny-1)
           else if(i==1.and.j==ny)then
              psi(i,j)=psi(nx-1,2)
           else if(i==nx.and.j==1)then
              psi(i,j)=psi(2,ny-1)
           else if(i==nx.and.j==ny)then
              psi(i,j)=psi(2,2)
           end if

        case (8)  ! 両方フラックス一定で左下角か右上角
           if(i==1.and.j==1)then  ! -- 評価 1
              psi(i,j)=psi(i+1,j+1)-0.5*(bnd(i+1,j)*dy(j)+bnd(i,j+1)*dx(i))

           else if(i==nx.and.j==ny)then  ! -- 評価 2
              psi(i,j)=psi(i-1,j-1)+0.5*(bnd(i-1,j)*dy(j)+bnd(i,j-1)*dx(i))

           else if(ib(i-1,j)==10.or.ib(i,j-1)==10)then
              ! -- 評価 1 と同じ
              psi(i,j)=psi(i+1,j+1)-0.5*(bnd(i+1,j)*dy(j)+bnd(i,j+1)*dx(i))

           else if(ib(i+1,j)==10.or.ib(i,j+1)==10)then
              ! -- 評価 2 と同じ
              psi(i,j)=psi(i-1,j-1)+0.5*(bnd(i-1,j)*dy(j)+bnd(i,j-1)*dx(i))

           end if

        case (-8)  ! 両方フラックス一定で右下角か左上角
           if(i==1.and.j==ny)then  ! -- 評価 1
              psi(i,j)=psi(i+1,j+1)+0.5*(bnd(i+1,j)*dy(j)-bnd(i,j-1)*dx(i))

           else if(i==nx.and.j==1)then  ! -- 評価 2
              psi(i,j)=psi(i-1,j-1)+0.5*(-bnd(i-1,j)*dy(j)+bnd(i,j+1)*dx(i))

           else if(ib(i-1,j)==10.or.ib(i,j+1)==10)then
              ! -- 評価 1 と同じ
              psi(i,j)=psi(i+1,j+1)+0.5*(bnd(i+1,j)*dy(j)-bnd(i,j-1)*dx(i))

           else if(ib(i+1,j)==10.or.ib(i,j-1)==10)then
              ! -- 評価 2 と同じ
              psi(i,j)=psi(i-1,j-1)+0.5*(-bnd(i-1,j)*dy(j)+bnd(i,j+1)*dx(i))

           end if
        end select 
     end do
  end do

end subroutine calculate_bound
Subroutine :
bound :character(4), intent(in)
: 領域境界のフラグ
ib(:,:) :integer, intent(inout)
: 全境界の判別整数
inner_flag(size(ib,1),size(ib,2)) :logical, intent(inout)
: 全領域境界フラグ
inner_bound(size(ib,1),size(ib,2)) :integer, intent(in), optional
: 内部領域境界の判別整数

各反復法ルーチンにおいて設定される境界条件のフラグをチェック, 設定する.

[Source]

subroutine set_bound( bound, ib, inner_flag, inner_bound )
! 各反復法ルーチンにおいて設定される境界条件のフラグをチェック, 設定する.
  implicit none
  character(4), intent(in) :: bound   ! 領域境界のフラグ
  integer, intent(inout) :: ib(:,:)   ! 全境界の判別整数
  logical, intent(inout) :: inner_flag(size(ib,1),size(ib,2))
                            ! 全領域境界フラグ
  integer, intent(in), optional :: inner_bound(size(ib,1),size(ib,2))
                            ! 内部領域境界の判別整数
  integer :: i, j, nx, ny

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

!-- 出力変数の初期化

  ib=0
  inner_flag=.false.

!-- 周期境界の設定確認.
!-- 周期境界なので, 両端とも 3 が設定されていないといけない.
  if(bound(1:1)=='3')then
     if(bound(3:3)/='3')then
        write(*,*) "### ERROR ###"
        write(*,*) "if bound = 3, bound(1:1)==bound(3:3). STOP."
        stop
     end if
  end if

  if(bound(3:3)=='3')then
     if(bound(1:1)/='3')then
        write(*,*) "### ERROR ###"
        write(*,*) "if bound = 3, bound(1:1)==bound(3:3). STOP."
        stop
     end if
  end if

  if(bound(2:2)=='3')then
     if(bound(4:4)/='3')then
        write(*,*) "### ERROR ###"
        write(*,*) "if bound = 3, bound(2:2)==bound(4:4). STOP."
        stop
     end if
  end if

  if(bound(4:4)=='3')then
     if(bound(2:2)/='3')then
        write(*,*) "### ERROR ###"
        write(*,*) "if bound = 3, bound(2:2)==bound(4:4). STOP."
        stop
     end if
  end if

  select case (bound(1:1))
  case ('1')
     do i=2,nx-1
        ib(i,1)=1
     end do

  case ('2')
     do i=2,nx-1
        ib(i,1)=4  ! y 方向のフラックスで参照値は上側
     end do

  case ('3')
     do i=2,nx-1
        ib(i,1)=3  ! y 方向 (x 軸) 周期境界
     end do
  end select

  select case (bound(2:2))
  case ('1')
     do j=2,ny-1
        ib(1,j)=1
     end do

  case ('2')
     do j=2,ny-1
        ib(1,j)=2  ! x 方向のフラックスで参照値は右側
     end do

  case ('3')
     do j=2,ny-1
        ib(1,j)=-3  ! x 方向 (y 軸) 周期境界
     end do
  end select

  select case (bound(3:3))
  case ('1')
     do i=2,nx-1
        ib(i,ny)=1
     end do

  case ('2')
     do i=2,nx-1
        ib(i,ny)=-4  ! y 方向のフラックスで参照値は下側
     end do

  case ('3')
     do i=2,nx-1
        ib(i,ny)=3  ! y 方向 (x 軸) 周期境界
     end do
  end select

  select case (bound(4:4))
  case ('1')
     do j=2,ny-1
        ib(nx,j)=1
     end do

  case ('2')
     do j=2,ny-1
        ib(nx,j)=-2  ! x 方向のフラックスで参照値は左側
     end do

  case ('3')
     do j=2,ny-1
        ib(nx,j)=-3  ! x 方向 (y 軸) 周期境界
     end do
  end select

!-- 領域隅境界は 2 辺が重なるので, 境界での強制方法は以下の順に優先する.
!-- (1) どちらかが 1 のとき, 隅領域は 1,
!-- (2) (1) 以外でどちらかが 3 のとき, 隅領域は 3.
!-- (3) (1), (2) 以外で 2. つまり, 両軸とも 2 なら 2 となる.

  ib(1,1)=0
  ib(1,ny)=0
  ib(nx,1)=0
  ib(nx,ny)=0

  if(bound(1:1)=='1')then
     ib(1,1)=1
     ib(nx,1)=1
  end if
  if(bound(2:2)=='1')then
     ib(1,1)=1
     ib(1,ny)=1
  end if
  if(bound(3:3)=='1')then
     ib(1,ny)=1
     ib(nx,ny)=1
  end if
  if(bound(4:4)=='1')then
     ib(nx,1)=1
     ib(nx,ny)=1
  end if

!-- 4 隅とも周期境界の場合
  if(bound(1:2)=='33')then
     ib(1,1)=-3
     ib(1,ny)=-3
     ib(nx,1)=-3
     ib(nx,ny)=-3
  end if

  if(bound(1:1)=='3')then
     if(ib(1,1)==0)then
        ib(1,1)=3
     end if
     if(ib(nx,1)==0)then
        ib(nx,1)=3
     end if
  end if
  if(bound(2:2)=='3')then
     if(ib(1,1)==0)then
        ib(1,1)=3
     end if
     if(ib(1,ny)==0)then
        ib(1,ny)=3
     end if
  end if
  if(bound(3:3)=='3')then
     if(ib(1,ny)==0)then
        ib(1,ny)=3
     end if
     if(ib(nx,ny)==0)then
        ib(nx,ny)=3
     end if
  end if
  if(bound(4:4)=='3')then
     if(ib(nx,1)==0)then
        ib(nx,1)=3
     end if
     if(ib(nx,ny)==0)then
        ib(nx,ny)=3
     end if
  end if

!-- 領域の隅であることを 8, -8 で示す.
!-- 8 の場合, 左下か右上であることを, -8 の場合, 左上か右下であることを示す.

  if(ib(1,1)==0)then
     ib(1,1)=8
  end if
  if(ib(nx,1)==0)then
     ib(nx,1)=-8
  end if
  if(ib(1,ny)==0)then
     ib(1,ny)=-8
  end if
  if(ib(nx,ny)==0)then
     ib(nx,ny)=8
  end if

!-- 内部境界の設定

  if(present(inner_bound))then
     do j=2,ny-1
        do i=2,nx-1
           ib(i,j)=inner_bound(i,j)
           if(inner_bound(i,j)/=0)then
              inner_flag(i,j)=.true.
           end if
        end do
     end do
  end if

  do i=1,nx  ! inner_bound が設定されているに関わらず, 領域端は計算しない.
     inner_flag(i,1)=.true.
     inner_flag(i,ny)=.true.
  end do

  do j=1,ny
     inner_flag(1,j)=.true.
     inner_flag(nx,j)=.true.
  end do

end subroutine set_bound
Subroutine :
coe(:,:) :real, intent(inout)
: 代入される配列
ext(size(coe,1),size(coe,2)) :real, intent(in), optional
: 代入する配列
def :real, intent(in), optional
: 代入する一定値

2 次元配列に ext で指定された値もしくは def で指定された一定値を代入する. ext, def どちらも optional であるが, 必ずどちらかは指定されていないといけない.

[Source]

subroutine set_coe( coe, ext, def )
! 2 次元配列に ext で指定された値もしくは def で指定された一定値を代入する.
! ext, def どちらも optional であるが, 必ずどちらかは指定されていないといけない.
  implicit none
  real, intent(inout) :: coe(:,:)  ! 代入される配列
  real, intent(in), optional :: ext(size(coe,1),size(coe,2))  ! 代入する配列
  real, intent(in), optional :: def  ! 代入する一定値
  integer :: i, j, nx, ny

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

  if(present(ext))then
     do j=1,ny
        do i=1,nx
           coe(i,j)=ext(i,j)
        end do
     end do
  else if(present(def))then
     do j=1,ny
        do i=1,nx
           coe(i,j)=def
        end do
     end do
  else
     write(*,*) "### ERROR ###"
     write(*,*) "subroutine set_coe must be set optional argument 'ext' or 'def'"
     write(*,*) "STOP."
     stop
  end if

end subroutine set_coe
Subroutine :
ib(:,:) :integer, intent(in)
: 境界条件判別整数
bnd(size(ib,1),size(ib,2)) :real, intent(inout)
: 境界での値
psi(size(ib,1),size(ib,2)) :real, intent(inout)
: 応答
bound_opt(size(ib,1),size(ib,2)) :real, intent(in), optional
: 境界での値

境界値を境界条件判別整数をもとに設定する.

[Source]

subroutine setval_bound( ib, bnd, psi, bound_opt )
! 境界値を境界条件判別整数をもとに設定する.
  implicit none
  integer, intent(in) :: ib(:,:)  ! 境界条件判別整数
  real, intent(inout) :: bnd(size(ib,1),size(ib,2))  ! 境界での値
  real, intent(inout) :: psi(size(ib,1),size(ib,2))  ! 応答
  real, intent(in), optional :: bound_opt(size(ib,1),size(ib,2))  ! 境界での値
  integer :: i, j, nx, ny

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

!-- 境界値の設定
  if(present(bound_opt))then
     do j=1,ny
        do i=1,nx
           bnd(i,j)=bound_opt(i,j)
        end do
     end do
  else
     do j=1,ny
        do i=1,nx
           bnd(i,j)=0.0
        end do
     end do
  end if

!-- 境界条件の代入 "ib(i,j)==1 の場合のみ"
!-- 内部領域についてもここで代入してしまう.

  do j=1,ny
     do i=1,nx
        if(ib(i,j)==1)then
           psi(i,j)=bnd(i,j)
        end if
     end do
  end do

end subroutine setval_bound