program unit_test1

  use dc_types, only: &
    & DP, STRING
	
	use dc_message, only: &
		& MessageNotify
		
  use LocalProcRect_Manager
  
  use IcGrid2D_FVM_Manager, only: &
  	& IcGrid2D_FVM, &
  	& IcGrid2D_FVM_Init, IcGrid2D_FVM_Final, &
  	& malloc_GP_array, malloc_CV_array, &
  	& get_LocalProcRectInfo
  
  use IcGrid2D_FVM_Builder, only: &
    & construct_icosahedral_grid, & 
    & GTYPE_STDGCGRID

  use Field_IcGrid2D_Manager
  
  use mpi

  implicit none

  type(IcGrid2D_FVM) :: icgrid2D, icgrid2D_mpi
  type(LocalProcRectInfo), pointer :: lcpcRCInfo
  type(Field_IcGrid2D) :: field2D
    
  integer :: glevel = 5
  real(DP) :: ic_radius = 1.0d0
  integer :: procNum, my_rank, ierr
  integer :: rcID, i, j
  
  real(DP) :: global_mean_val
  
  character(STRING) :: nc_filename = 'parallel_IO_file.nc'

	! MPI の初期化 
	! Initialization of MPI
  call MPI_Init(ierr)
  
  ! MPI のプロセス数とランクを取得する. 
  call MPI_Comm_size(MPI_COMM_WORLD, procNum, ierr)
  call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)

	! unit test 1-1
	!
  call unit_test_1_1()
	
	! unit test 1-2
	!
  call unit_test_1_2()

	! unit test 1-3
	!
	if ( my_rank == 0 ) then
	  write(*,*) '= Unit Test 1.3 [ Field_IcGrid2D Parallel Mode ]'
	  write(*,*) '* Initialize and finalize a variable of dderived type Field_IcGrid2D (MPI Parallel Mode).'
	  write(*,*) '* Test some procedures in Field_IcGrid2D_Manager module.'
	end if
	 
	!
	call Field_IcGrid2D_Init(field2D, icgrid2D_mpi, 'test_field', 1, 'This is a test field', 'm/s')
  lcpcRCInfo => get_LocalProcRectInfo(icgrid2D_mpi)
	field2D%val(:,:,:,:) = lcpcRCInfo%lcproc_MPIRank
	
	field2D%val(:,1,1,:) = -1d0
  field2D%val(:,2,1,:) = lcpcRCInfo%lcproc_MPIRank 
  field2D%val(:,1,2,:) = mod(lcpcRCInfo%lcproc_MPIRank+1, 5)
  field2D%val(:,2,2,:) = lcpcRCInfo%lcproc_MPIRank 

	call MessageNotify('M', 'Field_IcGrid2D(Parllel)', &
     & 'Test paste_field2D_margin(MPI parallel mode)', &
     & rank_mpi=0 )

	call MPI_Barrier(MPI_COMM_WORLD, ierr)

	call paste_field2D_margin(field2D)

	if ( my_rank == 0 .or. my_rank == 5 * 4**lcpcRCInfo%rlevel ) then
	  write(*,*) '* RCID=', lcpcRCInfo%globalRCID
	  lcpcRCInfo => get_LocalProcRectInfo(icgrid2D_mpi)

	  write(*,'(3f5.1)') field2D%val(1, 1:3, lcpcRCInfo%jIdMin, lcpcRCInfo%globalRCID)
	  write(*,'(3f5.1)') field2D%val(1, lcpcRCInfo%iEMax-2:lcpcRCInfo%iEMax, lcpcRCInfo%jIdMin, lcpcRCInfo%globalRCID)
	  
	  write(*,'(3f5.1)') field2D%val(1, lcpcRCInfo%iIdMax, 1:3, lcpcRCInfo%globalRCID)
	  write(*,'(3f5.1)') field2D%val(1, lcpcRCInfo%iIdMax, lcpcRCInfo%jEMax-2:lcpcRCInfo%jEMax, lcpcRCInfo%globalRCID)
    
    write(*,'(3f5.1)') field2D%val(1, 1:3, lcpcRCInfo%jIdMax, lcpcRCInfo%globalRCID)
	  write(*,'(3f5.1)') field2D%val(1, lcpcRCInfo%iEMax-2:lcpcRCInfo%iEMax, lcpcRCInfo%jIdMax, lcpcRCInfo%globalRCID)

    write(*,'(3f5.1)') field2D%val(1, lcpcRCInfo%iIdMin, 1:3, lcpcRCInfo%globalRCID)
	  write(*,'(3f5.1)') field2D%val(1, lcpcRCInfo%iIdMin, lcpcRCInfo%jEMax-2:lcpcRCInfo%jEMax, lcpcRCInfo%globalRCID)

	end if

	call MPI_Barrier(MPI_COMM_WORLD, ierr)

	call MessageNotify('M', 'Field_IcGrid2D(MPI parllel mode)', &
     & 'Check value of array storing data near both poles. ', &
     & rank_mpi=0 )
  
  if ( my_rank == 0 ) then
    do rcID=1,5
      write(*,'(i2,a,20f5.1)') rcID, ':', field2D%dataNearNPole(1, 1:2,1:2, rcID)
    end do
  end if
  
  if ( my_rank == 6 * 4**lcpcRCInfo%rlevel - 1) then
    write(*,'(20f5.1)') field2D%dataNearSPole(1, 32:33,32:33,6:10)
  end if

	call MPI_Barrier(MPI_COMM_WORLD, ierr)

  !
  field2D%val(:,:,:,:) = lcpcRCInfo%globalRCID
  global_mean_val = global_mean(field2D)
  
  call MessageNotify('M', 'Field_IcGrid2D(Parllel)', &
     & 'Test global_mean(MPI parallel mode)', &
     & rank_mpi=0 )
	
	call MessageNotify('M', 'Field_IcGrid2D::global_mean', &
     & 'calculated val=%f, answer=%f',                   &
     & d=(/ global_mean_val, 10 * 11 / 2.0d0 / 10d0 /),  &
     & rank_mpi=0 )


	call MPI_Barrier(MPI_COMM_WORLD, ierr)  
  call output_test_field(field2D, icgrid2D_mpi)

	call MPI_Barrier(MPI_COMM_WORLD, ierr)  
  call read_test_field(field2D, icgrid2D_mpi)
  
	! 後片付け
	!
	call Field_IcGrid2D_Final(field2D)
  call IcGrid2D_FVM_Final(icgrid2D_mpi)

	call MPI_Barrier(MPI_COMM_WORLD, ierr)  
  call MPI_Finalize(ierr)

