program test_TerrainFCoordinate

  ! モジュール引用 ; Use statements
  !
  use dc_types, only: &
    & DP, TOKEN, STRING

  use dc_test, only: AssertLessThan, AssertEqual

  use igmcore_math, only: &
    & PI, toRadians
    
  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos

  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & IcGrid2D_FVM_Init, IcGrid2D_FVM_Final, &
    & get_EffSize_Min, get_EffSize_Max, &
    & RC_REGIONS_NUM

  use IcGrid2D_FVM_Builder, only: &
    & construct_icosahedral_grid, &
    & GTYPE_STDGRID, GTYPE_STDGCGRID

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

  ! 1 次元物理場データの管理
  ! Manager of 1D physical field data.
  !
  use Field_Grid1D_Manager, only: &
    & Field_Grid1D, &
    & Field_Grid1D_Init, Field_Grid1D_Final
  
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & paste_field2D_margin

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

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

  use TFCoord_IcGrid3D_helper, only: &
    & TFCoord_IcGrid3D, &
    & TFCoord_IcGrid3D_Init, TFCoord_IcGrid3D_Final, &
    & TF_to_Z_IcGrid3D, calc_metric_G_1div2, calc_metric_G_Z, &
    & global_mean_Field_TFIcGrid3Dr, &
    & interpolate_FieldZ1D_into_FieldTF3D

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

  use Earth_Constants, only: &
    & earth_radius

  use Field_Pattern_Builder, only: &
    & create_cosine_bell_field

  use IcGrid_ncGridLoader, only: load_gridData

  implicit none

  type(IcGrid2D_FVM), pointer :: icgrid2D
  integer :: glevel = 5
  real(DP) :: ic_radius = earth_radius

  type(IcGrid3D_FVM) :: icgrid3D
  integer, parameter :: vertical_lv_num = 20

  type(TFCoord_IcGrid3D) :: TFCoordinate

  real(DP) :: vinterval_rate(vertical_lv_num-1)
  type(Field_IcGrid2D) :: Z_s
  real(DP) :: Z_T = 1.0e04

  type(Derivate_Field_IcGrid2D) :: derivate
  character(STRING) :: outputfile_name = 'metrics_data.nc'
  
  type(Field_IcGrid3D) :: test_field
  real(DP) :: gm_test_field(1)
  real(DP) :: gm_check
  
  type(Field_Grid1D) :: test_field1D
  real(DP) :: grid1D(vertical_lv_num)
  integer :: LvID, rcID, i, j
  integer :: glEMin, glEMax
  real(DP) :: z, ip_error, s_pos(3)

  ! Unit test 4.0
  !
  ! 準備: icgrid3D の生成, 地表高度場の生成, 微分演算オブジェクトの生成
  !
  write(*,*) '= Unit Test4.0 [ Prepairation ]'
  write(*,*) '* Constructs IcGrid3D(grid type: STD-GC-grid)'
  write(*,*) '* Generates surface height'
  write(*,*) '* Generates the object for differntial operations'

  ! IcGrid 3D の生成
  !
  allocate(icgrid2D)
  call construct_icosahedral_grid(GTYPE_STDGCGRID, glevel, ic_radius, icgrid2D)
 !call load_gridData('../../../IGModel_tool/src/icgridgen/STD-grid_glevel5.nc', icgrid2D)

  call IcGrid3D_FVM_Init(icgrid3D, icgrid2D, vertical_lv_num)

  ! 地表高度場の生成
  !

  call Field_IcGrid2D_Init( &
    & self=Z_s, &
    & icgrid=get_icgrid2D(icgrid3D), name='hs', rank=1, &
    & long_name='surface height', units='m' &
    & )

  ! コサインベル型の高度場を生成する.
  ! * 中心 : (東経 0 度, 北緯 0 度)
  ! * 半径 : 球の半径の 1/2
  ! * 最大高度 : 大気上端高度の 1/10

  write(*,*) 'location=(', toRadians(0d0), ',', toRadians(0d0), ')'
  write(*,*) 'cbR=', ic_radius/4.0d0
  write(*,*) 'h0=', Z_T/10.0d0

  call create_cosine_bell_field( &
    & field=Z_s, &
    & lambda_c=toRadians(0.0d0), theta_c=toRadians(0.0d0), &
    & cbR=ic_radius/2.0d0, h0=Z_T/10.0d0 &
    & )

  call paste_field2D_margin(Z_s)


  write(*,*) '--------------------------'

  ! 微分演算オブジェクトの生成
  call Derivate_Field_IcGrid2D_Init(derivate, get_icgrid2D(icgrid3D))


  ! Unit test 4.1
  !
  ! Terrain following coordinate の生成, メトリックの計算
  !
  write(*,*) '= Unit Test4.1 [ TFCoord_IcGrid3D_helper ]'
  write(*,*) '* Constructs terrain following coordinate'
  write(*,*) '* Calculates metrics G^1/2, G^Z'

  vinterval_rate(:) = 1.0d0
  call redistribute_verticalGrid(icgrid3D, 0.0d0, Z_T, vinterval_rate(:))
  write(*,*) '** terrain following coordinate'
  write(*,*) icgrid3D%verticalGrid

  call TFCoord_IcGrid3D_Init(TFCoordinate, icgrid3D, Z_s, Z_T)
  call calc_metric_G_1div2(TFCoordinate)
  call calc_metric_G_Z(TFCoordinate, derivate)

  write(*,*) '--------------------------'

  ! Unit test 4.2
  !
  ! 計算したメトリックのデータを NetCDF に出力する
  !
  write(*,*) '= Unit Test4.2 [ TFCoord_IcGrid3D_helper ]'
  write(*,*) '* Ouputs the calculated metrics data intot NetCDF file.'

  call output_metrics(TFCoordinate%G_1div2, TFCoordinate%G_Z)

  write(*,*) '--------------------------'


  ! Unit test 4.3
  !
  ! 
  !
  write(*,*) '= Unit Test4.3 [ TFCoord_IcGrid3D_helper ]'
  write(*,*) '* Calculate the global mean of the physical field on 3D icosahedral grid.'

  call Field_IcGrid3D_Init(test_field, icgrid3D, 'test_field', 1)

	call fill_Value(test_field, 2.0d0)
	gm_test_field =  global_mean_Field_TFIcGrid3Dr(TFCoordinate, test_field, 1)
	
	write(*,*) 'gm', gm_test_field
	gm_check = ( gm_test_field(1) - 2.0d0 ) / 2.0d0
      
	call AssertLessThan( &
     & message='Check the global mean value of physical field on 3D icosahedral grid' &
      	  & // '(applying the terrain following coordinate as vertical grid).', &
     & answer=1.0d-12, check=abs(gm_check) &
    & )
      	
  write(*,*) '--------------------------'

  ! Unit test 4.4
  !
  ! 
  !
  write(*,*) '= Unit Test4.4 [ TFCoord_IcGrid3D_helper ]'
  write(*,*) '* Check a subroutine to interpolate 1D field data on Z coordinate into 3D field data on terrain following coordinate.'

  !
  do LvID=1, vertical_lv_num
    grid1D(lvID) = Z_T / real( vertical_lv_num - 1 ) * ( LvID - 1 )
  end do
  
  call Field_Grid1D_Init(test_field1D, grid1D, 'test_field1D', 1)
  
  do LvID=1, vertical_lv_num
    test_field1D%val(LvID, 1) = sin(2d0 * PI * grid1D(lvID) / Z_T)
  end do
  
  !
  call interpolate_FieldZ1D_into_FieldTF3D(TFCoordinate, test_field1D, test_field)
  
  !
  glEMin = get_EffSize_Min(icgrid2D)
  glEMax = get_EffSize_Max(icgrid2D)
  
  ip_error = 0d0
  
  do rcID=1, RC_REGIONS_NUM
    do j=glEMin, glEMax
      do i=glEMin, glEMax
        s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:, i,j,rcID)) 

        do LvID=1, vertical_lv_num

          z = TF_to_Z_IcGrid3D(TFCoordinate, grid1D(LvID), rcID, i, j)
          ip_error = ip_error + abs( test_field%field2Ds(lvID)%val(1, i,j,rcID) - sin(2d0 * PI * z / Z_T) ) / 1.0d0
          if (i==1 .and. j==10 .and. rcID == 2) then
          !  write(*,*) LvID, ':', test_field%field2Ds(lvID)%val(1, i,j,rcID),  sin(2d0 * PI * z / Z_T)
          end if
          
          if ( abs(s_pos(1)) < 1.0e-02 .and. abs(s_pos(2)) < 1.0e-02 ) then
            write(*,*) LvID, ':', test_field%field2Ds(lvID)%val(1, i,j,rcID),  sin(2d0 * PI * z / Z_T)
          end if
          
        end do
      end do
    end do
  end do
  
  ip_error = ip_error / ( vertical_lv_num * RC_REGIONS_NUM * ( glEMax - glEMin + 1 )**2 ) * 100d0
  write(*,*) 'Error of interpolated field data(%): ', ip_error * 100d0
	call AssertLessThan( &
     & message='Check the accuracy of interpolated field on terrain following coordinate.' &
      	  & // '(applying the terrain following coordinate as vertical grid).', &
     & answer=0.01d0, check=abs(ip_error) &
    & )
      
  call Field_Grid1D_Final(test_field1D)

  ! 後片付け
  !
      
	call Field_IcGrid3D_Final(test_field)
  call TFCoord_IcGrid3D_Final(TFCoordinate)
  call Derivate_Field_IcGrid2D_Final(derivate)
  call Field_IcGrid2D_Final(Z_s)
  call IcGrid3D_FVM_Final(icgrid3D)
  call IcGrid2D_FVM_Final(icgrid2D)

