program TestDerivateManupulation

  use dc_test, only: AssertLessThan, AssertEqual

  use dc_types, only: &
    & DP

  use igmcore_math, only: &
    & PI

  use Derivate_Field_UniformGrid1D, only: &
    & delx_2ord_3pts_ufGrid, &
    & deldelx_2ord_3pts_ufGrid

  use Derivate_Field_UnuniGrid1D_Manager, only: &
    & Derivate_Field_UnuniGrid1D, &
    & Derivate_Field_UnuniGrid1D_Init, &
    & Derivate_Field_UnuniGrid1D_Final, &
    & delx_2ord_3pts_unufGrid, &
    & deldelx_2ord_3pts_unufGrid

  use VDerivate_Field_IcGrid3D_Manager, only: &
    & VDerivate_Field_IcGrid3D, &
    & VDerivate_Field_IcGrid3D_Init, &
    & VDerivate_Field_IcGrid3D_Final, &
    & d_dVCoord, dd_dVCoord

  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & IcGrid2D_FVM_Init, IcGrid2D_FVM_Final


  use IcGrid2D_FVM_Builder, only: &
    & construct_icosahedral_grid, GTYPE_STDGCGRID

  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & IcGrid3D_FVM_init, IcGrid3D_FVM_Final, &
    & get_icgrid2D, get_vertical_level_num, &
    & redistribute_verticalGrid

  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final

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

  implicit none

  integer, parameter :: gridNum = 41
  integer, parameter :: fieldRank = 1
  real(DP) :: grid1D(gridNum)
  real(DP) :: field1D(gridNum, fieldRank)
  real(DP) :: dfield1Ddz(gridNum, fieldRank)
  real(DP) :: ddfield1Ddz(gridNum, fieldRank)

  type(Derivate_Field_UnuniGrid1D) :: unufGridDelTool
  integer :: k
  real(DP) :: delz
  real(DP) :: delzAve
  real(DP) :: error(gridNum)

  type(IcGrid2D_FVM) :: icgrid2D
  type(IcGrid3D_FVM) :: icgrid3D

  type(VDerivate_Field_IcGrid3D) :: vderivate
  type(Field_IcGrid3D) :: field3Drank3
  type(Field_IcGrid2D), pointer :: field3Drank3_k
  type(Field_IcGrid3D) :: Delfield3Drank3
  type(Field_IcGrid2D), pointer :: Delfield3Drank3_k
  ! Unit test 5.0
  !
  ! 準備: 1 次元格子, 物理場(sin) の生成
  !
  write(*,*) '= Unit Test 5.0 [ Prepairation ]'
  write(*,*) '* Constructs 1D Grid'
  write(*,*) '* Generate 1D field described sine curve'

  grid1D(1) = 0.0d0
  delz = PI / real( gridNum - 1 )

  do k= 2, gridNum
    grid1D(k) = grid1D(k-1) + delz
  end do

!  write(*,*) 'grid vertex---'
!  write(*,*) grid1D(1:gridNum)

  do k=1, gridNum
    field1D(k,1) = sin(grid1D(k))
  end do