contains

subroutine output_LocalProcRectInfo(icgrid)

	type(IcGrid2D_FVM), intent(inout) :: icgrid

  lcpcRCInfo => get_LocalProcRectInfo(icgrid)

  call MessageNotify('M', 'LocalProcRectInfo', &
    & 'MPIRank=%d, rlevel=%d, ' &
    & // 'iEMin=%d, iEMax=%d, jEMin=%d, jEMax=%d, '           &
    & // 'iIdMin=%d, iIdMax=%d, jIdMin=%d, jIdMax=%d, '       &
    & // 'globalRCID=%d, globalRCIDMin=%d, globalRCIDMax=%d', &
    & i=(/ lcpcRCInfo%lcproc_MPIRank, lcpcRCInfo%rlevel, &
     &     lcpcRCInfo%iEMin, lcpcRCInfo%iEMax,    &
     &     lcpcRCInfo%jEMin, lcpcRCInfo%jEMax,    &
     &     lcpcRCInfo%iIdMin, lcpcRCInfo%iIdMax,    &
     &     lcpcRCInfo%jIdMin, lcpcRCInfo%jIdMax,  &
     &     lcpcRCInfo%globalRCID, lcpcRCInfo%globalRCIDMin, lcpcRCInfo%globalRCIDMax &
     &     /), &
    & rank_mpi=-1 )

