!= Module DynFunc
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu 
! Version::   $Id: dynamic_vely.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_vely
  !
  !陽開放を用いた力学過程の各項の計算モジュール. 
  !具体的には以下の項を計算するための関数を格納する.  
  !  * 移流項
  !  * 浮力項
  !  * 気圧傾度力項
  !

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

  use gridset, only:  imin,           &! x 方向の配列の下限
    &                 imax,           &! x 方向の配列の上限
    &                 jmin,           &! y 方向の配列の下限
    &                 jmax,           &! y 方向の配列の上限
    &                 kmin,           &! z 方向の配列の下限
    &                 kmax             ! z 方向の配列の上限
  use basicset, only: xyz_PotTempBasicZ   !基本場の温位
  use constants,only: CpDry,             &!乾燥成分の比熱
    &                 Grav,              &!重力加速度
    &                 Nu,                &!
    &                 NuH2,              &!
    &                 NuV2
  use axesset, only: &
    &                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_VelY_Long
  public Tendency_Dynamic_VelY_Short

contains

  subroutine Tendency_Dynamic_VelY_Long(                          &
    &            xqz_VelYBl, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, &
    &            xqz_Dyn)
    
    implicit none

    real(DP), intent(in)  :: xqz_VelYBl(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(out) :: xqz_Dyn(imin:imax, jmin:jmax, kmin:kmax)

    real(DP) :: xqz_Adv(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xqz_Diff(imin:imax, jmin:jmax, kmin:kmax)

    if (jmin == jmax) return

    xqz_Adv  = xqz_AdvVelY(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
    xqz_Diff = xqz_NumDiffVelY(xqz_VelYBl)
!    xqz_Turb = xqz_TurbVelY(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl) 
    
    xqz_Dyn = xqz_Adv + xqz_Diff 
    
  end subroutine Tendency_Dynamic_VelY_Long
  
  
  subroutine Tendency_Dynamic_VelY_Short(                          &
    &            pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, xyz_ExnerNs, &
    &            xqz_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(out) :: xqz_Dyn(imin:imax, jmin:jmax, kmin:kmax)

!    real(DP) :: pyz_PGrd(imin:imax, jmin:jmax, kmin:kmax)

    if (jmin == jmax) return
    
    xqz_Dyn = - xqz_GradPi(xyz_ExnerNs, pyz_VelXNs, xqz_VelYNs, xyr_VelZNs)
    
  end subroutine Tendency_Dynamic_VelY_Short
  
!!!------------------------------------------------------------------------!!!
  function xqz_AdvVelY(pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! x 方向に半格子ずれた点における移流を計算
    !
    
    !モジュール読み込み
    use xyz_deriv_c4_module, only: pqz_dx_xqz, xyz_dy_xqz, xqr_dz_xqz
    
    !暗黙の型宣言禁止
    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)             :: xqz_AdvVelY &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                                        !鉛直風の移流
  
!    xr_AdvVelZ = 0.0d0  !初期化
    xqz_AdvVelY =                                                 &
      & - xqz_avr_pqz( pqz_avr_pyz( pyz_VelX ) * pqz_dx_xqz( xqz_VelY ) ) &
      & - xqz_VelY * xqz_avr_xyz( xyz_dy_xqz( xqz_VelY ) ) &
      & - xqz_avr_xqr( xqr_avr_xyr( xyr_VelZ ) * xqr_dz_xqz( xqz_VelY ) )
    
  end function xqz_AdvVelY

 
!!!------------------------------------------------------------------------!!!
  function xqz_GradPi(xyz_Exner, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! y 方向に半格子ずれた点での圧力傾度力項の計算. 
    ! 音波減衰項を含めた形式で定式化してあることに注意.
    !
    
    !モジュール読み込み
    use xyz_deriv_module, only: xyz_dx_pyz, xyz_dy_xqz, xyz_dz_xyr, xqz_dy_xyz
        
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)  :: xyz_Exner &
      &                    (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)              :: xqz_GradPi &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                               !圧力傾度力
    real(DP)              :: xyz_DivVel &
      &                    (imin:imax,jmin:jmax,kmin:kmax)
                                               !速度の収束

    !速度の収束
    xyz_DivVel = xyz_dx_pyz( pyz_VelX ) + xyz_dy_xqz( xqz_VelY ) &
      &          + xyz_dz_xyr( xyr_VelZ )
    
    !圧力傾度
!    xqz_GradPi = 0.0d0
    xqz_GradPi =  &
      & xqz_avr_xyz( CpDry * xyz_PotTempBasicZ )    &
      &   * ( xqz_dy_xyz( xyz_Exner )               &
      &       - xqz_dy_xyz( Nu * xyz_DivVel ) )  
    
  end function xqz_GradPi


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

    !モジュール読み込み
    use xyz_deriv_module, only: xqz_dx_pqz, pqz_dx_xqz, xqz_dy_xyz, xyz_dy_xqz, xqz_dz_xqr, xqr_dz_xqz
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: xqz_VarY &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                    !物理量
    real(DP)             :: xqz_NumDiffVelY &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                    !数値拡散
    
    xqz_NumDiffVelY = &
      &   NuH2 * ( xqz_dx_pqz( pqz_dx_xqz( xqz_VarY ) ) ) &
      & + NuH2 * ( xqz_dy_xyz( xyz_dy_xqz( xqz_VarY ) ) ) &
      & + NuV2 * ( xqz_dz_xqr( xqr_dz_xqz( xqz_VarY ) ) )
    
  end function xqz_NumDiffVelY

  
end module Dynamic_vely

