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

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

  use gridset, only:  imin,           &! x 方向の配列の下限
    &                 imax,           &! x 方向の配列の上限
    &                 jmin,           &! y 方向の配列の下限
    &                 jmax,           &! y 方向の配列の上限
    &                 kmin,           &! z 方向の配列の下限
    &                 kmax,           &! z 方向の配列の上限
    &                 nx,ny,nz,nf      !
  use basicset, only: xyza_MixRtBasicZ    !基本場の混合比
  use constants,only: CpDry,             &!乾燥成分の比熱
    &                 Grav,              &!重力加速度
    &                 NuV, NuH,          &!数値拡散
    &                 Nu,                &!音波の減衰係数
    &                 Grav                !重力加速度
  use moistset, only: SpcWetSymbol
  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 Dynamic_MixRt_init
  public Tendency_Dynamic_MixRt


contains

  subroutine Dynamic_MixRt_init

    use gtool_historyauto, only: HistoryAutoAddVariable

    implicit none

    integer :: l

    do l = 1, nf
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Adv', &
        & dims=(/'x','y','z','t'/),     &
        & longname='Advection term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='double')
      
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Diff', & 
        & dims=(/'x','y','z','t'/),     &
        & longname='Diffusion term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='double')
    end do

  end subroutine Dynamic_MixRt_init

!!!------------------------------------------------------------------------!!!

  subroutine Tendency_Dynamic_MixRt(    &
    & Time, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, xyza_MixRtBl, xyza_MixRtNl, xyza_Dyn )
    ! 温位の移流拡散の計算 
    ! Advection and diffusion of potential temperature.
    !

    implicit none

    real(DP), intent(in)  :: Time
    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)  :: xyza_MixRtBl(imin:imax, jmin:jmax, kmin:kmax, nf)
    real(DP), intent(in)  :: xyza_MixRtNl(imin:imax, jmin:jmax, kmin:kmax, nf)
    real(DP), intent(out) :: xyza_Dyn(imin:imax, jmin:jmax, kmin:kmax, nf)

    real(DP) :: xyza_Adv(imin:imax, jmin:jmax, kmin:kmax, nf)
    real(DP) :: xyza_Diff(imin:imax, jmin:jmax, kmin:kmax, nf)
    integer :: l

    xyza_Adv  = xyza_AdvScalar(xyza_MixRtNl + xyza_MixRtBasicZ, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
    xyza_Diff = xyza_NumDiffScalar(xyza_MixRtBl) 
    
    xyza_Dyn = xyza_Adv + xyza_Diff 
    
    ! Output
    do l = 1, nf
      call HistoryAutoPut(Time, trim(SpcWetSymbol(l))//'_Adv',  xyza_Adv(1:nx,1:ny,1:nz,l))
      call HistoryAutoPut(Time, trim(SpcWetSymbol(l))//'_Diff', xyza_Diff(1:nx,1:ny,1:nz,l))
    end do
    
  end subroutine Tendency_Dynamic_MixRt

!!!-------------------------------------------------------------------!!!
  function xyza_NumDiffScalar(xyza_Scalar)
    !
    ! x, z 方向に半格子ずれた点での数値拡散項を評価
    !
    use xyz_deriv_module, only: xyz_dx_pyz, pyz_dx_xyz, xyz_dy_xqz, xqz_dy_xyz, xyz_dz_xyr, xyr_dz_xyz
    
    !暗黙の型宣言禁止
    implicit none    

    !変数定義
    real(DP), intent(in) :: xyza_Scalar &
      &                     (imin:imax,jmin:jmax,kmin:kmax,nf)
                                                    !スカラー量
    real(DP)             :: xyza_NumDiffScalar &
      &                     (imin:imax,jmin:jmax,kmin:kmax,nf)
                                                    !水平方向の数値拡散
    integer             :: s
    

    do s = 1, nf
      xyza_NumDiffScalar(:,:,:,s) =   &
        &    NuH * (xyz_dx_pyz(pyz_dx_xyz( xyza_Scalar(:,:,:,s) ))) &
        &  + NuH * (xyz_dy_xqz(xqz_dy_xyz( xyza_Scalar(:,:,:,s) ))) &
        &  + NuV * (xyz_dz_xyr(xyr_dz_xyz( xyza_Scalar(:,:,:,s) ))) 
    end do

  end function xyza_NumDiffScalar


!!!------------------------------------------------------------------------!!!
  function xyza_AdvScalar(xyza_Var, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! x, z 方向に半格子ずれた点における移流を計算
    !
    
    !モジュール読み込み
    use xyz_deriv_c4_module, only: pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz
!    use xyz_deriv_module, only: pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz

    !暗黙の型宣言禁止
    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), intent(in) :: xyza_Var &
      &                    (imin:imax,jmin:jmax,kmin:kmax,nf)
                                                        !スカラー量
    real(DP)             :: xyza_AdvScalar &
      &                    (imin:imax,jmin:jmax,kmin:kmax,nf)
                                                        !スカラー量の水平移流
    integer :: s ! ループ変数

    do s = 1, nf        
      xyza_AdvScalar(:,:,:,s) =                                  &
        & - xyz_avr_pyz(pyz_VelX * pyz_dx_xyz(xyza_Var(:,:,:,s))) &
        & - xyz_avr_xqz(xqz_VelY * xqz_dy_xyz(xyza_Var(:,:,:,s))) &
        & - xyz_avr_xyr(xyr_VelZ * xyr_dz_xyz(xyza_Var(:,:,:,s)))    
    end do

  end function xyza_AdvScalar
 
end module Dynamic_mixrt
