| Class | Derivation |
| In: |
derivation.f90
|
微分演算計算モジュール
| Subroutine : | |||
| signal : | character(1)
| ||
| x(:) : | real, intent(in)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| phi(size(x),size(y),size(z)) : | real, intent(in)
| ||
| rho(size(x),size(y),size(z)) : | real, intent(in)
| ||
| 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
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| sfcphi(size(x),size(y)) : | real, intent(in), optional
|
スカラー量の乱流フラックスを計算する.
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)
| ||
| x(:) : | real, intent(in)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| u(size(x),size(y),size(z)) : | real, intent(in)
| ||
| v(size(x),size(y),size(z)) : | real, intent(in)
| ||
| w(size(x),size(y),size(z)) : | real, intent(in)
| ||
| rho(size(x),size(y),size(z)) : | real, intent(in)
| ||
| 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
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| sfctau(size(x),size(y)) : | real, intent(in), optional
|
デカルト座標系におけるレイノルズ応力テンソルを計算する.
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)
| ||
| y(:) : | real, intent(in)
| ||
| u(size(x),size(y)) : | real, intent(in)
| ||
| v(size(x),size(y)) : | real, intent(in)
| ||
| val(size(x),size(y)) : | real, intent(inout)
| ||
| undeff : | real, intent(in), optional | ||
| hx(size(x),size(y)) : | real, intent(in), optional
| ||
| hy(size(x),size(y)) : | real, intent(in), optional
| ||
| ord : | logical, intent(in), optional
|
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 が微分計算を担当するので, 境界の計算も自動で行う.
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)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| u(size(x),size(y),size(z)) : | real, intent(in)
| ||
| v(size(x),size(y),size(z)) : | real, intent(in)
| ||
| w(size(x),size(y),size(z)) : | real, intent(in)
| ||
| zeta(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| eta(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| xi(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| undeff : | real, intent(in), optional | ||
| hx(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
|
3 次元渦度を計算する. 引数の順番は右手系で x, y, z のデカルト座標系, それに対応するベクトル成分 u, v, w を代入すると, それに対応した渦度ベクトル成分 zeta, eta, xi が計算される. 実質は grad_1d が計算を担当するので, 境界の処理も自動で行う.
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)
| ||
| x(:) : | real, intent(in)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| u(size(x),size(y),size(z)) : | real, intent(in)
| ||
| v(size(x),size(y),size(z)) : | real, intent(in)
| ||
| w(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
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
|
デカルト座標系における変形速度テンソルを計算する.
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)
| ||
| y(:) : | real, intent(in)
| ||
| u(size(x),size(y)) : | real, intent(in)
| ||
| v(size(x),size(y)) : | real, intent(in)
| ||
| val(size(x),size(y)) : | real, intent(inout)
| ||
| undeff : | real, intent(in), optional | ||
| hx(size(x),size(y)) : | real, intent(in), optional
| ||
| hy(size(x),size(y)) : | real, intent(in), optional
|
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 ) という形で計算を行えば, 境界も自動的に計算可能.
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)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| u(size(x),size(y),size(z)) : | real, intent(in)
| ||
| v(size(x),size(y),size(z)) : | real, intent(in)
| ||
| w(size(x),size(y),size(z)) : | real, intent(in)
| ||
| val(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| undeff : | real, intent(in), optional | ||
| hx(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
|
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 ルーチンが計算を行うので, 境界の処理も自動で行う.
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)
| ||
| val(size(x)) : | real, intent(inout)
| ||
| undef : | real, intent(in), optional | ||
| hx(size(x)) : | real, intent(in), optional
|
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 次の差分近似で計算するので, 少し精度が 落ちる.
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)
| ||
| y(:) : | real, intent(in)
| ||
| u(size(x),size(y)) : | real, intent(in)
| ||
| valx(size(x),size(y)) : | real, intent(inout)
| ||
| valy(size(x),size(y)) : | real, intent(inout)
| ||
| undeff : | real, intent(in), optional | ||
| hx(size(x),size(y)) : | real, intent(in), optional
| ||
| hy(size(x),size(y)) : | real, intent(in), optional
|
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 回呼び出すことにしている.
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)
| ||
| val(size(x)) : | real, intent(inout)
| ||
| undef : | real, intent(in), optional | ||
| hx(size(x)) : | real, intent(in), optional
|
1 次元のスカラー変数の勾配を計算する $frac{partial p}{partial x} $ を 2 次の中央差分近似で書き換えると, 点 $(i)$ での勾配は $frac{p_{i+1}-p_{i-1}}{2dx} $ とできる. これを用いて 1 次元勾配を計算. データ点が足りない隅の領域では, 1 次の差分近似で計算するので, 少し精度が 落ちる.
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)
| ||
| y(:) : | real, intent(in)
| ||
| u(size(x),size(y)) : | real, intent(in)
| ||
| valx(size(x),size(y)) : | real, intent(inout)
| ||
| valy(size(x),size(y)) : | real, intent(inout)
| ||
| undeff : | real, intent(in), optional | ||
| hx(size(x),size(y)) : | real, intent(in), optional
| ||
| hy(size(x),size(y)) : | real, intent(in), optional
|
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 回呼び出すことにしている.
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)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| u(size(x),size(y),size(z)) : | real, intent(in)
| ||
| valx(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| valy(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| valz(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| undeff : | real, intent(in), optional | ||
| hx(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
|
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 回呼び出すことにしている.
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)
| ||
| u(size(x)) : | real, intent(in)
| ||
| val(size(x)) : | real, intent(inout)
| ||
| 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 を代入する.
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)
| ||
| y(:) : | real, intent(in)
| ||
| u(size(x),size(y)) : | real, intent(in)
| ||
| val(size(x),size(y)) : | real, intent(inout)
| ||
| 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 を代入する. 一般直交座標系への対応はまだ.
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)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| u(size(x),size(y),size(z)) : | real, intent(in)
| ||
| val(size(x),size(y),size(z)) : | real, intent(inout)
| ||
| 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 を代入する. 一般直交座標系への対応はまだ.
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)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| phi(size(x),size(y),size(z)) : | real, intent(in)
| ||
| rho(size(x),size(y),size(z)) : | real, intent(in)
| ||
| 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
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| sfcphi(size(x),size(y)) : | real, intent(in), optional
|
乱流拡散項を計算する.
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)
| ||
| x(:) : | real, intent(in)
| ||
| y(:) : | real, intent(in)
| ||
| z(:) : | real, intent(in)
| ||
| u(size(x),size(y),size(z)) : | real, intent(in)
| ||
| v(size(x),size(y),size(z)) : | real, intent(in)
| ||
| w(size(x),size(y),size(z)) : | real, intent(in)
| ||
| rho(size(x),size(y),size(z)) : | real, intent(in)
| ||
| 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
| ||
| hy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| hz(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
| sfctau(size(x),size(y)) : | real, intent(in), optional
|
デカルト座標系における乱流粘性項を計算する.
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)
| ||
| zf(size(z,1),size(z,2)) : | real, intent(in)
| ||
| zt(size(z,1),size(z,2)) : | real, intent(in)
| ||
| zeta(size(z,1),size(z,2),size(z,3)) : | real, intent(inout)
|
幾何座標 z を terrain following 座標に変換する.
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)
| ||
| y : | real, dimension(:), intent(in)
| ||
| zeta : | real, dimension(size(x),size(y)), intent(in)
| ||
| zf : | real, dimension(size(x),size(y)), intent(in)
| ||
| zt : | real, dimension(size(x),size(y)), intent(in)
| ||
| u : | real, dimension(size(x),size(y)), intent(in)
| ||
| v : | real, dimension(size(x),size(y)), intent(in)
| ||
| w : | real, dimension(size(x),size(y)), intent(in)
| ||
| wh : | real, dimension(size(x),size(y)), intent(inout)
| ||
| undef : | real, intent(in), optional
|
terrain following 座標系で評価されている各風速成分をデカルト座標系での 各風速成分に変換する. ただし, terrain following 系の水平風速はデカルト座標系でも 大きさが変化しないため, このルーチンでは鉛直風速のみ変換を行う. また, 3 次元デカルト系の格子点上に変換するのではなく, terrain following 系の 各点で風速成分をデカルト系方向に変換するだけ. 現在, 水平方向にはデカルト系にのみ対応している. もし, 座標点も変換する場合は, vert_coord_conv モジュールを使用のこと.
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)
| ||
| y : | real, dimension(size(zeta,2)), intent(in)
| ||
| zeta : | real, dimension(:,:,:), intent(in)
| ||
| zf : | real, dimension(size(zeta,1),size(zeta,2)), intent(in)
| ||
| zt : | real, dimension(size(zeta,1),size(zeta,2)), intent(in)
| ||
| u : | real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)),
intent(in)
| ||
| v : | real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)),
intent(in)
| ||
| w : | real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)),
intent(in)
| ||
| wh : | real, dimension(size(zeta,1),size(zeta,2),size(zeta,3)),
intent(inout)
| ||
| undef : | real, intent(in), optional
|
terrain following 座標系で評価されている各風速成分をデカルト座標系での 各風速成分に変換する. ただし, terrain following 系の水平風速はデカルト座標系でも 大きさが変化しないため, このルーチンでは鉛直風速のみ変換を行う. また, 3 次元デカルト系の格子点上に変換するのではなく, terrain following 系の 各点で風速成分をデカルト系方向に変換するだけ. 現在, 水平方向にはデカルト系にのみ対応している. もし, 座標点も変換する場合は, vert_coord_conv モジュールを使用のこと.
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