!  write(*,*) 'sin field---'
!  write(*,*) field1D(1:gridNum, 1)


  ! Unit test 5.1
  !
  ! 準備: 1 次元空間微分のテスト(等間隔格子)
  !
  write(*,*) '= Unit Test 5.1 '
  write(*,*) '* Check the spatial first differntial operation on 1D uniform grid.'
  write(*,*) '* Check the spatial second differntial operation on 1D uniform grid.'

  write(*,*) 'Derivate sin field (= cosine )---'

  dfield1Ddz(:,:) = delx_2ord_3pts_ufgrid(field1D, delz)

  error(:) = 0.0d0

  write(*,*) '[grid x] [d sin /dx ] [error]'
  do k=1, gridNum
    if ( abs(cos(grid1D(k))) > 1.0e-15 ) then
      error(k) = abs(cos(grid1D(k)) - dfield1Ddz(k, 1))
     write(*,*) grid1D(k), dfield1Ddz(k, 1), error(k), error(k) / abs(cos(grid1D(k)))
    else
     write(*,*) grid1D(k), dfield1Ddz(k, 1), '---'
    end if
  end do

  write(*,*) 'delZ**2=', delz**2, 'average error=', sum(error) / gridNum

  call AssertLessThan( &
      & message='Check the accuracy of 2-order first differential operation.', &
      & answer=delz**2, check=sum(error) / gridNum )

    write(*,*) 'Derivate sin field (= cosine )---'

  ddfield1Ddz(:,:) = deldelx_2ord_3pts_ufgrid(field1D, delz)

  error(:) = 0.0d0

  write(*,*) '[grid x] [d^2 sin /dx^2 ] [error]'
  do k=1, gridNum
    if ( abs(-sin(grid1D(k))) > 1.0e-15 ) then
      error(k) = abs(-sin(grid1D(k)) - ddfield1Ddz(k, 1))
      write(*,*) grid1D(k), ddfield1Ddz(k, 1), error(k), error(k) / abs(sin(grid1D(k)))
    else
      write(*,*) grid1D(k), ddfield1Ddz(k, 1), '---'
    end if
  end do

  write(*,*) 'delZ**2=', delz**2, 'average error=', sum(error) / gridNum

  call AssertLessThan( &
      & message='Check the accuracy of 2-order second differential operation.', &
      & answer=delz**2, check=sum(error) / gridNum )


  ! Unit test 5.2
  !
  ! 準備: 1 次元空間微分のテスト(非等間隔格子)
  !
  write(*,*) '= Unit Test 5.2 '
  write(*,*) '* Check the spatial first differntial operation on 1D ununiform grid.'
  write(*,*) '* Check the spatial second differntial operation on 1D ununiform grid.'

  !
  !
  write(*,*) 'Derivate cos field (= -sine )---'


  Delz = (0.5d0 * PI) /  real( gridNum - 1 )
  grid1D(1) = 0.0d0
  do k=2, gridNum
    grid1D(k) = grid1D(k-1) + Delz
  end do

  grid1D(:) = PI /2.0d0 * ( sqrt(grid1D(:)) / sqrt(grid1D(gridNum)) )
  field1D(:,1) = cos(grid1D(:))


  call Derivate_Field_UnuniGrid1D_Init(unufGridDelTool, grid1D)

  dfield1Ddz(:,:) = delx_2ord_3pts_unufGrid(unufGridDelTool, field1D)

  error(:) = 0.0d0

  write(*,*) '[grid x] [d cos/dx ] [error]'
  do k=1, gridNum
    if ( abs(-sin(grid1D(k))) > 1.0e-15 ) then
      error(k) = abs(-sin(grid1D(k)) - dfield1Ddz(k, 1))
      write(*,*) grid1D(k), dfield1Ddz(k, 1), error(k)
    else
      write(*,*) grid1D(k), dfield1Ddz(k, 1), '---'
    end if
  end do

  write(*,*) 'delZ(Max)**2=', (grid1D(2) - grid1D(1))**2, &
    &         'average error=', sum(error) / gridNum

  call AssertLessThan( &
      & message='Check the accuracy of 2-order first differential operation.', &
      & answer=(grid1D(2) - grid1D(1))**2, check=sum(error) / gridNum )

  !
  !
  write(*,*) 'Second Derivate sin field (= - sine )---'


  ddfield1Ddz(:,:) = deldelx_2ord_3pts_unufGrid(unufGridDelTool, field1D)

  error(:) = 0.0d0

  write(*,*) '[grid x] [d^2 cos /dx^2 ] [error]'
  do k=1, gridNum
    if ( abs(-cos(grid1D(k))) > 1.0e-15 ) then
      error(k) = abs(-cos(grid1D(k)) - ddfield1Ddz(k, 1))
      write(*,*) grid1D(k), ddfield1Ddz(k, 1), error(k)
    else
      write(*,*) grid1D(k), ddfield1Ddz(k, 1), '---'
    end if
  end do

  write(*,*) 'delZ(Max)**2=', (grid1D(2) - grid1D(1))**2, &
    &         'average error=', sum(error) / gridNum

  call AssertLessThan( &
      & message='Check the accuracy of 2-order first differential operation.', &
      & answer=(grid1D(2) - grid1D(1))**2, check=sum(error) / gridNum )


  ddfield1Ddz(:,:) = deldelx_2ord_3pts_unufGrid(unufGridDelTool, field1D, dfield1Ddz)

  error(:) = 0.0d0

  write(*,*) '[grid x] [d^2 cos /dx^2 ] [error]'
  do k=1, gridNum
    if ( abs(-cos(grid1D(k))) > 1.0e-15 ) then
      error(k) = abs(-cos(grid1D(k)) - ddfield1Ddz(k, 1))
      write(*,*) grid1D(k), ddfield1Ddz(k, 1), error(k)
    else
      write(*,*) grid1D(k), ddfield1Ddz(k, 1), '---'
    end if
  end do

  write(*,*) 'delZ(Max)**2=', (grid1D(2) - grid1D(1))**2, &
    &         'average error=', sum(error) / gridNum

  call AssertLessThan( &
      & message='Check the accuracy of 2-order first differential operation.', &
      & answer=(grid1D(2) - grid1D(1))**2, check=sum(error) / gridNum )


  call Derivate_Field_UnuniGrid1D_Final(unufGridDelTool)

  ! Unit test 5.3
  !
  ! 準備: 3 次元正二十面格子上の物理場の鉛直微分のテスト(非等間隔格子)
  !
  write(*,*) '= Unit Test 5.3 '
  write(*,*) '* Check the vertical differential operation for the physical field on 3D icosahedral grid.'

  write(*,*) 'initialize 3D icosahedral grid and field....'

  call construct_icosahedral_grid( &
    & grid_type=GTYPE_STDGCGRID, glevel=5, ic_radius=1.0d0, &
    & icgrid=icgrid2D )

  call IcGrid3D_FVM_Init(icgrid3D, icgrid2D, gridNum)
  call Field_IcGrid3D_Init(field3Drank3, icgrid3D, 'field3Drank3', 3)
  call Field_IcGrid3D_Init(Delfield3Drank3, icgrid3D, 'Delfield3Drank3', 3)

  ! Generate uniform grid
  grid1D(1) = 0.0d0
  delz = 1.1d0 * PI / real( gridNum - 1 )

  do k= 2, gridNum
    grid1D(k) = grid1D(k-1) + delz
  end do
  icgrid3D%verticalGrid(:) = grid1D(:)

  ! Set the value to field.  
  do k=1, gridNum
    field3Drank3_k => get_field_layer(field3Drank3, k)
    field3Drank3_k%val(:,:,:,:) = cos(grid1D(k))
  end do

  write(*,*) 'Vertical derivative....'
  call VDerivate_Field_IcGrid3D_Init(vderivate, icgrid3D)

  ! field1D is cos(grid1D(:)).
  call d_dVCoord(vderivate, field3Drank3, Delfield3Drank3)

  write(*,*) '[cal value] [answer d cos(x) /dx ] -------------------------'
  do k=1, gridNum
    write(*,*) Delfield3Drank3%field2Ds(k)%val(1,3,3,5), -sin(grid1D(k)), &
      & abs(Delfield3Drank3%field2Ds(k)%val(1,3,3,5) - (-sin(grid1D(k)) )) / abs(sin(grid1D(k)))
  end do  

  ! field1D is cos(grid1D(:)).
  call dd_dVCoord(vderivate, field3Drank3, Delfield3Drank3)

  write(*,*) '[cal value] [answer d^2 cos(x) / dx^2] -------------------------'
  do k=1, gridNum
    write(*,*) Delfield3Drank3%field2Ds(k)%val(1,3,3,5), -cos(grid1D(k)), &
      & abs(Delfield3Drank3%field2Ds(k)%val(1,3,3,5) - (-cos(grid1D(k)) )) / abs(cos(grid1D(k)))
  end do

  call Field_IcGrid3D_Final(field3Drank3)
  call Field_IcGrid3D_Final(Delfield3Drank3)
  call VDerivate_Field_IcGrid3D_Final(vderivate)
  call IcGrid3D_FVM_Final(icgrid3D)
  call IcGrid2D_FVM_Final(icgrid2D)

end program TestDerivateManupulation
