module helmholtzEq_ForWDevAs_solver_mod


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

  !* gtool ***************
  !*

  ! 種類型パラメータ
  ! Kind type parameters
  !
  use dc_types, only: &
    & DP

  !* IGMBaseLib ***************
  !*
  
  ! 線形代数
  ! Linear algebra
  !
  use igmcore_SimLinearAlgebraEq, only: &
    & LUdecomp_TriDiagMat, sol_SimLinearAlEq_lutri

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

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

  ! 正二十面格子における広域矩形領域のサブドメインである局所矩形領域(各 MPI プロセスと対応)の情報管理
  ! Manager of information about a local rectangle region (corresponding to each MPI process)
  !
  use LocalProcRect_Manager, only: &
    & LocalProcRectInfo, &
    & get_InteriorRegion_IndexInfo
    
  ! 3 次元正二十面体格子の管理
  ! Manager of 3D icosahedral grid
  !
  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_vertical_level_num, get_icgrid2D


  !
  !
  !
  use IcosahedralLorentzGrid_Manager, only: &
    & IcosahedralLorentzGrid, &
    & get_TerrainFollowingCoordinate, &
    & LORENTZGRID_HARF_LEVEL, LORENTZGRID_FULL_LEVEL
  	
  !
  !
  !
  use Field_Grid1D_Manager, only: &
    & Field_Grid1D

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

  ! 3 次元正二十面格子上の物理場データの管理
  ! Manager of the physical field on 3D icosahedral grid
  !
  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & Field_IcGrid3D_Init, Field_IcGrid3D_Final, &
    & get_field_layer, get_icgrid3D

  use TFCoord_IcGrid3D_helper, only: &
    & TFCoord_IcGrid3D, &
    & TF_to_Z_IcGrid3D, calc_gamma

  use Derivate_Field_UnuniGrid1D_Manager, only: &
    & Derivate_Field_UnuniGrid1D, &
    & delx_2ord_3pts_unufGrid

  use VDerivate_Field_IcGrid3D_Manager, only: &
    & VDerivate_Field_IcGrid3D, & 
    & d_dVCoord, DELOPTRTYPE_FVM_CELL_DEL

  !* DIGMO ***********************
  !*

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private
  
  ! 公開手続き
  ! Public procedures
  !
  public :: solve_helmholtzEq_forWDevAs
	
	!
	!
	!
  
  real(DP), allocatable :: xyr_EConf(:,:,:,:,:)
    
