module ProgVar_Preparation
  
  ! モジュール引用 ; Use statements
  !
    
  use, intrinsic :: iso_c_binding

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

  use dc_types, only: &
    & DP

  !* IGMBaseLib ********************************
  !*
  
  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos
    
  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

  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_IcGrid2D, get_vertical_level_num

  use Field_Grid1D_Manager, only:  &
    & Field_Grid1D

  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D

  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & Field_IcGrid3D_Init, Field_IcGrid3D_Final, &
    & get_field_layer, get_icgrid3D, get_field3D_rank

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

  !* DIGMO ********************************
  !*
  
  use Vertical_Interpolation, only: &
    & interpolate_fieldr_to_z, &
    & interpolate_fieldz_to_r

  use Utility_GeoCoord, only: &
    & local_k

  ! 宣言文 ; Declaration statements
  !  
  implicit none
  private
  
  ! 公開手続き
  ! Public statements
  !
  public :: prepare_BasicStateVars, & 
    &       prepare_ProgVars1, prepare_ProgVars2, prepare_ProgVars3, &
    &        calc_initial_TE
    
contains
subroutine prepare_BasicStateVars( &
  & cptr_z_RhoBasic, cptr_z_PressBasic, &
  & cptr_xyz_RhoBasic, cptr_xyz_PressBasic, &
  & cptr_tfcoordz &
  & ) bind(C, name='prepare_BasicStateVars')

  type(c_ptr), value :: cptr_z_RhoBasic
  type(c_ptr), value :: cptr_z_PressBasic
  type(c_ptr), value :: cptr_xyz_RhoBasic
  type(c_ptr), value :: cptr_xyz_PressBasic
  type(c_ptr), value :: cptr_tfcoordz

  type(Field_Grid1D), pointer :: ptr_z_RhoBasic
  type(Field_Grid1D), pointer :: ptr_z_PressBasic
  type(Field_IcGrid3D), pointer :: ptr_xyz_RhoBasic
  type(Field_IcGrid3D), pointer :: ptr_xyz_PressBasic
  type(TFCoord_IcGrid3D), pointer :: ptr_tfcoordz

  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid2D), pointer :: xyz_RhoBasic_k
  type(Field_IcGrid2D), pointer :: xyz_PressBasic_k

  integer :: iEMin, iEMax, jEMin, jEMax, kzMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: rcID, i, j, k
  real(DP) :: gamma_, radius
  real(DP) :: G_1div2_rcijk

  call c_f_pointer(cptr_z_RhoBasic, ptr_z_RhoBasic)
  call c_f_pointer(cptr_z_PressBasic, ptr_z_PressBasic)
  call c_f_pointer(cptr_xyz_RhoBasic, ptr_xyz_RhoBasic)
  call c_f_pointer(cptr_xyz_PressBasic, ptr_xyz_PressBasic)
  call c_f_pointer(cptr_tfcoordz, ptr_tfcoordz)

  icgrid3Dz => get_icgrid3D(ptr_xyz_PressBasic)
  icgrid2D => get_icgrid2D(icgrid3Dz)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  kzMax = get_vertical_level_num(icgrid3Dz)
  radius = get_IcRadius(icgrid2D)

  call interpolate_FieldZ1D_into_FieldTF3D( &
    & ptr_tfcoordz, ptr_z_RhoBasic, ptr_xyz_RhoBasic &
    & )

  call interpolate_FieldZ1D_into_FieldTF3D( &
    & ptr_tfcoordz,  ptr_z_PressBasic, ptr_xyz_PressBasic &
    & )

  do k=1, kzMax

    xyz_RhoBasic_k => get_field_layer(ptr_xyz_RhoBasic, k)
    xyz_PressBasic_k => get_field_layer(ptr_xyz_PressBasic, k)

    do rcID=rcIDMin, rcIDMax
      do j=jEMin, jEMax
        do i=iEMin, iEMax
          gamma_ = calc_gamma(TF_to_Z_IcGrid3D(ptr_tfCoordz, icgrid3Dz%verticalGrid(k), rcID, i, j), radius)
          G_1div2_rcijk = ptr_tfcoordz%G_1div2%val(1, i,j,rcID)


          xyz_RhoBasic_k%val(1, i,j,rcID) = &
            & xyz_RhoBasic_k%val(1, i,j,rcID) * G_1div2_rcijk * gamma_**2

          xyz_PressBasic_k%val(1, i,j,rcID) = &
            & xyz_PressBasic_k%val(1, i,j,rcID) * G_1div2_rcijk * gamma_**2

        end do
      end do
    end do
  end do

