!= Module EnergyMonit
!
! Authors::   ODAKA Masatsugu 
! Version::   $Id: energymonit.f90,v 1.1 2012-12-06 08:02:21 odakker Exp $
! Tag Name::  $Name: arare5-20130129 $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!保存量の計算と出力を行うためのモジュール
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!
!== Future Plans
!
module EnergyMonit
  ! 保存量の計算と出力を行うためのモジュール

  !モジュール読み込み
  use dc_types,   only: DP, STRING
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify
  use gtool_historyauto, &
    &             only: HistoryAutoAddVariable, HistoryAutoPut

  use mpi_wrapper,only: myrank
  use gridset,  only: imin,         & !x 方向の配列の下限
    &                 imax,         & !x 方向の配列の上限
    &                 jmin,         & !y 方向の配列の下限
    &                 jmax,         & !y 方向の配列の上限
    &                 kmin,         & !z 方向の配列の下限
    &                 kmax,         & !z 方向の配列の上限
    &                 nx, ny, nz
  use axesset,  only: xyz_dx, xyz_dy, xyz_dz, &           !格子点間隔
    &                 xyz_Z,                  &           !座標
    &                 xmin, xmax, ymin, ymax, zmin, zmax,&!計算領域
    &                 xyz_avr_pyz, xyz_avr_xqz, xyz_avr_xyr
  use constants,only: CpDry,        & !定圧比熱
    &                 CvDry,        & !定積比熱
    &                 GasRDry,      & !期待定数
    &                 PressSfc,     & !地表気圧
    &                 Grav            !重力加速度
  use basicset, only: xyz_ExnerBZ,  & !基本場の圧力関数
    &                 xyz_PressBZ,  & !基本場の気圧
    &                 xyz_PTempBZ,  & !基本場の温位
    &                 xyz_TempBZ,   & !基本場の温度
    &                 xyzf_QMixBZ,  & !基本場の混合非
    &                 xyz_DensBZ,   & !基本場の密度
    &                 xyz_VelSoundBZ  !音速

  use timeset, only:  TimeN
  use namelist_util, only: namelist_filename

  ! 変数定義
  character(*), parameter:: module_name = 'energymonit'
                              ! モジュールの名称.

                              ! Module name
  

contains
!!!------------------------------------------------------------------------!!!
  subroutine EnergyMonit_init
    !
    !出力変数の設定を行う. 
    !

    !暗黙の型宣言禁止
    implicit none

    call HistoryAutoAddVariable(   &
      & varname='DensDev',         &
      & dims=(/'t'/),  &
      & longname='Total density deviation', &
      & units='kg',            &
      & xtype='float')

    call HistoryAutoAddVariable(   &
      & varname='KinEnrgy',        &
      & dims=(/'t'/),  &
      & longname='Total kinetic energy', &
      & units='J',             &
      & xtype='float')

    call HistoryAutoAddVariable(   &
      & varname='ElstEnrgy',       &
      & dims=(/'t'/),  &
      & longname='Total elastic energy', &
      & units='J',             &
      & xtype='float')

    call HistoryAutoAddVariable(   &
      & varname='PotEnrgy',        &
      & dims=(/'t'/),  &
      & longname='Total potential energy', &
      & units='J',             &
      & xtype='float')

  end subroutine EnergyMonit_init
!!!------------------------------------------------------------------------!!!
  subroutine EnergyMonit_exec( &
    & pyz_VelX, xqz_VelY, xyr_VelZ, xyz_Exner, xyz_PTemp &
    )
    !
    !保存量の計算と出力を行う. 計算する量は
    !質量, 運動エネルギー, 弾性エネルギー, ポテンシャルエネルギー
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP),intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)   !速度(x成分)
    real(DP),intent(in) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)   !速度(y成分)
    real(DP),intent(in) :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax)   !速度(z成分)
    real(DP),intent(in) :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)  !圧力関数
    real(DP),intent(in) :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)  !温位

    real(DP) :: xyz_DensDev(imin:imax,jmin:jmax,kmin:kmax) 
                                     !密度(偏差)
    real(DP) :: xyz_KinEnrgyDens(imin:imax,jmin:jmax,kmin:kmax) 
                                     !運動エネルギー密度
    real(DP) :: xyz_ElstEnrgyDens(imin:imax,jmin:jmax,kmin:kmax)
                                     !弾性エネルギー密度
    real(DP) :: xyz_PotEnrgyDens(imin:imax,jmin:jmax,kmin:kmax)
                                     !ポテンシャルエネルギー密度

    real(DP) :: MassTotalBZ          !全質量(基本場)
    real(DP) :: MassTotalDev         !全質量(偏差)
    real(DP) :: KinEnrgyTotal        !全運動エネルギー
    real(DP) :: ElstEnrgyTotal       !全弾性エネルギー
    real(DP) :: PotEnrgyTotal        !全ポテンシャルエネルギー

