!> \~japanese
!! \brief 浅水方程式ソルバを提供するモジュール.
!! \~english
!! \brief This module provides the shallow water equations solver.
!! \~
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
!!
module sw_equation_solver

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

  !* gtool **
  !*

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

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

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


  !* IGModel-SW **
  !*

  ! 物理場の管理
  ! Manager of the physical fields
  !
  use field_manager, only: &
    & xy_Vel_TL_list, xy_Height_TL_list, xy_Htopo, &
    & DVelDt_TL_list, DHeightDt_TL_list, &
    & TL_N, TL_Nplus1, &
    & TL_DDT_N, TL_DDT_Nminus1, TL_DDT_Nminus2, &
    & xy_Coli, diff_eval

  ! 運動方程式の評価
  ! Evaluate motion equation
  !
  use motion_equation, only: &
    & motion_equation_Init, &
    & calc_motion_eq_dvdt

  ! 連続の式の評価
  ! Evaluate continious equation
  !
  use continious_equation, only: &
    & continious_equation_Init, &
    & calc_continious_eq_dhdt

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

  ! 宣言文 ; Declaration statement
  !
  implicit none
  private
    
  ! 公開手続き
  ! Public procedure
  !
  public :: init_sw_equation_solver, temporal_integration

  ! 非公開変数
  ! Private variables
  !

contains

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

  ! 実行文 ; Executable statements
  !

  !write(*,*) 'init sw_equation_solver'

  ! 運動方程式・連続の式の評価用モジュールの初期化. 
  call motion_equation_Init(icgrid)
  call continious_equation_Init(icgrid)

end subroutine init_sw_equation_solver


