module EnergyAnalysis

  use dc_types, only: &
    & DP

  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_vertical_level_num, &
    & get_icgrid2d

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

  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & paste_field2D_margin, global_mean

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

  use Vertical_Interpolation, only: &
    & interpolate_fieldr_to_z

  use Utility_GeoCoord, only: &
    & local_k
  
  implicit none
  private

	public :: Analize_SystemEnergy
	
contains

subroutine Analize_SystemEnergy( &
	& gm_TotalE, gm_KineticE, gm_InternalE, gm_PotentialE, &
	& gm_Mass, &
	& xyz_TE, xyz_Vh, xyr_W, xyz_RhoBasic, xyz_R, xyz_E, &
	& tfcoordz, tfcoordr, Grav )
	
	real(DP), intent(out) :: gm_TotalE
	real(DP), intent(out) :: gm_KineticE
	real(DP), intent(out) :: gm_InternalE
	real(DP), intent(out) :: gm_PotentialE
	real(DP), intent(out) :: gm_Mass
	type(Field_IcGrid3D), intent(in) :: xyz_TE
	type(Field_IcGrid3D), intent(in) :: xyz_Vh
	type(Field_IcGrid3D), intent(in) :: xyr_W
	type(Field_IcGrid3D), intent(in) :: xyz_RhoBasic
	type(Field_IcGrid3D), intent(in) :: xyz_R
	type(Field_IcGrid3D), intent(in) :: xyz_E
	type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
	type(TFCoord_IcGrid3D), intent(in) :: tfcoordr
	real(DP), intent(in) :: Grav

	type(Field_IcGrid3D) :: TEField3D
	type(Field_IcGrid3D) :: KEField3D
	type(Field_IcGrid3D) :: PEField3D
	type(Field_IcGrid3D) :: IEField3D
	type(Field_IcGrid3D) :: DensField3D
	
	real(DP) :: gm_TE(1) 
	real(DP) :: gm_KE(1), gm_IE(1), gm_PE(1)
	real(DP) :: gm_dens(1)
	
	type(IcGrid3D_FVM), pointer :: icgrid3Dz
	
	icgrid3Dz => get_icgrid3D(xyz_Vh)
	
	call Field_IcGrid3D_Init(TEField3D, icgrid3Dz, 'TEField3D', 1)
	call Field_IcGrid3D_Init(KEField3D, icgrid3Dz, 'KEField3D', 1)
	call Field_IcGrid3D_Init(PEField3D, icgrid3Dz, 'PEField3D', 1)
	call Field_IcGrid3D_Init(IEField3D, icgrid3Dz, 'IEField3D', 1)
	call Field_IcGrid3D_Init(DensField3D, icgrid3Dz, 'DensField3D', 1)	
	
	call calc_TotalEnergyField( &
  	& TEField3D, &
  	& xyz_TE,  &
  	& tfcoordz &
  	& )
  	
	call calc_KineticEnergyField( &
  	& KEField3D, &
  	& xyz_Vh, xyr_W, xyz_RhoBasic, xyz_R, &
  	& tfcoordz &
  	& )
 	
	call calc_PotentialEnergyField( &
  	& PEField3D, &
  	& xyz_RhoBasic, xyz_R, &
  	& tfcoordz, Grav &
  	& )

	call calc_InternalEnergyField( &
  	& IEField3D, &
  	& xyz_E, &
  	& tfcoordz &
  	& ) 	

	call calc_DensityField( &
  	& DensField3D, &
  	& xyz_RhoBasic, xyz_R, &
  	& tfcoordz &
  	& )
  		
	gm_TE = global_mean_Field_TFIcGrid3Dz(tfcoordz, TEField3D, 1, tfcoordr)
	gm_KE = global_mean_Field_TFIcGrid3Dz(tfcoordz, KEField3D, 1, tfcoordr)
	gm_PE = global_mean_Field_TFIcGrid3Dz(tfcoordz, PEField3D, 1, tfcoordr)
	gm_IE = global_mean_Field_TFIcGrid3Dz(tfcoordz, IEField3D, 1, tfcoordr)
	gm_dens = global_mean_Field_TFIcGrid3Dz(tfcoordz, DensField3D, 1, tfcoordr)
	
	gm_TotalE = gm_TE(1)
	gm_KineticE = gm_KE(1)
	gm_PotentialE = gm_PE(1)
	gm_InternalE = gm_IE(1)
	gm_Mass = gm_dens(1)
	
	call Field_IcGrid3D_Final(TEField3D)
	call Field_IcGrid3D_Final(KEField3D)
	call Field_IcGrid3D_Final(PEField3D)
	call Field_IcGrid3D_Final(IEField3D)
	call Field_IcGrid3D_Final(DensField3D)
	
end subroutine Analize_SystemEnergy



