Class | Thermo_Advanced_Routine |
In: |
thermo_advanced_routine.f90
|
基礎ルーチン, 関数集を use して複雑な熱力学関数を計算するモジュール
Subroutine : | |||
z(:) : | real, intent(in)
| ||
pt(size(z)) : | real, intent(in)
| ||
BV(size(z)) : | real, intent(inout)
| ||
undeff : | real, intent(in), optional |
ブラントバイサラ振動数の 2 乗を計算する.
subroutine Brunt_Freq( z, pt, BV, undeff ) ! ブラントバイサラ振動数の 2 乗を計算する. use Derivation use Phys_Const implicit none real, intent(in) :: z(:) ! z 方向の座標変数 [m] real, intent(in) :: pt(size(z)) ! 温位 [K] real, intent(inout) :: BV(size(z)) ! ブラントバイサラ振動数 [1/s] real, intent(in), optional :: undeff integer :: i, j, k integer :: nz ! 第 3 要素数 real :: dz ! z 方向の格子間隔 [m] if(present(undeff))then call grad_1d( z, pt, BV, undeff ) else call grad_1d( z, pt, BV ) end if !-- 以下, grad_1d によって undef 点が増えているので, BV で判定. do k=1,nz if(present(undeff))then if(BV(k)==undeff)then BV(k)=undeff else BV(k)=(g/pt(k))*BV(k) end if else BV(k)=(g/pt(k))*BV(k) end if end do end subroutine
Subroutine : | |
rhop(:,:,:) : | real, intent(in) |
rhob(size(rhop,3)) : | real, intent(in) |
buo(size(rhop,1),size(rhop,2),size(rhob)) : | real, intent(inout) |
qall(size(rhop,1),size(rhop,2),size(rhob)) : | real, intent(in), optional |
subroutine Buoyanc( rhop, rhob, buo, qall ) use phys_const implicit none real, intent(in) :: rhop(:,:,:) real, intent(in) :: rhob(size(rhop,3)) real, intent(inout) :: buo(size(rhop,1),size(rhop,2),size(rhob)) real, intent(in), optional :: qall(size(rhop,1),size(rhop,2),size(rhob)) integer :: nx, ny, nz integer :: i, j, k nx=size(rhop,1) ny=size(rhop,2) nz=size(rhop,3) if(present(qall))then do k=1,nz do j=1,ny do i=1,nx buo(i,j,k)=g*(rhop(i,j,k)-rhob(k))/(rhob(k)) end do end do end do else do k=1,nz do j=1,ny do i=1,nx buo(i,j,k)=g*(rhop(i,j,k)-rhob(k))/(rhob(k)) end do end do end do end if end subroutine Buoyanc
Subroutine : | |||
types : | character(1), intent(in)
| ||
dt : | real, intent(in)
| ||
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)
| ||
pt(size(x),size(y),size(z)) : | real, intent(in)
| ||
nuth(size(x),size(y),size(z)) : | real, intent(inout)
| ||
nutv(size(x),size(y),size(z)) : | real, intent(inout)
| ||
nuhh(size(x),size(y),size(z)) : | real, intent(inout)
| ||
nuhv(size(x),size(y),size(z)) : | real, intent(inout)
| ||
undef : | real, intent(in), optional
|
subroutine EDC_SMA( types, dt, x, y, z, u, v, w, pt, nuth, nutv, nuhh, nuhv, undef ) use Phys_Const use Statistics use Derivation implicit none ! スマゴリンスキースキームによる渦粘性係数を計算する. character(1), intent(in) :: types ! 乱流の種類の指定. ! 等方性乱流の場合 = 'i', 非等方性乱流の場合 = 'r'. real, intent(in) :: dt ! 時間ステップ [s] 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 方向の速度 [m/s] real, intent(in) :: v(size(x),size(y),size(z)) ! y 方向の速度 [m/s] real, intent(in) :: w(size(x),size(y),size(z)) ! z 方向の速度 [m/s] real, intent(in) :: pt(size(x),size(y),size(z)) ! 基本場の温位 [K] real, intent(inout) :: nuth(size(x),size(y),size(z)) ! 水平渦粘性係数 [m^2/s] real, intent(inout) :: nutv(size(x),size(y),size(z)) ! 鉛直渦粘性係数 [m^2/s] real, intent(inout) :: nuhh(size(x),size(y),size(z)) ! 水平渦拡散係数 [m^2/s] real, intent(inout) :: nuhv(size(x),size(y),size(z)) ! 鉛直渦拡散係数 [m^2/s] real, intent(in), optional :: undef ! 未定義値 real :: BV(size(x),size(y),size(z)) real, parameter :: alpha=1.0e-6 integer :: i, j, k integer :: nx ! 第 1 要素数 integer :: ny ! 第 2 要素数 integer :: nz ! 第 3 要素数 real :: dx(size(x)) ! x 方向の格子間隔 [m] real :: dy(size(y)) ! y 方向の格子間隔 [m] real :: dz(size(z)) ! z 方向の格子間隔 [m] real, dimension(size(x),size(y),size(z)) :: dsh, dsv, def, tmp real :: ck, coea intrinsic :: min, max logical, dimension(size(x),size(y),size(z)) :: undeflag real, parameter :: kmin=0.125 ! nu の最大値は kmin / dt (By CReSS) real, parameter :: csnum=0.21 ! スマゴリンスキー定数 real, parameter :: Praiv=3.0 ! プラントル数 (定数) の逆数. nx=size(x) ny=size(y) nz=size(z) def=0.0 tmp=0.0 BV=0.0 ck=kmin/dt coea=2.0/3.0 undeflag=.false. ! undef オプションがないときは常に false. do i=2,nx-1 dx(i)=0.5*(x(i+1)-x(i-1)) end do do j=2,ny-1 dy(j)=0.5*(y(j+1)-y(j-1)) end do do k=2,nz-1 dz(k)=0.5*(z(k+1)-z(k-1)) end do dx(1)=x(2)-x(1) dx(nx)=x(nx)-x(nx-1) dy(1)=y(2)-y(1) dy(ny)=y(ny)-y(ny-1) dz(1)=z(2)-z(1) dz(nz)=z(nz)-z(nz-1) ! 鉛直方向に ptb のブラントバイサラ振動数を計算する. ! 以降, undef のあり, なしで分ける. if(present(undef))then !$omp parallel default(shared) !$omp do schedule(dynamic) private(i,j) do j=1,ny do i=1,nx call Brunt_Freq( z, pt(i,j,:), BV(i,j,:), undeff=undef ) end do end do !$omp end do !$omp end parallel call deform_tensor( '11', x, y, z, u, v, w, tmp, undef=undef ) call temporary_pow( tmp, undef=undef ) call temporary_add( tmp, def, undef=undef, coe=0.5 ) call deform_tensor( '22', x, y, z, u, v, w, tmp, undef=undef ) call temporary_pow( tmp, undef=undef ) call temporary_add( tmp, def, undef=undef, coe=0.5 ) call deform_tensor( '33', x, y, z, u, v, w, tmp, undef=undef ) call temporary_pow( tmp, undef=undef ) call temporary_add( tmp, def, undef=undef, coe=0.5 ) call deform_tensor( '12', x, y, z, u, v, w, tmp, undef=undef ) call temporary_pow( tmp, undef=undef ) call temporary_add( tmp, def, undef=undef ) call deform_tensor( '13', x, y, z, u, v, w, tmp, undef=undef ) call temporary_pow( tmp, undef=undef ) call temporary_add( tmp, def, undef=undef ) call deform_tensor( '23', x, y, z, u, v, w, tmp, undef=undef ) call temporary_pow( tmp, undef=undef ) call temporary_add( tmp, def, undef=undef ) call div_3d( x, y, z, u, v, w, tmp, undeff=undef ) call temporary_pow( tmp, undef=undef ) call temporary_add( tmp, def, undef=undef, coe=-coea ) !$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(BV(i,j,k)/=undef)then undeflag(i,j,k)=.false. else undeflag(i,j,k)=.true. end if end do end do end do !$omp end do !$omp end parallel else !$omp parallel default(shared) !$omp do schedule(dynamic) private(i,j) do j=1,ny do i=1,nx call Brunt_Freq( z, pt(i,j,:), BV(i,j,:) ) end do end do !$omp end do !$omp end parallel call deform_tensor( '11', x, y, z, u, v, w, tmp ) call temporary_pow( tmp ) call temporary_add( tmp, def, coe=0.5 ) call deform_tensor( '22', x, y, z, u, v, w, tmp ) call temporary_pow( tmp ) call temporary_add( tmp, def, coe=0.5 ) call deform_tensor( '33', x, y, z, u, v, w, tmp ) call temporary_pow( tmp ) call temporary_add( tmp, def, coe=0.5 ) call deform_tensor( '12', x, y, z, u, v, w, tmp ) call temporary_pow( tmp ) call temporary_add( tmp, def ) call deform_tensor( '13', x, y, z, u, v, w, tmp ) call temporary_pow( tmp ) call temporary_add( tmp, def ) call deform_tensor( '23', x, y, z, u, v, w, tmp ) call temporary_pow( tmp ) call temporary_add( tmp, def ) call div_3d( x, y, z, u, v, w, tmp ) call temporary_pow( tmp ) call temporary_add( tmp, def, coe=-coea ) end if !-- 以降, undeflag によって, undef オプションのあるなしに関わらず, !-- 統一の計算式で計算できる. if(types(1:1)=='i')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 dsh(i,j,k)=(dx(i)*dy(j)*dz(k))**(1.0/3.0) dsv(i,j,k)=dsh(i,j,k) end if end do end do end do !$omp end do !$omp end parallel else if(types(1:1)=='r')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 dsh(i,j,k)=sqrt(dx(i)*dy(j)) dsv(i,j,k)=dz(k) end if end do 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 nuth(i,j,k)=((csnum*dsh(i,j,k))**2)*(def(i,j,k)-BV(i,j,k)*Praiv) if(nuth(i,j,k)<0.0)then nuth(i,j,k)=0.0 end if !-- CReSS original nuth(i,j,k)=min(nuth(i,j,k),ck*sqrt(dx(i)*dy(j))) if(types(1:1)=='i')then ! 等方性乱流の場合 nutv(i,j,k)=nuth(i,j,k) nuhh(i,j,k)=Praiv*nuth(i,j,k) nuhv(i,j,k)=nuhh(i,j,k) else if(types(1:1)=='r')then ! 非等方性乱流の場合 nutv(i,j,k)=((csnum*dsv(i,j,k))**2)*(def(i,j,k)-BV(i,j,k)*Praiv) if(nutv(i,j,k)<0.0)then nutv(i,j,k)=0.0 end if !-- CReSS original nutv(i,j,k)=min(nutv(i,j,k),ck*sqrt(dx(i)*dy(j))) nuhh(i,j,k)=Praiv*nuth(i,j,k) nuhv(i,j,k)=Praiv*nutv(i,j,k) end if else nuth(i,j,k)=undef nutv(i,j,k)=undef nuhh(i,j,k)=undef nuhv(i,j,k)=undef end if end do end do end do !$omp end do !$omp end parallel !-- contains subroutine contains subroutine temporary_add( vali, valo, undef, coe ) ! calc. valo = valo + coe*vali implicit none real, intent(in) :: vali(:,:,:) ! adding value real, intent(inout) :: valo(size(vali,1),size(vali,2),size(vali,3)) ! orig. value real, intent(in), optional :: undef real, intent(in), optional :: coe integer :: i, j, k, nx, ny, nz real :: coef nx=size(vali,1) ny=size(vali,2) nz=size(vali,3) if(present(coe))then coef=coe else coef=1.0 end if if(present(undef))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(vali(i,j,k)/=undef.and.valo(i,j,k)/=undef)then valo(i,j,k)=valo(i,j,k)+coef*vali(i,j,k) else valo(i,j,k)=undef 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 valo(i,j,k)=valo(i,j,k)+coef*vali(i,j,k) end do end do end do !$omp end do !$omp end parallel end if end subroutine temporary_add subroutine temporary_pow( vali, undef, coe ) ! calc. vali = vali**coe ! coe = 2 (default) implicit none real, intent(inout) :: vali(:,:,:) ! powering value real, intent(in), optional :: undef integer, intent(in), optional :: coe integer :: i, j, k, nx, ny, nz integer :: coef nx=size(vali,1) ny=size(vali,2) nz=size(vali,3) if(present(coe))then coef=coe else coef=2 end if if(present(undef))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(vali(i,j,k)/=undef)then vali(i,j,k)=vali(i,j,k)**coef 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 vali(i,j,k)=vali(i,j,k)**coef end do end do end do !$omp end do !$omp end parallel end if end subroutine temporary_pow !-- contains subroutine end subroutine EDC_SMA
Subroutine : | |||
types : | character(1), intent(in)
| ||
dt : | real, intent(in)
| ||
x(:) : | real, intent(in)
| ||
y(:) : | real, intent(in)
| ||
z(:) : | real, intent(in)
| ||
pt(size(x),size(y),size(z)) : | real, intent(in)
| ||
tke(size(x),size(y),size(z)) : | real, intent(in)
| ||
nuth(size(x),size(y),size(z)) : | real, intent(inout)
| ||
nutv(size(x),size(y),size(z)) : | real, intent(inout)
| ||
nuhh(size(x),size(y),size(z)) : | real, intent(inout)
| ||
nuhv(size(x),size(y),size(z)) : | real, intent(inout)
| ||
undef : | real, intent(in), optional
|
subroutine EDC_TKE( types, dt, x, y, z, pt, tke, nuth, nutv, nuhh, nuhv, undef ) use Phys_Const use Statistics use Derivation implicit none ! 1.5 次の TKE を用いた渦粘性係数を計算する. character(1), intent(in) :: types ! 乱流の種類の指定. ! 等方性乱流の場合 = 'i', 非等方性乱流の場合 = 'r'. real, intent(in) :: dt ! 時間ステップ [s] real, intent(in) :: x(:) ! x 方向の座標変数 [m] real, intent(in) :: y(:) ! y 方向の座標変数 [m] real, intent(in) :: z(:) ! z 方向の座標変数 [m] real, intent(in) :: pt(size(x),size(y),size(z)) ! 基本場の温位 [K] real, intent(in) :: tke(size(x),size(y),size(z)) ! tke [J/kg] real, intent(inout) :: nuth(size(x),size(y),size(z)) ! 水平渦粘性係数 [m^2/s] real, intent(inout) :: nutv(size(x),size(y),size(z)) ! 鉛直渦粘性係数 [m^2/s] real, intent(inout) :: nuhh(size(x),size(y),size(z)) ! 水平渦拡散係数 [m^2/s] real, intent(inout) :: nuhv(size(x),size(y),size(z)) ! 鉛直渦拡散係数 [m^2/s] real, intent(in), optional :: undef ! 未定義値 real :: BV(size(x),size(y),size(z)) real, parameter :: alpha=1.0e-6 integer :: i, j, k integer :: nx ! 第 1 要素数 integer :: ny ! 第 2 要素数 integer :: nz ! 第 3 要素数 real :: dx(size(x)) ! x 方向の格子間隔 [m] real :: dy(size(y)) ! y 方向の格子間隔 [m] real :: dz(size(z)) ! z 方向の格子間隔 [m] real :: dsh(size(x),size(y),size(z)) real :: dsv(size(x),size(y),size(z)) real :: lh(size(x),size(y),size(z)) real, dimension(size(x),size(y),size(z)) :: lv, ls real :: tmp, ck intrinsic :: min, max logical, dimension(size(x),size(y),size(z)) :: undeflag real, parameter :: kmin=0.125 ! nu の最大値は kmin / dt (By CReSS) nx=size(x) ny=size(y) nz=size(z) ck=kmin/dt undeflag=.false. ! undef オプションがないときは常に false. do i=2,nx-1 dx(i)=0.5*(x(i+1)-x(i-1)) end do do j=2,ny-1 dy(j)=0.5*(y(j+1)-y(j-1)) end do do k=2,nz-1 dz(k)=0.5*(z(k+1)-z(k-1)) end do dx(1)=x(2)-x(1) dx(nx)=x(nx)-x(nx-1) dy(1)=y(2)-y(1) dy(ny)=y(ny)-y(ny-1) dz(1)=z(2)-z(1) dz(nz)=z(nz)-z(nz-1) ! 鉛直方向に ptb のブラントバイサラ振動数を計算する. ! 以降, undef のあり, なしで分ける. if(present(undef))then !$omp parallel default(shared) !$omp do schedule(dynamic) private(i,j) do j=1,ny do i=1,nx call Brunt_Freq( z, pt(i,j,:), BV(i,j,:), undeff=undef ) 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(BV(i,j,k)/=undef.and.tke(i,j,k)/=undef)then ls(i,j,k)=0.76*sqrt(abs(tke(i,j,k)/BV(i,j,k))) undeflag(i,j,k)=.false. else undeflag(i,j,k)=.true. end if end do end do end do !$omp end do !$omp end parallel else !$omp parallel default(shared) !$omp do schedule(dynamic) private(i,j) do j=1,ny do i=1,nx call Brunt_Freq( z, pt(i,j,:), BV(i,j,:) ) end do end do !$omp end do !$omp end parallel !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx ls(i,j,k)=0.76*sqrt(abs(tke(i,j,k)/BV(i,j,k))) end do end do end do !$omp end do !$omp end parallel end if !-- 以降, undeflag によって, undef オプションのあるなしに関わらず, !-- 統一の計算式で計算できる. if(types(1:1)=='i')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 dsh(i,j,k)=(dx(i)*dy(j)*dz(k))**(1.0/3.0) dsv(i,j,k)=dsh(i,j,k) if(BV(i,j,k)>0.0)then lh(i,j,k)=min(dsh(i,j,k),ls(i,j,k)) else lh(i,j,k)=dsh(i,j,k) end if lv(i,j,k)=lh(i,j,k) end if end do end do end do !$omp end do !$omp end parallel else if(types(1:1)=='r')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 dsh(i,j,k)=sqrt(dx(i)*dy(j)) dsv(i,j,k)=dz(k) lh(i,j,k)=dsh(i,j,k) if(BV(i,j,k)>0.0)then lv(i,j,k)=min(dsh(i,j,k),ls(i,j,k)) else lv(i,j,k)=dsv(i,j,k) end if end if end do 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 nuth(i,j,k)=max(0.1*sqrt(tke(i,j,k))*lh(i,j,k),alpha*(dsh(i,j,k)**2)) nutv(i,j,k)=max(0.1*sqrt(tke(i,j,k))*lv(i,j,k),alpha*(dsv(i,j,k)**2)) !-- CReSS original nuth(i,j,k)=min(nuth(i,j,k),ck*(dsh(i,j,k)**2)) nutv(i,j,k)=min(nutv(i,j,k),ck*(dsv(i,j,k)**2)) if(types(1:1)=='i')then ! 等方性乱流の場合 nuhh(i,j,k)=nuth(i,j,k)*(1.0+2.0*(lv(i,j,k)/dsv(i,j,k))) nuhv(i,j,k)=nutv(i,j,k)*(1.0+2.0*(lv(i,j,k)/dsv(i,j,k))) else if(types(1:1)=='r')then ! 非等方性乱流の場合 nuhh(i,j,k)=3.0*nuth(i,j,k) nuhv(i,j,k)=nutv(i,j,k)*(1.0+2.0*(lv(i,j,k)/dsv(i,j,k))) end if else nuth(i,j,k)=undef nutv(i,j,k)=undef nuhh(i,j,k)=undef nuhv(i,j,k)=undef end if end do end do end do !$omp end do !$omp end parallel end subroutine EDC_TKE
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)
| ||
rho(size(x),size(y),size(z)) : | real, intent(in)
| ||
pt(size(x),size(y),size(z)) : | real, intent(in)
| ||
cor(size(x),size(y)) : | real, intent(in)
| ||
pv(size(x),size(y),size(z)) : | real, intent(inout)
| ||
undeff : | real, intent(in), optional | ||
sx(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
sy(size(x),size(y),size(z)) : | real, intent(in), optional
| ||
sz(size(x),size(y),size(z)) : | real, intent(in), optional
|
エルテルのポテンシャル渦度を計算する
subroutine Ertel_PV( x, y, z, u, v, w, rho, pt, cor, pv, undeff, sx, sy, sz ) ! エルテルのポテンシャル渦度を計算する use Thermo_Function use Thermo_Routine use derivation 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 成分 [m/s] real, intent(in) :: v(size(x),size(y),size(z)) ! 速度場の y 成分 [m/s] real, intent(in) :: w(size(x),size(y),size(z)) ! 速度場の z 成分 [m/s] real, intent(in) :: rho(size(x),size(y),size(z)) ! 密度場 [kg/m^3] real, intent(in) :: pt(size(x),size(y),size(z)) ! 温位場 [K] real, intent(in) :: cor(size(x),size(y)) ! コリオリパラメータ [/s] real, intent(inout) :: pv(size(x),size(y),size(z)) ! PV [Km^2/kgs] real, intent(in), optional :: undeff real, intent(in), optional :: sx(size(x),size(y),size(z)) ! スケーリングファクター real, intent(in), optional :: sy(size(x),size(y),size(z)) ! スケーリングファクター real, intent(in), optional :: sz(size(x),size(y),size(z)) ! スケーリングファクター real :: tmp1(size(x),size(y),size(z)) real :: tmp2(size(x),size(y),size(z)) real :: tmp3(size(x),size(y),size(z)) real :: tmp4(size(x),size(y),size(z)) real :: tmp5(size(x),size(y),size(z)) real :: tmp6(size(x),size(y),size(z)) real :: tmp7(size(x),size(y),size(z)) integer :: i, j, k integer :: nx ! 第 1 要素数 integer :: ny ! 第 2 要素数 integer :: nz ! 第 3 要素数 real :: scalex(size(x),size(y),size(z)) real :: scaley(size(x),size(y),size(z)) real :: scalez(size(x),size(y),size(z)) nx=size(x) ny=size(y) nz=size(z) if(present(sx).and.present(sy).and.present(sz))then do k=1,nz do j=1,ny do i=1,nx scalex(i,j,k)=sx(i,j,k) scaley(i,j,k)=sy(i,j,k) scalez(i,j,k)=sz(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 ! 温位の空間勾配を計算. call grad_3d( x, y, z, pt, tmp1, tmp2, tmp3, undeff, scalex, scaley, scalez ) ! 3 次元 rotation を計算. call curl_3d( x, y, z, u, v, w, tmp4, tmp5, tmp6, undeff, scalex, scaley, scalez ) ! omega と grad pt の内積を計算 call dot_prod_3d( tmp4, tmp5, tmp6, tmp1, tmp2, tmp3, tmp7, undeff ) else ! 温位の空間勾配を計算. call grad_3d( x, y, z, pt, tmp1, tmp2, tmp3, hx=scalex, hy=scaley, hz=scalez ) ! 3 次元 rotation を計算. call curl_3d( x, y, z, u, v, w, tmp4, tmp5, tmp6, hx=scalex, hy=scaley, hz=scalez ) ! omega と grad pt の内積を計算 call dot_prod_3d( tmp4, tmp5, tmp6, tmp1, tmp2, tmp3, tmp7 ) end if 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(tmp7(i,j,k)==undeff.or.rho(i,j,k)==undeff)then pv(i,j,k)=undeff else pv(i,j,k)=(tmp7(i,j,k)+cor(i,j)*tmp3(i,j,k))/rho(i,j,k) ! 計算するのは絶対渦度なので, コリオリを足し合わせ. end if end do end do end do !$omp end do !$omp end parallel else !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx pv(i,j,k)=(tmp7(i,j,k)+cor(i,j)*tmp3(i,j,k))/rho(i,j,k) end do end do end do !$omp end do !$omp end parallel end if end subroutine
Subroutine : | |||
z : | real, intent(in)
| ||
z0m : | real, intent(in), dimension(:,:)
| ||
richard : | real, intent(in), dimension(size(z0m,1),size(z0m,2))
| ||
Lo : | real, intent(inout), dimension(size(z0m,1),size(z0m,2))
|
Louis(1980) で提案されている大気の不安定度を考慮したバルク係数の補正係数を計算する関数
subroutine Louis_horizon( z, z0m, richard, Lo ) ! Louis(1980) で提案されている大気の不安定度を考慮したバルク係数の補正係数を計算する関数 use Thermo_Advanced_Function implicit none real, intent(in) :: z ! cm を求める高度 [m] real, intent(in), dimension(:,:) :: z0m ! モデルで計算される粗度高度 [m] real, intent(in), dimension(size(z0m,1),size(z0m,2)) :: richard ! バルクリチャードソン数 real, intent(inout), dimension(size(z0m,1),size(z0m,2)) :: Lo ! 補正係数 real, parameter :: b=5.0, c=5.0 integer :: i, j, nx, ny nx=size(z0m,1) ny=size(z0m,2) !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j) do j=1,ny do i=1,nx Lo(i,j)=Louis( z, z0m(i,j), richard(i,j) ) end do end do !$omp end do !$omp end parallel end subroutine
Subroutine : | |||
za : | real, intent(in)
| ||
pta : | real, intent(in), dimension(:,:)
| ||
ptg : | real, intent(in), dimension(size(pta,1),size(pta,2))
| ||
va : | real, intent(in), dimension(size(pta,1),size(pta,2))
| ||
qva : | real, intent(in), dimension(size(pta,1),size(pta,2))
| ||
qvs : | real, intent(in), dimension(size(pta,1),size(pta,2))
| ||
Ri : | real, intent(inout), dimension(size(pta,1),size(pta,2))
|
バルクリチャードソン数を計算するルーチン
subroutine Rich_horizon( za, pta, ptg, va, qva, qvs, Ri ) ! バルクリチャードソン数を計算するルーチン use Thermo_Advanced_Function implicit none real, intent(in) :: za ! リチャードソン数を計算する高度 [m] real, intent(in), dimension(:,:) :: pta ! za での仮温位 [K] real, intent(in), dimension(size(pta,1),size(pta,2)) :: ptg ! 地表面での温位 [K] real, intent(in), dimension(size(pta,1),size(pta,2)) :: va ! 高度 za での水平風速の絶対値 [m/s] real, intent(in), dimension(size(pta,1),size(pta,2)) :: qva ! za での混合比 [kg/kg] real, intent(in), dimension(size(pta,1),size(pta,2)) :: qvs ! 地表面での飽和混合比 [kg/kg] real, intent(inout), dimension(size(pta,1),size(pta,2)) :: Ri ! 求めるリチャードソン数 integer :: i, j, nx, ny nx=size(pta,1) ny=size(pta,2) !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j) do j=1,ny do i=1,nx Ri(i,j)=Rich( za, pta(i,j), ptg(i,j), va(i,j), qva(i,j), qvs(i,j) ) end do end do !$omp end do !$omp end parallel end subroutine
Subroutine : | |||
z : | real, intent(in)
| ||
z0m : | real, intent(in), dimension(:,:)
| ||
coem : | real, intent(inout), dimension(size(z0m,1),size(z0m,2))
| ||
richard : | real, intent(in), dimension(size(z0m,1),size(z0m,2)), optional
|
運動量に関するバルク係数を計算するルーチン
subroutine cm_horizon( z, z0m, coem, richard ) ! 運動量に関するバルク係数を計算するルーチン use Thermo_Advanced_Function implicit none real, intent(in) :: z ! cm を求める高度 [m] real, intent(in), dimension(:,:) :: z0m ! モデルで計算される粗度高度 [m] real, intent(inout), dimension(size(z0m,1),size(z0m,2)) :: coem ! バルク係数 real, intent(in), dimension(size(z0m,1),size(z0m,2)), optional :: richard ! Louis (1980) のスキームで計算する場合のバルクリチャードソン数 integer :: i, j, nx, ny nx=size(z0m,1) ny=size(z0m,2) if(present(richard))then !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j) do j=1,ny do i=1,nx coem(i,j)=cm( z, z0m(i,j), richard(i,j) ) end do end do !$omp end do !$omp end parallel else !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j) do j=1,ny do i=1,nx coem(i,j)=cm( z, z0m(i,j), richard(i,j) ) end do end do !$omp end do !$omp end parallel end if end subroutine
Subroutine : | |||
cmd : | real, intent(in), dimension(:,:)
| ||
va : | real, intent(in), dimension(size(cmd,1),size(cmd,2))
| ||
velst : | real, intent(inout), dimension(size(cmd,1),size(cmd,2))
|
バルク係数, 速度から摩擦速度 u_* を計算するルーチン
subroutine cmdva_2_ust_horizon( cmd, va, velst ) ! バルク係数, 速度から摩擦速度 u_* を計算するルーチン use Thermo_Advanced_Function implicit none real, intent(in), dimension(:,:) :: cmd ! 高度 za でのバルク係数 real, intent(in), dimension(size(cmd,1),size(cmd,2)) :: va ! 高度 za での水平風の絶対値 [m/s] real, intent(inout), dimension(size(cmd,1),size(cmd,2)) :: velst ! 摩擦速度 [m/s] integer :: i, j, nx, ny nx=size(cmd,1) ny=size(cmd,2) !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j) do j=1,ny do i=1,nx velst(i,j)=cmdva_2_ust( cmd(i,j), va(i,j) ) end do end do !$omp end do !$omp end parallel end subroutine
Subroutine : | |||
taux : | real, intent(in), dimension(:,:)
| ||
tauy : | real, intent(in), dimension(size(taux,1),size(taux,2))
| ||
rho : | real, intent(in), dimension(size(taux,1),size(taux,2))
| ||
velst : | real, intent(inout), dimension(size(taux,1),size(taux,2))
|
風応力のデカルト水平 2 成分とその高度での密度から摩擦速度 u_* を計算するルーチン
subroutine taurho_2_ust_horizon( taux, tauy, rho, velst ) ! 風応力のデカルト水平 2 成分とその高度での密度から摩擦速度 u_* を計算するルーチン use Thermo_Advanced_Function implicit none real, intent(in), dimension(:,:) :: taux ! 基準高度での風応力のデカルト x 成分 [N/m] real, intent(in), dimension(size(taux,1),size(taux,2)) :: tauy ! 基準高度での風応力のデカルト y 成分 [N/m] real, intent(in), dimension(size(taux,1),size(taux,2)) :: rho ! 基準高度での密度 [kg/m^3] real, intent(inout), dimension(size(taux,1),size(taux,2)) :: velst ! 摩擦速度 [m/s] integer :: i, j, nx, ny nx=size(taux,1) ny=size(taux,2) !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j) do j=1,ny do i=1,nx velst(i,j)=taurho_2_ust( (/ taux(i,j), tauy(i,j) /), rho(i,j) ) end do end do !$omp end do !$omp end parallel end subroutine