!
!> \~japanese
!! \brief 浅水方程式系の時間積分を行う.
!!
!! @param[in] tstep       時間スッテプ数
!! @param[in] dt          時間刻み幅.
!! @param[in] icgrid  構造型 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Performs time integration of the shallow water equations.
!!
!! @param[in] tstep       The number of time step.
!! @param[in] dt          Time step.
!! @param[in] icgrid  The variable of derived type IcGrid2D_FVM.
!!
subroutine temporal_integration( &
  & tstep, dt,                   &  ! (in)
  & icgrid                       &  ! (inout)
  )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: tstep
  real(DP), intent(in) :: dt
  type(IcGrid2D_FVM), intent(in) :: icgrid

  ! 作業変数
  ! Work variables
  !
  integer :: idMin, EMin, EMax

  ! 実行文 ; Executable statement
  !
  

  ! 最初の 2 Step は, 4 次の Runge=Kutta 法を使って時間積分する. 
  ! それ以降は, 3 次の Adams=Bashforth 法を使って時間積分する. 
  !
  if( tstep <= 2 ) then
    idMin = get_IdMin(icgrid)

    ! 4 次の Runge=Kutta 法により, 浅水方程式を時間積分する. 
    ! 各時間スッテプにおいて, dvdt, dhdt を記憶しておく. 
    call RungeKutta_fourth_order( &
      & xy_Vel_TL_list(TL_Nplus1)%val, xy_Vel_TL_list(TL_N)%val, DVelDt_TL_list(TL_DDT_N)%val, &
      & xy_Height_TL_list(TL_Nplus1)%val, xy_Height_TL_list(TL_N)%val, DHeightDt_TL_list(TL_DDT_N)%val, &
      & dt, icgrid, idMin )

  else

    ! 3 次の Adams=Bashforth 法により, 浅水方程式を時間積分する. 
    ! 運動方程式を時間積分して, タイムレベル n+1 の速度場を求める. 
    ! motionEq_flag が偽のときは, 運動方程式の時間積分は行われない. 
    !
    
    EMin = get_EffSize_Min(icgrid)
    EMax = get_EffSize_Max(icgrid)

    if ( motionEq_flag ) then
      ! 速度場の時間変化率を計算する. 
      call calc_motion_eq_dvdt( &
        & DVelDtN    = DVelDt_TL_list(TL_DDT_N), &
        & xy_VelN    = xy_Vel_TL_list(TL_N),     &
        & xy_HeightN = xy_Height_TL_list(TL_N), &
        & xy_Coli    = xy_Coli, &
        & diff_eval = diff_eval )
    end if
    
    ! 高度場の時間変化率を計算する.
    call calc_continious_eq_dhdt( &
      & DHeightDtN = DHeightDt_TL_list(TL_DDT_N), &
      & xy_VelN    = xy_Vel_TL_list(TL_N), &
      & xy_HeightN = xy_Height_TL_list(TL_N), &
      & xy_Htopo   = xy_Htopo, &
      & diff_eval=diff_eval )
    
    !write(*,*) tstep, ':dh=', dhdt_TL_list(TL_DDT_N)%val(10,EMax,EMax,1)*dt

    !$omp sections
    
    !$omp section
    if ( motionEq_flag ) then
      ! 速度場に対する時間積分を行う. 
      call Adams_Bashforth_third_order( &
        & U_nplus1     = xy_Vel_TL_list(TL_Nplus1)%val(:,EMin:EMax,EMin:EMax,:), &
        & U_n          = xy_Vel_TL_list(TL_N)%val(:,EMin:EMax,EMin:EMax,:),      &
        & dUdt_n       = DVelDt_TL_list(TL_DDT_N)%val(:,EMin:EMax,EMin:EMax,:), &
        & dUdt_nminus1 = DVelDt_TL_list(TL_DDT_Nminus1)%val(:,EMin:EMax,EMin:EMax,:), &
        & dUdt_nminus2 = DVelDt_TL_list(TL_DDT_Nminus2)%val(:,EMin:EMax,EMin:EMax,:), &
        & dt           = dt )

    end if
    
    !$omp section
    ! 高度場に対する時間積分を行う. 
    call Adams_Bashforth_third_order(&
      & U_nplus1     = xy_Height_TL_list(TL_Nplus1)%val(:,EMin:EMax,EMin:EMax,:), &
      & U_n          = xy_Height_TL_list(TL_N)%val(:,EMin:EMax,EMin:EMax,:), &
      & dUdt_n       = DHeightDt_TL_list(TL_DDT_N)%val(:,EMin:EMax,EMin:EMax,:), &
      & dUdt_nminus1 = DHeightDt_TL_list(TL_DDT_Nminus1)%val(:,EMin:EMax,EMin:EMax,:), &
      & dUdt_nminus2 = DHeightDt_TL_list(TL_DDT_Nminus2)%val(:,EMin:EMax,EMin:EMax,:), &
      & dt           = dt )
    
    !$omp end sections

  end if

end subroutine temporal_integration


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

