!> \~japanese
!! \brief 表面高度場の時間変化量を計算する(連続の式を解く)手続きを提供するモジュール.
!! \~english
!! \brief This module provides some subroutines which calculate the rate of change of surface height over time
!! (i.e solve the continious equation).
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module continious_equation

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

  !* gtool ****
  !*

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

  !* IGMBaseLib ****
  !

  ! 正二十面体格子上の物理場の管理
  ! Manager of a physical field on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, &
    & get_icgrid, paste_field2D_margin

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

  ! 正二十面体格子上の物理場に作用する微分演算子の提供
  ! Providing some differential operators acting on the physical field on icosahedral grid
  !
  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & divergence_op

  !* IGModel-SW ****
  !*

  ! シミュレーションパラメータの管理
  ! Manager of simulation parameter
  !
  use param_manager, only: &
    & Grav, Omega, alpha
  
  ! OpenMP
  !
  !
  use omp_lib

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public statements
  !
  public :: continious_equation_Init, calc_continious_eq_dhdt

  ! 非公開変数
  ! Pricate statements
  !

  !> \~japanese 質量フラックスの場を管理する構造型 Field_IcGrid2D の変数.
  !! \~english The variable of derived type Field_IcGrid2D to manage the field data of mass flux.
  type(Field_IcGrid2D), save :: mass_flux

contains

!
!> \~japanese
!! \brief continious_equation モジュールの初期化.
!!
!! @param[in,out] icgrid 構造型 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Initializes the continious_equation module.
!!
!! @param[in,out] icgrid The variable of derived type IcGrid2D_FVM.
!!
subroutine continious_equation_Init( &
  & icgrid                           &  ! (inout)
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid
 
  ! 実行文 ; Executable statement
  !

  !  質量フラックスを管理する構造型 Field_IcGrid2D の変数を初期化する.
  call Field_IcGrid2D_Init(self=mass_flux, icgrid=icgrid, name='mass_flux', rank=3)
 
  
end subroutine continious_equation_Init

!
!> \~japanese
!! \brief 連続の式の半離散式の右辺を評価し, 全格子点に対して, 表面高度の時間変化率を計算する.
!!
!! 具体的には, このサブルーチンでは, 以下の Tomita, etal(2001) の式(2) の半離散式の右辺を評価する.
!!
!! \~english
!! \brief Evaluates the right-hand side of the semi-discrete continious equation,
!! and calculates the rate of change of surface height at time level n for all grid points.
!!
!! Specifically, this subroutine evaluates the right-hand side of the following demi-discrete equation
!! as described in Eq (2) in Tomita et al(2001).
!!
!! \~
!! \f[ \left(\DP{h}{t}\right)_n = - \nabla \cdot ( h_n - h_s ) \Dvect{v}_n  \f]
!!
!! \~japanese
!! ここで, \f$h_s\f$は下部境界における地形の高度場であり, 添字の\f$n\f$はタイムレベル n の物理場であることを表す.
!!
!! @param[in,out] DHeightDtN    タイムレベル n の表面高度場の時間変化率を管理する構造型 Field_IcGrid2D の変数.
!! @param[in]     xy_VelN    タイムレベル n の速度場を管理する構造型 Field_IcGrid2D の変数.
!! @param[in]     xy_HeightN タイムレベル n の表面高度場を管理する構造型 Field_IcGrid2D の変数.
!! @param[in]     xy_Htopo   下部境界の高度場を管理する構造型 Field_IcGrid2D の変数.
!! @param[in]     diff_eval 正二十面格子上の物理場の微分演算に対するデータを管理する構造型 Derivate_Field_IcGrid2D の変数.
!!
!! \~english
!! where \f$h_s\f$ is the height field of underlying mountain and
!! the subscript \f$n\f$ denotes the physical field at time level n.
!!
!! @param[in,out] DHeightDtN    The variable of derived type Field_IcGrid2D which manages the rate of change of surface height with time at time level n.
!! @param[in]     xy_VelN       The variable of derived type Field_IcGrid2D which manages the velocity field at time level n.
!! @param[in]     xy_HeightN    The variable of derived type Field_IcGrid2D which manages the surface height field at time level n.
!! @param[in]     xy_Htopo      The variable of derived type Field_IcGrid2D which manages the height field of underlying mountain.
!! @param[in]     diff_eval The variable of derived type IcGridHDiffOptr which manages some data for the differential operators acting on a physical field on icosahedral grid.
!!
subroutine calc_continious_eq_dhdt( &
  & DHeightDtN,                    &  ! (inout)
  & xy_VelN, xy_HeightN, xy_Htopo, &  ! (in)
  & diff_eval                      &  ! (inout)
  & )
 
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(inout) :: DHeightDtN
  type(Field_IcGrid2D), intent(in) :: xy_VelN
  type(Field_IcGrid2D), intent(in) :: xy_HeightN
  type(Field_IcGrid2D), intent(in) :: xy_Htopo
  type(Derivate_Field_IcGrid2D), intent(inout) :: diff_eval


  ! 作業変数
  ! Work variables
  !
  integer :: EMax, EMin
  integer :: rcID,i,j
  real(DP) :: s_pos(3)
  real(DP) :: f
  type(IcGrid2D_FVM), pointer :: icgrid

  ! 実行文 ; Executable statements
  !

  icgrid => get_icgrid(xy_velN)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  ! 質量フラックスを計算する.
  ! ただし, ここではマイナス符号を付加した \f$ -h\Dvect{v}\f$ を計算することにする. 
  ! Calculates the mass flux.
  ! Note that \f$ -h\Dvect{v}\f$ which is mass flux muitiplied by -1 is actually calculated here.

  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do
    do j=EMin, EMax
      do i=EMin, EMax

        mass_flux%val(rcID,i,j,1:3) = &
          & - ( xy_HeightN%val(rcID,i,j,1) - xy_Htopo%val(rcID,i,j, 1) ) * xy_VelN%val(rcID,i,j,1:3)

      end do
    end do
  end do

  ! ノリシロの貼りあわせ
  call paste_field2D_margin(mass_flux)
  
  ! 質量フラックスの収束を計算し, その結果を高度の時間変化率とする.
  ! Calculates the convergence of mass flux,
  ! and the result becomes the rate of change of height wit respect to time.

  call divergence_op(diff_eval, mass_flux, DHeightDtN)
  
end subroutine calc_continious_eq_dhdt

end module continious_equation