subroutine calc_KineticEnergyField( &
  & xyz_KE, xyz_Vh, xyr_W, xyz_RhoBasic, xyz_R, &
  & tfcoordz &
  & )

  type(Field_IcGrid3D), intent(inout) :: xyz_KE
  type(Field_IcGrid3D), intent(in) :: xyz_Vh
  type(Field_IcGrid3D), intent(in) :: xyr_W
  type(Field_IcGrid3D), intent(in) :: xyz_RhoBasic
  type(Field_IcGrid3D), intent(in) :: xyz_R
	type(TFCoord_IcGrid3D), intent(in) :: tfcoordz

	type(Field_IcGrid3D) :: xyz_W
  type(Field_IcGrid2D), pointer :: xyz_KE_k
  type(Field_IcGrid2D), pointer :: xyz_Vh_k
  type(Field_IcGrid2D), pointer :: xyz_W_k
  type(Field_IcGrid2D), pointer :: xyz_RhoBasic_k
  type(Field_IcGrid2D), pointer :: xyz_R_k

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

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

  real(DP) :: xi, gamma_rijk, radius
  
  icgrid3Dz => get_icgrid3D(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(xyr_W, xyz_W)
	
	do k=1, KzMax
		xyz_KE_k => get_field_layer(xyz_KE, k)
		xyz_Vh_k => get_field_layer(xyz_Vh, k)
		xyz_W_k => get_field_layer(xyz_W, k)
		xyz_RhoBasic_k => get_field_layer(xyz_RhoBasic, k)
		xyz_R_k => get_field_layer(xyz_R, k)

		xi = icgrid3Dz%verticalGrid(k)
		
		do rcID=rcIDMin, rcIDMax
			do j=jEMin, jEMax
				do i=iEMin, iEMax
					
					gamma_rijk = &
						& calc_gamma( TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j), radius)
					
					xyz_KE_k%val(1, i,j,rcID) = &
						& calc_KE_rcijk( &
							&	xyz_Vh_k%val(1:3, i,j,rcID), xyz_W_k%val(1, i,j,rcID), &
							& xyz_RhoBasic_k%val(1, i,j,rcID), xyz_R_k%val(1, i,j,rcID), &
							& icgrid2D%rcs_AGrid(1:3, i,j,rcID), &
							& tfcoordz%G_1div2%val(1, i,j,rcID), gamma_rijk &
							& )
							  
				end do
			end do
		end do
	end do
	
	call Field_Icgrid3D_Final(xyz_W)
	
end subroutine calc_KineticEnergyField

subroutine calc_PotentialEnergyField( &
	& xyz_PE, xyz_RhoBasic, xyz_R, &
	& tfcoordz, Grav &
	& )
	
	type(Field_IcGrid3D), intent(inout) :: xyz_PE
	type(Field_IcGrid3D), intent(in) :: xyz_RhoBasic
	type(Field_IcGrid3D), intent(in) :: xyz_R
	type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
	real(DP), intent(in) :: Grav
	
  type(Field_IcGrid2D), pointer :: xyz_PE_k
  type(Field_IcGrid2D), pointer :: xyz_RhoBasic_k
  type(Field_IcGrid2D), pointer :: xyz_R_k

  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(IcGrid2D_FVM), pointer :: icgrid2D
	real(DP) :: rho_rijk
  real(DP) :: xi, gamma_rijk
  real(DP) :: radius

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

  icgrid3Dz => get_icgrid3D(xyz_R)
  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)
	
	do k=1, KzMax

		xyz_PE_k => get_field_layer(xyz_PE, k)
		xyz_RhoBasic_k => get_field_layer(xyz_RhoBasic, k)
		xyz_R_k => get_field_layer(xyz_R, k)
		xi = icgrid3Dz%verticalGrid(k)
		
		do rcID=rcIDMin, rcIDMax
			do j=jEMin, jEMax
				do i=iEMin, iEMax
					
					gamma_rijk = &
						& calc_gamma( TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j), radius)
			
					rho_rijk = ( xyz_RhoBasic_k%val(1, i,j,rcID) + xyz_R_k%val(1, i,j,rcID) ) &
												&  / ( tfcoordz%G_1div2%val(1, i,j,rcID) * gamma_rijk**2 )
												
					xyz_PE_k%val(1, i,j,rcID) = &
							&  rho_rijk 	* ( Grav * TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j) )
						 
				end do
			end do
		end do
	end do
	
end subroutine calc_PotentialEnergyField


subroutine calc_InternalEnergyField( &
	& xyz_IE, xyz_E, &
	& tfcoordz &
	& )
	
	type(Field_IcGrid3D), intent(inout) :: xyz_IE
	type(Field_IcGrid3D), intent(in) :: xyz_E
	type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
	
  type(Field_IcGrid2D), pointer :: xyz_IE_k
  type(Field_IcGrid2D), pointer :: xyz_E_k

  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(IcGrid2D_FVM), pointer :: icgrid2D
	real(DP) :: rho_rijk
  real(DP) :: xi, gamma_rijk
  real(DP) :: radius

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

  icgrid3Dz => get_icgrid3D(xyz_E)
  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)
	
	do k=1, KzMax

		xyz_IE_k => get_field_layer(xyz_IE, k)
		xyz_E_k => get_field_layer(xyz_E, k)
		xi = icgrid3Dz%verticalGrid(k)
		
		do rcID=rcIDMin, rcIDMax
			do j=jEMin, jEMax
				do i=iEMin, iEMax
					
					gamma_rijk = &
						& calc_gamma( TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j), radius)
			
					xyz_IE_k%val(1, i,j,rcID) = &
							&  xyz_E_k%val(1, i,j,rcID) / ( tfcoordz%G_1div2%val(1, i,j,rcID) * gamma_rijk**2 )
						 
				end do
			end do
		end do
	end do
	