!
!> \~japanese
!! \brief 4 次の Runge=Kutta 法を用いて, 浅水方程式系を時間積分する.
!!
!! \~english
!! \brief Performs time integration of the shallow water equations using the forth-order Runge=Kutta method.
!!
!! \~japanese
!! このサブルーチンでは, 時刻\f$t\f$(タイムレベル n )の物理場の値から未来の時刻\f$t+\Delta t\f$(タイムレベル n+1 )の物理場の値を予報する.
!! 今, 速度場と高度場の変数の組を\f$ \Dvect{V} = \{\Dvect{v}, h\} \f$と書くことにする.
!! また, \f$ L_{motion}, L_{continious} \f$をそれぞれ運動方程式と連続の式の半離散式の右辺を表す関数とし,
!! \f$\Dvect{L}=\{ L_{motion}, L_{continious} \} \f$とする.
!! このとき, 以下の式を計算し, 未来の物理場を予報する.
!!
!! \~english
!! In this subroutine, the future physical fields at \f$t+\Delta t\f$(time level n+1) will be predicted in terms of the fields
!! at current time level n.
!!
!! The future physical field will be predicted using the following equation as
!!
!! \~
!! \f[ \Dvect{V}_{n+1} = \Dvect{V}_n + \dfrac{\Dvect{k}_1 + 2\Dvect{k}_2 + 2\Dvect{k}_3 + \Dvect{k}_4}{6} \f]
!!
!! \~japanese
!! ここで, \f$\Dvect{k}_i\f$は次のように定義される.
!! \~english
!! Where \f$ \Dvect{V} = \{\Dvect{v}, h\} \f$ is a set of variables of velocity and surface height field.
!! \f$ L_{motion}, L_{continious} \f$ are functions which represent righ-hand sides of semi-discrete equation of motion, and continious equation, respectively.
!! Then, a set of the functions is defined as \f$\Dvect{L}=\{ L_{motion}, L_{continious} \} \f$.
!! \f$\Dvect{k}_i\f$ are defined as
!!
!! \~
!! \f{eqnarray*}{
!!  \Dvect{k}_1 &=& \Delta t \; \Dvect{L}(t_n, \Dvect{v}, h),  \\\\
!!  \Dvect{k}_2 &=& \Delta t \; \Dvect{L}(t_n + \dfrac{\Delta t}{2}, \Dvect{v}+\dfrac{\Dvect{k}_{1v}}{2}, h+\dfrac{k_{1h}}{2}) , \\\\
!!  \Dvect{k}_3 &=& \Delta t \; \Dvect{L}(t_n + \dfrac{\Delta t}{2}, \Dvect{v}+\dfrac{\Dvect{k}_{2v}}{2}, h+\dfrac{k_{2h}}{2}) , \\\\
!!  \Dvect{k}_4 &=& \Delta t \; \Dvect{L}(t_n+\Delta t, \Dvect{v}+\Dvect{k}_{3v}, h+k_{3h}).
!! \f}
!!
!! \~japanese
!! @param[in,out] xy_VelA   タイムレベル n+1 の速度場データを保持する配列.
!! @param[in]     xy_VelN   タイムレベル n の速度場データを保持する配列.
!! @param[in,out] DVelDtN    タイムレベル n の速度場データの時間変化率を保持する配列.
!! @param[in,out] h_nplus1  タイムレベル n+1 の表面高度場データを保持する配列.
!! @param[in]     h_n       タイムレベル n の表面高度場データを保持する配列.
!! @param[in,out] dhdt_n    タイムレベル n の表面高度場の時間変化率のデータを保持する配列.
!! @param[in]     dt        時間刻み幅 [s].
!! @param[in]     icgrid    構造型 IcGrid2D_FVM の変数.
!! @param[in]     idMin     各種物理場を保持する配列の下限境界インデックス.
!! \~english
!! @param[in,out] xy_VelA   The array in which the velocity field data at time level n+1 is stored.
!! @param[in]     xy_VelN   The array in which the velocity field data at time level n is stored.
!! @param[in,out] DVelDtN    The array in which the time rate of change of velocity field data at time level n is stored.
!! @param[in,out] h_nplus1  The array in which the surface height field data at time level n+1 is stored.
!! @param[in]     h_n       The array in which the surface height field data at time level n is stored.
!! @param[in,out] dhdt_n    The array in which the time rate of change of surface height field data at time level n is stored.
!! @param[in]     dt        Time step [s].
!! @param[in]     icgrid    The variable of derived type IcGrid2D_FVM.
!! @param[in]     idMin     The lower bound index of arrays in which various physical fields data is stored.
!!
subroutine RungeKutta_fourth_order( &
  & xy_VelA,                          &  ! (inout)
  & xy_VelN, DVelDtN,                 &  ! (in)
  & xy_HeightA,                       &  ! (inout)
  & xy_HeightN, DHeightDtN,           &  ! (in)
  & dt, icgrid, idMin                 &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  integer, intent(in)    :: idMin
  real(DP), intent(inout) :: xy_VelA(:,idMin:,idMin:,:)
  real(DP), intent(in)    :: xy_VelN(:,idMin:,idMin:,:)
  real(DP), intent(inout) :: DVelDtN(:,idMin:,idMin:,:)
  real(DP), intent(inout) :: xy_HeightA(:,idMin:,idMin:,:)
  real(DP), intent(in)    :: xy_HeightN(:,idMin:,idMin:,:)
  real(DP), intent(inout) :: DHeightDtN(:,idMin:,idMin:,:)
  real(DP), intent(in)    :: dt
  type(IcGrid2D_FVM), intent(in) :: icgrid ! IcGrid2D_FVM クラスのインスタンスの参照.

  ! 作業変数 
  ! Work variables
  !
  type(Field_IcGrid2D) :: k_v(4)
  type(Field_IcGrid2D) :: v_predictor
  type(Field_IcGrid2D) :: k_h(4)
  type(Field_IcGrid2D) :: h_predictor
  
  integer, parameter :: K_NUM = 4
  real(DP), parameter :: k_coef(3) = (/ 0.5d0, 0.5d0, 1.0d0 /)
  integer :: m,n
 
  ! 実行文 ; Executable statements
  !
  
  ! 作業配列を初期化する. 
  ! Initializes the work variables. 

  do m=1, K_NUM
    call Field_IcGrid2D_Init( &
      & k_v(m), &
      & icgrid=icgrid, name='k_v', rank=3 )

    call Field_IcGrid2D_Init( &
      & k_h(m), &
      & icgrid=icgrid, name='k_h', rank=1 )

  end do
  
  call Field_IcGrid2D_Init( &
    & self=v_predictor, &
    & icgrid=icgrid, name='v_predictor', rank=3 )

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

  ! k_v(1), k_h(1) をまず求めるために, 半離散式 L の右辺の関数の引数 v_predictor, h_predictor として, 
  ! v_n, h_n を与える.  
  v_predictor%val(:,:,:,:) = xy_VelN(:,:,:,:)
  h_predictor%val(:,:,:,:) = xy_HeightN(:,:,:,:)

  do m=1, K_NUM
    call paste_field2D_margin(v_predictor)
    call paste_field2D_margin(h_predictor)

    ! (t_n, v_predictor, h_predictor) を引数とする運動方程式の半離散式を評価し, 
    ! k_v(m) を求める. 
    
    call calc_motion_eq_dvdt( &
      & DVelDtN=k_v(m), &
      & xy_VelN=v_predictor, xy_HeightN=h_predictor, xy_Coli=xy_Coli, &
      & diff_eval=diff_eval )

    ! (t_n, v_predictor, h_predictor, h_s) を引数とする運動方程式の半離散式を評価し, 
    ! k_h(m) を求める. 
    call calc_continious_eq_dhdt( &
      & DHeightDtN=k_h(m), &
      & xy_VelN=v_predictor, xy_HeightN=h_predictor, xy_Htopo=xy_Htopo, &
      & diff_eval=diff_eval )
    

    ! k_v(m+1), k_h(m+1) (ただし, 2 <= m < 4) を求めるために,
    ! 半離散式 L の右辺の関数の引数 v_predictor, h_predictor を更新する.

    if( m < K_NUM )then

      ! 速度場に対して
      if ( motionEq_flag ) then
        v_predictor%val(:,:,:,:) = xy_VelN(:,:,:,:) + k_coef(m) * dt * k_v(m)%val(:,:,:,:)
      else
        ! 運動方程式を使用しない場合は, 予測子にタイムレベル n 速度場を代入する.
        v_predictor%val(:,:,:,:) = xy_VelN(:,:,:,:)
      end if

      ! 表面高度場に対して
      h_predictor%val(:,:,:,:) = xy_HeightN(:,:,:,:) + k_coef(m) * dt * k_h(m)%val(:,:,:,:)

    end if

  end do
  
  ! 求めた k_v(1:4), k_h(1:4) を使って, タイムレベル n の局所時間変化率 dh/dt, dv/dt を求める. 
  !

  DVelDtN(:,:,:,:) = ( k_v(1)%val(:,:,:,:) + 2.0d0 * k_v(2)%val(:,:,:,:) &
    &                  + 2.0d0 * k_v(3)%val(:,:,:,:) + k_v(4)%val(:,:,:,:) &
    &                 ) / 6.0d0

  DHeightDtN(:,:,:,:) = ( k_h(1)%val(:,:,:,:) + 2.0d0 * k_h(2)%val(:,:,:,:) &
    &                     + 2.0d0 * k_h(3)%val(:,:,:,:) + k_h(4)%val(:,:,:,:) &
    &                    ) / 6.0d0

  ! タイムレベル n+1 の速度場, 高度場を求める. 
  !

  xy_VelA(:,:,:,:) = xy_VelN(:,:,:,:) + DVelDtN(:,:,:,:) * dt
  xy_HeightA(:,:,:,:) = xy_HeightN(:,:,:,:) + DHeightDtN(:,:,:,:) * dt


  ! 作業配列を破棄する. 
  !

  do m=1, K_NUM
    call Field_IcGrid2D_Final(k_v(m))
    call Field_IcGrid2D_Final(k_h(m))
  end do
   
end subroutine RungeKutta_fourth_order

!
!> \~japanese
!! \brief 3 次の Adams=Bashforth 法を用いて, 浅水方程式系を時間積分する.
!!
!! \~english
!! \brief Performs time integration of the shallow water equations using the third-order Adams-Bashforth method.
!! 
!! \~japanese
!! このサブルーチンでは, 現在の時刻\f$t\f$(タイムレベル n )の物理場の値と, タイムレベル n,n-1,n-2 の物理場の局所時間微分を用いて,
!! 時刻\f$t+\Delta t\f$(タイムレベル n+1 )の未来の物理場の値を予報する.
!!
!! 全格子点におけるタイムレベル n+1 の物理場\f$\Dvect{U}\f$の値は, 以下の式によって予測される.
!!
!! \~english
!! In this subroutine, the future physical fields at \f$t+\Delta t\f$(time level n+1) will be predicted in terms of
!! the fields at current time level n and the local time derivatives of fields at current and previous time levels(n, n-1, n-2).
!!
!! The value of fields at time level n+1 at all grid-points will be predicted using the equation as
!!
!! \~
!! \f[
!!   \Dvect{U}_{n+1} = \Dvect{U}_n + 
!!                      \dfrac{1}{12} \Delta t \left[
!!                        23 \left(\DP{\Dvect{U}}{t}\right)_{n} - 16 \left(\DP{\Dvect{U}}{t}\right)_{n-1} + 5 \left(\DP{\Dvect{U}}{t}\right)_{n-2}
!!                      \right]
!! \f]
!!
!! \~japanese
!! @param[in,out] U_nplus1     タイムレベル n+1 の物理場データを保持する配列.
!! @param[in]     dUdt_n       タイムレベル n の物理場の局所時間微分データを保持する配列.
!! @param[in]     dUdt_nminus1 タイムレベル n-1 の物理場の局所時間微分データを保持する配列.
!! @param[in]     dUdt_nminus2 タイムレベル n-2 の物理場の局所時間微分データを保持する配列.
!! @param[in]     dt           時間刻み幅 [s].
!! \~english
!! @param[in,out] U_nplus1     The array in which the physical field data at time level n+1 is stored.
!! @param[in]     dUdt_n       The array in which the local time derivative data of the physical field at time level n is stored.
!! @param[in]     dUdt_nminus1 The array in which the local time derivative data of the physical field at time level n-1 is stored.
!! @param[in]     dUdt_nminus2 The array in which the local time derivative data of the physical field at time level n-2 is stored.
!! @param[in]     dt           Time step [s].
!!
subroutine Adams_Bashforth_third_order( &
  & U_nplus1,                                &  ! (inout)
  & U_n, dUdt_n, dUdt_nminus1, dUdt_nminus2, &  ! (in)
  & dt )                                        ! (in)

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(inout) :: U_nplus1(:,:,:,:)
  real(DP), intent(in) :: U_n(:,:,:,:)
  real(DP), intent(in) :: dUdt_n(:,:,:,:)
  real(DP), intent(in) :: dUdt_nminus1(:,:,:,:) 
  real(DP), intent(in) :: dUdt_nminus2(:,:,:,:)
  real(DP), intent(in) :: dt

  ! 実行文 ; Executable statements
  !

  U_nplus1(:,:,:,:) = &
    & U_n(:,:,:,:) &
    & + dt * ( 23.0d0 * dUdt_n(:,:,:,:) - 16.0d0 * dUdt_nminus1(:,:,:,:) + 5.0d0 * dUdt_nminus2(:,:,:,:) ) / 12.0d0
 
end subroutine Adams_Bashforth_third_order

end module sw_equation_solver
