!> \~japanese
!! \brief 速度場の時間変化率を計算する(運動方程式を解く)手続きを提供するモジュール.
!!
!! \~english
!! \brief This module provides some subroutines which calculate the rate of change of velocity over time(i.e solve the equation of motion).
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module motion_equation
  
  ! モジュール引用 ; Use statements
  !

  !* gtool ****
  !*

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

  !* IGMBaseLib ****
  !*

  ! 数学
  ! Mathematica
  !
  use igmcore_math, only: PI

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

  ! 座標変換
  ! Coordinate conversion
  !
  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos, geo_to_orth_vec


  ! 正二十面体格子上の物理場の管理
  ! 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, check_pole, &
    & RC_REGIONS_NUM, NOT_POLE_FLAG

  ! 正二十面体格子上の物理場に作用する微分演算を提供するクラス
  ! Class providing some differential operators acting on a physcial field on icosahedral grid.  
  !
  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & vertical_curl_op, gradient_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 :: motion_equation_Init, calc_motion_eq_dvdt

  ! 非公開変数
  ! Private statements
  !

  !> \~japanese 相対渦度場データを管理する構造型 Field_IcGrid2D の変数.
  !! \~english The variable of derived type Field_IcGrid2D to manage the field data of relative vorticity.
  type(Field_IcGrid2D), save :: xy_Zeta

  !> \~japanese 力学的エネルギーの物理場データを管理する構造型 Field_IcGrid2D の変数.
  !! \~english The variable of derived type Field_IcGrid2D to manage the field data of mechanical energy.
  type(Field_IcGrid2D), save :: mecha_energy

  !> \~japanese ジオポテンシャルと運動エネルギーの和の勾配場を管理する構造型 Field_IcGrid2D の変数.
  !! \~english The variable of derived type Field_IcGrid2D to manage the gradient field data of the sum of geopotential and kinetic energy.
  type(Field_IcGrid2D), save :: grad_mecha_energy

contains

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

  ! 渦度場, 運動エネルギーの場, 運動エネルギーの勾配場を管理する, 構造型 Field_IcGrid2D の変数を初期化する.
  !

  call Field_IcGrid2D_Init( &
    & self=xy_Zeta, &
    & icgrid=icgrid, name='relative vorticity field', rank=1 )

  call Field_IcGrid2D_Init( &
    & self=mecha_energy, &
    & icgrid=icgrid, name='mechanical_energy_field', rank=1 )

  call Field_IcGrid2D_Init( &
    & self=grad_mecha_energy, &
    & icgrid=icgrid, name='gradient of mechanical energy field', rank=3 )
  
end subroutine motion_equation_Init