contains
subroutine solve_helmholtzEq_forWDevAs( &
  & WDevAs, WDevNs, &
  & RDevNs, PDevNs, xyz_VhDevAs, &
  & xyr_hNl, gtildeNl, &
  & xyr_WNl, SR, SP, SW, &
  & IcosaLorentzGrid, tfcoordr, vrderivate, &
  & Grav, Rd, Cp, Cv, DelTimeShort, &
  & hydrostatic_opt &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: WDevAs
  type(Field_IcGrid3D), intent(in) :: WDevNs
  type(Field_IcGrid3D), intent(in) :: RDevNs
  type(Field_IcGrid3D), intent(in) :: PDevNs
  type(Field_IcGrid3D), intent(in) :: xyz_VhDevAs
  type(Field_IcGrid3D), intent(in) :: xyr_WNl
  type(Field_IcGrid3D), intent(in) :: xyr_hNl
  type(Field_IcGrid3D), intent(in) :: gtildeNl
  type(Field_IcGrid3D), intent(in) :: SR
  type(Field_IcGrid3D), intent(in) :: SP
  type(Field_IcGrid3D), intent(in) :: SW
  type(IcosahedralLorentzGrid), intent(in) :: IcosaLorentzGrid
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordr
  type(VDerivate_Field_IcGrid3D), intent(in) :: vrderivate
  real(DP), intent(in) :: Grav
  real(DP), intent(in) :: Rd
  real(DP), intent(in) :: Cp
  real(DP), intent(in) :: Cv
  real(DP), intent(in) :: DelTimeShort
  logical, intent(in) :: hydrostatic_opt

  ! 作業変数
  ! Local variables
  !
  type(IcGrid3D_FVM), pointer :: icgrid3Dr
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid2D), pointer :: xyr_WDevAs_k
  type(Field_IcGrid2D), pointer :: xyr_RHSVec_k

  real(DP) :: h_rijk
  real(DP) :: radius
  real(DP), allocatable :: C(:,:)
  real(DP), allocatable :: RHSVec(:,:,:,:)
  real(DP), allocatable :: AnswerVec(:,:,:,:)
  real(DP), allocatable :: Fl(:), Fd(:), Fu(:)

  integer :: iEMin, iEMax, jEMin, jEMax, KrMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  
  integer :: rcID, i, j, k
  real(DP) :: s_pos(3)

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3Dr => get_icgrid3D(WDevNs)
  icgrid2D => get_icgrid2D(icgrid3Dr)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  KrMax = get_vertical_level_num(icgrid3Dr)
  radius = get_IcRadius(icgrid2D)

  allocate(RHSVec(krMax, iEMin:iEMax,jEMin:jEMax, rcIDMin:rcIDMax))
  allocate(AnswerVec(KrMax, iEMin:iEMax,jEMin:jEMax, rcIDMin:rcIDMax))

  call prepare_RHSVec( &
    & RHSVec(:,:,:,:), &
    & WDevNs, RDevNs, PDevNs, &
    & xyr_WNl, SR, SP, SW, &
    & IcosaLorentzGrid, tfcoordr, vrderivate, &
    & icgrid3Dr%verticalGrid(:), &
    & radius, Grav, DelTimeShort, &
    & hydrostatic_opt, &
    & iEMin, jEMin, rcIDMin   &
    & )
  
  AnswerVec(:,:,:,:) = 0.0d0
  
  !$omp parallel  private(C, Fl, Fd, Fu)
  allocate(Fl(krMax))
  allocate(Fd(krMax))
  allocate(Fu(krMax))
  allocate(C(5,krMax))
  !$omp end parallel
  
  do rcID=rcIDMin, rcIDMax
    !$omp parallel do
    do j=jEMin, jEMax
      do i=iEMin, iEMax

         ! 鉛直運動量 W に関する連立一次方程式
        !  F * W = RHSVec
         ! の係数行列(三重対角行列) F を構築する.
         !

         ! 係数行列 F を求めるために用いる作業変数の値を求める.
        call prepare_C_rcij( &
          & rcID, i, j, &
          & C(:,:), &
          & xyr_hNl, gtildeNl, &
          & tfcoordr, icgrid3Dr%verticalGrid(:), &
          & radius, Grav, Rd, Cp, Cv, DelTimeShort, &
          & hydrostatic_opt &
           & )

         ! 上で求めた, A, B を使って, 三重対角行列を構築する.
        call construct_tridiagonal_matrix( &
          & Fl(:), Fd(:), Fu(:), &
          & C(:,:), &
          & icgrid3Dr%verticalGrid(:) &
           & )

         ! 連立一次方程式 F * W = RHSVec を LU 分解を使って解く.
         !

         ! 右辺のベクトルに対して, 下部境界条件を考慮する.
        RHSVec(2, i,j,rcID) = RHSVec(2, i,j,rcID) - Fl(2) * WDevAs%field2Ds(1)%val(1, i,j,rcID)
        !RHSVec(KrMax-1, i,j,rcID) = RHSVec(KrMax-1, i,j,rcID) - Fu(KrMax-1) * WDevAs%field2Ds(KrMax)%val(1, i,j,rcID)
       
         ! 三重対角行列 F を LU 分解する.
        call LUdecomp_TriDiagMat(Fl(2:KrMax-1), Fd(2:KrMax-1), Fu(2:KrMax-1))
        ! LU 分解法による連立一次方程式ソルバを使って, W を求める.
        AnswerVec(2:KrMax-1, i,j,rcID) = &
          &  sol_SimLinearAlEq_lutri(Fl(2:KrMax-1), Fd(2:KrMax-1), Fu(2:KrMax-1), RHSVec(2:krMax-1, i,j,rcID))

                  
      end do
    end do
  end do
  
  !$omp parallel do
  do k=2, KrMax-1
  !  xyr_WDevAs_k => get_field_layer(WDevAs, k)
    WDevAs%field2Ds(k)%val(1, iEMin:iEMax,jEMin:jEMax,:) = AnswerVec(k,:,:,:)
  end do
  
