Class Derivation
In: derivation.f90

微分演算計算モジュール

Methods

Public Instance methods

Subroutine :
signal :character(1)
: デカルト座標系の何番目の乱流フラックス成分かを判定する.
1
= デカルト座標右手系における x 座標成分
2
= デカルト座標右手系における y 座標成分
3
= デカルト座標右手系における z 座標成分
x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
z(:) :real, intent(in)
: z 方向の空間座標 [m]
phi(size(x),size(y),size(z)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
rho(size(x),size(y),size(z)) :real, intent(in)
: 基本場の密度 [kg/m^3]
nuh(size(x),size(y),size(z)) :real, intent(in)
: 水平渦拡散係数
nuv(size(x),size(y),size(z)) :real, intent(in)
: 鉛直渦拡散係数
val(size(x),size(y),size(z)) :real, intent(inout)
: 計算されたテンソル成分 現在, 以下のオプションは使用していない.
undef :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子
sfcphi(size(x),size(y)) :real, intent(in), optional
: モデル最下層でのフラックスが与えられていれば, その値を代入. この値は何もせず, 単に val の最下層に代入されるだけ.

スカラー量の乱流フラックスを計算する.

[Source]

subroutine Reynolds_scal( signal, x, y, z, phi, rho, nuh, nuv, val, undef, hx, hy, hz, sfcphi )
! スカラー量の乱流フラックスを計算する.
  implicit none
  character(1) :: signal  ! デカルト座標系の何番目の乱流フラックス成分かを判定する.
                          ! [1] = デカルト座標右手系における x 座標成分
                          ! [2] = デカルト座標右手系における y 座標成分
                          ! [3] = デカルト座標右手系における z 座標成分
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: z(:)  ! z 方向の空間座標 [m]
  real, intent(in) :: phi(size(x),size(y),size(z))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: rho(size(x),size(y),size(z))  ! 基本場の密度 [kg/m^3]
  real, intent(in) :: nuh(size(x),size(y),size(z))  ! 水平渦拡散係数
  real, intent(in) :: nuv(size(x),size(y),size(z))  ! 鉛直渦拡散係数
  real, intent(inout) :: val(size(x),size(y),size(z))  ! 計算されたテンソル成分
! 現在, 以下のオプションは使用していない.
  real, intent(in), optional :: undef
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  real, intent(in), optional :: sfcphi(size(x),size(y))  ! モデル最下層でのフラックスが与えられていれば, その値を代入.
        ! この値は何もせず, 単に val の最下層に代入されるだけ.
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real, dimension(size(x),size(y)) :: stau
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez
  logical, dimension(size(x),size(y),size(z)) :: undeflag

  nx=size(x)
  ny=size(y)
  nz=size(z)
  undeflag=.false.

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

  val=0.0

  if(present(sfcphi))then
     if(signal(1:1)=='3')then
        stau(:,:)=sfcphi(:,:)
     end if
  end if

  select case (signal(1:1))
  case ('1')
     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, phi(:,j,k), val(:,j,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if(val(i,j,k)==undef.or.nuh(i,j,k)==undef.or.nuv(i,j,k)==undef)then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, phi(:,j,k), val(:,j,k) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=rho(i,j,k)*nuh(i,j,k)*val(i,j,k)/scalex(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  case ('2')
     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, phi(i,:,k), val(i,:,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if(val(i,j,k)==undef.or.nuh(i,j,k)==undef.or.nuv(i,j,k)==undef)then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, phi(i,:,k), val(i,:,k) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=rho(i,j,k)*nuh(i,j,k)*val(i,j,k)/scaley(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  case ('3')
     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, phi(i,j,:), val(i,j,:), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do j=1,ny
           do i=1,nx
              do k=1,nz
                 if(val(i,j,k)==undef.or.nuh(i,j,k)==undef.or.nuv(i,j,k)==undef)then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, phi(i,j,:), val(i,j,:) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=rho(i,j,k)*nuv(i,j,k)*val(i,j,k)/scalez(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

     do j=1,ny
        do i=1,nx
           if(present(sfcphi))then
              if(undeflag(i,j,1).eqv..false.)then
                 val(i,j,1)=stau(i,j)
              else
                 val(i,j,1)=undef
              end if
           end if
        end do
     end do

  end select

end subroutine
Subroutine :
signal :character(2)
: 計算するテンソル成分.
‘11’, ‘22’, ‘33‘
= それぞれ対角テンソル成分
‘12’, ‘13’, ‘21’, ‘23’, ‘31’, ‘32‘
= それぞれ非対角

テンソル成分. ただし, 対称テンソルであるため, ‘12’=’21’ を 計算していることに注意.

x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
z(:) :real, intent(in)
: z 方向の空間座標 [m]
u(size(x),size(y),size(z)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
v(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
w(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
rho(size(x),size(y),size(z)) :real, intent(in)
: 基本場の密度 [kg/m^3]
nuh(size(x),size(y),size(z)) :real, intent(in)
: 水平渦粘性係数
nuv(size(x),size(y),size(z)) :real, intent(in)
: 鉛直渦粘性係数
val(size(x),size(y),size(z)) :real, intent(inout)
: 計算されたテンソル成分 現在, 以下のオプションは使用していない.
undef :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子
sfctau(size(x),size(y)) :real, intent(in), optional
: モデル最下層での応力が与えられていれば, その値を代入. この値は何もせず, 単に val の最下層に代入されるだけ.

デカルト座標系におけるレイノルズ応力テンソルを計算する.

[Source]

subroutine Reynolds_stress( signal, x, y, z, u, v, w, rho, nuh, nuv, val, undef, hx, hy, hz, sfctau )
! デカルト座標系におけるレイノルズ応力テンソルを計算する.
  implicit none
  character(2) :: signal  ! 計算するテンソル成分.
                  ! ['11', '22', '33'] = それぞれ対角テンソル成分
                  ! ['12', '13', '21', '23', '31', '32'] = それぞれ非対角
                  ! テンソル成分. ただし, 対称テンソルであるため, '12'='21' を
                  ! 計算していることに注意.
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: z(:)  ! z 方向の空間座標 [m]
  real, intent(in) :: u(size(x),size(y),size(z))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: v(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: w(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: rho(size(x),size(y),size(z))  ! 基本場の密度 [kg/m^3]
  real, intent(in) :: nuh(size(x),size(y),size(z))  ! 水平渦粘性係数
  real, intent(in) :: nuv(size(x),size(y),size(z))  ! 鉛直渦粘性係数
  real, intent(inout) :: val(size(x),size(y),size(z))  ! 計算されたテンソル成分
! 現在, 以下のオプションは使用していない.
  real, intent(in), optional :: undef
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  real, intent(in), optional :: sfctau(size(x),size(y))  ! モデル最下層での応力が与えられていれば, その値を代入.
        ! この値は何もせず, 単に val の最下層に代入されるだけ.
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real :: stau(size(x),size(y))
  real :: sxx(size(x),size(y),size(z)), nu(size(x),size(y),size(z))
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez
  logical, dimension(size(x),size(y),size(z)) :: undeflag

  nx=size(x)
  ny=size(y)
  nz=size(z)
  undeflag=.false.

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

  val=0.0
  stau=0.0

!-- 地表面フラックスの値を代入.

  if(present(sfctau))then
     if(signal(2:2)=='3'.and.signal(1:1)/='3')then
        stau(:,:)=sfctau(:,:)
     end if
  end if

!-- [NOTE]
!-- 以下, 文字で case の or ができないため, 
!-- if 文の入れ子ではなく, if 文の並列表記で case と同じように見せかける.
!-- これはもちろん, 上から順に if をたどるが, どの場合も 2 種類以上の if に
!-- 合致しないことが既知であるために可能となる書き方であり,
!-- 並列表記した if の 2 パターン以上に合致してしまうような条件文では,
!-- case の代用には用いることができないことに注意.
!-- 本ライブラリでこのような紛らわしい表記をしている場合は必ず NOTE が入る.

!-- 応力テンソルの計算
  if(present(undef))then
     call deform_tensor( signal, x, y, z, u, v, w, sxx, hx=scalex, hy=scaley, hz=scalez, undef=undef )
  else
     call deform_tensor( signal, x, y, z, u, v, w, sxx, hx=scalex, hy=scaley, hz=scalez )
  end if

  if(signal(1:2)=='12'.or.signal(1:2)=='21')then

     do k=1,nz
        do j=1,ny
           do i=1,nx
              nu(i,j,k)=nuh(i,j,k)
           end do
        end do
     end do

  else if(signal(1:2)=='23'.or.signal(1:2)=='32')then

     if(signal(2:2)=='3')then
        do k=1,nz
           do j=1,ny
              do i=1,nx
                 nu(i,j,k)=nuv(i,j,k)
              end do
           end do
        end do
     else
        do k=1,nz
           do j=1,ny
              do i=1,nx
                 nu(i,j,k)=nuh(i,j,k)
              end do
           end do
        end do
     end if

  else if(signal(1:2)=='13'.or.signal(1:2)=='31')then

     if(signal(2:2)=='3')then
        do k=1,nz
           do j=1,ny
              do i=1,nx
                 nu(i,j,k)=nuv(i,j,k)
              end do
           end do
        end do
     else
        do k=1,nz
           do j=1,ny
              do i=1,nx
                 nu(i,j,k)=nuh(i,j,k)
              end do
           end do
        end do
     end if

  else if(signal(1:2)=='11')then

     if(present(undef))then
        call div_3d( x, y, z, u, v, w, val, hx=scalex, hy=scaley, hz=scalez, undeff=undef )
     else
        call div_3d( x, y, z, u, v, w, val, hx=scalex, hy=scaley, hz=scalez )
     end if

     do k=1,nz
        do j=1,ny
           do i=1,nx
              nu(i,j,k)=nuh(i,j,k)
           end do
        end do
     end do

  else if(signal(1:2)=='22')then

     if(present(undef))then
        call div_3d( x, y, z, u, v, w, val, hx=scalex, hy=scaley, hz=scalez, undeff=undef )
     else
        call div_3d( x, y, z, u, v, w, val, hx=scalex, hy=scaley, hz=scalez )
     end if

     do k=1,nz
        do j=1,ny
           do i=1,nx
              nu(i,j,k)=nuh(i,j,k)
           end do
        end do
     end do

  else if(signal(1:2)=='33')then

     if(present(undef))then
        call div_3d( x, y, z, u, v, w, val, hx=scalex, hy=scaley, hz=scalez, undeff=undef )
     else
        call div_3d( x, y, z, u, v, w, val, hx=scalex, hy=scaley, hz=scalez )
     end if

     do k=1,nz
        do j=1,ny
           do i=1,nx
              nu(i,j,k)=nuv(i,j,k)
           end do
        end do
     end do
  end if

  if(present(undef))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              if(val(i,j,k)==undef.or.nu(i,j,k)==undef.or.sxx(i,j,k)==undef)then
                 undeflag(i,j,k)=.true.
              end if
           end do
        end do
     end do
  end if

!-- 以下, 最下層は地表面フラックスを代入するかどうかのオプションのため, 別ループ

  if(present(sfctau))then
     do j=1,ny
        do i=1,nx
           if(undeflag(i,j,1).eqv..false.)then
              val(i,j,1)=stau(i,j)
           else
              val(i,j,1)=undef
           end if
        end do
     end do
  else
!$omp parallel default(shared)
!$omp do schedule(dynamic) private(i,j)
     do j=1,ny
        do i=1,nx
           if(undeflag(i,j,1).eqv..false.)then
              val(i,j,1)=rho(i,j,1)*nu(i,j,1)*(sxx(i,j,1)-(2.0/3.0)*val(i,j,1))
           else
              val(i,j,1)=undef
           end if
        end do
     end do
!$omp end do
!$omp end parallel
  end if

!$omp parallel default(shared)
!$omp do schedule(dynamic) private(i,j,k)
  do k=2,nz
     do j=1,ny
        do i=1,nx
           if(undeflag(i,j,k).eqv..false.)then
              val(i,j,k)=rho(i,j,k)*nu(i,j,k)*(sxx(i,j,k)-(2.0/3.0)*val(i,j,k))
           else
              val(i,j,k)=undef
           end if
        end do
     end do
  end do
!$omp end do
!$omp end parallel

end subroutine
Subroutine :
x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
u(size(x),size(y)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
v(size(x),size(y)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
val(size(x),size(y)) :real, intent(inout)
: 2次元渦度
undeff :real, intent(in), optional
hx(size(x),size(y)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y)) :real, intent(in), optional
: y 方向のスケール因子
ord :logical, intent(in), optional
: 微分の順番を入れ替えるオプション. true なら, 入れ替えない, false なら, 入れ替える. デフォルトは true, du/dz-dw/dx を計算するときのみ用いる.

2 次元渦度計算ルーチン

x, y は配列の次元の若い順に必ず並べること.

u, v は偶置換の向きに配置すれば, 任意の2次元渦度が計算可能 ただし, du/dz-dw/dx を計算するときのみ, (x,z,u,w) の順番で, ord オプション false.

$frac{partial v}{partial x} -frac{partial u}{partial y} $ を 2 次の中央差分近似で書き換えると, 点 $(i,j)$ での発散は $frac{v_{i,j+1}-v_{i,j-1}}{2dx} -frac{u_{i+1,j}-u_{i-1,j}}{2dy} $ とできる. これを用いて2次元発散を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる. 実質的には grad_1d が微分計算を担当するので, 境界の計算も自動で行う.

[Source]

subroutine curl( x, y, u, v, val, undeff, hx, hy, ord )
! 2 次元渦度計算ルーチン
!
! x, y は配列の次元の若い順に必ず並べること.
!
! u, v は偶置換の向きに配置すれば, 任意の2次元渦度が計算可能
! ただし, du/dz-dw/dx を計算するときのみ, (x,z,u,w) の順番で, ord オプション false.
!
! $\frac{\partial v}{\partial x} -\frac{\partial u}{\partial y} $ を
! 2 次の中央差分近似で書き換えると, 点 $(i,j)$ での発散は
! $\frac{v_{i,j+1}-v_{i,j-1}}{2dx} -\frac{u_{i+1,j}-u_{i-1,j}}{2dy} $
! とできる. これを用いて2次元発散を計算.
! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
! 落ちる.
! 実質的には grad_1d が微分計算を担当するので, 境界の計算も自動で行う.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: u(size(x),size(y))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: v(size(x),size(y))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(inout) :: val(size(x),size(y))  ! 2次元渦度
  real, intent(in), optional :: undeff
  real, intent(in), optional :: hx(size(x),size(y))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y))  ! y 方向のスケール因子
  logical, intent(in),  optional :: ord  ! 微分の順番を入れ替えるオプション.
                 ! true なら, 入れ替えない, false なら, 入れ替える.
                 ! デフォルトは true, du/dz-dw/dx を計算するときのみ用いる.
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  logical :: order
  real :: scalex(size(x),size(y)), scaley(size(x),size(y))
  real :: dvdx(size(x),size(y)), dudy(size(x),size(y))
  real :: tmpu(size(x),size(y)), tmpv(size(x),size(y))

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

!-- スケーリング変数の設定.

  if(present(hx).and.present(hy))then
     do j=1,ny
        do i=1,nx
           scalex(i,j)=hx(i,j)
           scaley(i,j)=hy(i,j)
        end do
     end do
  else
     do j=1,ny
        do i=1,nx
           scalex(i,j)=1.0
           scaley(i,j)=1.0
        end do
     end do
  end if

!-- 方針は, x 方向に dvdy を計算し, y 方向に dudx を計算する.

  if(present(undeff))then
     do j=1,ny
        do i=1,nx
           if(u(i,j)/=undeff.and.v(i,j)/=undeff)then
              tmpu(i,j)=scalex(i,j)*u(i,j)
              tmpv(i,j)=scaley(i,j)*v(i,j)
           else
              tmpu(i,j)=undeff
              tmpv(i,j)=undeff
           end if
        end do
     end do

     do i=1,nx
        call grad_1d( y, tmpu(i,:), dudy(i,:), undeff )
     end do
     do j=1,ny
        call grad_1d( x, tmpv(:,j), dvdx(:,j), undeff )
     end do

     do j=1,ny
        do i=1,nx
           if(dudy(i,j)/=undeff.and.dvdx(i,j)/=undeff)then
              if(scalex(i,j)/=0.0.and.scaley(i,j)/=0.0)then
                 val(i,j)=(dvdx(i,j)-dudy(i,j))/(scalex(i,j)*scaley(i,j))
              else
                 val(i,j)=0.0
              end if
           else
              val(i,j)=undeff
           end if
        end do
     end do

  else

     do j=1,ny
        do i=1,nx
           tmpu(i,j)=scalex(i,j)*u(i,j)
           tmpv(i,j)=scaley(i,j)*v(i,j)
        end do
     end do

     do i=1,nx
        call grad_1d( y, tmpu(i,:), dudy(i,:) )
     end do
     do j=1,ny
        call grad_1d( x, tmpv(:,j), dvdx(:,j) )
     end do

     do j=1,ny
        do i=1,nx
           if(scalex(i,j)/=0.0.and.scaley(i,j)/=0.0)then
              val(i,j)=(dvdx(i,j)-dudy(i,j))/(scalex(i,j)*scaley(i,j))
           else
              val(i,j)=0.0
           end if
        end do
     end do
  end if

!-- 回転の順番を逆にするオプション
!-- false なら, 順番を入れ替えて出力する.

  if(present(ord))then
     order=ord
  else
     order=.true.
  end if

  if(order.eqv..false.)then
     do j=1,ny
        do i=1,nx
           val(i,j)=-val(i,j)
        end do
     end do
  end if

end subroutine curl
Subroutine :
x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
z(:) :real, intent(in)
: z 方向の空間座標 [m]
u(size(x),size(y),size(z)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
v(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
w(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
zeta(size(x),size(y),size(z)) :real, intent(inout)
: 渦度ベクトル x 成分
eta(size(x),size(y),size(z)) :real, intent(inout)
: 渦度ベクトル y 成分
xi(size(x),size(y),size(z)) :real, intent(inout)
: 渦度ベクトル z 成分
undeff :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子

3 次元渦度を計算する. 引数の順番は右手系で x, y, z のデカルト座標系, それに対応するベクトル成分 u, v, w を代入すると, それに対応した渦度ベクトル成分 zeta, eta, xi が計算される. 実質は grad_1d が計算を担当するので, 境界の処理も自動で行う.

[Source]

subroutine curl_3d( x, y, z, u, v, w, zeta, eta, xi, undeff, hx, hy, hz )
! 3 次元渦度を計算する.
! 引数の順番は右手系で x, y, z のデカルト座標系,
! それに対応するベクトル成分 u, v, w を代入すると,
! それに対応した渦度ベクトル成分 zeta, eta, xi が計算される.
! 実質は grad_1d が計算を担当するので, 境界の処理も自動で行う.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: z(:)  ! z 方向の空間座標 [m]
  real, intent(in) :: u(size(x),size(y),size(z))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: v(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: w(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(inout) :: zeta(size(x),size(y),size(z))  ! 渦度ベクトル x 成分
  real, intent(inout) :: eta(size(x),size(y),size(z))  ! 渦度ベクトル y 成分
  real, intent(inout) :: xi(size(x),size(y),size(z))  ! 渦度ベクトル z 成分
  real, intent(in), optional :: undeff
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez

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

!-- スケーリング変数の設定.

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

  if(present(undeff))then

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i)
! x 軸に垂直な面上での x 方向の渦度ベクトルを 3 次元全域で計算.
     do i=1,nx
        call curl( y, z, v(i,:,:), w(i,:,:), zeta(i,:,:), undeff, hx=scaley(i,:,:), hy=scalez(i,:,:) )
     end do
!$omp end do
!$omp end parallel

  else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i)
! x 軸に垂直な面上での x 方向の渦度ベクトルを 3 次元全域で計算.
     do i=1,nx
        call curl( y, z, v(i,:,:), w(i,:,:), zeta(i,:,:), hx=scaley(i,:,:), hy=scalez(i,:,:) )
     end do
!$omp end do
!$omp end parallel

  end if

! y 軸に垂直な面上での y 方向の渦度ベクトルを 3 次元全域で計算.

  if(present(undeff))then

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j)
     do j=1,ny
        call curl( x, z, u(:,j,:), w(:,j,:), eta(:,j,:), undeff, ord=.false., hx=scalex(:,j,:), hy=scalez(:,j,:) )
     end do
!$omp end do
!$omp end parallel

  else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j)
     do j=1,ny
        call curl( x, z, u(:,j,:), w(:,j,:), eta(:,j,:), ord=.false., hx=scalex(:,j,:), hy=scalez(:,j,:) )
     end do
!$omp end do
!$omp end parallel

  end if

! z 軸に垂直な面上での z 方向の渦度ベクトルを 3 次元全域で計算.

  if(present(undeff))then

!$omp parallel default(shared)
!$omp do schedule(runtime) private(k)
     do k=1,nz
        call curl( x, y, u(:,:,k), v(:,:,k), xi(:,:,k), undeff, hx=scalex(:,:,k), hy=scaley(:,:,k) )
     end do
!$omp end do
!$omp end parallel

  else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(k)
     do k=1,nz
        call curl( x, y, u(:,:,k), v(:,:,k), xi(:,:,k), hx=scalex(:,:,k), hy=scaley(:,:,k) )
     end do
!$omp end do
!$omp end parallel

  end if

end subroutine
Subroutine :
signal :character(2)
: 計算するテンソル成分.
‘11’, ‘22’, ‘33‘
= それぞれ対角テンソル成分
‘12’, ‘13’, ‘21’, ‘23’, ‘31’, ‘32‘
= それぞれ非対角

テンソル成分. ただし, 対称テンソルであるため, ‘12’=’21’ を 計算していることに注意.

x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
z(:) :real, intent(in)
: z 方向の空間座標 [m]
u(size(x),size(y),size(z)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
v(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
w(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
val(size(x),size(y),size(z)) :real, intent(inout)
: 計算されたテンソル成分 現在, 以下のオプションは使用していない.
undef :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子

デカルト座標系における変形速度テンソルを計算する.

[Source]

subroutine deform_tensor( signal, x, y, z, u, v, w, val, undef, hx, hy, hz )
! デカルト座標系における変形速度テンソルを計算する.
  implicit none
  character(2) :: signal  ! 計算するテンソル成分.
                  ! ['11', '22', '33'] = それぞれ対角テンソル成分
                  ! ['12', '13', '21', '23', '31', '32'] = それぞれ非対角
                  ! テンソル成分. ただし, 対称テンソルであるため, '12'='21' を
                  ! 計算していることに注意.
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: z(:)  ! z 方向の空間座標 [m]
  real, intent(in) :: u(size(x),size(y),size(z))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: v(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: w(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(inout) :: val(size(x),size(y),size(z))  ! 計算されたテンソル成分
! 現在, 以下のオプションは使用していない.
  real, intent(in), optional :: undef
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real, dimension(size(x),size(y),size(z)) :: tmpu, tmpv, tmpw
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez
  real, dimension(size(x),size(y),size(z)) :: ddx, ddy, ddz
  logical, dimension(size(x),size(y),size(z)) :: undeflag

  nx=size(x)
  ny=size(y)
  nz=size(z)
  undeflag=.false.

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

!-- [NOTE]
!-- 以下, 文字で case の or ができないため, 
!-- if 文の入れ子ではなく, if 文の並列表記で case と同じように見せかける.
!-- これはもちろん, 上から順に if をたどるが, どの場合も 2 種類以上の if に
!-- 合致しないことが既知であるために可能となる書き方であり,
!-- 並列表記した if の 2 パターン以上に合致してしまうような条件文では,
!-- case の代用には用いることができないことに注意.
!-- 本ライブラリでこのような紛らわしい表記をしている場合は必ず NOTE が入る.

  if(present(undef))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              if(u(i,j,k)==undef.or.v(i,j,k)==undef.or.w(i,j,k)==undef)then
                 undeflag(i,j,k)=.true.
              end if
           end do
        end do
     end do
  end if

  if(signal(1:2)=='12'.or.signal(1:2)=='21')then

     do k=1,nz
        do j=1,ny
           do i=1,nx
              if(undeflag(i,j,k).eqv..false.)then
                 tmpu(i,j,k)=u(i,j,k)/scalex(i,j,k)
                 tmpv(i,j,k)=v(i,j,k)/scaley(i,j,k)
              else
                 tmpu(i,j,k)=undef
                 tmpv(i,j,k)=undef
              end if
           end do
        end do
     end do

     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, tmpv(:,j,k), ddx(:,j,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, tmpu(i,:,k), ddy(i,:,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if((ddx(i,j,k)==undef).or.(ddy(i,j,k)==undef))then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, tmpv(:,j,k), ddx(:,j,k) )
           end do
        end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, tmpu(i,:,k), ddy(i,:,k) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=ddx(i,j,k)*scaley(i,j,k)/scalex(i,j,k)+ ddy(i,j,k)*scalex(i,j,k)/scaley(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  else if(signal(1:2)=='23'.or.signal(1:2)=='32')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(undeflag(i,j,k).eqv..false.)then
                 tmpv(i,j,k)=v(i,j,k)/scaley(i,j,k)
                 tmpw(i,j,k)=w(i,j,k)/scalez(i,j,k)
              else
                 tmpv(i,j,k)=undef
                 tmpw(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, tmpw(i,:,k), ddy(i,:,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, tmpv(i,j,:), ddz(i,j,:), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if((ddy(i,j,k)==undef).or.(ddz(i,j,k)==undef))then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, tmpw(i,:,k), ddy(i,:,k) )
           end do
        end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, tmpv(i,j,:), ddz(i,j,:) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=ddy(i,j,k)*scalez(i,j,k)/scaley(i,j,k)+ ddz(i,j,k)*scaley(i,j,k)/scalez(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  else if(signal(1:2)=='13'.or.signal(1:2)=='31')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(undeflag(i,j,k).eqv..false.)then
                 tmpu(i,j,k)=u(i,j,k)/scalex(i,j,k)
                 tmpw(i,j,k)=w(i,j,k)/scalez(i,j,k)
              else
                 tmpu(i,j,k)=undef
                 tmpw(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, tmpu(i,j,:), ddz(i,j,:), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, tmpw(:,j,k), ddx(:,j,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if((ddz(i,j,k)==undef).or.(ddx(i,j,k)==undef))then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, tmpu(i,j,:), ddz(i,j,:) )
           end do
        end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, tmpw(:,j,k), ddx(:,j,k) )
           end do
        end do
!$omp end do
!$omp end parallel
     end if

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=ddz(i,j,k)*scalex(i,j,k)/scalez(i,j,k)+ ddx(i,j,k)*scalez(i,j,k)/scalex(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  else if(signal(1:2)=='11')then

     !-- scale は undef 指定していないので, undef 計算はしない.
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
     do k=1,nz
        do i=1,nx
           call grad_1d( y, scalex(i,:,k), ddy(i,:,k) )
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do j=1,ny
        do i=1,nx
           call grad_1d( z, scalex(i,j,:), ddz(i,j,:) )
        end do
     end do
!$omp end do
!$omp end parallel

     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, u(:,j,k), ddx(:,j,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if(ddx(i,j,k)==undef)then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
        do k=1,nz
           do j=1,ny
              call grad_1d( x, u(:,j,k), ddx(:,j,k) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=2.0*(ddx(i,j,k)+ddy(i,j,k)*v(i,j,k)/scaley(i,j,k)+ ddz(i,j,k)*w(i,j,k)/scalez(i,j,k))/scalex(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  else if(signal(1:2)=='22')then

     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, v(i,:,k), ddy(i,:,k), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if(ddy(i,j,k)==undef)then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
        do k=1,nz
           do i=1,nx
              call grad_1d( y, v(i,:,k), ddy(i,:,k) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

     !-- scale は undef 指定していないので, undef 計算はしない.
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do j=1,ny
        do i=1,nx
           call grad_1d( z, scaley(i,j,:), ddz(i,j,:) )
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
     do k=1,nz
        do j=1,ny
           call grad_1d( x, scaley(:,j,k), ddx(:,j,k) )
        end do
     end do
!$omp end do
!$omp end parallel

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=2.0*(ddy(i,j,k)+ddz(i,j,k)*w(i,j,k)/scalez(i,j,k)+ ddx(i,j,k)*u(i,j,k)/scalex(i,j,k))/scaley(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  else if(signal(1:2)=='33')then

     !-- scale は undef 指定していないので, undef 計算はしない.
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
     do k=1,nz
        do i=1,nx
           call grad_1d( y, scalez(i,:,k), ddy(i,:,k) )
        end do
     end do
!$omp end do
!$omp end parallel

     if(present(undef))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, w(i,j,:), ddz(i,j,:), undef=undef )
           end do
        end do
!$omp end do
!$omp end parallel

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 if(ddz(i,j,k)==undef)then
                    undeflag(i,j,k)=.true.
                 end if
              end do
           end do
        end do

     else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
        do j=1,ny
           do i=1,nx
              call grad_1d( z, w(i,j,:), ddz(i,j,:) )
           end do
        end do
!$omp end do
!$omp end parallel

     end if

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
     do k=1,nz
        do j=1,ny
           call grad_1d( x, scalez(:,j,k), ddx(:,j,k) )
        end do
     end do
!$omp end do
!$omp end parallel

!$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(undeflag(i,j,k).eqv..false.)then
                 val(i,j,k)=2.0*(ddz(i,j,k)+ddx(i,j,k)*u(i,j,k)/scalex(i,j,k)+ ddy(i,j,k)*v(i,j,k)/scaley(i,j,k))/scalez(i,j,k)
              else
                 val(i,j,k)=undef
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  end if

end subroutine
Subroutine :
x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
u(size(x),size(y)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
v(size(x),size(y)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
val(size(x),size(y)) :real, intent(inout)
: 2次元発散値
undeff :real, intent(in), optional
hx(size(x),size(y)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y)) :real, intent(in), optional
: y 方向のスケール因子

2次元発散計算ルーチン

$frac{partial u}{partial x} +frac{partial v}{partial y} $ を 2 次の中央差分近似で書き換えると, 点 $(i,j)$ での発散は $frac{u_{i+1,j}-u_{i-1,j}}{2dx} + frac{v_{i,j+1}-v_{i,j-1}}{2dy} $ とできる. これを用いて2次元発散を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる. 実質的には, grad_1d ルーチンを組み合わせることで計算を行う. du/dx = grad_1d( x, u, dudx ), dv/dy = grad_1d( y, v, dvdy ) という形で計算を行えば, 境界も自動的に計算可能.

[Source]

subroutine div( x, y, u, v, val, undeff, hx, hy )
! 2次元発散計算ルーチン
! 
! $\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y} $ を
! 2 次の中央差分近似で書き換えると, 点 $(i,j)$ での発散は
! $\frac{u_{i+1,j}-u_{i-1,j}}{2dx} + \frac{v_{i,j+1}-v_{i,j-1}}{2dy} $
! とできる. これを用いて2次元発散を計算.
! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
! 落ちる.
! 実質的には, grad_1d ルーチンを組み合わせることで計算を行う.
! du/dx = grad_1d( x, u, dudx ), dv/dy = grad_1d( y, v, dvdy ) という形で計算を行えば,
! 境界も自動的に計算可能.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: u(size(x),size(y))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: v(size(x),size(y))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(inout) :: val(size(x),size(y))  ! 2次元発散値
  real, intent(in), optional :: undeff
  real, intent(in), optional :: hx(size(x),size(y))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y))  ! y 方向のスケール因子
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  real :: scalex(size(x),size(y)), scaley(size(x),size(y))
  real :: dudx(size(x),size(y)), dvdy(size(x),size(y))
  real :: tmpu(size(x),size(y)), tmpv(size(x),size(y))

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

!-- スケーリング変数の設定.

  if(present(hx).and.present(hy))then
     do j=1,ny
        do i=1,nx
           scalex(i,j)=hx(i,j)
           scaley(i,j)=hy(i,j)
        end do
     end do
  else
     do j=1,ny
        do i=1,nx
           scalex(i,j)=1.0
           scaley(i,j)=1.0
        end do
     end do
  end if

!-- 方針は, x 方向に dvdy を計算し, y 方向に dudx を計算する.

  if(present(undeff))then
     do j=1,ny
        do i=1,nx
           if(u(i,j)/=undeff.and.v(i,j)/=undeff)then
              tmpu(i,j)=scaley(i,j)*u(i,j)
              tmpv(i,j)=scalex(i,j)*v(i,j)
           else
              tmpu(i,j)=undeff
              tmpv(i,j)=undeff
           end if
        end do
     end do

     do i=1,nx
        call grad_1d( y, tmpv(i,:), dvdy(i,:), undeff )
     end do
     do j=1,ny
        call grad_1d( x, tmpu(:,j), dudx(:,j), undeff )
     end do

     do j=1,ny
        do i=1,nx
           if(dudx(i,j)/=undeff.and.dvdy(i,j)/=undeff)then
              if(scalex(i,j)/=0.0.and.scaley(i,j)/=0.0)then
                 val(i,j)=(dudx(i,j)+dvdy(i,j))/(scalex(i,j)*scaley(i,j))
              else
                 val(i,j)=0.0
              end if
           else
              val(i,j)=undeff
           end if
        end do
     end do

  else

     do j=1,ny
        do i=1,nx
           tmpu(i,j)=scaley(i,j)*u(i,j)
           tmpv(i,j)=scalex(i,j)*v(i,j)
        end do
     end do

     do i=1,nx
        call grad_1d( y, tmpv(i,:), dvdy(i,:) )
     end do
     do j=1,ny
        call grad_1d( x, tmpu(:,j), dudx(:,j) )
     end do

     do j=1,ny
        do i=1,nx
           if(scalex(i,j)/=0.0.and.scaley(i,j)/=0.0)then
              val(i,j)=(dudx(i,j)+dvdy(i,j))/(scalex(i,j)*scaley(i,j))
           else
              val(i,j)=0.0
           end if
        end do
     end do
  end if

end subroutine div
Subroutine :
x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
z(:) :real, intent(in)
: z 方向の空間座標 [m]
u(size(x),size(y),size(z)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
v(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
w(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
val(size(x),size(y),size(z)) :real, intent(inout)
: 3 次元発散値
undeff :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子

3次元発散計算ルーチン

$frac{partial u}{partial x} +frac{partial v}{partial y} +frac{partial w}{partial z} $ を 2 次の中央差分近似で書き換えると, 点 $(i,j,k)$ での発散は $frac{u_{i+1,j,k}-u_{i-1,j,k}}{2dx} + frac{v_{i,j+1,k}-v_{i,j-1,k}}{2dy} + frac{w_{i,j,k+1}-w_{i,j,k-1}}{2dz} $ とできる. これを用いて 3 次元発散を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる. 実質的には, grad_1d ルーチンが計算を行うので, 境界の処理も自動で行う.

[Source]

subroutine div_3d( x, y, z, u, v, w, val, undeff, hx, hy, hz )
! 3次元発散計算ルーチン
!
! $\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y} +\frac{\partial w}{\partial z} $ を
! 2 次の中央差分近似で書き換えると, 点 $(i,j,k)$ での発散は
! $\frac{u_{i+1,j,k}-u_{i-1,j,k}}{2dx} + \frac{v_{i,j+1,k}-v_{i,j-1,k}}{2dy} + \frac{w_{i,j,k+1}-w_{i,j,k-1}}{2dz} $
! とできる. これを用いて 3 次元発散を計算.
! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
! 落ちる.
! 実質的には, grad_1d ルーチンが計算を行うので, 境界の処理も自動で行う.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: z(:)  ! z 方向の空間座標 [m]
  real, intent(in) :: u(size(x),size(y),size(z))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: v(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: w(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(inout) :: val(size(x),size(y),size(z))  ! 3 次元発散値
  real, intent(in), optional :: undeff
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez
  real, dimension(size(x),size(y),size(z)) :: dudx, dvdy, dwdz
  real, dimension(size(x),size(y),size(z)) :: tmpu, tmpv, tmpw

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

!-- スケーリング変数の設定.

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

!-- 方針は, x 方向に dvdy を計算し, y 方向に dudx を計算する.

  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(u(i,j,k)/=undeff.and.v(i,j,k)/=undeff.and.w(i,j,k)/=undeff)then
                 tmpu(i,j,k)=scaley(i,j,k)*scalez(i,j,k)*u(i,j,k)
                 tmpv(i,j,k)=scalez(i,j,k)*scalex(i,j,k)*v(i,j,k)
                 tmpw(i,j,k)=scalex(i,j,k)*scaley(i,j,k)*w(i,j,k)
              else
                 tmpu(i,j,k)=undeff
                 tmpv(i,j,k)=undeff
                 tmpw(i,j,k)=undeff
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
     do k=1,nz
        do i=1,nx
           call grad_1d( y, tmpv(i,:,k), dvdy(i,:,k), undeff )
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
     do k=1,nz
        do j=1,ny
           call grad_1d( x, tmpu(:,j,k), dudx(:,j,k), undeff )
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do j=1,ny
        do i=1,nx
           call grad_1d( z, tmpw(i,j,:), dwdz(i,j,:), undeff )
        end do
     end do
!$omp end do
!$omp end parallel

!$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(dudx(i,j,k)/=undeff.and.dvdy(i,j,k)/=undeff.and. dwdz(i,j,k)/=undeff)then
                 if(scalex(i,j,k)/=0.0.and.scaley(i,j,k)/=0.0.and. scalez(i,j,k)/=0.0)then
                    val(i,j,k)=(dudx(i,j,k)+dvdy(i,j,k)+dwdz(i,j,k))/ (scalex(i,j,k)*scaley(i,j,k)*scalez(i,j,k))
                 else
                    val(i,j,k)=0.0
                 end if
              else
                 val(i,j,k)=undeff
              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
              tmpu(i,j,k)=scaley(i,j,k)*scalez(i,j,k)*u(i,j,k)
              tmpv(i,j,k)=scalez(i,j,k)*scalex(i,j,k)*v(i,j,k)
              tmpw(i,j,k)=scalex(i,j,k)*scaley(i,j,k)*w(i,j,k)
           end do
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,k)
     do k=1,nz
        do i=1,nx
           call grad_1d( y, tmpv(i,:,k), dvdy(i,:,k) )
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(j,k)
     do k=1,nz
        do j=1,ny
           call grad_1d( x, tmpu(:,j,k), dudx(:,j,k) )
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do j=1,ny
        do i=1,nx
           call grad_1d( z, tmpw(i,j,:), dwdz(i,j,:) )
        end do
     end do
!$omp end do
!$omp end parallel

!$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(scalex(i,j,k)/=0.0.and.scaley(i,j,k)/=0.0.and.scalez(i,j,k)/=0.0)then
                 val(i,j,k)=(dudx(i,j,k)+dvdy(i,j,k)+dwdz(i,j,k))/ (scalex(i,j,k)*scaley(i,j,k)*scalez(i,j,k))
              else
                 val(i,j,k)=0.0
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  end if

end subroutine div_3d
Subroutine :
x(:) :real, intent(in)
: 空間座標
u(size(x)) :real, intent(in)
: 上の空間配列に対応する 1 次元スカラー値
val(size(x)) :real, intent(inout)
: スカラー値の x 方向の勾配
undef :real, intent(in), optional
hx(size(x)) :real, intent(in), optional
: x のスケール因子

1 次元のスカラー変数の勾配を計算する $frac{partial p}{partial x} $ を 4 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は $(2/3)*frac{p_{i+1}-p_{i-1}}{dx} -(p_{i+2}-p_{i-2})/(12dx)$ とできる. これを用いて 1 次元勾配を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる.

[Source]

subroutine grad4_1d( x, u, val, undef, hx )
! 1 次元のスカラー変数の勾配を計算する
! $\frac{\partial p}{\partial x} $ を
! 4 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は
! $(2/3)*\frac{p_{i+1}-p_{i-1}}{dx} -(p_{i+2}-p_{i-2})/(12dx)$
! とできる. これを用いて 1 次元勾配を計算.
! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
! 落ちる.
  implicit none
  real, intent(in) :: x(:)  ! 空間座標
  real, intent(in) :: u(size(x))  ! 上の空間配列に対応する 1 次元スカラー値
  real, intent(inout) :: val(size(x))  ! スカラー値の x 方向の勾配
  real, intent(in), optional :: undef
  real, intent(in), optional :: hx(size(x))  ! x のスケール因子
  integer :: i  ! イタレーション用添字
  integer :: nx  ! 配列要素数
  real :: scalex(size(x)), dx(size(x)), coe23, coe112

  nx=size(x)
  do i=2,nx-1
     dx(i)=0.5*(x(i+1)-x(i-1))
  end do
  dx(1)=x(2)-x(1)
  dx(nx)=x(nx)-x(nx-1)
  coe23=2.0/3.0
  coe112=1.0/12.0

  if(present(hx))then
     do i=1,nx
        scalex(i)=hx(i)
     end do
  else
     do i=1,nx
        scalex(i)=1.0
     end do
  end if

  if(present(undef))then
     do i=3,nx-2
        if(u(i+1)==undef.or.u(i-1)==undef.or. u(i+2)==undef.or.u(i-2)==undef)then
           val(i)=undef
        else
           if(scalex(i)/=0.0)then
              val(i)=(coe23*(u(i+1)-u(i-1))-coe112*(u(i+2)-u(i-2))) /(scalex(i)*dx(i))
           else
              val(i)=0.0
           end if
        end if
     end do

!-- データ数のない両端の処理 (両端の 1 つ内側) ---
     if(u(1)==undef.or.u(3)==undef)then
        val(2)=undef
     else
        if(scalex(2)/=0.0)then
           val(2)=0.5*(u(3)-u(1))/(scalex(2)*dx(2))
        else
           val(2)=0.0
        end if
     end if

     if(u(nx)==undef.or.u(nx-2)==undef)then
        val(nx-1)=undef
     else
        if(scalex(nx-1)/=0.0)then
           val(nx-1)=0.5*(u(nx)-u(nx-2))/(scalex(nx-1)*dx(nx-1))
        else
           val(nx-1)=0.0
        end if
     end if

!-- データ数のない両端の処理 ---
     if(u(1)==undef.or.u(2)==undef)then
        val(1)=undef
     else
        if(scalex(1)/=0.0)then
           val(1)=(u(2)-u(1))/(scalex(1)*(x(2)-x(1)))
        else
           val(1)=0.0
        end if
     end if
     if(u(nx)==undef.or.u(nx-1)==undef)then
        val(nx)=undef
     else
        if(scalex(nx)/=0.0)then
           val(nx)=(u(nx)-u(nx-1))/(scalex(nx)*(x(nx)-x(nx-1)))
        else
           val(nx)=0.0
        end if
     end if
  else
     do i=3,nx-2
        if(scalex(i)/=0.0)then
           val(i)=(coe23*(u(i+1)-u(i-1))-coe112*(u(i+2)-u(i-2))) /(scalex(i)*dx(i))
        else
           val(i)=0.0
        end if
     end do

!-- データ数のない両端の処理 (両端の 1 つ内側) ---
     if(scalex(2)/=0.0)then
        val(2)=0.5*(u(3)-u(1))/(scalex(2)*dx(2))
     else
        val(2)=0.0
     end if

     if(scalex(nx-1)/=0.0)then
        val(nx-1)=0.5*(u(nx)-u(nx-2))/(scalex(nx-1)*dx(nx-1))
     else
        val(nx-1)=0.0
     end if

!-- データ数のない両端の処理 ---
     if(scalex(1)/=0.0)then
        val(1)=(u(2)-u(1))/(scalex(1)*(x(2)-x(1)))
     else
        val(1)=0.0
     end if
     if(scalex(nx)/=0.0)then
        val(nx)=(u(nx)-u(nx-1))/(scalex(nx)*(x(nx)-x(nx-1)))
     else
        val(nx)=0.0
     end if
  end if

end subroutine grad4_1d
Subroutine :
x(:) :real, intent(in)
: x 方向の座標変数 [m]
y(:) :real, intent(in)
: y 方向の座標変数 [m]
u(size(x),size(y)) :real, intent(in)
: 勾配をとる 2 次元スカラー成分
valx(size(x),size(y)) :real, intent(inout)
: 計算された y 方向の 2 次元勾配ベクトル
valy(size(x),size(y)) :real, intent(inout)
: 計算された y 方向の 2 次元勾配ベクトル
undeff :real, intent(in), optional
hx(size(x),size(y)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y)) :real, intent(in), optional
: y 方向のスケール因子

1 次元スカラー勾配のルーチンを用いて 2 次元勾配のベクトルを計算 $nabla _hp =left(frac{partial p}{partial x} ,\; frac{partial p}{partial y} right) $ を 4 次の中央差分近似で書き換えると, 点 $(i,j)$ での勾配は $left(frac{p_{i+1,j}-p_{i-1,j}}{2dx} ,\; frac{p_{i,j+1}-p_{i,j-1}}{2dy} right) $ $(2/3)*frac{p_{i+1,j}-p_{i-1,j}}{dx} -(p_{i+2,j}-p_{i-2,j})/(12dx),\; (2/3)*frac{p_{i,j+1}-p_{i,j-1}}{dy} -(p_{i,j+2}-p_{i,j-2})/(12dy)$ とできる. これを用いて2次元勾配を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる. 先に用いた 1 次元勾配計算ルーチンを 2 回呼び出すことにしている.

[Source]

subroutine grad4_2d( x, y, u, valx, valy, undeff, hx, hy )
  ! 1 次元スカラー勾配のルーチンを用いて 2 次元勾配のベクトルを計算
  ! $\nabla _hp =\left(\frac{\partial p}{\partial x} ,\; \frac{\partial p}{\partial y} \right) $ を
  ! 4 次の中央差分近似で書き換えると, 点 $(i,j)$ での勾配は
  ! $\left(\frac{p_{i+1,j}-p_{i-1,j}}{2dx} ,\; \frac{p_{i,j+1}-p_{i,j-1}}{2dy} \right) $
  ! $(2/3)*\frac{p_{i+1,j}-p_{i-1,j}}{dx} -(p_{i+2,j}-p_{i-2,j})/(12dx),\; 
  ! (2/3)*\frac{p_{i,j+1}-p_{i,j-1}}{dy} -(p_{i,j+2}-p_{i,j-2})/(12dy)$
  ! とできる. これを用いて2次元勾配を計算.
  ! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
  ! 落ちる.
  ! 先に用いた 1 次元勾配計算ルーチンを 2 回呼び出すことにしている.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の座標変数 [m]
  real, intent(in) :: y(:)  ! y 方向の座標変数 [m]
  real, intent(in) :: u(size(x),size(y))  ! 勾配をとる 2 次元スカラー成分
  real, intent(inout) :: valx(size(x),size(y))  ! 計算された y 方向の 2 次元勾配ベクトル
  real, intent(inout) :: valy(size(x),size(y))  ! 計算された y 方向の 2 次元勾配ベクトル
  real, intent(in), optional :: undeff
  real, intent(in), optional :: hx(size(x),size(y))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y))  ! y 方向のスケール因子
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: nx  ! x の配列要素数(size 関数で自動的に計算)
  integer :: ny  ! y の配列要素数(size 関数で自動的に計算)
  real :: scalex(size(x),size(y)), scaley(size(x),size(y))

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

  if(present(hx).or.present(hy))then
     do j=1,ny
        do i=1,nx
           scalex(i,j)=hx(i,j)
           scaley(i,j)=hy(i,j)
        end do
     end do
  else
     do j=1,ny
        do i=1,nx
           scalex(i,j)=1.0
           scaley(i,j)=1.0
        end do
     end do
  end if

  if(present(undeff))then
     do i=1,ny
        call grad4_1d(x, u(:,i), valx(:,i), undeff)
     end do

     do i=1,nx
        call grad4_1d(y, u(i,:), valy(i,:), undeff)
     end do

     do j=1,ny
        do i=1,nx
           if(u(i,j)/=undeff)then
              if(scalex(i,j)/=0.0)then
                 valx(i,j)=valx(i,j)/scalex(i,j)
              else
                 valx(i,j)=0.0
              end if
              if(scaley(i,j)/=0.0)then
                 valy(i,j)=valy(i,j)/scaley(i,j)
              else
                 valy(i,j)=0.0
              end if
!-- ここで, else しないのは, grad_1d ルーチンですでに undeff が入っているから同じ作業に
!-- なるので, 割愛.
           end if
        end do
     end do

  else

     do i=1,ny
        call grad4_1d(x, u(:,i), valx(:,i) )
     end do

     do i=1,nx
        call grad4_1d(y, u(i,:), valy(i,:) )
     end do

     do j=1,ny
        do i=1,nx
           if(scalex(i,j)/=0.0)then
              valx(i,j)=valx(i,j)/scalex(i,j)
           else
              valx(i,j)=0.0
           end if
           if(scaley(i,j)/=0.0)then
              valy(i,j)=valy(i,j)/scaley(i,j)
           else
              valy(i,j)=0.0
           end if
        end do
     end do

  end if

end subroutine grad4_2d
Subroutine :
x(:) :real, intent(in)
: 空間座標
u(size(x)) :real, intent(in)
: 上の空間配列に対応する 1 次元スカラー値
val(size(x)) :real, intent(inout)
: スカラー値の x 方向の勾配
undef :real, intent(in), optional
hx(size(x)) :real, intent(in), optional
: x のスケール因子

1 次元のスカラー変数の勾配を計算する $frac{partial p}{partial x} $ を 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は $frac{p_{i+1}-p_{i-1}}{2dx} $ とできる. これを用いて 1 次元勾配を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる.

[Source]

subroutine grad_1d( x, u, val, undef, hx )
! 1 次元のスカラー変数の勾配を計算する
! $\frac{\partial p}{\partial x} $ を
! 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は
! $\frac{p_{i+1}-p_{i-1}}{2dx} $
! とできる. これを用いて 1 次元勾配を計算.
! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
! 落ちる.
  implicit none
  real, intent(in) :: x(:)  ! 空間座標
  real, intent(in) :: u(size(x))  ! 上の空間配列に対応する 1 次元スカラー値
  real, intent(inout) :: val(size(x))  ! スカラー値の x 方向の勾配
  real, intent(in), optional :: undef
  real, intent(in), optional :: hx(size(x))  ! x のスケール因子
  integer :: i  ! イタレーション用添字
  integer :: nx  ! 配列要素数
  real :: scalex(size(x))

  nx=size(x)

  if(present(hx))then
     do i=1,nx
        scalex(i)=hx(i)
     end do
  else
     do i=1,nx
        scalex(i)=1.0
     end do
  end if

  if(present(undef))then
     do i=2,nx-1
        if(u(i+1)==undef.or.u(i-1)==undef)then
           val(i)=undef
        else
           if(scalex(i)/=0.0)then
              val(i)=(u(i+1)-u(i-1))/(scalex(i)*(x(i+1)-x(i-1)))
           else
              val(i)=0.0
           end if
        end if
     end do
!-- データ数のない両端の処理 ---
     if(u(1)==undef.or.u(2)==undef)then
        val(1)=undef
     else
        if(scalex(1)/=0.0)then
           val(1)=(u(2)-u(1))/(scalex(1)*(x(2)-x(1)))
        else
           val(1)=0.0
        end if
     end if
     if(u(nx)==undef.or.u(nx-1)==undef)then
        val(nx)=undef
     else
        if(scalex(nx)/=0.0)then
           val(nx)=(u(nx)-u(nx-1))/(scalex(nx)*(x(nx)-x(nx-1)))
        else
           val(nx)=0.0
        end if
     end if
  else
     do i=2,nx-1
        if(scalex(i)/=0.0)then
           val(i)=(u(i+1)-u(i-1))/(scalex(i)*(x(i+1)-x(i-1)))
        else
           val(i)=0.0
        end if
     end do
!-- データ数のない両端の処理 ---
     if(scalex(1)/=0.0)then
        val(1)=(u(2)-u(1))/(scalex(1)*(x(2)-x(1)))
     else
        val(1)=0.0
     end if
     if(scalex(nx)/=0.0)then
        val(nx)=(u(nx)-u(nx-1))/(scalex(nx)*(x(nx)-x(nx-1)))
     else
        val(nx)=0.0
     end if
  end if

end subroutine grad_1d
Subroutine :
x(:) :real, intent(in)
: x 方向の座標変数 [m]
y(:) :real, intent(in)
: y 方向の座標変数 [m]
u(size(x),size(y)) :real, intent(in)
: 勾配をとる 2 次元スカラー成分
valx(size(x),size(y)) :real, intent(inout)
: 計算された y 方向の 2 次元勾配ベクトル
valy(size(x),size(y)) :real, intent(inout)
: 計算された y 方向の 2 次元勾配ベクトル
undeff :real, intent(in), optional
hx(size(x),size(y)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y)) :real, intent(in), optional
: y 方向のスケール因子

1 次元スカラー勾配のルーチンを用いて 2 次元勾配のベクトルを計算 $nabla _hp =left(frac{partial p}{partial x} ,\; frac{partial p}{partial y} right) $ を 2 次の中央差分近似で書き換えると, 点 $(i,j)$ での勾配は $left(frac{p_{i+1,j}-p_{i-1,j}}{2dx} ,\; frac{p_{i,j+1}-p_{i,j-1}}{2dy} right) $ とできる. これを用いて2次元勾配を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる. 先に用いた 1 次元勾配計算ルーチンを 2 回呼び出すことにしている.

[Source]

subroutine grad_2d( x, y, u, valx, valy, undeff, hx, hy )
  ! 1 次元スカラー勾配のルーチンを用いて 2 次元勾配のベクトルを計算
  ! $\nabla _hp =\left(\frac{\partial p}{\partial x} ,\; \frac{\partial p}{\partial y} \right) $ を
  ! 2 次の中央差分近似で書き換えると, 点 $(i,j)$ での勾配は
  ! $\left(\frac{p_{i+1,j}-p_{i-1,j}}{2dx} ,\; \frac{p_{i,j+1}-p_{i,j-1}}{2dy} \right) $
  ! とできる. これを用いて2次元勾配を計算.
  ! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
  ! 落ちる.
  ! 先に用いた 1 次元勾配計算ルーチンを 2 回呼び出すことにしている.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の座標変数 [m]
  real, intent(in) :: y(:)  ! y 方向の座標変数 [m]
  real, intent(in) :: u(size(x),size(y))  ! 勾配をとる 2 次元スカラー成分
  real, intent(inout) :: valx(size(x),size(y))  ! 計算された y 方向の 2 次元勾配ベクトル
  real, intent(inout) :: valy(size(x),size(y))  ! 計算された y 方向の 2 次元勾配ベクトル
  real, intent(in), optional :: undeff
  real, intent(in), optional :: hx(size(x),size(y))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y))  ! y 方向のスケール因子
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: nx  ! x の配列要素数(size 関数で自動的に計算)
  integer :: ny  ! y の配列要素数(size 関数で自動的に計算)
  real :: scalex(size(x),size(y)), scaley(size(x),size(y))

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

  if(present(hx).or.present(hy))then
     do j=1,ny
        do i=1,nx
           scalex(i,j)=hx(i,j)
           scaley(i,j)=hy(i,j)
        end do
     end do
  else
     do j=1,ny
        do i=1,nx
           scalex(i,j)=1.0
           scaley(i,j)=1.0
        end do
     end do
  end if

  if(present(undeff))then
     do i=1,ny
        call grad_1d(x, u(:,i), valx(:,i), undeff)
     end do

     do i=1,nx
        call grad_1d(y, u(i,:), valy(i,:), undeff)
     end do

     do j=1,ny
        do i=1,nx
           if(u(i,j)/=undeff)then
              if(scalex(i,j)/=0.0)then
                 valx(i,j)=valx(i,j)/scalex(i,j)
              else
                 valx(i,j)=0.0
              end if
              if(scaley(i,j)/=0.0)then
                 valy(i,j)=valy(i,j)/scaley(i,j)
              else
                 valy(i,j)=0.0
              end if
!-- ここで, else しないのは, grad_1d ルーチンですでに undeff が入っているから同じ作業に
!-- なるので, 割愛.
           end if
        end do
     end do

  else

     do i=1,ny
        call grad_1d(x, u(:,i), valx(:,i) )
     end do

     do i=1,nx
        call grad_1d(y, u(i,:), valy(i,:) )
     end do

     do j=1,ny
        do i=1,nx
           if(scalex(i,j)/=0.0)then
              valx(i,j)=valx(i,j)/scalex(i,j)
           else
              valx(i,j)=0.0
           end if
           if(scaley(i,j)/=0.0)then
              valy(i,j)=valy(i,j)/scaley(i,j)
           else
              valy(i,j)=0.0
           end if
        end do
     end do

  end if

end subroutine grad_2d
Subroutine :
x(:) :real, intent(in)
: x 方向の座標変数 [m]
y(:) :real, intent(in)
: y 方向の座標変数 [m]
z(:) :real, intent(in)
: z 方向の座標変数 [m]
u(size(x),size(y),size(z)) :real, intent(in)
: 勾配をとる 2 次元スカラー成分
valx(size(x),size(y),size(z)) :real, intent(inout)
: 計算された y 方向の 2 次元勾配ベクトル
valy(size(x),size(y),size(z)) :real, intent(inout)
: 計算された y 方向の 2 次元勾配ベクトル
valz(size(x),size(y),size(z)) :real, intent(inout)
: 計算された z 方向の 2 次元勾配ベクトル
undeff :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子

1 次元スカラー勾配のルーチンを用いて 3 次元勾配のベクトルを計算. $nabla p =left(frac{partial p}{partial x} ,\; frac{partial p}{partial y} ,\; frac{partial p}{partial z} right) $ を 2 次の中央差分近似で書き換えると, 点 $(i,j,k)$ での勾配は $left(frac{p_{i+1,j,k}-p_{i-1,j,k}}{2dx} ,\; frac{p_{i,j+1,k}-p_{i,j-1,k}}{2dy} ,\; frac{p_{i,j,k+1}-p_{i,j,k-1}}{2dz} right) $ とできる. これを用いて 3 次元勾配を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる. 先に用いた 1 次元勾配計算ルーチンを 3 回呼び出すことにしている.

[Source]

subroutine grad_3d( x, y, z, u, valx, valy, valz, undeff, hx, hy, hz )
  ! 1 次元スカラー勾配のルーチンを用いて 3 次元勾配のベクトルを計算.
  ! $\nabla p =\left(\frac{\partial p}{\partial x} ,\; \frac{\partial p}{\partial y} ,\; \frac{\partial p}{\partial z} \right) $ を
  ! 2 次の中央差分近似で書き換えると, 点 $(i,j,k)$ での勾配は
  ! $\left(\frac{p_{i+1,j,k}-p_{i-1,j,k}}{2dx} ,\; \frac{p_{i,j+1,k}-p_{i,j-1,k}}{2dy} ,\; \frac{p_{i,j,k+1}-p_{i,j,k-1}}{2dz} \right) $
  ! とできる. これを用いて 3 次元勾配を計算.
  ! データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が
  ! 落ちる.
  ! 先に用いた 1 次元勾配計算ルーチンを 3 回呼び出すことにしている.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の座標変数 [m]
  real, intent(in) :: y(:)  ! y 方向の座標変数 [m]
  real, intent(in) :: z(:)  ! z 方向の座標変数 [m]
  real, intent(in) :: u(size(x),size(y),size(z))  ! 勾配をとる 2 次元スカラー成分
  real, intent(inout) :: valx(size(x),size(y),size(z))  ! 計算された y 方向の 2 次元勾配ベクトル
  real, intent(inout) :: valy(size(x),size(y),size(z))  ! 計算された y 方向の 2 次元勾配ベクトル
  real, intent(inout) :: valz(size(x),size(y),size(z))  ! 計算された z 方向の 2 次元勾配ベクトル
  real, intent(in), optional :: undeff
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez

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

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

  if(present(undeff))then
!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do i=1,nz
        do j=1,ny
           call grad_1d(x, u(:,j,i), valx(:,j,i),undeff)
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do i=1,nz
        do j=1,nx
           call grad_1d(y, u(j,:,i), valy(j,:,i),undeff)
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do i=1,ny
        do j=1,nx
           call grad_1d(z, u(j,i,:), valz(j,i,:),undeff)
        end do
     end do
!$omp end do
!$omp end parallel

!$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(u(i,j,k)/=undeff)then
                 if(scalex(i,j,k)/=0.0)then
                    valx(i,j,k)=valx(i,j,k)/scalex(i,j,k)
                 else
                    valx(i,j,k)=0.0
                 end if
                 if(scaley(i,j,k)/=0.0)then
                    valy(i,j,k)=valy(i,j,k)/scaley(i,j,k)
                 else
                    valy(i,j,k)=0.0
                 end if
                 if(scalez(i,j,k)/=0.0)then
                    valz(i,j,k)=valz(i,j,k)/scalez(i,j,k)
                 else
                    valz(i,j,k)=0.0
                 end if
              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)
     do i=1,nz
        do j=1,ny
           call grad_1d(x, u(:,j,i), valx(:,j,i))
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do i=1,nz
        do j=1,nx
           call grad_1d(y, u(j,:,i), valy(j,:,i))
        end do
     end do
!$omp end do
!$omp end parallel

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j)
     do i=1,ny
        do j=1,nx
           call grad_1d(z, u(j,i,:), valz(j,i,:))
        end do
     end do
!$omp end do
!$omp end parallel

!$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(scalex(i,j,k)/=0.0)then
                 valx(i,j,k)=valx(i,j,k)/scalex(i,j,k)
              else
                 valx(i,j,k)=0.0
              end if
              if(scaley(i,j,k)/=0.0)then
                 valy(i,j,k)=valy(i,j,k)/scaley(i,j,k)
              else
                 valy(i,j,k)=0.0
              end if
              if(scalez(i,j,k)/=0.0)then
                 valz(i,j,k)=valz(i,j,k)/scalez(i,j,k)
              else
                 valz(i,j,k)=0.0
              end if
           end do
        end do
     end do
!$omp end do
!$omp end parallel

  end if

end subroutine grad_3d
Subroutine :
x(:) :real, intent(in)
: x 方向の座標変数 [m]
u(size(x)) :real, intent(in)
: 上の空間配列に対応する 1 次元スカラー値
val(size(x)) :real, intent(inout)
: スカラー値の x 方向の勾配
undef :real, intent(in), optional

1 次元のスカラー変数のラプラシアンを計算する $frac{partial ^2p}{partial x^2} $ を 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は $frac{p_{i+1}+p_{i-1}-2p_i}{dx^2} $ とできる. これを用いて 1 次元ラプラシアンを計算. データ点が足りない隅の領域では, undef を定義する. option で undef が定義されていない場合は, 0.0 を代入する.

[Source]

subroutine laplacian_1d( x, u, val, undef )
! 1 次元のスカラー変数のラプラシアンを計算する
! $\frac{\partial ^2p}{\partial x^2} $ を
! 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は
! $\frac{p_{i+1}+p_{i-1}-2p_i}{dx^2} $
! とできる. これを用いて 1 次元ラプラシアンを計算.
! データ点が足りない隅の領域では, undef を定義する.
! option で undef が定義されていない場合は, 0.0 を代入する.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の座標変数 [m]
  real, intent(in) :: u(size(x))  ! 上の空間配列に対応する 1 次元スカラー値
  real, intent(inout) :: val(size(x))  ! スカラー値の x 方向の勾配
  real, intent(in), optional :: undef
  integer :: i  ! イタレーション用添字
  integer :: nx  ! 配列要素数

  nx=size(x)

  if(present(undef))then
     do i=2,nx-1
        if(u(i+1)==undef.or.u(i-1)==undef.or.u(i)==undef)then
           val(i)=undef
        else
           val(i)=4.0*(u(i+1)+u(i-1)-2.0*u(i))/((x(i+1)-x(i-1))**2)
        end if
     end do
!-- データ数のない両端の処理 ---
!     if(u(1)==undef.or.u(2)==undef)then
!        val(1)=undef
!     else
!        val(1)=(u(2)-u(1))/dx
!     end if
!     if(u(nx)==undef.or.u(nx-1)==undef)then
!        val(nx)=undef
!     else
!        val(nx)=(u(nx)-u(nx-1))/dx
!     end if
     val(1)=undef
     val(nx)=undef
  else
     do i=2,nx-1
        val(i)=4.0*(u(i+1)+u(i-1)-2.0*u(i))/((x(i+1)-x(i-1))**2)
     end do
!-- データ数のない両端の処理 ---
!     val(1)=(u(2)-u(1))/dx
!     val(nx)=(u(nx)-u(nx-1))/dx
     val(1)=0.0
     val(nx)=0.0
  end if

end subroutine laplacian_1d
Subroutine :
x(:) :real, intent(in)
: x 方向の座標変数 [m]
y(:) :real, intent(in)
: y 方向の座標変数 [m]
u(size(x),size(y)) :real, intent(in)
: 上の空間配列に対応する 2 次元スカラー値
val(size(x),size(y)) :real, intent(inout)
: スカラー値の 2 次元方向のラプラシアン
undef :real, intent(in), optional

2 次元のスカラー変数のラプラシアンを計算する $frac{partial ^2p}{partial x^2} $ を 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は $frac{p_{i+1}+p_{i-1}-2p_i}{dx^2} $ とできる. これを用いて 1 次元ラプラシアンを計算. データ点が足りない隅の領域では, undef を定義する. option で undef が定義されていない場合は, 0.0 を代入する. 一般直交座標系への対応はまだ.

[Source]

subroutine laplacian_2d( x, y, u, val, undef )
! 2 次元のスカラー変数のラプラシアンを計算する
! $\frac{\partial ^2p}{\partial x^2} $ を
! 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は
! $\frac{p_{i+1}+p_{i-1}-2p_i}{dx^2} $
! とできる. これを用いて 1 次元ラプラシアンを計算.
! データ点が足りない隅の領域では, undef を定義する.
! option で undef が定義されていない場合は, 0.0 を代入する.
! 一般直交座標系への対応はまだ.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の座標変数 [m]
  real, intent(in) :: y(:)  ! y 方向の座標変数 [m]
  real, intent(in) :: u(size(x),size(y))  ! 上の空間配列に対応する 2 次元スカラー値
  real, intent(inout) :: val(size(x),size(y))  ! スカラー値の 2 次元方向のラプラシアン
  real, intent(in), optional :: undef
  integer :: i, j  ! イタレーション用添字
  integer :: nx, ny  ! 配列要素数
!  real :: scalex(size(x),size(y))
  real :: tmpu(size(x),size(y)), tmpv(size(x),size(y))

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

  if(present(undef))then

     val=undef

     do j=1,ny
        call laplacian_1d( x, u(1:nx,j), tmpu(1:nx,j), undef )
     end do
     do i=1,nx
        call laplacian_1d( y, u(i,1:ny), tmpv(i,1:ny), undef )
     end do

     do j=1,ny
        do i=1,nx
           if(tmpu(i,j)/=undef.and.tmpv(i,j)/=undef)then
              val(i,j)=tmpu(i,j)+tmpv(i,j)
           end if
        end do
     end do

  else

     do j=1,ny
        call laplacian_1d( x, u(1:nx,j), tmpu(1:nx,j) )
     end do
     do i=1,nx
        call laplacian_1d( y, u(i,1:ny), tmpv(i,1:ny) )
     end do

     do j=1,ny
        do i=1,nx
           val(i,j)=tmpu(i,j)+tmpv(i,j)
        end do
     end do

  end if

end subroutine laplacian_2d
Subroutine :
x(:) :real, intent(in)
: x 方向の座標変数 [m]
y(:) :real, intent(in)
: y 方向の座標変数 [m]
z(:) :real, intent(in)
: z 方向の座標変数 [m]
u(size(x),size(y),size(z)) :real, intent(in)
: 上の空間配列に対応する 3 次元スカラー値
val(size(x),size(y),size(z)) :real, intent(inout)
: スカラー値の 3 次元方向のラプラシアン
undef :real, intent(in), optional

3 次元のスカラー変数のラプラシアンを計算する $frac{partial ^2p}{partial x^2} $ を 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は $frac{p_{i+1}+p_{i-1}-2p_i}{dx^2} $ とできる. これを用いて 1 次元ラプラシアンを計算. データ点が足りない隅の領域では, undef を定義する. option で undef が定義されていない場合は, 0.0 を代入する. 一般直交座標系への対応はまだ.

[Source]

subroutine laplacian_3d( x, y, z, u, val, undef )
! 3 次元のスカラー変数のラプラシアンを計算する
! $\frac{\partial ^2p}{\partial x^2} $ を
! 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は
! $\frac{p_{i+1}+p_{i-1}-2p_i}{dx^2} $
! とできる. これを用いて 1 次元ラプラシアンを計算.
! データ点が足りない隅の領域では, undef を定義する.
! option で undef が定義されていない場合は, 0.0 を代入する.
! 一般直交座標系への対応はまだ.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の座標変数 [m]
  real, intent(in) :: y(:)  ! y 方向の座標変数 [m]
  real, intent(in) :: z(:)  ! z 方向の座標変数 [m]
  real, intent(in) :: u(size(x),size(y),size(z))  ! 上の空間配列に対応する 3 次元スカラー値
  real, intent(inout) :: val(size(x),size(y),size(z))  ! スカラー値の 3 次元方向のラプラシアン
  real, intent(in), optional :: undef
  integer :: i, j, k  ! イタレーション用添字
  integer :: nx, ny, nz  ! 配列要素数
!  real :: scalex(size(x),size(y))
  real :: tmpu(size(x),size(y),size(z)), tmpv(size(x),size(y),size(z))

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

  if(present(undef))then

     val=undef

!$omp parallel default(shared)
!$omp do schedule(dynamic) private(k)

     do k=1,nz
        call laplacian_2d( x, y, u(1:nx,1:ny,k), tmpu(1:nx,1:ny,k), undef )
     end do

!$omp end do
!$omp end parallel

     do j=1,ny
        do i=1,nx
           call laplacian_1d( z, u(i,j,1:nz), tmpv(i,j,1:nz), undef )
        end do
     end do

!$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(tmpu(i,j,k)/=undef.and.tmpv(i,j,k)/=undef)then
                 val(i,j,k)=tmpu(i,j,k)+tmpv(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(dynamic) private(k)

     do k=1,nz
        call laplacian_2d( x, y, u(1:nx,1:ny,k), tmpu(1:nx,1:ny,k) )
     end do

!$omp end do
!$omp end parallel

     do j=1,ny
        do i=1,nx
           call laplacian_1d( z, u(i,j,1:nz), tmpv(i,j,1:nz) )
        end do
     end do

!$omp parallel default(shared)
!$omp do schedule(runtime) private(i,j,k)

     do k=1,nz
        do j=1,ny
           do i=1,nx
              val(i,j,k)=tmpu(i,j,k)+tmpv(i,j,k)
           end do
        end do
     end do

!$omp end do
!$omp end parallel

  end if

end subroutine laplacian_3d
Subroutine :
x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
z(:) :real, intent(in)
: z 方向の空間座標 [m]
phi(size(x),size(y),size(z)) :real, intent(in)
: 乱流拡散を計算するスカラー量
rho(size(x),size(y),size(z)) :real, intent(in)
: 基本場の密度 [kg/m^3]
nuh(size(x),size(y),size(z)) :real, intent(in)
: 水平渦粘性係数
nuv(size(x),size(y),size(z)) :real, intent(in)
: 鉛直渦粘性係数
val(size(x),size(y),size(z)) :real, intent(inout)
: 3 次元発散値
undef :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子
sfcphi(size(x),size(y)) :real, intent(in), optional
: 地表面からのフラックス

乱流拡散項を計算する.

[Source]

subroutine turb_diff( x, y, z, phi, rho, nuh, nuv, val, undef, hx, hy, hz, sfcphi )
! 乱流拡散項を計算する.
  implicit none
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: z(:)  ! z 方向の空間座標 [m]
  real, intent(in) :: phi(size(x),size(y),size(z))  ! 乱流拡散を計算するスカラー量
  real, intent(in) :: rho(size(x),size(y),size(z))  ! 基本場の密度 [kg/m^3]
  real, intent(in) :: nuh(size(x),size(y),size(z))  ! 水平渦粘性係数
  real, intent(in) :: nuv(size(x),size(y),size(z))  ! 鉛直渦粘性係数
  real, intent(inout) :: val(size(x),size(y),size(z))  ! 3 次元発散値
  real, intent(in), optional :: undef
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  real, intent(in), optional :: sfcphi(size(x),size(y))      ! 地表面からのフラックス
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real, dimension(size(x),size(y),size(z),3) :: tau  ! signal 方向に
              ! 作用する 1,2,3 面に垂直な応力
  character(1) :: signaltau(3)
  integer :: id
  real, dimension(size(x),size(y)) :: stau
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez

  signaltau=(/ '1', '2', '3' /)

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

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

  val=0.0

  do id=1,3
     if(present(sfcphi))then
        stau(:,:)=sfcphi(:,:)
        if(present(undef))then
           call Reynolds_scal( signaltau(id), x, y, z, phi, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez, sfcphi=stau, undef=undef )
        else
           call Reynolds_scal( signaltau(id), x, y, z, phi, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez, sfcphi=stau )
        end if
     else
        if(present(undef))then
           call Reynolds_scal( signaltau(id), x, y, z, phi, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez, undef=undef )
        else
           call Reynolds_scal( signaltau(id), x, y, z, phi, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez )
        end if
     end if
  end do

!-- 乱流項は 3 次元発散と同じ形となるので, div_3d でまとめる.
  if(present(undef))then
     call div_3d( x, y, z, tau(:,:,:,1), tau(:,:,:,2), tau(:,:,:,3), val, hx=scalex, hy=scaley, hz=scalez, undeff=undef )
  else
     call div_3d( x, y, z, tau(:,:,:,1), tau(:,:,:,2), tau(:,:,:,3), val, hx=scalex, hy=scaley, hz=scalez )
  end if

end subroutine
Subroutine :
signal :character(1)
: デカルト座標系の何番目の乱流成分かを判定する.
1
= デカルト座標右手系における x 座標成分 (方程式 u 成分)
2
= デカルト座標右手系における y 座標成分 (方程式 v 成分)
3
= デカルト座標右手系における z 座標成分 (方程式 w 成分)
x(:) :real, intent(in)
: x 方向の空間座標 [m]
y(:) :real, intent(in)
: y 方向の空間座標 [m]
z(:) :real, intent(in)
: z 方向の空間座標 [m]
u(size(x),size(y),size(z)) :real, intent(in)
: x に対応する方向の 2 次元ベクトル成分
v(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
w(size(x),size(y),size(z)) :real, intent(in)
: y に対応する方向の 2 次元ベクトル成分
rho(size(x),size(y),size(z)) :real, intent(in)
: 基本場の密度 [kg/m^3]
nuh(size(x),size(y),size(z)) :real, intent(in)
: 水平渦粘性係数
nuv(size(x),size(y),size(z)) :real, intent(in)
: 鉛直渦粘性係数
val(size(x),size(y),size(z)) :real, intent(inout)
: 乱流フラックス
undef :real, intent(in), optional
hx(size(x),size(y),size(z)) :real, intent(in), optional
: x 方向のスケール因子
hy(size(x),size(y),size(z)) :real, intent(in), optional
: y 方向のスケール因子
hz(size(x),size(y),size(z)) :real, intent(in), optional
: z 方向のスケール因子
sfctau(size(x),size(y)) :real, intent(in), optional
: 地表面からのフラックス

デカルト座標系における乱流粘性項を計算する.

[Source]

subroutine turb_visc( signal, x, y, z, u, v, w, rho, nuh, nuv, val, undef, hx, hy, hz, sfctau )
! デカルト座標系における乱流粘性項を計算する.
  implicit none
  character(1) :: signal  ! デカルト座標系の何番目の乱流成分かを判定する.
                  ! [1] = デカルト座標右手系における x 座標成分 (方程式 u 成分)
                  ! [2] = デカルト座標右手系における y 座標成分 (方程式 v 成分)
                  ! [3] = デカルト座標右手系における z 座標成分 (方程式 w 成分)
  real, intent(in) :: x(:)  ! x 方向の空間座標 [m]
  real, intent(in) :: y(:)  ! y 方向の空間座標 [m]
  real, intent(in) :: z(:)  ! z 方向の空間座標 [m]
  real, intent(in) :: u(size(x),size(y),size(z))  ! x に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: v(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: w(size(x),size(y),size(z))  ! y に対応する方向の 2 次元ベクトル成分
  real, intent(in) :: rho(size(x),size(y),size(z))  ! 基本場の密度 [kg/m^3]
  real, intent(in) :: nuh(size(x),size(y),size(z))  ! 水平渦粘性係数
  real, intent(in) :: nuv(size(x),size(y),size(z))  ! 鉛直渦粘性係数
  real, intent(inout) :: val(size(x),size(y),size(z))  ! 乱流フラックス
  real, intent(in), optional :: undef
  real, intent(in), optional :: hx(size(x),size(y),size(z))  ! x 方向のスケール因子
  real, intent(in), optional :: hy(size(x),size(y),size(z))  ! y 方向のスケール因子
  real, intent(in), optional :: hz(size(x),size(y),size(z))  ! z 方向のスケール因子
  real, intent(in), optional :: sfctau(size(x),size(y))  ! 地表面からのフラックス
  integer :: i   ! イタレーション用添字
  integer :: j   ! イタレーション用添字
  integer :: k   ! イタレーション用添字
  integer :: nx  ! 空間配列要素数 1 次元目
  integer :: ny  ! 空間配列要素数 2 次元目
  integer :: nz  ! 空間配列要素数 3 次元目
  real, dimension(size(x),size(y),size(z),3) :: tau  ! signal 方向に
              ! 作用する 1,2,3 面に垂直な応力
  character(1) :: signaltau(3)
  integer :: id
  real, dimension(size(x),size(y)) :: stau
  real, dimension(size(x),size(y),size(z)) :: scalex, scaley, scalez

  signaltau=(/ '1', '2', '3' /)
  nx=size(x)
  ny=size(y)
  nz=size(z)

  if(present(hx).and.present(hy).and.present(hz))then
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=hx(i,j,k)
              scaley(i,j,k)=hy(i,j,k)
              scalez(i,j,k)=hz(i,j,k)
           end do
        end do
     end do
  else
     do k=1,nz
        do j=1,ny
           do i=1,nx
              scalex(i,j,k)=1.0
              scaley(i,j,k)=1.0
              scalez(i,j,k)=1.0
           end do
        end do
     end do
  end if

  val=0.0

  do id=1,3
     if(present(sfctau))then
        stau(:,:)=sfctau(:,:)
        if(present(undef))then
           call Reynolds_stress( signal//signaltau(id), x, y, z, u, v, w, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez, sfctau=stau, undef=undef )
        else
           call Reynolds_stress( signal//signaltau(id), x, y, z, u, v, w, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez, sfctau=stau )
        end if
     else
        if(present(undef))then
           call Reynolds_stress( signal//signaltau(id), x, y, z, u, v, w, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez, undef=undef )
        else
           call Reynolds_stress( signal//signaltau(id), x, y, z, u, v, w, rho, nuh, nuv, tau(:,:,:,id), hx=scalex, hy=scaley, hz=scalez )
        end if
     end if
  end do

!-- 乱流項の計算は, 3 次元発散を行うことと等価であるため,
!-- 以下では, div_3d を用いて計算を行う.

  if(present(undef))then
     call div_3d( x, y, z, tau(:,:,:,1), tau(:,:,:,2), tau(:,:,:,3), val, hx=scalex, hy=scaley, hz=scalez, undeff=undef )
  else
     call div_3d( x, y, z, tau(:,:,:,1), tau(:,:,:,2), tau(:,:,:,3), val, hx=scalex, hy=scaley, hz=scalez )
  end if

end subroutine
Subroutine :
z(:,:,:) :real, intent(in)
: デカルト幾何座標系鉛直座標 [m]
zf(size(z,1),size(z,2)) :real, intent(in)
: 地表面高度 [m]
zt(size(z,1),size(z,2)) :real, intent(in)
: 定義域最高度 [m]
zeta(size(z,1),size(z,2),size(z,3)) :real, intent(inout)
: terrain following 座標 [m]

幾何座標 z を terrain following 座標に変換する.

[Source]

subroutine z_2_zeta( z, zf, zt, zeta )
! 幾何座標 z を terrain following 座標に変換する.
  implicit none
  real, intent(in) :: z(:,:,:)    ! デカルト幾何座標系鉛直座標 [m]
  real, intent(in) :: zf(size(z,1),size(z,2))   ! 地表面高度 [m]
  real, intent(in) :: zt(size(z,1),size(z,2))   ! 定義域最高度 [m]
  real, intent(inout) :: zeta(size(z,1),size(z,2),size(z,3))  ! terrain following 座標 [m]

  integer :: i, j, k
  integer :: nx, ny, nz

  nx=size(z,1)
  ny=size(z,2)
  nz=size(z,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
           zeta(i,j,k)=zt(i,j)*(z(i,j,k)-zf(i,j))/(zt(i,j)-zf(i,j))
        end do
     end do
  end do

!$omp end do
!$omp end parallel

end subroutine z_2_zeta
Subroutine :
x :real, dimension(:), intent(in)
: 東西方向の座標 [m]
y :real, dimension(:), intent(in)
: 東西方向の座標 [m]
zeta :real, dimension(size(x),size(y)), intent(in)
: terrain 系鉛直座標 [m]
zf :real, dimension(size(x),size(y)), intent(in)
: 地表面高度座標 [m]
zt :real, dimension(size(x),size(y)), intent(in)
: 座標最上端 [m]
u :real, dimension(size(x),size(y)), intent(in)
: zeta における東西風 [m/s]
v :real, dimension(size(x),size(y)), intent(in)
: zeta における南北風 [m/s]
w :real, dimension(size(x),size(y)), intent(in)
: zeta における鉛直風 [m/s]
wh :real, dimension(size(x),size(y)), intent(inout)
: デカルト系における鉛直風 [m/s]
undef :real, intent(in), optional
: 欠損値

terrain following 座標系で評価されている各風速成分をデカルト座標系での 各風速成分に変換する. ただし, terrain following 系の水平風速はデカルト座標系でも 大きさが変化しないため, このルーチンでは鉛直風速のみ変換を行う. また, 3 次元デカルト系の格子点上に変換するのではなく, terrain following 系の 各点で風速成分をデカルト系方向に変換するだけ. 現在, 水平方向にはデカルト系にのみ対応している. もし, 座標点も変換する場合は, vert_coord_conv モジュールを使用のこと.

[Source]

subroutine zast_2_w_2d( x, y, zeta, zf, zt, u, v, w, wh, undef )
! terrain following 座標系で評価されている各風速成分をデカルト座標系での
! 各風速成分に変換する.
! ただし, terrain following 系の水平風速はデカルト座標系でも
! 大きさが変化しないため, このルーチンでは鉛直風速のみ変換を行う.
! また, 3 次元デカルト系の格子点上に変換するのではなく, terrain following 系の
! 各点で風速成分をデカルト系方向に変換するだけ.
! 現在, 水平方向にはデカルト系にのみ対応している.
! もし, 座標点も変換する場合は, vert_coord_conv モジュールを使用のこと.
  implicit none
  real, dimension(:), intent(in) :: x  ! 東西方向の座標 [m]
  real, dimension(:), intent(in) :: y  ! 東西方向の座標 [m]
  real, dimension(size(x),size(y)), intent(in) :: zeta  ! terrain 系鉛直座標 [m]
  real, dimension(size(x),size(y)), intent(in) :: zf  ! 地表面高度座標 [m]
  real, dimension(size(x),size(y)), intent(in) :: zt  ! 座標最上端 [m]
  real, dimension(size(x),size(y)), intent(in) :: u  ! zeta における東西風 [m/s]
  real, dimension(size(x),size(y)), intent(in) :: v  ! zeta における南北風 [m/s]
  real, dimension(size(x),size(y)), intent(in) :: w  ! zeta における鉛直風 [m/s]
  real, dimension(size(x),size(y)), intent(inout) :: wh  ! デカルト系における鉛直風 [m/s]
  real, intent(in), optional :: undef  ! 欠損値
  integer :: i, j, nx, ny
  real, dimension(size(x),size(y)) :: dx, dy
  real :: j31, j32, jd, coe

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

  call grad_2d( x, y, zf, dx, dy )

  if(present(undef))then
     do j=1,ny
        do i=1,nx
           if(u(i,j)/=undef.and.v(i,j)/=undef.and.w(i,j)/=undef)then
              jd=1.0-zf(i,j)/zt(i,j)
              coe=zeta(i,j)/zt(i,j)-1.0
              j31=coe*dx(i,j)
              j32=coe*dy(i,j)
              wh(i,j)=(u(i,j)*j31+v(i,j)*j32+w(i,j))/jd
           else
              wh(i,j)=undef
           end if
        end do
     end do
  else
     do j=1,ny
        do i=1,nx
           jd=1.0-zf(i,j)/zt(i,j)
           coe=zeta(i,j)/zt(i,j)-1.0
           j31=coe*dx(i,j)
           j32=coe*dy(i,j)
           wh(i,j)=(u(i,j)*j31+v(i,j)*j32+w(i,j))/jd
        end do
     end do
  end if

end subroutine zast_2_w_2d
Subroutine :
x :real, dimension(size(zeta,1)), intent(in)
: 東西方向の座標 [m]
y :real, dimension(size(zeta,2)), intent(in)
: 東西方向の座標 [m]
zeta :real, dimension(:,:,:), intent(in)
: terrain 系の鉛直座標 [m]
zf :real, dimension(size(zeta,1),size(zeta,2)), intent(in)
: 地表面高度座標 [m]
zt :real, dimension(size(zeta,1),size(zeta,2)), intent(in)
: 座標最上端 [m]
u :real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(in)
: zeta における東西風 [m/s]
v :real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(in)
: zeta における南北風 [m/s]
w :real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(in)
: zeta における鉛直風 [m/s]
wh :real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(inout)
: デカルト系における鉛直風 [m/s]
undef :real, intent(in), optional
: 欠損値

terrain following 座標系で評価されている各風速成分をデカルト座標系での 各風速成分に変換する. ただし, terrain following 系の水平風速はデカルト座標系でも 大きさが変化しないため, このルーチンでは鉛直風速のみ変換を行う. また, 3 次元デカルト系の格子点上に変換するのではなく, terrain following 系の 各点で風速成分をデカルト系方向に変換するだけ. 現在, 水平方向にはデカルト系にのみ対応している. もし, 座標点も変換する場合は, vert_coord_conv モジュールを使用のこと.

[Source]

subroutine zast_2_w_3d( x, y, zeta, zf, zt, u, v, w, wh, undef )
! terrain following 座標系で評価されている各風速成分をデカルト座標系での
! 各風速成分に変換する.
! ただし, terrain following 系の水平風速はデカルト座標系でも
! 大きさが変化しないため, このルーチンでは鉛直風速のみ変換を行う.
! また, 3 次元デカルト系の格子点上に変換するのではなく, terrain following 系の
! 各点で風速成分をデカルト系方向に変換するだけ.
! 現在, 水平方向にはデカルト系にのみ対応している.
! もし, 座標点も変換する場合は, vert_coord_conv モジュールを使用のこと.
  implicit none
  real, dimension(:,:,:), intent(in) :: zeta  ! terrain 系の鉛直座標 [m]
  real, dimension(size(zeta,1)), intent(in) :: x  ! 東西方向の座標 [m]
  real, dimension(size(zeta,2)), intent(in) :: y  ! 東西方向の座標 [m]
  real, dimension(size(zeta,1),size(zeta,2)), intent(in) :: zf  ! 地表面高度座標 [m]
  real, dimension(size(zeta,1),size(zeta,2)), intent(in) :: zt  ! 座標最上端 [m]
  real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(in) :: u  ! zeta における東西風 [m/s]
  real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(in) :: v  ! zeta における南北風 [m/s]
  real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(in) :: w  ! zeta における鉛直風 [m/s]
  real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)), intent(inout) :: wh  ! デカルト系における鉛直風 [m/s]
  real, intent(in), optional :: undef  ! 欠損値
  integer :: k, nx, ny, nz

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

  if(present(undef))then

!$omp parallel default(shared)
!$omp do schedule(runtime) private(k)
     do k=1,nz
        call zast_2_w_2d( x, y, zeta(:,:,k), zf, zt, u(:,:,k), v(:,:,k), w(:,:,k), wh(:,:,k), undef )
     end do
!$omp end do
!$omp end parallel

  else

!$omp parallel default(shared)
!$omp do schedule(runtime) private(k)
     do k=1,nz
        call zast_2_w_2d( x, y, zeta(:,:,k), zf, zt, u(:,:,k), v(:,:,k), w(:,:,k), wh(:,:,k) )
     end do
!$omp end do
!$omp end parallel

  end if

end subroutine zast_2_w_3d