!= 
!
! Authors::   Yuta Kawai
! Version::   $Id: grid_test.f90,v 1.2 2011-05-17 15:34:20 ykawai Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
program main
  use dc_types

  use dc_test, only: AssertLessThan, AssertEqual

  use igmcore_math

  ! 球面三角法
  ! Spherical trigonometry
  !
  use igmcore_spherical_trigonometry, only: &
    & geodesic_arc_length
  
  use IcGrid2D_FVM_Manager

  use IcGrid2D_FVM_Builder

  use IcGrid_ncGridLoader, only: load_gridData

  implicit none
  
  type(IcGrid2D_FVM) :: ic_grid
  integer :: div_level = 5
  real(DP) :: ic_radius = 1.0d0

  !
  !
  write(*,*) '= grid_test 1'
  write(*,*) '* construct_icosahedral_grid(1)'
  call construct_icosahedral_grid(GTYPE_STDGRID, div_level, ic_radius, ic_grid)

  call AssertEqual( message='Check glevel', &
    & answer=div_level, check=get_glevel(ic_grid)  )

  call AssertEqual( message='Check the mean distance between neighboring grid points.', &
    & answer=2.0d0 * PI * ic_radius / ( 10d0 * 2.0d0**(div_level-1) ), &
    & check=get_IcGridMeanInterval(ic_grid)  )


  call IcGrid2D_FVM_Final(ic_grid)

  !
  !
  write(*,*) '= grid_test 2'
  write(*,*) '* construct_icosahedral_grid(2)'
  call IcGrid2D_FVM_Init(ic_grid, div_level, 1.0d0 )
  call malloc_GP_array(ic_grid)
  call malloc_CV_array(ic_grid)
  call construct_icosahedral_grid(GTYPE_STDGRID, ic_grid)

  
  ! netCDF 書き込みテスト
  !
  write(*,*) '= grid_test 3'
  write(*,*) '* Writes icosahedral grid data into NetCDF file.'
  call grid_test3(ic_grid)

  ! netCDF 読み込みテスト
  !
  write(*,*) '= grid_test 4'
  write(*,*) '* Reads icosahedral grid data from NetCDF file.'
  call grid_test4(ic_grid)

  ! 
  write(*,*) '= grid_test 5'
  write(*,*) '* Calculates the vertex and area of control volumes'

  call calc_CV_vertex(ic_grid)
  call calc_CV_Area(ic_grid)

  ! 
  write(*,*) '= grid_test 6'
  write(*,*) '* Check control volume shape'

!  call construct_icosahedral_grid(GTYPE_STDGRID, ic_grid)
  call load_gridData( '../diffOptr_test/SPR-GC-grid_glevel5.nc', ic_grid )
  call grid_test6(ic_grid)


  call IcGrid2D_FVM_Final(ic_grid)

contains
subroutine grid_test3(icgrid)
  use IcGrid_ncWriter_mod

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

  type(IcGrid_ncWriter) :: writer
  integer :: EMin, EMax, EffSize

  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  EffSize = EMax - EMin + 1

  call IcGrid_ncWriter_Init(writer, icgrid)
  call open_ncFile(writer, 'icgrid.nc')
  call ncdef_GridData(writer)
  call end_ncdef(writer)

  call write_GridData(writer)
  call close_ncFile(writer)

  write(*,*) 'succeed to write netCDF file'

end subroutine grid_test3

subroutine grid_test4(icgrid)
  use IcGrid_ncReader_mod

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

  type(IcGrid_ncReader) :: reader
  integer :: EMin, EMax, EffSize

  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  EffSize = EMax - EMin + 1

  call IcGrid_ncReader_Init(reader, icgrid)
  call open_ncFile(reader, 'icgrid.nc')
  call read_GridData(reader)
  call close_ncFile(reader)

  write(*,*) 'succeed to read netCDF file'

end subroutine grid_test4