end subroutine solve_helmholtzEq_forWDevAs

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! モジュール外非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine prepare_RHSVec( &
  & RHSVec, &
  & WDevNs, RDevNs, PDevNs, &
  & WNl, SR, SP, SW, &
  & IcosaLorentzGrid, tfcoordr, vrderivate, &
  & Xi, radius, Grav, DelTimeShort, &
  & hydrostatic_opt, &
  & iEMin_, jEMin_, rcIDMin_ &
  & )
use igmcore_math  
  ! 宣言文 ; Declaration statements
  !
  integer, intent(in) :: iEMin_, jEMin_, rcIDMin_
  real(DP), intent(out) :: RHSVec(:, iEMin_:,jEMin_:,rcIDMin_:)
  type(Field_IcGrid3D), intent(in) :: WDevNs
  type(Field_IcGrid3D), intent(in) :: RDevNs
  type(Field_IcGrid3D), intent(in) :: PDevNs
  type(Field_IcGrid3D), intent(in) :: WNl
  type(Field_IcGrid3D), intent(in) :: SR
  type(Field_IcGrid3D), intent(in) :: SP
  type(Field_IcGrid3D), intent(in) :: SW
  type(IcosahedralLorentzGrid), intent(in) :: IcosaLorentzGrid
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordr
  type(VDerivate_Field_IcGrid3D), intent(in) :: vrderivate
  real(DP), intent(in) :: Xi(:)
  real(DP), intent(in) :: radius
  real(DP), intent(in) :: Grav
  real(DP), intent(in) :: DelTimeShort
  logical, intent(in) :: hydrostatic_opt

  ! 作業変数
  ! Local variables
  !
  integer :: rcID, i, j, k
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KrMax, KzMax, rcIDMin, rcIDMax

  real(DP) :: G_1div2
  real(DP) :: gamma_, gammaPowM2_
  real(DP) :: gammaPowM2(size(Xi))

  type(Field_IcGrid3D) :: tmp_RHSvderivateTerm
  type(Field_IcGrid3D) :: RHSvderivateTerm
  type(TFCoord_IcGrid3D), pointer :: tfcoordz
  
  type(IcGrid3D_FVM), pointer :: icgrid3Dr, icgrid3Dz
  type(Field_IcGrid2D), pointer :: WDevNs_k
  type(Field_IcGrid2D), pointer :: RDevNs_k
  type(Field_IcGrid2D), pointer :: PDevNs_k
  type(Field_IcGrid2D), pointer :: SR_k
  type(Field_IcGrid2D), pointer :: SP_k
  type(Field_IcGrid2D), pointer :: SW_k
  type(Field_IcGrid2D), pointer :: xyr_RHS_k
  type(Field_IcGrid2D), pointer :: tmp_RHSvderivateTerm_k
  type(Field_IcGrid2D), pointer :: RHSvderivateTerm_k
  real(DP) :: alpha, s_pos(3), dumpCoef
  type(IcGrid2D_FVM), pointer :: icgrid2D

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3Dr => get_icgrid3D(WDevNs)
  icgrid3Dz => get_icgrid3D(PDevNs)
  icgrid2D => get_icgrid2D(icgrid3Dr)
  tfcoordz => get_TerrainFollowingCoordinate(IcosaLorentzGrid, LORENTZGRID_FULL_LEVEL)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  KrMax = get_vertical_level_num(icgrid3Dr)
  KzMax = get_vertical_level_num(icgrid3Dz)

  if ( hydrostatic_opt ) then
    alpha = 0.0d0
  else
    alpha = 1.0d0
  end if

  call Field_IcGrid3D_Init(tmp_RHSvderivateTerm, icgrid3Dr, 'tmp_RHSvderivateTerm', 1)
  call Field_IcGrid3D_Init(RHSvderivateTerm, icgrid3Dr, 'RHSvderivateTerm', 1)

  do k=1, krMax
    tmp_RHSvderivateTerm_k => get_field_layer(tmp_RHSvderivateTerm, k)
    PDevNs_k => get_field_layer(PDevNs, k)
    SP_k => get_field_layer(SP, k)

    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private(i, rcID, gamma_, G_1div2)
      do j=jEMin, jEMax
        do i=iEMin, iEMax


          gamma_ = calc_gamma(TF_to_Z_IcGrid3D(tfcoordr, icgrid3Dr%verticalGrid(k), rcID, i, j), radius)
          G_1div2 = tfcoordr%G_1div2%val(1, i,j,rcID)

          tmp_RHSvderivateTerm_k%val(1, i,j,rcID) =  &
            & ( PDevNs_k%val(1, i,j,rcID) + DelTimeShort * SP_k%val(1, i,j,rcID) ) / ( G_1div2 * gamma_**2 )

        end do
      end do
    end do
  end do

  call d_dVCoord( &
    & vrderivate, tmp_RHSvderivateTerm, RHSvderivateTerm)

  do k=1, krMax
    WDevNs_k => get_field_layer(WDevNs, k)
    SW_k => get_field_layer(SW, k)
    SR_k => get_field_layer(SR, k)
    RDevNs_k => get_field_layer(RDevNs, k)
    RHSvderivateTerm_k => get_field_layer(RHSvderivateTerm, k)