end subroutine prepare_BasicStateVars

subroutine prepare_ProgVars1( &
  & cptr_xyz_Phi, &
  & cptr_tfcoord &
  & ) bind(C, name='prepare_ProgVars1')

  type(c_ptr), value :: cptr_xyz_Phi
  type(c_ptr), value :: cptr_tfcoord

  type(Field_IcGrid3D), pointer :: ptr_xyz_Phi
  type(TFCoord_IcGrid3D), pointer :: ptr_tfcoordz

  integer :: iEMin, iEMax, jEMin, jEMax, kMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: rcID, i, j, k

  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid2D), pointer :: xyz_Phi_k

  real(DP) :: gamma_, radius
  real(DP) :: G_1div2_rcijk
  real(DP) :: xi

  call c_f_pointer(cptr_xyz_Phi, ptr_xyz_Phi)
  call c_f_pointer(cptr_tfcoord, ptr_tfcoordz)

  icgrid3D => get_icgrid3D(ptr_xyz_Phi)
  icgrid2D => get_icgrid2D(icgrid3D)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  kMax = get_vertical_level_num(icgrid3D)
  radius = get_IcRadius(icgrid2D)

  do k=1, kMax
    xi = icgrid3D%verticalGrid(k)
    xyz_Phi_k => get_field_layer(ptr_xyz_Phi, k)

    do rcID=rcIDMin, rcIDMax
      do j=jEMin, jEMax
        do i=iEMin, iEMax
         gamma_ = calc_gamma(TF_to_Z_IcGrid3D(ptr_tfcoordz, xi, rcID, i, j), radius)
         G_1div2_rcijk = ptr_tfcoordz%G_1div2%val(1, i,j,rcID)

         xyz_Phi_k%val(:, i,j,rcID) = &
          & xyz_Phi_k%val(:, i,j,rcID) * G_1div2_rcijk * gamma_**2

        end do
      end do
    end do

  end do

end subroutine prepare_ProgVars1

subroutine prepare_ProgVars2( &
  & cptr_xyz_Phi, &
  & cptr_xyz_RhoBasic, cptr_xyz_R &
  & ) bind(C, name='prepare_ProgVars2')

  type(c_ptr), value :: cptr_xyz_Phi
  type(c_ptr), value :: cptr_xyz_RhoBasic
  type(c_ptr), value :: cptr_xyz_R

  type(Field_IcGrid3D), pointer :: ptr_xyz_Phi
  type(Field_IcGrid3D), pointer :: ptr_xyz_RhoBasic
  type(Field_IcGrid3D), pointer :: ptr_xyz_R

  integer :: iEMin, iEMax, jEMin, jEMax, kMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: rcID, i, j, k
  integer :: rank
  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid2D), pointer :: xyz_Phi_k
  type(Field_IcGrid2D), pointer :: xyz_RhoBasic_k, xyz_R_k

  call c_f_pointer(cptr_xyz_Phi, ptr_xyz_Phi)
  call c_f_pointer(cptr_xyz_RhoBasic, ptr_xyz_RhoBasic)
  call c_f_pointer(cptr_xyz_R, ptr_xyz_R)

  icgrid3D => get_icgrid3D(ptr_xyz_Phi)
  icgrid2D => get_icgrid2D(icgrid3D)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  kMax = get_vertical_level_num(icgrid3D)

  rank = get_field3D_rank(ptr_xyz_Phi)

  do k=1, kMax
     xyz_Phi_k => get_field_layer(ptr_xyz_Phi, k)
     xyz_RhoBasic_k => get_field_layer(ptr_xyz_RhoBasic, k)
     xyz_R_k => get_field_layer(ptr_xyz_R, k)

    do rcID=rcIDMin, rcIDMax
      do j=jEMin, jEMax
        do i=iEMin, iEMax
          xyz_Phi_k%val(1:rank, i,j,rcID) = xyz_Phi_k%val(1:rank, i,j,rcID) &
                &          * ( xyz_RhoBasic_k%val(1, i,j,rcID) + xyz_R_k%val(1, i,j,rcID) )
         end do
       end do
     end do
  end do

