!= Module DynFunc
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu 
! Version::   $Id: dynamic_velz.f90,v 1.1 2011-04-26 04:08:10 sugiyama Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!モデルの力学過程を計算するために必要となる関数群を束ねたモジュール
!具体的には以下の項を計算するための関数を格納する.  
!  * 移流項
!  * 浮力項
!  * 気圧傾度力項
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * エクスナー関数の空間方向の離散化において, 2 次精度の離散化を陽に利用しているため, 気圧傾度力項の計算プログラムにおいて differentiate_center4 モジュールを指定することはできないので注意.
!
!== Future Plans
!

module dynamic_velz
  !
  !陽開放を用いた力学過程の各項の計算モジュール. 
  !具体的には以下の項を計算するための関数を格納する.  
  !  * 移流項
  !  * 浮力項
  !  * 気圧傾度力項
  !

  !モジュール読み込み
  use dc_types, only: DP

  use gridset, only:  imin,           &! x 方向の配列の下限
    &                 imax,           &! x 方向の配列の上限
    &                 jmin,           &! y 方向の配列の下限
    &                 jmax,           &! y 方向の配列の上限
    &                 kmin,           &! z 方向の配列の下限
    &                 kmax,           &! z 方向の配列の上限
    &                 nf              !
  use basicset, only: xyz_PotTempBasicZ   !基本場の温位
  use constants,only: CpDry,             &!乾燥成分の比熱
    &                 Grav,              &!重力加速度
    &                 MolWtDry,          &!乾燥成分の分子量
    &                 Nu, NuH2, NuV2 
  use moistset, only: CondNum,          &!凝結過程の数
    &                 IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,            &!凝結過程(雲)の配列添え字
    &                 GasNum,           &!気体の数
    &                 IdxG,             &!気体の配列添え字
    &                 IdxNH3,           &!NH3(蒸気)の配列添え字
    &                 IdxH2S,           &!H2S(蒸気)の配列添え字
    &                 SpcWetID,         &!凝縮成分のID
    &                 MolWtWet           !凝縮成分の分子量
  use basicset, only:  &
    &                 xyz_PotTempBasicZ,  &!基本場の温位
    &                 xyz_EffMolWtBasicZ, &!基本場の分子量効果
    &                 xyz_ExnerBasicZ,    &!基本場の無次元圧力
    &                 xyza_MixRtBasicZ     !基本場の混合比
  use axesset, only: x_dx, y_dy, z_dz, &
    &                xyz_avr_pyz, xyr_avr_pyr, xqz_avr_pqz, &
    &                pyz_avr_xyz, pyr_avr_xyr, pqz_avr_xqz, &
    &                xyz_avr_xqz, pyz_avr_pqz, xyr_avr_xqr, &
    &                xqz_avr_xyz, pqz_avr_pyz, xqr_avr_xyr, &
    &                xyz_avr_xyr, pyz_avr_pyr, xqz_avr_xqr, &
    &                xyr_avr_xyz, pyr_avr_pyz, xqr_avr_xqz 

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !移流項計算のための関数を public にする
  public Tendency_Dynamic_VelZ_Long
  public Tendency_Dynamic_VelZ_Long
  public Dynamic_VelZ_init

  real(DP), allocatable :: xyr_MixRtBasicZPerMolWt(:,:,:)
                              !基本場の混合比 / 分子量 の総和
  real(DP), allocatable :: xyr_MixRtBasicZ(:,:,:)
                              !基本場の混合比の総和
  save xyr_MixRtBasicZPerMolWt
  save xyr_MixRtBasicZ 