contains

subroutine output_metrics(G_1div2, G_Z)

  !
  !
  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_vec, orth_to_geo_pos

  use IcGrid_ncWriter_mod, only: &
    & IcGrid_ncWriter, &
    & IcGrid_ncWriter_Init, &
    & open_ncFile, close_ncFile, &
    & end_ncdef, ncdef_GridData, ncdef_FieldData, &
    & write_GridData, write_FieldData, &
    & ncdef_Simulation_Parameter

  !
  !
  type(Field_IcGrid2D), intent(in) :: G_1div2
  type(Field_IcGrid3D), intent(in) :: G_Z

  !
  !
  !
  type(IcGrid_ncWriter) :: writer

  integer :: Z_s_ID
  integer :: G_1div2_ID
  integer :: G_Zonal_ID
  integer :: G_Merio_ID
  integer :: G_Rad_ID
  integer :: G_Zonal_k_ID

  type(Field_IcGrid3D) :: G_Z_Zonal
  type(Field_IcGrid3D) :: G_Z_Merio
  type(Field_IcGrid3D) :: G_Z_Rad

  type(Field_IcGrid2D), pointer :: G_Z_Zonal_k
  type(Field_IcGrid2D), pointer :: G_Z_Merio_k
  type(Field_IcGrid2D), pointer :: G_Z_Rad_k
  type(Field_IcGrid2D), pointer :: G_Z_k

  integer :: vt_num_max
  integer :: verticalLvID
  real(DP) :: s_pos(3)
  integer :: EMin, EMax
  integer :: rcID, i, j
  real(DP) :: G_Z_krij(3)
  !
  !

  call Field_IcGrid3D_Init(G_Z_Zonal, icgrid3D, 'G_Z_Zonal', 1, 'Zonal component of G_Z', '1')
  call Field_IcGrid3D_Init(G_Z_Merio, icgrid3D, 'G_Z_Merio', 1, 'Meridional component of G_Z', '1')
  call Field_IcGrid3D_Init(G_Z_Rad, icgrid3D, 'G_Z_Rad', 1, 'Vertical component of G_Z', '1')

  vt_num_max = get_vertical_level_num(icgrid3D)
  EMin = get_EffSize_Min(icgrid2D)
  EMax = get_EffSize_Max(icgrid2D)

  do verticalLvID=1, vt_num_max
    G_Z_k => get_field_layer(G_Z, verticalLvID)
    G_Z_Zonal_k => get_field_layer(G_Z_Zonal, verticalLvID)
    G_Z_Merio_k => get_field_layer(G_Z_Merio, verticalLvID)
    G_Z_Rad_k => get_field_layer(G_Z_Rad, verticalLvID)

    do rcID= 1, RC_REGIONS_NUM
      do j=EMin, EMax
        do i=EMin, EMax
          s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:, i,j,rcID))

          G_Z_krij(:) = orth_to_geo_vec(G_Z_k%val(:,i,j,rcID), s_pos)
          G_Z_Zonal_k%val(1, i,j,rcID) =  G_Z_krij(1)
          G_Z_Merio_k%val(1, i,j,rcID) =  G_Z_krij(2)
          G_Z_Rad_k%val(1, i,j,rcID) =  G_Z_krij(3)

        end do
      end do
    end do
  end do

  G_Z_Zonal_k => get_field_layer(G_Z_Zonal,1)
  !
  !

  call IcGrid_ncWriter_Init(writer, icgrid3D)

  call open_ncFile(writer, outputfile_name)

  call ncdef_GridData(writer)
  call ncdef_Simulation_Parameter(writer, 1.0d0, 1.0d0, 1.0d0)

  Z_s_ID = ncdef_FieldData(writer, Z_s)
  G_1div2_ID = ncdef_FieldData(writer, G_1div2)
  G_Zonal_ID = ncdef_FieldData(writer, G_Z_Zonal)
  G_Merio_ID = ncdef_FieldData(writer, G_Z_Merio)
  G_Rad_ID = ncdef_FieldData(writer, G_Z_Rad)
  G_Zonal_k_ID = ncdef_FieldData(writer, G_Z_Zonal_k)

  call end_ncdef(writer)

  call write_GridData(writer)
  call write_FieldData(writer, Z_s_ID, Z_s)
  call write_FieldData(writer, G_1div2_ID, G_1div2)
  call write_FieldData(writer, G_Zonal_ID, G_Z_Zonal)
  call write_FieldData(writer, G_Zonal_k_ID, G_Z_Zonal_k)
  call write_FieldData(writer, G_Merio_ID, G_Z_Merio)
  call write_FieldData(writer, G_Rad_ID, G_Z_Rad)

  call close_ncFile(writer)

end subroutine output_metrics

end program test_TerrainFCoordinate