end subroutine prepare_ProgVars2

subroutine prepare_ProgVars3( &
  & cptr_xyr_Phi, &
  & cptr_xyz_RhoBasic, cptr_xyz_R &
  & ) bind(C, name='prepare_ProgVars3')

  type(c_ptr), value :: cptr_xyr_Phi
  type(c_ptr), value :: cptr_xyz_RhoBasic
  type(c_ptr), value :: cptr_xyz_R

  type(Field_IcGrid3D), pointer :: ptr_xyr_Phi
  type(Field_IcGrid3D), pointer :: ptr_xyz_RhoBasic
  type(Field_IcGrid3D), pointer :: ptr_xyz_R

  integer :: kMax
  integer :: k
  integer :: rank

  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid3D) :: xyr_RhoBasic
  type(Field_IcGrid3D) :: xyr_R

  type(Field_IcGrid2D), pointer :: xyr_Phi_k
  type(Field_IcGrid2D), pointer :: xyr_RhoBasic_k, xyr_R_k

  call c_f_pointer(cptr_xyr_Phi, ptr_xyr_Phi)
  call c_f_pointer(cptr_xyz_RhoBasic, ptr_xyz_RhoBasic)
  call c_f_pointer(cptr_xyz_R, ptr_xyz_R)

  icgrid3D => get_icgrid3D(ptr_xyr_Phi)
  icgrid2D => get_icgrid2D(icgrid3D)
  kMax = get_vertical_level_num(icgrid3D)

  rank = get_field3D_rank(ptr_xyr_Phi)

  call Field_IcGrid3D_Init(xyr_RhoBasic, icgrid3D, 'xyr_RhoBasic', 1)
  call Field_IcGrid3D_Init(xyr_R, icgrid3D, 'xyr_R', 1)

  call interpolate_fieldz_to_r(ptr_xyz_RhoBasic, xyr_RhoBasic)
  call interpolate_fieldz_to_r(ptr_xyz_R, xyr_R)

  do k=1, kMax
     xyr_Phi_k => get_field_layer(ptr_xyr_Phi, k)
     xyr_RhoBasic_k => get_field_layer(xyr_RhoBasic, k)
     xyr_R_k => get_field_layer(xyr_R, k)

     xyr_Phi_k%val(1:rank, :,:,:) = xyr_Phi_k%val(1:rank, :,:,:) &
      &                         * ( xyr_RhoBasic_k%val(1:1, :,:,:) + xyr_R_k%val(1:1, :,:,:) )

  end do

  call Field_IcGrid3D_Final(xyr_RhoBasic)
  call Field_IcGrid3D_Final(xyr_R)


end subroutine prepare_ProgVars3

