!> \~japanese
!! \brief 系の保存量に対する解析を行うための手続きを提供するモジュール.
!! \~english
!! \brief This module provides some subroutines to perform the analysis for the conservative quantities of a system.
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module conservation_analysis

  ! モジュール引用 ; Use statement
  !

  !* gtool **
  !*

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP ! 倍精度実数型. Double precision.

  !* IGMBaseLib **
  !*

  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: &
    & dot

  ! 正二十面体格子の管理
  ! Manager of icosahedral grid
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, &
    & RC_REGIONS_NUM

  ! 正二十面格子上の物理場データの管理
  ! Manager of the physical field data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & get_IcGrid, global_mean

  ! 正二十面体格子上の物理場に作用する微分演算
  ! Some differntial operators acting on a physical field on icosahedral grid
  !
  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & vertical_curl_op

  !* IGModel-SW **
  !*

  ! シミュレーションパラメータの管理
  ! Manager of the simulation parameters
  !
  use param_manager, only: Grav

  ! 宣言文 ; Declaration statement
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: init_conservation_analysis, calc_total_energy, calc_potential_enstrophy

  ! 非公開変数
  ! Private variables

  !> \~japanese 初期状態におけるポテンシャルエネルギー.
  !! \~english The potential energy in the initial state.
  real(DP) :: Ep0

contains

!
!> \~japanese
!! \brief conservation_analysis モジュールの初期化を行う.
!!
!! 基準となるポテンシャルエネルギーを決定するために, 初期のポテンシャルエネルギーを計算する.
!!
!! @param[in] init_h 初期の流体表面の高度場を管理する構造体 Field_IcGrid2D の変数.
!! @param[in] hs     下部地形の高度場を管理する構造体 Field_IcGrid2D の変数.
!!
!! \~english
!! \brief Initializes the conservation_analysis module.
!!
!! 基準となるポテンシャルエネルギーを決定するために, 初期のポテンシャルエネルギーを計算する. 
!!
!! @param[in] init_h The variable of derived type Field_IcGrid2D which manages the initial height field of fluid surface.
!! @param[in] hs     The variable of derived type Field_IcGrid2D which manages the initial height field of underlying mountain.
!!
subroutine init_conservation_analysis(init_h, hs)
  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(in) :: init_h
  type(Field_IcGrid2D), intent(in) :: hs

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid2D) :: tmp_Ep_field
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: rcID, i, j
  integer :: EMax, EMin

  ! 実行文 ; Executable statements
  !
  icgrid => get_IcGrid(init_h)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  call Field_IcGrid2D_Init(tmp_Ep_field, icgrid, 'tmp_PE_field', 1)

  do rcID=1, RC_REGIONS_NUM
    do j= EMin, EMax
      do i= EMin, EMax
        tmp_Ep_field%val(rcID,i,j,1) = &
          & 0.5d0 * Grav * &
          &   ( init_h%val(rcID,i,j,1)**2 - hs%val(rcID,i,j,1)**2 )

      end do
    end do
  end do

  Ep0 = global_mean(tmp_Ep_field)
  !write(*,*) 'Ep0=',Ep0
  call Field_IcGrid2D_Final(tmp_Ep_field)

end subroutine init_conservation_analysis