end subroutine calc_InternalEnergyField


subroutine calc_TotalEnergyField( &
	& xyz_TotalE, xyz_TE, &
	& tfcoordz &
	& )
	
	type(Field_IcGrid3D), intent(inout) :: xyz_TotalE
	type(Field_IcGrid3D), intent(in) :: xyz_TE
	type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
	
  type(Field_IcGrid2D), pointer :: xyz_TotalE_k
  type(Field_IcGrid2D), pointer :: xyz_TE_k

  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(IcGrid2D_FVM), pointer :: icgrid2D
	real(DP) :: rho_rijk
  real(DP) :: xi, gamma_rijk
  real(DP) :: radius

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

  icgrid3Dz => get_icgrid3D(xyz_TE)
  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)
	
	do k=1, KzMax

		xyz_TotalE_k => get_field_layer(xyz_TotalE, k)
		xyz_TE_k => get_field_layer(xyz_TE, k)
		xi = icgrid3Dz%verticalGrid(k)
		
		do rcID=rcIDMin, rcIDMax
			do j=jEMin, jEMax
				do i=iEMin, iEMax
					
					gamma_rijk = &
						& calc_gamma( TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j), radius)
			
					xyz_TotalE_k%val(1, i,j,rcID) = &
							&  xyz_TE_k%val(1, i,j,rcID) / ( tfcoordz%G_1div2%val(1, i,j,rcID) * gamma_rijk**2 )
						 
				end do
			end do
		end do
	end do
	
end subroutine calc_TotalEnergyField

subroutine calc_DensityField( &
	& xyz_dens, xyz_RhoBasic, xyz_R, &
	& tfcoordz &
	& )
	
	type(Field_IcGrid3D), intent(inout) :: xyz_dens
	type(Field_IcGrid3D), intent(in) :: xyz_RhoBasic
	type(Field_IcGrid3D), intent(in) :: xyz_R
	type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
	
  type(Field_IcGrid2D), pointer :: xyz_dens_k
  type(Field_IcGrid2D), pointer :: xyz_RhoBasic_k
  type(Field_IcGrid2D), pointer :: xyz_R_k

  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(IcGrid2D_FVM), pointer :: icgrid2D
	real(DP) :: rho_rijk
  real(DP) :: xi, gamma_rijk
  real(DP) :: radius

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

  icgrid3Dz => get_icgrid3D(xyz_R)
  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)
	
	do k=1, KzMax

		xyz_dens_k => get_field_layer(xyz_dens, k)
		xyz_RhoBasic_k => get_field_layer(xyz_RhoBasic, k)
		xyz_R_k => get_field_layer(xyz_R, k)
		xi = icgrid3Dz%verticalGrid(k)
		
		do rcID=rcIDMin, rcIDMax
			do j=jEMin, jEMax
				do i=iEMin, iEMax
					
					gamma_rijk = &
						& calc_gamma( TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j), radius)
			
					xyz_dens_k%val(1, i,j,rcID) = ( xyz_RhoBasic_k%val(1, i,j,rcID) + xyz_R_k%val(1, i,j,rcID) ) &
												&  / ( tfcoordz%G_1div2%val(1, i,j,rcID) * gamma_rijk**2 )
												
				end do
			end do
		end do
	end do
	
end subroutine calc_DensityField



!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function calc_KE_rcijk( &
	& Vh_rijk, W_rijk, &
	& RhoBasic_rijk, R_rijk, &
	& gp_OrthPos, G_1div2_rij, gamma_rijk &
	& ) result(KE)
	
	real(DP), intent(in) :: Vh_rijk(3)
	real(DP), intent(in) :: W_rijk
	real(DP), intent(in) :: RhoBasic_rijk
	real(DP), intent(in) :: R_rijk
	real(DP), intent(in) :: gp_OrthPos(3)
	real(DP), intent(in) :: G_1div2_rij
	real(DP), intent(in) :: gamma_rijk
	real(DP) KE
	
	real(DP) :: s_pos(3)
	real(DP) :: univec_k(3)
	real(DP) :: Vel_rijk(3)
	real(DP) :: Rho_rijk
	
	s_pos(:) = orth_to_geo_pos(gp_OrthPos(:))
	univec_k(:) = local_k(s_pos(1), s_pos(2))
	
	Vel_rijk(:) = ( Vh_rijk(:) + W_rijk * univec_k(:) ) / ( RhoBasic_rijk + R_rijk )
	Rho_rijk = ( RhoBasic_rijk + R_rijk ) / ( G_1div2_rij * gamma_rijk**2 )
	   
	KE = 0.5d0 * Rho_rijk * dot_product( Vel_rijk(:), Vel_rijk(:))
	
end function calc_KE_rcijk

end module EnergyAnalysis