!if ( xi(k) > 25d03 ) then
!DumpCoef = 0.5d0 * (1d0-cos(PI *(Xi(k)-25d03)/(40d03-25d03))) / 1800d0
!else
!DumpCoef = 0d0
!end if

    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private(i, gamma_, gammaPowM2_)
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          gamma_ = calc_gamma(TF_to_Z_IcGrid3D(tfcoordr, Xi(k), rcID, i, j), radius)
          gammaPowM2_ = 1.0d0 / gamma_**2

          RHSVec(k, i,j, rcID) = &
            &   gammaPowM2_ * ( alpha * WDevNs_k%val(1, i,j,rcID) + DelTimeShort * SW_k%val(1, i,j,rcID) ) &
!                & -dumpCoef * (WNl%field2Ds(k)%val(1,i,j,rcID)+WDevNs%field2Ds(k)%val(1,i,j,rcID)))) &
            & - DelTimeShort * RHSvderivateTerm_k%val(1, i,j,rcID) &
            & - DelTimeShort * gammaPowM2_ * Grav *                &
            &              ( RDevNs_k%val(1, i,j,rcID) + DelTimeShort * SR_k%val(1, i,j,rcID) )

        end do
      end do
    end do
  end do

  call Field_IcGrid3D_Final(tmp_RHSvderivateTerm)
  call Field_IcGrid3D_Final(RHSvderivateTerm)

end subroutine prepare_RHSVec