! 追加すべき量(要検討)
! * 全運動量
! * エンタルピー
! * トレーサーの総量
! * 潜熱

    ! 各格子点の大気質量密度偏差
    ! 
    xyz_DensDev = PressSfc * ( & 
      & (xyz_ExnerBZ + xyz_Exner)**( CvDry/GasRDry ) &
      &          - (xyz_ExnerBZ )**( CvDry/GasRDry ) &
      & )/xyz_PTempBZ 

    ! 領域全体の大気質量
    !
    MassTotalBZ  = (xyz_PressBz(1,1,1) - xyz_PressBz(1,1,nz))/Grav &
      &              *(xmax - xmin)*(ymax - ymin)
    MassTotalDev = sum(xyz_dx(1:nx,1:ny,1:nz)* &
      &                xyz_dy(1:nx,1:ny,1:nz)* &
      &                xyz_dz(1:nx,1:ny,1:nz)* &
      &                xyz_DensDev(1:nx,1:ny,1:nz))

   
    ! 各格子点の運動エネルギー密度
    ! * 密度は基本場の値で評価
    ! 
    xyz_KinEnrgyDens = 0.5d0 * xyz_DensBZ * ( &
      &             xyz_avr_pyz(pyz_VelX)**2  &
      &           + xyz_avr_xqz(xqz_VelY)**2  &
      &           + xyz_avr_xyr(xyr_VelZ)**2  &
      & )

    !各格子点の弾性エネルギー密度
    ! * 密度は基本場の値で評価
    !  
    xyz_ElstEnrgyDens = 0.5d0 * xyz_DensBZ     &
      &  * (CpDry * xyz_PTempBZ * xyz_Exner / xyz_VelSoundBZ)**2 
      
    ! 領域全体の運動/弾性エネルギー
    !
    KinEnrgyTotal = sum(xyz_dx(1:nx,1:ny,1:nz)* &
      &                  xyz_dy(1:nx,1:ny,1:nz)* &
      &                  xyz_dz(1:nx,1:ny,1:nz)* &
      &                  xyz_KinEnrgyDens(1:nx,1:ny,1:nz))

    ElstEnrgyTotal = sum(xyz_dx(1:nx,1:ny,1:nz)* &
      &                   xyz_dy(1:nx,1:ny,1:nz)* &
      &                   xyz_dz(1:nx,1:ny,1:nz)* &
      &                   xyz_ElstEnrgyDens(1:nx,1:ny,1:nz))


    !各格子点のポテンシャルエネルギー密度
    ! 
    xyz_PotEnrgyDens = &
      &  - Grav * xyz_DensBZ * xyz_PTemp * xyz_Z / xyz_PTempBZ

    ! 領域全体のポテンシャルエネルギー
    !
    PotEnrgyTotal = sum(xyz_dx(1:nx,1:ny,1:nz)* &
      &                  xyz_dy(1:nx,1:ny,1:nz)* &
      &                  xyz_dz(1:nx,1:ny,1:nz)* &
      &                  xyz_PotEnrgyDens(1:nx,1:ny,1:nz))


    !ファイルへの出力
    !
    call HistoryAutoPut(TimeN, 'DensDev',   MassTotalDev  )
    call HistoryAutoPut(TimeN, 'KinEnrgy',  KinEnrgyTotal )
    call HistoryAutoPut(TimeN, 'ElstEnrgy', ElstEnrgyTotal)
    call HistoryAutoPut(TimeN, 'PotEnrgy',  PotEnrgyTotal )

  end subroutine EnergyMonit_exec

end module EnergyMonit
