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