!>
!!
!!
!
program unit_test1
  
  use dc_types, only: DP, TOKEN, STRING

  use igmcore_math, only: PI

  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & IcGrid2D_FVM_Init, &
    & get_EffSize_Min, get_EffSize_Max, RC_REGIONS_NUM
  
  use IcGrid2D_FVM_Builder, only: &
    & construct_icosahedral_grid, GTYPE_STDGCGRID

  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, &
    & paste_field2D_margin, get_field2D_name

  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & Derivate_Field_IcGrid2D_Init, &
    & divergence_op

  use Field_Pattern_Builder, only: create_solid_rotation_field
  
  use rcregion_vtkwriter


  implicit none

  type(IcGrid2D_FVM) :: icgrid
  integer, parameter :: glevel = 5
  real(DP), parameter :: ic_radius = 1.0d0

  type(Field_IcGrid2D) :: v_vec
  type(Field_IcGrid2D) :: div_v
  type(Field_IcGrid2D) :: true_div_v
  
  type(Derivate_Field_IcGrid2D) :: diff_eval

  real(DP) :: angular_speed
  real(DP), parameter :: alpha = 0.0d0!PI/2.0d0 


  !***************
  !
  !***************

  !
  write(*,*) 'Constructs STD-GC-grid(glevel=5, ic_radius=1.0d0)'
  call construct_icosahedral_grid(GTYPE_STDGCGRID, glevel, ic_radius, icgrid)

  !
  !
  call Field_IcGrid2D_Init(v_vec, icgrid, 'solid_rotation_v', 3)
  call Field_IcGrid2D_Init(div_v, icgrid, 'div_v', 1)

  angular_speed = 2.0d0 * PI / ( 3600d0 * 24d0 * 12d0 )
  call create_solid_rotation_field(v_vec, angular_speed, alpha)
  call paste_field2D_margin(v_vec)
 
  !
  call Derivate_Field_IcGrid2D_Init(diff_eval, icgrid)
  call divergence_op(diff_eval, v_vec, div_v)

  !
  call write_field_as_vtkfile()

contains 
subroutine write_field_as_vtkfile()
  integer :: EMin, EMax
  integer :: rcID, i,j
  integer ::deviceID = 8

  real(DP), allocatable :: rc_CVArea(:,:) 

  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  allocate(rc_CVArea(EMin:EMax,EMin:EMax))
  
  do rcID=1,RC_REGIONS_NUM
     call init_rcregion_vtkwriter(deviceID, rcID, './', '00')
     call begin_rcregion_grid_data((EMax-EMin+1)**2)
     call write_rcregion_grid_data(icgrid%rcs_AGrid(:,EMin:EMax,EMin:EMax,rcID))

     call begin_point_data((EMax-EMin+1)**2)
     call write_point_data_scalars(trim(get_field2D_name(div_v)), div_v%val(1,EMin:EMax,EMin:EMax,rcID))
     
     do j=EMin, EMax
       do i=EMin, EMax
         rc_CVArea(i,j) = icgrid%rcs_CVInfo(i,j,rcID)%cv_area
       end do
     end do

     call write_point_data_scalars('CVArea', rc_CVArea(:,:))

     call close_rcregion_vtkwriter()
  end do

end subroutine write_field_as_vtkfile

end program unit_test1