subroutine calc_initial_TE( &
  & cptr_xyz_TE, &
  & cptr_xyz_RhoBasic, cptr_xyz_R, cptr_xyz_E, &
  & cptr_xyz_Vh, cptr_xyr_W, &
  & cptr_tfcoordz, &
  & Grav &
  & ) bind(C, name='calc_initial_TE')

  type(c_ptr), value :: cptr_xyz_TE
  type(c_ptr), value :: cptr_xyz_RhoBasic
  type(c_ptr), value :: cptr_xyz_R
  type(c_ptr), value :: cptr_xyz_E
  type(c_ptr), value :: cptr_xyz_Vh
  type(c_ptr), value :: cptr_xyr_W
  type(c_ptr), value :: cptr_tfcoordz
  real(DP), value :: Grav

  type(Field_IcGrid3D), pointer :: ptr_xyz_TE
  type(Field_IcGrid3D), pointer :: ptr_xyz_RhoBasic
  type(Field_IcGrid3D), pointer :: ptr_xyz_R
  type(Field_IcGrid3D), pointer :: ptr_xyz_E
  type(Field_IcGrid3D), pointer :: ptr_xyz_Vh
  type(Field_IcGrid3D), pointer :: ptr_xyr_W
  type(TFCoord_IcGrid3D), pointer :: ptr_tfcoordz

  type(Field_IcGrid3D) :: xyz_W
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz

  type(Field_IcGrid2D), pointer :: xyz_Vptr_k
  type(Field_IcGrid2D), pointer :: xyz_W_k
  type(Field_IcGrid2D), pointer :: xyz_E_k
  type(Field_IcGrid2D), pointer :: xyz_RhoBasic_k
  type(Field_IcGrid2D), pointer :: xyz_R_k
  type(Field_IcGrid2D), pointer :: xyz_TE_k

  real(DP) :: z_rijk, G_1div2_rijk
  real(DP) :: Vel_rijk(3), Rho_rijk
  real(DP) :: s_pos(3), univec_k(3)
  real(DP) :: xi, gamma_
  real(DP) :: radius

  integer :: iEMin, iEMax, jEMin, jEMax, KzMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: rcID, i, j, k

  call c_f_pointer(cptr_xyz_TE, ptr_xyz_TE)
  call c_f_pointer(cptr_xyz_RhoBasic, ptr_xyz_RhoBasic)
  call c_f_pointer(cptr_xyz_R, ptr_xyz_R)
  call c_f_pointer(cptr_xyz_E, ptr_xyz_E)
  call c_f_pointer(cptr_xyz_Vh, ptr_xyz_Vh)
  call c_f_pointer(cptr_xyr_W, ptr_xyr_W)
  call c_f_pointer(cptr_tfcoordz, ptr_tfcoordz)

  icgrid3Dz => get_icgrid3D(ptr_xyz_Vh)
  icgrid2D => get_icgrid2D(icgrid3Dz)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  KzMax = get_vertical_level_num(icgrid3Dz)
  radius = get_IcRadius(icgrid2D)

  call Field_IcGrid3D_Init(xyz_W, icgrid3Dz, 'xyz_W', 1)
  call interpolate_fieldr_to_z(ptr_xyr_W, xyz_W)

  do k=1, KzMax
    xi = icgrid3Dz%verticalGrid(k)

    xyz_E_k => get_field_layer(ptr_xyz_E, k)
    xyz_Vptr_k => get_field_layer(ptr_xyz_Vh, k)
    xyz_W_k => get_field_layer(xyz_W, k)
    xyz_TE_k => get_field_layer(ptr_xyz_TE, k)
    xyz_RhoBasic_k => get_field_layer(ptr_xyz_RhoBasic, k)
    xyz_R_k => get_field_layer(ptr_xyz_R, k)

    do rcID=rcIDMin, rcIDMax
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:, i,j,rcID))
          univec_k(:) = local_k(s_pos(1), s_pos(2))

          z_rijk = TF_to_Z_IcGrid3D(ptr_tfCoordz, xi, rcID, i, j)
          gamma_ = calc_gamma(z_rijk, radius)
          G_1div2_rijk = ptr_tfcoordz%G_1div2%val(1, i,j,rcID)

          Rho_rijk = xyz_RhoBasic_k%val(1, i,j,rcID) + xyz_R_k%val(1, i,j,rcID)
          Vel_rijk(1:3) = ( xyz_vptr_k%val(1:3, i,j,rcID) + xyz_W_k%val(1, i,j,rcID) * univec_k(1:3) ) &
            &                  / Rho_rijk

          xyz_TE_k%val(1, i,j,rcID) =   xyz_E_k%val(1, i,j,rcID)      &
            &                        + Rho_rijk * (                 &
            &                           + Grav * z_rijk             &
            &                           + 0.5d0 * dot_product(Vel_rijk(:), Vel_rijk(:)) &
             &                              )

        end do
      end do
    end do

  end do

  call Field_IcGrid3D_Final(xyz_W)

end subroutine calc_initial_TE

end module ProgVar_Preparation