subroutine grid_test6(icgrid)
  use IcGrid_ncWriter_mod
  use Field_IcGrid2D_Manager
  
  type(IcGrid2D_FVM), intent(inout) :: icgrid
  
  type(Field_IcGrid2D) :: l_field, q_field
  integer :: l_fieldID, q_fieldID
  type(IcGrid_ncWriter) :: writer
  
  integer :: EMin, EMax, EffSize
  integer :: rcID, i, j, cvID, cvIDup
  
  integer :: cv_num
  real(DP) :: l(3,6), theta(6), q
  real(DP) :: q_max
  integer :: rcIDMax, iMax, jMax 
  real(DP) :: cv_pts(3,6) 
  real(DP) :: ic_radius
  
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  EffSize = EMax - EMin + 1
  
  call Field_IcGrid2D_Init(l_field, icgrid, 'l', 1)
  call Field_IcGrid2D_Init(q_field, icgrid, 'q', 1)
  
  ic_radius = get_IcRadius(icgrid)
  q_max = 0d0
  
  l_field%val(:,:,:,:) = -5000d0
  
  do rcID=1, RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax

        l_field%val(1, i,j,rcID) = sqrt(icgrid%rcs_CVInfo(i,j,rcID)%cv_area)

        if ( check_pole(icgrid, rcID, i, j) /= NOT_POLE_FLAG ) then
        else
        
        if( ( rcID > 5 .and. i==EMin .and. j==EMin ) .or. ( i==EMin .and. j==EMax ) .or. &
            & ( i==EMax .and. j==EMin ) .or. ( rcID <= 5 .and. i==EMax .and. j==EMax )     &
            &       ) then
            Cv_num = 5
        else
            cv_num = 6
        end if

     
        cv_pts(1:3, 1:cv_num) = icgrid%rcs_CVInfo(i,j,rcID)%cv_pts(1:3, 1:cv_num)
        do cvID=1, cv_num
           cvIDup = 1 + mod(cvID, cv_num)        
           l(1, cvID) = geodesic_arc_length(ic_radius, icgrid%rcs_AGrid(1:3, i,j,rcID), cv_pts(1:3, cvID)) / ic_radius
           l(2, cvID) = geodesic_arc_length(ic_radius, cv_pts(1:3, cvIDup), cv_pts(1:3, cvID)) / ic_radius          
           l(3, cvID) = geodesic_arc_length(ic_radius, icgrid%rcs_AGrid(1:3, i,j,rcID), cv_pts(1:3, cvIDup)) / ic_radius

           theta(cvID) = dacos( &
            &        ( dcos(l(2,cvID)) - dcos(l(1,cvID)) * dcos(l(3,cvID)) ) / ( dsin(l(1,cvID)) * dsin(l(3,cvID)) ) &
              &         )
    
        end do
        
        q = maxval(theta(1:cv_num)) / minval(theta(1:cv_num)) - 1d0
        q_field%val(1, i,j,rcID) = q
        
        !write(*,*) rcID,i,j, q
        if ( q > q_max ) then
          q_max = q
        !  write(*,*) rcID,i,j, q
        !  write(*,*) (theta(cvID), cvID=1,cv_num)
        !  write(*,*) (l(2, cvID), cvID=1,cv_num)
            
          rcIDMax = rcID; iMax = i; jMax = j;
        end if
        
        end if
      end do
    end do
  end do


 write(*,*) 'Check control volume shape'
 write(*,*) 'Max distortion location:', q_max, '-',rcIDMax, iMax, jMax

  call IcGrid_ncWriter_Init(writer, icgrid)
  
  call open_ncFile(writer, 'shapeData.nc')
  call ncdef_GridData(writer)
  call ncdef_simulation_parameter(writer, 1d0, 1d0, 1d0)
  l_fieldID = ncdef_FieldData(writer, l_field)
  q_fieldID = ncdef_FieldData(writer, q_field)
  call end_ncdef(writer)
  
  call write_GridData(writer)
  call write_FieldData(writer, l_fieldID, l_field)
  call write_FieldData(writer, q_fieldID, q_field)

  call close_ncFile(writer)
  
  call Field_IcGrid2D_Final(l_field)
  call Field_IcGrid2D_Final(q_field)
  call IcGrid_ncWriter_Final(writer)
  
end subroutine grid_test6

end program
 