subroutine prepare_C_rcij( &
  & rcID, i, j, &
  & C, &
  & xyr_hNl, gtildeNl, &
  & tfcoordr, &
  & Xi, radius, Grav, Rd, Cp, Cv, DelTimeShort, &
  & hydrostatic_opt &
  & )
  
  ! 宣言文 ; Declaration statements
  !
  integer, intent(in) :: rcID
  integer, intent(in) :: i
  integer, intent(in) :: j
  real(DP), intent(out) :: C(:,:)
  type(Field_IcGrid3D), intent(in) :: xyr_hNl
  type(Field_IcGrid3D), intent(in) :: gtildeNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordr
  real(DP), intent(in) :: Xi(:)
  real(DP), intent(in) :: radius
  real(DP), intent(in) :: Grav
  real(DP), intent(in) :: Rd
  real(DP), intent(in) :: Cp
  real(DP), intent(in) :: Cv
  real(DP), intent(in) :: DelTimeShort
  logical, intent(in) :: hydrostatic_opt

  ! 作業変数
  ! Local variables
  !
  integer k, krMax
  real(DP) :: G_1div2
  real(DP) :: gamma_

  type(IcGrid3D_FVM), pointer :: icgrid3Dr
  type(Field_IcGrid2D), pointer :: xyr_hNl_k
  type(Field_IcGrid2D), pointer :: gtildeNl_k
  real(DP) :: constCoef
  real(DP) :: alpha
  real(DP) :: DelTimeShortPow2, gammaPow2
  
  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  krMax = size(Xi)
  icgrid3Dr => get_icgrid3D(xyr_hNl)

  G_1div2 = tfcoordr%G_1div2%val(1, i,j,rcID)
  DelTimeShortPow2 = DelTimeShort**2
  constCoef = DelTimeShortPow2 * Rd / ( Cv * G_1div2 )

  if ( hydrostatic_opt ) then
    alpha = 0.0d0
  else
    alpha = 1.0d0
  end if

  do k=1, krMax
    xyr_hNl_k => get_field_layer(xyr_hNl, k)
    gtildeNl_k => get_field_layer(gtildeNl, k)

    gamma_ = calc_gamma(TF_to_Z_IcGrid3D(tfcoordr, Xi(k), rcID, i, j), radius)
    gammaPow2 = gamma_ * gamma_
    
    C(1,k) = 1d0 / ( G_1div2 * gammaPow2 )
    C(2,k) = constCoef * xyr_hNl_k%val(1, i,j,rcID)
    C(3,k) = constCoef * gtildeNl_k%val(1, i,j,rcID) / gammaPow2
    C(4,k) = DelTimeShortPow2 * Grav * C(1,k)     
    C(5,k) = alpha / gamma_**2
  end do

end subroutine prepare_C_rcij

subroutine construct_tridiagonal_matrix( &
  & Fl, Fd, Fu, &
  & C, Xi &
  & )
  
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(out) :: Fl(:)
  real(DP), intent(out) :: Fd(:)
  real(DP), intent(out) :: Fu(:)
  real(DP), intent(in) :: C(:,:)
  real(DP), intent(in) :: Xi(:)

  ! 作業変数
  ! Local variables
  !
  integer :: k
  integer :: kmax
  real(DP) :: DelXipm1
  real(DP) :: DelXiU, DelXiL
  real(DP) :: DelXipm05
  real(DP) :: C1Up, C1Low
  ! 実行文 ; Executable statements
  !
  
  kmax = size(Fl)

  do k=2, kmax-1
    DelXipm1 = Xi(k+1) - Xi(k-1)
    DelXipm05 = 0.5d0 * DelXipm1
    DelXiU = Xi(k+1) - Xi(k)
    DelXiL = Xi(k) - Xi(k-1)
    
    C1Up = 0.5d0 * ( C(1,k+1) + C(1,k) )
    C1Low = 0.5d0 * ( C(1,k) + C(1,k-1) )
    
    Fl(k) = - C1Low * C(2,k-1)  / ( DelXipm05 * DelXiL ) + ( C(3,k-1) + C(4,k) ) / DelXipm1 
    Fd(k) = C(5,k) + ( C1up / DelXiU + C1Low / DelXiL ) * C(2,k) / DelXipm05
    Fu(k) = - C1Up * C(2,k+1)  / ( DelXipm05 * DelXiU ) - ( C(3,k+1) + C(4,k) ) / DelXipm1
  end do

end subroutine construct_tridiagonal_matrix

end module helmholtzEq_ForWDevAs_solver_mod