end subroutine output_LocalProcRectInfo

subroutine unit_test_1_1()

	! IcGrid2D_FVM の変数を初期化する. 
	! 格子点, コントロールボリュームのメモリを確保する. 
		
	if ( my_rank == 0 ) then
	  write(*,*) '= Unit Test 1.1 [ IcGrid2D_FVM Sequencial Mode ]'

	  call MessageNotify('M', 'IcGrid2D_FVM_Seq', &
     & 'Constructs STD-GC-grid(glevel=%d, ic_radius=%f)', &
     & i=(/ glevel /), d=(/ ic_radius /), rank_mpi=-1 )

	  call IcGrid2D_FVM_Init(icgrid2D, glevel, ic_radius)
		! 局所矩形領域の情報を出力する. 
		call output_LocalProcRectInfo(icgrid2D)
		call IcGrid2D_FVM_Final(icgrid2D)
		
	end if
	
end subroutine unit_test_1_1

subroutine unit_test_1_2()

	if ( my_rank == 0 ) then
	  write(*,*) '= Unit Test 1.2 [ IcGrid2D_FVM Parallel Mode ]'
	end if
	
	! IcGrid2D_FVM の変数を初期化する. 
	! 格子点, コントロールボリュームのメモリを確保する. 
		
	call MessageNotify('M', 'IcGrid2D_FVM_Parllel', &
     & 'Constructs STD-GC-grid(glevel=%d, ic_radius=%f)', &
     & i=(/ glevel /), d=(/ ic_radius /), rank_mpi=0 )

	call MPI_Barrier(MPI_COMM_WORLD, ierr)
	
  call IcGrid2D_FVM_Init(icgrid2D_mpi, glevel, ic_radius, mpi_flag=.true.)
  call malloc_GP_array(icgrid2D_mpi)
  call malloc_CV_array(icgrid2D_mpi)
  
  ! 局所矩形領域の情報を出力する. 
	call output_LocalProcRectInfo(icgrid2D_mpi)
  
  ! STD-GC-grid を構築する. 
  call construct_icosahedral_grid(GTYPE_STDGCGRID, icgrid2D_mpi)

end subroutine unit_test_1_2

subroutine output_test_field(field, icgrid_mpi)
  use IcGrid_ncWriter_mod

  type(Field_IcGrid2D), intent(in) :: field
  type(IcGrid2D_FVM), intent(in) :: icgrid_mpi
  
  type(IcGrid_ncWriter) :: writer
  integer :: ncfieldID
  
  call IcGrid_ncWriter_Init(writer, icgrid_mpi, parallel_IO_flag=.true.)
  
  call open_ncFile(writer, nc_filename)
  
  call ncdef_GridData(writer)
  ncfieldID = ncdef_FieldData(writer, field)
  call ncdef_Simulation_Parameter(writer, 1d0, 0d0, 1d0)
  call end_ncdef(writer)
  
  call write_GridData(writer)
  call write_FieldData(writer, ncfieldID, field)
  
  call close_ncFile(writer)

  call IcGrid_ncWriter_Final(writer)
  
end subroutine output_test_field

subroutine read_test_field(field, icgrid_mpi)
  use IcGrid_ncReader_mod

  type(Field_IcGrid2D), intent(inout) :: field
  type(IcGrid2D_FVM), intent(inout) :: icgrid_mpi

  type(IcGrid_ncReader) :: reader
  integer :: ncfieldID

  call IcGrid_ncReader_Init(reader, icgrid_mpi, parallel_IO_flag=.true.)

  call open_ncFile(reader, nc_filename)

  call read_GridData(reader)
  call read_FieldData(reader, get_field2D_name(field), field, time_level=1)

  call close_ncFile(reader)

  call IcGrid_ncReader_Final(reader)

end subroutine read_test_field

end program unit_test1
