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 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 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 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 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 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 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 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 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 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