!
!> \~japanese
!! \brief 運動方程式の半離散式の右辺を評価し, 全格子点に対してタイムレベル n の速度変化率を計算する.
!!
!! 具体的には, このサブルーチンでは, 以下の Tomita, etal(2001) の式(6) の半離散式の右辺を評価する.
!!
!! \~english
!! \brief Evaluates the right-hand side of the semi-discrete equation of motion,
!! and calculates the rate of change of velocity at time level n for all grid points.
!!
!! Specifically, this subroutine evaluates the right-hand side of the following semi-discrete equation
!! as described in Eq (6) in Tomita et al(2001).
!!
!! \~
!! \f[ \left( \DP{\Dvect{v}}{t} \right)_n =
!!     - ( \zeta_n + f ) \Dvect{k} \times \Dvect{v}_n
!!     - \nabla \left( gh + \dfrac{\Dvect{v} \cdot \Dvect{v}}{2} \right) \f]
!!
!! \~japanese
!! ここで, 添字の\f$n\f$はタイムレベル n の物理場であることを表す.
!!
!! @param[in,out] DVelDtN    タイムレベル n の速度場の時間変化率を管理する構造型 Field_IcGrid2D の変数.
!! @param[in]     xy_VelN    タイムレベル n の速度場を管理する構造型 Field_IcGrid2D の変数.
!! @param[in]     xy_HeightN タイムレベル n の表面高度場を管理する構造型 Field_IcGrid2D の変数.
!! @param[in]     xy_Coli
!! @param[in]     diff_eval 正二十面格子上の物理場の微分演算に対するデータを管理する構造型 Derivate_Field_IcGrid2D の変数.
!!
!! \~english
!! where the subscript \f$n\f$ denotes the physical field at time level n.
!!
!! @param[in,out] DVelDtN    The variable of derived type Field_IcGrid2D which manages the rate of change of velocity field 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_Coli
!! @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_motion_eq_dvdt( &
  & DVelDtN,                       &  ! (inout)
  & xy_VelN, xy_HeightN, xy_Coli,  &  ! (in)
  & diff_eval                      &  ! (inout)
  & )
 
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(inout) :: DVelDtN
  type(Field_IcGrid2D), intent(in) :: xy_VelN
  type(Field_IcGrid2D), intent(in) :: xy_HeightN
  type(Field_IcGrid2D), intent(in) :: xy_Coli
  type(Derivate_Field_IcGrid2D), intent(inout) :: diff_eval


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

  ! 実行文 ; Executable statements
  !

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

  
  ! 速度場に回転演算子を作用させて, 相対渦度場を計算する. 
  ! Curl operator acts on the velocity field to calculate the relative vorticity field.  
  call vertical_curl_op(diff_eval, xy_VelN, xy_Zeta )
  
  ! ジオポテンシャルと運動エネルギーの和の空間勾配場を計算する.
  ! Calculates the spatial gradient geopotential and kinetic energy.
  !

  ! まず, 各格子点において力学的エネルギーを計算する.
  ! First, calculates the mechanical energy at each grid point.
  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do
    do j=EMin, EMax
      do i=EMin, EMax
        mecha_energy%val(rcID,i,j,:) = &
          &   Grav * xy_HeightN%val(rcID,i,j,:) &
          & + 0.5d0 * dot( xy_VelN%val(rcID,i,j,:), xy_VelN%val(rcID,i,j,:) )
      end do
    end do
  end do

  ! 袖領域の貼りあわせ
  call paste_field2D_margin(mecha_energy)

  ! 力学的エネルギーの場に, 勾配演算子を作用させる. 
  ! Gradient operator acts on the mechanical energy field. 
  call gradient_op(diff_eval, mecha_energy, grad_mecha_energy)
  

  ! 運動方程式の半離散式の右辺を評価し, タイムレベル n における速度変化率を計算する. 
  ! Evaluates the right-hand side of the semi-discreted equation of motion,
  ! and calculates the rate of change of velocity at time level n with respect to time.
  !
  
  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do private(s_pos, unit_k)
    do j=EMin, EMax
      do i=EMin, EMax
        
        ! 格子点位置を極座標として取得する. 
        s_pos(:) = orth_to_geo_pos( icgrid%rcs_AGrid(rcID,i,j,:) )
        
        ! 格子点上におけるコリオリパラメータと鉛直方向の単位ベクトル(デカルト座標系表記)を求める. 
        if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
          unit_k(:) = local_k(s_pos(:))
        else
          ! 両極
          if ( rcID <= 5 ) then
            unit_k(:) = (/ 0.0d0, 0.0d0, 1.0d0 /)
            s_pos(2) = PI / 2.0d0
          else
            unit_k(:) = (/ 0.0d0, 0.0d0, -1.0d0 /)
            s_pos(2) = - PI / 2.0d0
          end if
        end if

       ! 速度変化率を計算する.
       DVelDtN%val(rcID,i,j,:) = &
         & - ( xy_Zeta%val(rcID,i,j,1) + xy_Coli%val(rcID,i,j,1) ) * cross( unit_k(:), xy_VelN%val(rcID,i,j,:) ) &
         & - grad_mecha_energy%val(rcID,i,j,:)
        
      end do
    end do
  end do
  
end subroutine calc_motion_eq_dvdt

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 非公開手続き
! Private procedure
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!> \~japanese
!! \brief 地理座標系の鉛直方向の単位ベクトルを直交座標系の単位ベクトルを使って表す.
!!
!! @param[in] s_pos  地理座標系における位置ベクトル\f$(\phi, \theta, r)\f$.
!! @return           引数で与えられた経度・緯度における鉛直方向の単位ベクトル(成分は直交座標系表示).
!!
!! \~english
!! \brief Represents the vertical unit vector in geographic coordinate system in terms of cartesian unit vectors.
!!
!! @param[in] s_pos  地理座標系における位置ベクトル\f$(\phi, \theta, r)\f$.  
!! @return           引数で与えられた経度・緯度における鉛直方向の単位ベクトル(成分は直交座標系表示). 
!! 
function local_k( &
    & s_pos        &  ! (in)
    & ) result(vec)

  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: s_pos(3)
  real(DP) :: vec(3)
  
  ! 作業変数
  ! Work variables
  !
  
  ! 地理座標系における鉛直方向の単位ベクトル.
  real(DP), parameter :: s_k(3) = (/ 0.0d0, 0.0d0, 1.0d0 /) 

  ! 実行文
  !

  ! 地理座標系における鉛直方向の単位ベクトルを, 直交座標系へと変換する. 
  vec = geo_to_orth_vec( s_k, s_pos)

end function local_k

end module motion_equation
