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

  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