contains

  subroutine Dynamic_VelZ_init( )

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer               :: s
    real(DP), allocatable :: xyza_MixRtBasicZPerMolWt(:,:,:,:)
                                  !基本場の混合比 / 分子量
    real(DP), allocatable :: xyz_MixRtBasicZPerMolWt(:,:,:)
                              !基本場の混合比 / 分子量 の総和
    real(DP), allocatable :: xyz_MixRtBasicZ(:,:,:)
                              !基本場の混合比の総和

    !-----------------------------------------------------------
    ! 作業配列の初期化. 
    !
    allocate( &
      & xyza_MixRtBasicZPerMolWt(imin:imax,jmin:jmax, kmin:kmax, GasNum),  &
      & xyr_MixRtBasicZPerMolWt(imin:imax,jmin:jmax,kmin:kmax),           &
      & xyz_MixRtBasicZPerMolWt(imin:imax,jmin:jmax,kmin:kmax),           &
      & xyr_MixRtBasicZ(imin:imax,jmin:jmax,kmin:kmax),                   &
      & xyz_MixRtBasicZ(imin:imax,jmin:jmax,kmin:kmax)                    &
      & )
    
    xyza_MixRtBasicZPerMolWt = 0.0d0
    xyr_MixRtBasicZPerMolWt = 0.0d0
    xyz_MixRtBasicZPerMolWt = 0.0d0
    xyr_MixRtBasicZ = 0.0d0
    xyz_MixRtBasicZ = 0.0d0

    do s = 1, GasNum
      xyza_MixRtBasicZPerMolWt(:,:,:,s) = &
        & xyza_MixRtBasicZ(:,:,:,IdxG(s)) / MolWtWet(IdxG(s))
    end do

    xyz_MixRtBasicZPerMolWt = sum(xyza_MixRtBasicZPerMolWt, 4) 
    xyr_MixRtBasicZPerMolWt = xyr_avr_xyz( sum(xyza_MixRtBasicZPerMolWt, 4) )
    xyz_MixRtBasicZ         = sum(xyza_MixRtBasicZ, 4) 
    xyr_MixRtBasicZ         = xyr_avr_xyz( sum(xyza_MixRtBasicZ, 4) )     
    
  end subroutine Dynamic_VelZ_init


  subroutine Tendency_Dynamic_VelZ_Long(                         &
    &            xyr_VelZBl, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, &
    &            xyz_PotTempNl, xyza_MixRtNl,                    &
    &            sw_moist, xyr_Dyn)
    
    implicit none

    real(DP), intent(in)  :: xyr_VelZBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: pyz_VelXNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyz_PotTempNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyza_MixRtNl(imin:imax, jmin:jmax, kmin:kmax, nf)
    logical,  intent(in)  :: sw_moist
    real(DP), intent(out) :: xyr_Dyn(imin:imax, jmin:jmax, kmin:kmax)
    
    real(DP) :: xyr_Adv(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyr_Diff(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyr_BuoyT(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyr_BuoyM(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyr_BuoyD(imin:imax, jmin:jmax, kmin:kmax)
    
    xyr_Adv  = xyr_AdvVelZ(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
    xyr_Diff = xyr_NumDiffVelZ(xyr_VelZBl)
!    xyr_Turb = xyr_TurbVelZ(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl) 

    xyr_BuoyT = xyr_Buoy(xyz_PotTempNl)
    if (sw_moist) then 
      xyr_BuoyM = xyr_BuoyMolWt(xyza_MixRtNl)
      xyr_BuoyD = xyr_BuoyDrag(xyza_MixRtNl) 
    else
      xyr_BuoyM = 0.0d0
      xyr_BuoyD = 0.0d0
    end if
    xyr_Dyn = xyr_Adv + xyr_Diff + xyr_BuoyT + xyr_BuoyM + xyr_BuoyD

  end subroutine Tendency_Dynamic_VelZ_Long

  
  subroutine Tendency_Dynamic_VelZ_Short(              &
    &            pyz_VelXNs, xqz_VelYNs, xyr_VelZNs,  &
    &            xyz_ExnerNs, xyz_ExnerAs,            &
    &            xyr_Dyn)
    
    implicit none

    real(DP), intent(in)  :: pyz_VelXNs(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNs(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNs(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyz_ExnerNs(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyz_ExnerAs(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(out) :: xyr_Dyn(imin:imax, jmin:jmax, kmin:kmax)

    xyr_Dyn = - xyr_GradPi(xyz_ExnerAs,xyz_ExnerNs, pyz_VelXNs, xqz_VelYNs, xyr_VelZNs)

  end subroutine Tendency_Dynamic_VelZ_Short
  

!!!------------------------------------------------------------------------!!!
  function xyr_AdvVelZ(pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! x 方向に半格子ずれた点における移流を計算
    !
    
    !モジュール読み込み
    use xyz_deriv_c4_module, only: pyr_dx_xyr, xqr_dy_xyr, xyz_dz_xyr
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: pyz_VelX &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                                        !水平風速
    real(DP), intent(in) :: xqz_VelY &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                                        !水平風速
    real(DP), intent(in) :: xyr_VelZ &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                                        !鉛直風速
    real(DP)             :: xyr_AdvVelZ &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                                        !鉛直風の移流
  
!    xr_AdvVelZ = 0.0d0  !初期化
    xyr_AdvVelZ =                                                 &
      & - xyr_avr_pyr( pyr_avr_pyz( pyz_VelX ) * pyr_dx_xyr( xyr_VelZ ) ) &
      & - xyr_avr_xqr( xqr_avr_xqz( xqz_VelY ) * xqr_dy_xyr( xyr_VelZ ) ) &
      & - xyr_VelZ * xyr_avr_xyz( xyz_dz_xyr( xyr_VelZ ) )
    
  end function xyr_AdvVelZ

 
!!!------------------------------------------------------------------------!!!
  function xyr_Buoy(xyz_PotTemp)
    !
    ! 鉛直方向の運動方程式に現れる浮力項を計算

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)  :: xyz_PotTemp(imin:imax,jmin:jmax,kmin:kmax)
                                               !温位擾乱
    real(DP)              :: xyr_Buoy &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                               !浮力項

!    !初期化
!    xr_Buoy = 0.0d0

    !浮力項の計算
    xyr_Buoy = Grav * xyr_avr_xyz(xyz_PotTemp / xyz_PotTempBasicZ)

  end function xyr_Buoy

!!!------------------------------------------------------------------------!!!
  function xyr_BuoyMolWt(xyza_MixRt)
    !
    ! 鉛直方向の運動方程式に現れる浮力項のうち, 
    ! 分子量の効果だけを求める
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)  :: xyza_MixRt(imin:imax,jmin:jmax,kmin:kmax, nf)
                                               !凝縮成分の混合比
    real(DP) :: xyr_BuoyMolWt(imin:imax,jmin:jmax,kmin:kmax)
                                               !浮力項(分子量効果)
    real(DP) :: xyza_MixRtPerMolWt(imin:imax,jmin:jmax, kmin:kmax, GasNum)
                                               !混合比/分子量
    integer              :: s
    
    !初期化
    xyr_BuoyMolWt           = 0.0d0
    
    !作業配列の初期化. 気体のみ利用
    do s = 1, GasNum
      xyza_MixRtPerMolWt(:,:,:,s) = xyza_MixRt(:,:,:,IdxG(s)) / MolWtWet(IdxG(s))
    end do
    
    !浮力項の計算
    xyr_BuoyMolWt =                                               &
      & + Grav * xyr_avr_xyz( sum(xyza_MixRtPerMolWt, 4) )        &
      &    / ( 1.0d0 / MolWtDry + xyr_MixRtBasicZPerMolWt )       &
      & - Grav * xyr_avr_xyz( sum(xyza_MixRt(:,:,:,1:GasNum), 4) )&
      &    / ( 1.0d0 + xyr_MixRtBasicZ )

  end function xyr_BuoyMolWt

!!!------------------------------------------------------------------------!!!
  function xyr_BuoyDrag(xyza_MixRt) 
    !
    ! 鉛直方向の運動方程式に現れる浮力項のうち, 
    ! 引きづりのの効果だけを求める
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)  :: xyza_MixRt(imin:imax,jmin:jmax,kmin:kmax, nf)
                                               !凝縮成分の混合比
    real(DP) :: xyr_BuoyDrag(imin:imax,jmin:jmax,kmin:kmax)
                                               !浮力項(分子量効果)
    !初期化
    xyr_BuoyDrag           = 0.0d0

    !浮力項の計算
    xyr_BuoyDrag =                                                     &
      & - Grav * xyr_avr_xyz( sum(xyza_MixRt(:,:,:,GasNum+1:nf), 4) ) &
      &    / ( 1.0d0 + xyr_MixRtBasicZ )

  end function xyr_BuoyDrag

!!!--------------------------------------------------------------------!!!
  function xyr_GradPi(xyz_ExnerA, xyz_ExnerN, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! x 方向に半格子ずれた点での圧力傾度力項の計算. 
    ! クランク・ニコルソン法を用いるために, 時刻 τ と τ+Δτでの
    ! エクスナー関数の値を引数として取る.
    ! 音波減衰項を含めた形式で定式化してあることに注意.
    !

    ! モジュールの読み込み
    use xyz_deriv_module, only: xyz_dx_pyz, xyz_dy_xqz, xyz_dz_xyr, &
      &                         xyr_dz_xyz, xyr_dz_xyz, xyr_dz_xyz


    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)  :: xyz_ExnerA &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !エクスナー関数の擾乱
    real(DP), intent(in)  :: xyz_ExnerN &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !エクスナー関数の擾乱
    real(DP), intent(in)  :: pyz_VelX &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !水平速度
    real(DP), intent(in)  :: xqz_VelY &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !水平速度
    real(DP), intent(in)  :: xyr_VelZ &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !鉛直速度
    real(DP)              :: xyr_GradPi &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !圧力傾度力
    real(DP)              :: xyz_DivVel &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !速度の収束
    real(DP)               :: beta  = 5.0d-1   !クランクニコルソン法なら 0.5
                                               !完全陰解法なら 1

    !速度の収束
    xyz_DivVel =  xyz_dx_pyz( pyz_VelX )   &
      &           + xyz_dy_xqz( xqz_VelY ) &
      &           + xyz_dz_xyr( xyr_VelZ )
    
    !速度 w の圧力勾配
!    xyr_GradPi = 0.0d0
    xyr_GradPi =   &
      & xyr_avr_xyz(CpDry * xyz_PotTempBasicZ )   &
      &   * (                                                      &
      &       beta * xyr_dz_xyz( xyz_ExnerA )                         &
      &       + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerN )             &
      &       - xyr_dz_xyz( Nu * xyz_DivVel )                  &
      &     )                                              
    
  end function xyr_GradPi
  

!!!-------------------------------------------------------------------!!!  
  function xyr_NumDiffVelZ(xyr_VarZ)
    !
    ! x 方向に半格子ずれた点での数値拡散項を評価
    !

    ! モジュールの読み込み
    use xyz_deriv_module, only: xyr_dx_pyr, pyr_dx_xyr, xyr_dy_xqr, xqr_dy_xyr, xyr_dz_xyz, xyz_dz_xyr
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: xyr_VarZ(imin:imax,jmin:jmax,kmin:kmax)
                                                    !物理量
    real(DP)             :: xyr_NumDiffVelZ(imin:imax,jmin:jmax,kmin:kmax)
                                                    !数値拡散

    xyr_NumDiffVelZ = &
      &   NuH2 * ( xyr_dx_pyr( pyr_dx_xyr( xyr_VarZ ) ) ) &
      & + NuH2 * ( xyr_dy_xqr( xqr_dy_xyr( xyr_VarZ ) ) ) & 
      & + NuV2 * ( xyr_dz_xyz( xyz_dz_xyr( xyr_VarZ ) ) )
    
  end function xyr_NumDiffVelZ
 
end module Dynamic_velz

