!> IGMBaseLib の IcGrid_ncWriter クラスを用いて netCDF ファイルに書き出された正二十面体格子(非構造格子)上の物理量データを, 
!! 地理座標格子(構造格子)上で補間し, netCDF ファイルとして書きだすプログラム.   
!!
!! 地理座標格子への補間アルゴリズムには, Willmott, et al(1985) を用いる. 
!!
!!
!> igmconvert のメインプログラム. 
!!
program igmconvert_main

  ! モジュール引用 ; Use statement
  ! 

  !* gtool5 ****
  !*

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP,    & ! 倍精度実数型. Double precision. 
    &                 TOKEN, &
    &                 STRING

  ! メッセージ出力
  ! Message dump
  !
  use dc_message, only: MessageNotify

  !* IGMBaseLib ****
  !*

  ! 地球の物理定数. 
  !
  !
  use Earth_Constants, only: earth_radius

  ! 正二十面体格子データを管理するクラス
  ! Manager of Icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & IcGrid2D_FVM_Init, malloc_GP_array

  ! 正二十面体格子上の物理場データの管理
  ! Manager of physical file data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, get_field2D_name

  ! 正二十面体格子データを管理するクラス
  ! Manager of Icosahedral grid data
  !
  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & IcGrid3D_FVM_Init, &
    & get_IcGrid3D_VerticalCoord_long_name, &
    & get_IcGrid3D_VerticalCoord_units

  ! 正二十面体格子上の物理場データの管理
  ! Manager of physical file data on icosahedral grid
  !
  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & Field_IcGrid3D_Init, get_field3D_name

  ! 正二十面体格子上の物理場データを地理座標系の構造格子上のデータとして補間するモジュール
  !
  !
  use geographic_mesh_interpolation, only: &
    & geographic_mesh_interpolation_Init, &
    & interpolate_Field_IcGrid2D, interpolate_Field_IcGrid3D, &
    & get_longitudes, get_latitudes
  
  !* igmconvert ****
  !*

  !
  !
  !
  use IcGrid_Field_Reader, only: &
    & IcGrid_Field_Reader_Init, IcGrid_Field_Reader_Final, &
    & read_NcHeaderIcGridData, read_Field_IcGrid2d_Data, read_Field_IcGrid3d_Data

  !
  !
  !
  use Interpolated_Field_Writer, only: &
    & prepair_output_ncfile, &
    & write_interpolated_fields2D, &
    & write_interpolated_fields3D

  ! 宣言文 ; Declaration statements
  !
  implicit none

  !> 正二十面体格子データの水平分割レベル. 
  integer :: glevel=5

  !>
  integer :: verticalLv_num = 1

  !> 正二十面体格子を内包する球の半径. 
  real(DP) :: ic_radius = earth_radius

  !>
  integer :: lon_num =  145

  !>
  integer :: lat_num =  73

  !>
  integer :: time_level_num = 145


  !> 読み込まれる netCDF ファイルの名前.
  !! The name of netCDF file to read from.
  character(STRING) :: read_ncFileName = 'output_data.nc'

  !> 補間完了後に物理場データが書き込まれる netCDF ファイルの名前. 
  !! The name of netCDF file to write to after interpolating. 
  character(STRING) :: written_ncFileName = 'interpolated_data.nc'
  
  !> 構造型 IcGrid2D_FVM の変数.
  type(IcGrid2D_FVM) :: icgrid2D
  
  !> 構造型 IcGrid3D_FVM の変数.
  type(IcGrid3D_FVM) :: icgrid3D
  type(IcGrid2D_FVM), pointer :: p_icgrid2D


  !> 変換対象となる物理場の種類の総数. 
  integer :: field2D_num = 0
  character(TOKEN), pointer :: field2Ds_name(:)
  type(Field_IcGrid2D), pointer :: fields_icgrid2D(:)

  integer :: field3D_num = 0
  character(TOKEN), pointer :: field3Ds_name(:)
  type(Field_IcGrid3D), pointer :: fields_icgrid3D(:)

  !
  real(DP) :: start_iptime = 0.0d0
  integer :: start_otstep
  real(DP) :: end_iptime = -1.0d0
  integer :: end_otstep
  integer :: iptstep_interval = 1.0d0
  real(DP) :: integration_time = 600d0
  real(DP) :: output_tick = 3000d0
  real(DP) :: time_step = 300d0

  ! igmconvert の実行モード
  !

  integer, parameter :: ICGRID2D_MODE = 1
  integer, parameter :: ICGRID3D_MODE = 2
  integer :: interpolation_mode = ICGRID2D_MODE

  ! 作業変数
  ! Work variables
  !
  integer :: fieldID
  integer :: tenpercent_timelevel
  integer :: timelevelID

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! 実行文 ; Executable statement
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !
  ! コマンドライン引数を解析し, パラメータを設定する.
  !

  call analize_cmdline_arg( &
    & read_ncFileName, written_ncFileName, &
    & glevel, &
    & verticalLv_num, &
    & interpolation_mode, field2D_num, field3D_num, &
    & start_iptime, end_iptime, iptstep_interval,   &
    & lon_num, lat_num )

  !
  !


  ! 構造型 IcGrid2D_FVM の変数を初期化する.
  ! 正二十面体格子の座標を保持する配列のメモリを確保するメンバー関数を呼ぶ. 
  call IcGrid2D_FVM_Init(icgrid2D, glevel, ic_radius)
  call malloc_GP_array(icgrid2D)

  if ( interpolation_mode == ICGRID2D_MODE ) then

    ! icgrid2D を引数として, IcGrid_Field_Reader モジュールを初期化する.
    call IcGrid_Field_Reader_Init(read_ncFileName, icgrid2D)


    ! NetCDF ファイルから, 正二十面体格子の格子点座標(icgrid2D)を読み取る.
    ! また, 出力間隔・積分時間も取得する.
    !write(*,*) 'Reading icosahedral grid data from netCDF file `', &
    !  & trim(read_ncFileName), '`.. '
    call read_NcHeaderIcGridData(output_tick, integration_time, time_step)


    ! icgrid2D を引数として, grgeographic_mesh_interpolation モジュールを初期化する.
    call geographic_mesh_interpolation_Init(lon_num, lat_num, field2D_num, icgrid2D)

  else if ( interpolation_mode == ICGRID3D_MODE ) then

    !
    call IcGrid3D_FVM_Init(icgrid3D, icgrid2D, verticalLv_num)

    !
    !

    write(*,*) 'Reading icosahedral grid 3D data from netCDF file `', &
      & trim(read_ncFileName), '`.. '


    ! icgrid3D を引数として, IcGrid_Field_Reader モジュールを初期化する.
    call IcGrid_Field_Reader_Init(read_ncFileName, icgrid3D)

    ! NetCDF ファイルから, 正二十面体格子の格子点座標(icgrid3D)を読み取る.
    ! また, 出力間隔・積分時間も取得する.
    call read_NcHeaderIcGridData(output_tick, integration_time, time_step)

    ! icgrid3D を引数として, grgeographic_mesh_interpolation モジュールを初期化する.
    call geographic_mesh_interpolation_Init( &
      & lon_num, lat_num, &
      & field2D_num, field3D_num, &
      & icgrid3D &
      & )

  end if

  ! コマンドライン引数と NetCDF から読み込んだシミュレーションパラメータをもとにして,
  ! 作成された設定を出力する.
  !
  
  ! 補間終了時刻が設定されていない場合は, 補間終了時刻として積分時間を指定する.
  if ( end_iptime == -1.0d0 ) then
    end_iptime = integration_time
  end if

  ! 時間情報をスッテプ数に変換する.
  start_otstep = int(start_iptime / output_tick)
  end_otstep = int(end_iptime / output_tick)
  time_level_num = int( ( end_otstep - start_otstep + 1 ) / iptstep_interval )

  call MessageNotify('M', 'igmconvert', &
    & 'The following configuration will be used.' &
    & )

  write(*,*) '* Horizontal division level(glevel)=', glevel
  write(*,*) '* Number of interpolated fields=', field2D_num
  write(*,*) '* Start time=', start_iptime, ', Start time level=', start_otstep
  write(*,*) '* End time=', end_iptime, ', End time level=',end_otstep
  write(*,*) '* Interpolated interval of time level=', iptstep_interval
  write(*,*) '* Number of interpolated time level=', time_level_num

  ! NetCDF から読み込んだ正二十面格子上の物理場データを保持する構造型 Field_IcGrid2D の配列変数
  ! を初期化する.

  !
  if ( field2D_num > 0 ) then
    write(*,*) 'field2D_num=', field2D_num!, time_level_num
    write(*,*) field2Ds_name(:)

    allocate(fields_icgrid2D(field2D_num))

    do fieldID=1,field2D_num
      call Field_IcGrid2D_Init(fields_icgrid2D(fieldID), icgrid2D, trim(field2Ds_name(fieldID)), 1)
    end do

    ! 補完した物理場データの書き出しに必要な情報を取得するために,
    ! 一度初期のタイムレベルのデータを NetCDF から読み込み, Field2D の物理場情報を収集してくる,
    call read_Field_IcGrid2d_Data( &
      & time_level=1, fields=fields_icgrid2D)
  else
      allocate(fields_icgrid2D(0))
  end if

  !
  if ( field3D_num > 0 ) then
    write(*,*) 'field3D_num=', field3D_num!, time_level_num
    write(*,*) field3Ds_name(:)

    allocate(fields_icgrid3D(field3D_num))

    do fieldID=1,field3D_num
      call Field_IcGrid3D_Init(fields_icgrid3D(fieldID), icgrid3D, trim(field3Ds_name(fieldID)), 1)
    end do

    ! 補完した物理場データの書き出しに必要な情報を取得するために,
    ! 一度初期のタイムレベルのデータを NetCDF から読み込み, Field2D の物理場情報を収集してくる,
    call read_Field_IcGrid3d_Data( &
      & time_level=1, fields=fields_icgrid3D )
  else
    allocate(fields_icgrid3D(0))
  end if

  
  ! 補間データを NetCDF ファイルとして, 書きだすための前処理を行う.
  !

  ! 補完データを gtool を使って, NetCDF に書き出すための準備を行う.
  if ( interpolation_mode == ICGRID2D_MODE ) then
    call prepair_output_ncfile( &
      & file_name=written_ncFileName, &
      & lons=get_longitudes(), lats=get_latitudes(), &
      & field2Ds=fields_icgrid2D,     &
      & tinit=start_iptime, output_tick=output_tick*iptstep_interval &
      & )

   else if ( interpolation_mode == ICGRID3D_MODE ) then

      call prepair_output_ncfile( &
        & file_name=written_ncFileName, &
        & lons=get_longitudes(), lats=get_latitudes(), &
        & levels=icgrid3D%verticalGrid(:), &
        & field2Ds=fields_icgrid2D, field3Ds=fields_icgrid3D, &
        & tinit=start_iptime, output_tick=output_tick*iptstep_interval, &
        & vertical_axis_name=get_IcGrid3D_VerticalCoord_long_name(icgrid3D), &
        & vertical_unit=get_IcGrid3D_VerticalCoord_units(icgrid3D) &
        & )

   end if
  
  tenpercent_timelevel = time_level_num * 0.1
  
  ! 実際に補完処理によって, 正二十面格子上の物理場データを緯度経度格子上のデータに変換する.
  !

  do timelevelID = start_otstep+1, end_otstep+1, iptstep_interval

    ! 進捗状況の表示
    !
    if( tenpercent_timelevel /= 0 .and. &
      & mod(timelevelID, tenpercent_timelevel) == 0 ) then
 
      write(*,*) '...', (timelevelID-1)*iptstep_interval, '/', time_level_num
    end if
  
    ! 時間レベル timelevelID の対象となる正二十面格子上の物理場データを読み込み,
    ! 緯度経度格子上への補間処理を行う.
    !
    
    if ( interpolation_mode == ICGRID2D_MODE ) then

      !
      call read_Field_IcGrid2d_Data(timelevelID, fields_icgrid2D)
      !
      call write_interpolated_fields2D( &
        & fields_icgrid2D, interpolate_Field_IcGrid2D(fields_icgrid2D) &
        & )

    else if ( interpolation_mode == ICGRID3D_MODE ) then

      if ( field2D_num > 0 ) then
        !
        !
        call read_Field_IcGrid2d_Data(timelevelID, fields_icgrid2D)

        call write_interpolated_fields2D( &
          & fields_icgrid2D, interpolate_Field_IcGrid2D(fields_icgrid2D) &
          & )

      end if

      if ( field3D_num > 0 ) then
        !
        !
        !
        call read_Field_IcGrid3d_Data(timelevelID, fields_icgrid3D)

        call write_interpolated_fields3D( &
          & fields_icgrid3D, interpolate_Field_IcGrid3D(fields_icgrid3D) &
          & )

      end if

    end if

  end do

  call IcGrid_Field_Reader_Final()