!
!> \~japanese
!! \brief 系の全エネルギー((正二十面格子の全表面積で規格化))を計算する.
!!
!! @param[in] h  流体表面の高度場を管理する構造体 Field_IcGrid2D の変数.
!! @param[in] hs 下部地形のの高度場を管理する構造体 Field_IcGrid2D の変数.
!! @param[in] v  速度場を管理する構造体 Field_IcGrid2D の変数.
!!
!! \~english
!! \brief Calculates the total energy (normalized by the total surface area of icosahedral grid) in a system.
!!
!! @param[in] h  The variable of derived type Field_IcGrid2D which manages initial height field of fluid surface.
!! @param[in] hs The variable of derived type Field_IcGrid2D which manages the height field of underlying mountain.
!! @param[in] v  The variable of derived type Field_IcGrid2D which manages the velocity field.
!!
function calc_total_energy(h, hs, v) result(total_energy)
  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(in) :: h
  type(Field_IcGrid2D), intent(in) :: hs
  type(Field_IcGrid2D), intent(in) :: v
  real(DP) total_energy

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid2D) :: tmp_TE_field
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: rcID, i, j
  integer :: EMax, EMin

  ! 実行文
  !
  icgrid => get_IcGrid(h)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  call Field_IcGrid2D_Init(tmp_TE_field, icgrid, 'tmp_TE_field', 1)

  do rcID=1, RC_REGIONS_NUM
    do j= EMin, EMax
      do i= EMin, EMax
        tmp_TE_field%val(rcID,i,j,1) = &
          &   0.5d0 * ( h%val(rcID,i,j,1) - hs%val(rcID,i,j,1) ) * dot(v%val(rcID,i,j,:), v%val(rcID,i,j,:)) &
          & + 0.5d0 * Grav * ( h%val(rcID,i,j,1)**2 - hs%val(rcID,i,j,1)**2 ) &
          & - Ep0
      end do
    end do
  end do

  total_energy = global_mean(tmp_TE_field)
  call Field_IcGrid2D_Final(tmp_TE_field)

end function calc_total_energy

!
!> \~japanese
!! \brief 系のポテンシャルエンストロフィー(正二十面格子の全表面積で規格化)を計算する.
!!
!! @param[in] h  流体表面の高度場を管理する構造体 Field_IcGrid2D の変数.
!! @param[in] hs 下部地形のの高度場を管理する構造体 Field_IcGrid2D の変数.
!! @param[in] v  速度場を管理する構造体 Field_IcGrid2D の変数.
!! @param[in] f  惑星渦度場を管理する構造体 Field_IcGrid2D の変数.
!! @param[in,out] diff_eval 正二十面格子上の物理場の微分演算に対するデータを管理する構造型 Derivate_Field_IcGrid2D の変数.
!!
!! \~english
!! \brief Calculates the potential energy (normalized by the total surface area of icosahedral grid) in a system.
!!
!! @param[in] h  The variable of derived type Field_IcGrid2D which manages initial height field of fluid surface.
!! @param[in] hs The variable of derived type Field_IcGrid2D which manages the height field of underlying mountain.
!! @param[in] v  The variable of derived type Field_IcGrid2D which manages the velocity field.
!! @param[in] f  The variable of derived type Field_IcGrid2D which manages the planetary vorticity field.
!! @param[in,out] diff_eval The variable of derived type IcGridHDiffOptr which manages some data for the differential operators acting on a physical field on icosahedral grid.
!!
function calc_potential_enstrophy(h, hs, v, f, diff_eval) result(potential_enstrophy)
  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(in) :: h
  type(Field_IcGrid2D), intent(in) :: hs
  type(Field_IcGrid2D), intent(in) :: v
  type(Field_IcGrid2D), intent(in) :: f
  type(Derivate_Field_IcGrid2D), intent(inout) :: diff_eval
  real(DP) potential_enstrophy

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid2D) :: tmp_PE_field
  type(Field_IcGrid2D) :: vorticity
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: rcID, i, j
  integer :: EMax, EMin

  ! 実行文
  !
  icgrid => get_IcGrid(h)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  call Field_IcGrid2D_Init(tmp_PE_field, icgrid, 'tmp_PE_field', 1)
  call Field_IcGrid2D_Init(vorticity, icgrid, 'vorticity', 1)

  ! 速度場に回転演算子を作用させて, 渦度場を計算する.
  !
  call vertical_curl_op(diff_eval, v, vorticity )

  do rcID=1, RC_REGIONS_NUM
    do j= EMin, EMax
      do i= EMin, EMax
        tmp_PE_field%val(rcID,i,j,1) = &
          & 0.5d0 * ( vorticity%val(rcID,i,j,1) + f%val(rcID,i,j,1) )**2 &
          &   / ( h%val(rcID,i,j,1) - hs%val(rcID,i,j,1) )
      end do
    end do
  end do

  potential_enstrophy = global_mean(tmp_PE_field)
  call Field_IcGrid2D_Final(tmp_PE_field)

end function calc_potential_enstrophy

end module conservation_analysis