contains

!
!>
!!
!!
subroutine analize_cmdline_arg( &
  & read_ncFileName, output_ncFileName, &
  & g_level, vertical_level_num, &
  & interpolationMode, & ! (inout)
  & field2D_num, &
  & field3D_num, &
  & start_iptime, end_iptime, iptstep_interval, & ! (inout)
  & nlons, nlats )  ! (inout)

  ! モジュール引用 ; Use statement
  !

  ! 文字列操作
  ! String manupulation
  !
  use dc_string, only: &
    & StoA, StoI, StoD, toChar, Split, Replace

  ! コマンドライン引数の解析
  !
  !
  use dc_args, only: ARGS, &
    & DCArgsOpen, DCArgsClose, DCArgsOption, &
    & DCArgsPutLine, DCArgsDebug, DCArgsHelp, DCArgsStrict 

  ! 宣言文 ; Declaration statement
  !
  character(*), intent(inout) :: read_ncFileName
  character(*), intent(inout) :: output_ncFileName
  integer, intent(inout) :: g_level
  integer, intent(inout) :: vertical_level_num
  integer, intent(inout) :: interpolationMode
  integer, intent(inout) :: field2D_num
  integer, intent(inout) :: field3D_num
  real(DP), intent(inout) :: start_iptime
  real(DP), intent(inout) :: end_iptime
  integer, intent(inout) :: iptstep_interval
  integer, intent(inout) :: nlons
  integer, intent(inout) :: nlats

  ! 作業変数
  ! Work variables
  !

  !
  type(ARGS) :: arg

  logical :: OPT_read_ncfname
  character(STRING) :: VAL_read_ncfname

  logical :: OPT_output_ncfname
  character(STRING) :: VAL_output_ncfname

  logical :: OPT_interpolationMode
  character(TOKEN) :: VAL_interpolationMode

  logical :: OPT_glevel
  character(TOKEN) :: VAL_glevel

  logical :: OPT_verticalLvNum
  character(TOKEN) :: VAL_verticalLvNum

  logical :: OPT_field2Ds_name
  character(STRING) :: VAL_field2Ds_name

  logical :: OPT_field3Ds_name
  character(STRING) :: VAL_field3Ds_name

  character(TOKEN), pointer :: p_fields_name(:)
  integer :: fieldID

  logical :: OPT_output_tl_info
  character(TOKEN) :: VAL_output_tl_info
  character(TOKEN), pointer :: p_tl_info(:)
  integer :: tl_info_num

  logical :: OPT_nlonlats
  character(STRING) :: VAL_nlonlats
  character(TOKEN), pointer :: p_nlonlats(:)


  ! 実行文; Executable statement
  !

  !
  call DCArgsOpen(arg)
  

  ! 読み込みファイルの名前を指定するためのオプション
  call DCArgsOption(arg, StoA('-R', '--Read_NetCDF'), OPT_read_ncfname, VAL_read_ncfname, &
    &    help="Specify NetCDF file output by IGModel. [defalut value is " // trim(read_ncFileName) // "]." )
  
  ! 書き出しファイルの名前を指定するためのオプション
  call DCArgsOption(arg, StoA('-O', '--Output_NetCDF'), OPT_output_ncfname, VAL_output_ncfname, &
    &    help="Specify NetCDF filename output by igmconvert. [defalut value is " // trim(output_ncFileName) // "]." )
  
  ! 補間モード指定のオプション
  call DCArgsOption(arg, StoA('-M', '--IcGrid_Mode'), OPT_interpolationMode, VAL_interpolationMode, &
    &    help="Specify the interpolation mode etheir 'IcGrid3D' or 'IcGrid2D'. [defalut value is 'IcGrid2D']." )

  ! glevel を指定するためのオプション
  call DCArgsOption(arg, StoA('-G', '--glevel'), OPT_glevel, VAL_glevel, &
    &    help="Specify division level of icosahedral grid. [defalut value is " // trim(toChar(glevel)) // "]." )

  ! 鉛直分解レベルを指定するためのオプション
  call DCArgsOption(arg, StoA('-V', '--vertical_level_num'), OPT_verticalLvNum, VAL_verticalLvNum, &
    &    help="Specify the number of vertical level. [defalut value is " // trim(toChar(vertical_level_num)) // "]." )

  ! 
  call DCArgsOption(arg, StoA('-F', '--Field2Ds_Name'), OPT_field2Ds_name, VAL_field2Ds_name, &
    &    help="Specify the name of physical 2D fields converting to  geographic mesh format. [example: 'temp, u, v']." )

  ! 
  call DCArgsOption(arg, StoA('-f', '--Field3Ds_Name'), OPT_field3Ds_name, VAL_field3Ds_name, &
    &    help="Specify the name of physical 3D fields converting to  geographic mesh format. [example: 'temp, u, v']." )

  !
  call DCArgsOption(arg, StoA('-T', '--Time_levels'), OPT_output_tl_info, VAL_output_tl_info, &
    & help="Specify the time at which the fields will be interpolated. &
    & For example, -T='<start time>,<end time>,<interpolated interval of output time level>'. &
    & [default value is -T='0,integration_time,1']." )

  ! 
  call DCArgsOption(arg, StoA('-N', '--Num_lon_lat'), OPT_nlonlats, VAL_nlonlats, &
    &    help="Specify the number of grid point in latitude and longitude. &
    & [default value is '"// trim(toChar(nlats)) // ',' // trim(toChar(nlons)) // "']." )


  ! デバッグ, ヘルプ, 厳密な引数指定のための設定を行う.
  call DCArgsDebug(arg)
  call DCArgsHelp(arg)
  call DCArgsStrict(arg)

  if( OPT_read_ncfname ) then
    read_ncFileName = VAL_read_ncfname
  end if

  if( OPT_output_ncfname ) then
    output_ncFileName = VAL_output_ncfname
  end if

  if ( OPT_interpolationMode ) then
    if ( VAL_interpolationMode == 'IcGrid3D' ) then
      interpolationMode = ICGRID3D_MODE

      call MessageNotify('M', 'igmconvert', &
        &  'Switch to IcGrid3D interpolation mode.' )

    end if
  end if

  if( OPT_glevel ) then
    g_level = StoI(VAL_glevel)
  end if

  if( OPT_verticalLvNum ) then
    vertical_level_num = StoI(VAL_verticalLvNum)
  end if

  ! 2D 物理場
  if ( OPT_field2Ds_name ) then
    call Split(Replace(VAL_field2Ds_name, ' ', '', .true.), p_fields_name, ',')

    field2D_num = size(p_fields_name)
    allocate(field2Ds_name(field2D_num))
    
    do fieldID=1, field2D_num
      field2Ds_name(fieldID) = trim(p_fields_name(fieldID))
    end do

    deallocate(p_fields_name)
  end if
  
  ! 3D 物理場
  if ( OPT_field3Ds_name ) then
    call Split(Replace(VAL_field3Ds_name, ' ', '', .true.), p_fields_name, ',')

    field3D_num = size(p_fields_name)
    allocate(field3Ds_name(field3D_num))

    do fieldID=1, field3D_num
      field3Ds_name(fieldID) = trim(p_fields_name(fieldID))
    end do

    deallocate(p_fields_name)
  end if

  if ( OPT_output_tl_info ) then
    call Split(Replace(VAL_output_tl_info, ' ', '', .true.), p_tl_info, ',')

    tl_info_num = size(p_tl_info)

    start_iptime = StoD(p_tl_info(1))
    end_iptime = StoD(p_tl_info(2))

    if ( tl_info_num == 3 ) then
      iptstep_interval = StoI(p_tl_info(3))
    else
      iptstep_interval = 1
    end if

  else
    start_iptime = 0.0d0
    end_iptime = -1.0d0
    iptstep_interval = 1
  end if

  if ( OPT_nlonlats ) then
    call Split(Replace(VAL_nlonlats, ' ', '', .true.), p_nlonlats, ',')
    
    if (size(p_nlonlats) == 2 ) then
      nlons = StoI(p_nlonlats(1))
      nlats = StoI(p_nlonlats(2))
    else
      write(*,*) 'Specifing option is invalid' 
    end if
    deallocate(p_nlonlats)
  end if

end subroutine analize_cmdline_arg


end program igmconvert_main
