program merge

  use vtype_module
  use ni3_module
  use fi_module

  implicit none

  character(extstr)     :: ctlfn
  integer               :: ctlfu
  character(extstr)     :: mode
  integer               :: ios

  integer               :: nprocs
  integer               :: jmax_global
  integer , allocatable :: a_jmax(:)
  real(DP), allocatable :: aa_Lat(:,:)
  integer , allocatable :: aa_Lat_index(:,:)
  real(DP), allocatable :: y_Lat_global(:)
  real(DP), allocatable :: a_Time (:)

  real(DP), allocatable :: y_Lat_weight_global(:)
  real(DP), allocatable :: ya_Lat_weight(:,:)

  integer , allocatable :: a_inncid(:)
  integer               :: outncid

  character(extstr)  :: varname
  character(extstr)  :: inncfn
  character(extstr)  :: inncfn_b
  character(extstr)  :: outncfn
  character(extstr)  :: outncfn_b
  character(extstr)  :: ncfn

  integer            :: ndims
  integer            :: ndims_exc_time
  integer            :: xtype

  character(extstr), allocatable :: a_dimnames(:)
  integer          , allocatable :: a_dimlen(:)

  integer            :: jmax_max
  integer            :: tmax
  integer , allocatable :: xy_ivar_global (:,:)
  integer , allocatable :: xyz_ivar_global(:,:,:)
  real(SP), allocatable :: xy_fvar_global (:,:)
  real(SP), allocatable :: xyz_fvar_global(:,:,:)
  real(DP), allocatable :: xy_dvar_global (:,:)
  real(DP), allocatable :: xyz_dvar_global(:,:,:)
  integer , allocatable :: xy_ivar (:,:)
  integer , allocatable :: xyz_ivar(:,:,:)
  real(SP), allocatable :: xy_fvar (:,:)
  real(SP), allocatable :: xyz_fvar(:,:,:)
  real(DP), allocatable :: xy_dvar (:,:)
  real(DP), allocatable :: xyz_dvar(:,:,:)

  logical            :: flag_inc_time_dim
  logical            :: flag_output_file_opened


  integer            :: j, jj, k, l
  integer            :: n
  integer            :: t

  namelist /proc/ nprocs
  namelist /item/ inncfn, outncfn, varname


  ctlfn = 'merge.nml'
  mode  = 'read'
  call fi_open( ctlfn, mode, ctlfu )


  nprocs = 1
  rewind( ctlfu )
  read( ctlfu, nml = proc, iostat = ios )
  write( 6, * ) nprocs
  if ( ios /= 0 ) stop 'Unable to read namelist file.'
  write( 6, proc )


  allocate( a_inncid( 0:nprocs-1 ) )


  rewind( ctlfu )


  varname = 'ZZZ'
  inncfn  = '-----'
  outncfn = '-----'
  read( ctlfu, nml = item, iostat = ios )
  write( 6, item )

  if ( inncfn  == '-----' ) inncfn  = trim( varname ) // '.nc'
  if ( outncfn == '-----' ) outncfn = inncfn


  flag_output_file_opened = .false.

  loop_namelist : do

    do n = 0, nprocs-1
      write( ncfn, '(a,a,i6.6,a)' ) inncfn(1:len_trim(inncfn)-3), &
        & '_rank', n, '.nc'
      mode = "read"
      call ni3_open( ncfn, mode, a_inncid(n) )
    end do


    ! Inquire variable
    !
    n = 0
    call ni3_inq_var( a_inncid(n), varname, ndims = ndims, xtype = xtype )
    allocate( a_dimnames( ndims ) )
    call ni3_inq_vardimnames( a_inncid(n), varname, ndims, a_dimnames )
    allocate( a_dimlen( ndims ) )
    do l = 1, ndims
      call ni3_inq_dimlen( a_inncid(n), a_dimnames(l), a_dimlen(l) )
    end do


    ! check dimensions
    !
    if ( a_dimnames(1) /= 'lon' ) then
      write( 6, * ) '1st dimension is not longitude, but ', trim( a_dimnames(1) )
      stop
    end if
    if ( a_dimnames(2) /= 'lat' ) then
      write( 6, * ) '1st dimension is not longitude, but ', trim( a_dimnames(1) )
      stop
    end if
    if ( a_dimnames(ndims) == 'time' ) then
      flag_inc_time_dim = .true.
    else
      flag_inc_time_dim = .false.
    end if

    if ( flag_inc_time_dim ) then
      ndims_exc_time = ndims - 1
    else
      ndims_exc_time = ndims
    end if


    ! read time
    !
    if ( flag_inc_time_dim ) then
      l = ndims
      allocate( a_Time( a_dimlen(l) ) )
      n = 0
      call ni3_get_var( a_inncid(n), a_dimnames(l), a_Time )
    end if


    ! set length of time dimension
    !
    if ( flag_inc_time_dim ) then
      tmax = a_dimlen(ndims)
    else
      tmax = 1
    end if


    ! set latitude on whole globe
    !
    jmax_global = a_dimlen(2)



    ! Check number of latitudinal grids on the globe
    !
    allocate( a_jmax( 0:nprocs-1 ) )
    !
    jmax_global = 0
    jmax_max    = 0
    do n = 0, nprocs-1
      call ni3_inq_dimlen( a_inncid(n), 'lat', a_jmax(n) )
      jmax_global = jmax_global + a_jmax(n)
      if( a_jmax(n) > jmax_max ) jmax_max = a_jmax(n)
    end do
    !
    ! Store latitude
    !
    allocate( aa_Lat( jmax_max, 0:nprocs-1 ) )
    !
    do n = 0, nprocs-1
      call ni3_get_var( a_inncid(n), 'lat', aa_Lat(1:a_jmax(n), n ) )
    end do
    !
    ! Set latitudinal grid on whole globe by sorting
    !
    allocate( y_Lat_global( jmax_global ) )
    allocate( aa_Lat_index( jmax_max, 0:nprocs-1 ) )
    jj = 1
    do n = nprocs-1, 0, -1
      do j = 1, a_jmax(n)/2
        y_Lat_global(jj) = aa_Lat(j,n)
        aa_Lat_index(j,n) = jj
        jj = jj + 1
      end do
    end do
    do n = 0, nprocs-1
      do j = a_jmax(n)/2+1, a_jmax(n)
        y_Lat_global(jj) = aa_Lat(j,n)
        aa_Lat_index(j,n) = jj
        jj = jj + 1
      end do
    end do



    if ( .not. flag_output_file_opened ) then
      n = 0
      call setupoutputfile( &
        & a_inncid(n), outncfn, &
        & jmax_global, y_Lat_global, &
        & outncid &
        & )
      flag_output_file_opened = .true.

      n = 0
      if ( ni3_chk_var( a_inncid(n), 'lat_weight' ) ) then
        allocate( y_Lat_weight_global( jmax_global ) )

        allocate( ya_Lat_weight( jmax_max, 0:nprocs-1 ) )
        do n = 0, nprocs-1
          call ni3_get_var( a_inncid(n), 'lat_weight', ya_Lat_weight( 1:a_jmax(n),n ) )
        end do
        jj = 1
        do n = nprocs-1, 0, -1
          do j = 1, a_jmax(n)/2
            y_Lat_weight_global(jj) = ya_Lat_weight(j,n)
            jj = jj + 1
          end do
        end do
        do n = 0, nprocs-1
          do j = a_jmax(n)/2+1, a_jmax(n)
            y_Lat_weight_global(jj) = ya_Lat_weight(j,n)
            jj = jj + 1
          end do
        end do
        deallocate( ya_Lat_weight )

        n = 0
        call putweight( &
          & a_inncid(n), outncid, &
          & jmax_global, y_Lat_weight_global &
          & )

        deallocate( y_Lat_weight_global )
      end if
    end if
    !
    ! define variable
    !
    call ni3_def_var( outncid, varname, xtype, ndims, a_dimnames )
    n = 0
    call ni3_cp_atts( a_inncid(n), outncid, varname )

    select case ( xtype )
    case ( NI3_INT )
      select case ( ndims_exc_time )
      case ( 2 )
        allocate( xy_ivar_global ( a_dimlen(1), jmax_global ) )
      case ( 3 )
        allocate( xyz_ivar_global( a_dimlen(1), jmax_global, a_dimlen(3) ) )
      end select
    case ( NI3_REAL )
      select case ( ndims_exc_time )
      case ( 2 )
        allocate( xy_fvar_global ( a_dimlen(1), jmax_global ) )
      case ( 3 )
        allocate( xyz_fvar_global( a_dimlen(1), jmax_global, a_dimlen(3) ) )
      end select
    case ( NI3_DOUBLE )
      select case ( ndims_exc_time )
      case ( 2 )
        allocate( xy_dvar_global ( a_dimlen(1), jmax_global ) )
      case ( 3 )
        allocate( xyz_dvar_global( a_dimlen(1), jmax_global, a_dimlen(3) ) )
      end select
    end select


    loop_time : do t = 1, tmax

      loop_proc : do n = 0, nprocs-1


        select case ( xtype )
        case ( NI3_INT )
          select case ( ndims_exc_time )
          case ( 2 )
            allocate( xy_ivar( a_dimlen(1), a_jmax(n) ) )
            if ( flag_inc_time_dim ) then
              call ni3_get_varss( a_inncid(n), varname, t, xy_ivar )
            else
              call ni3_get_var  ( a_inncid(n), varname, xy_ivar )
            end if
            do j = 1, a_jmax(n)
              xy_ivar_global(:,aa_Lat_index(j,n)) = xy_ivar(:,j)
            end do
            deallocate( xy_ivar )
          case ( 3 )
            allocate( xyz_ivar( a_dimlen(1), a_jmax(n), a_dimlen(3) ) )
            if ( flag_inc_time_dim ) then
              call ni3_get_varss( a_inncid(n), varname, t, xyz_ivar )
            else
              call ni3_get_var  ( a_inncid(n), varname, xyz_ivar )
            end if
            do k = 1, a_dimlen(3)
              do j = 1, a_jmax(n)
                xyz_ivar_global(:,aa_Lat_index(j,n),k) = xyz_ivar(:,j,k)
              end do
            end do
            deallocate( xyz_ivar )
          end select
        case ( NI3_REAL )
          select case ( ndims_exc_time )
          case ( 2 )
            allocate( xy_fvar( a_dimlen(1), a_jmax(n) ) )
            if ( flag_inc_time_dim ) then
              call ni3_get_varss( a_inncid(n), varname, t, xy_fvar )
            else
              call ni3_get_var  ( a_inncid(n), varname, xy_fvar )
            end if
            do j = 1, a_jmax(n)
              xy_fvar_global(:,aa_Lat_index(j,n)) = xy_fvar(:,j)
            end do
            deallocate( xy_fvar )
          case ( 3 )
            allocate( xyz_fvar( a_dimlen(1), a_jmax(n), a_dimlen(3) ) )
            if ( flag_inc_time_dim ) then
              call ni3_get_varss( a_inncid(n), varname, t, xyz_fvar )
            else
              call ni3_get_var  ( a_inncid(n), varname, xyz_fvar )
            end if
            do k = 1, a_dimlen(3)
              do j = 1, a_jmax(n)
                xyz_fvar_global(:,aa_Lat_index(j,n),k) = xyz_fvar(:,j,k)
              end do
            end do
            deallocate( xyz_fvar )
          end select
        case ( NI3_DOUBLE )
          select case ( ndims_exc_time )
          case ( 2 )
            allocate( xy_dvar( a_dimlen(1), a_jmax(n) ) )
            if ( flag_inc_time_dim ) then
              call ni3_get_varss( a_inncid(n), varname, t, xy_dvar )
            else
              call ni3_get_var  ( a_inncid(n), varname, xy_dvar )
            end if
            do j = 1, a_jmax(n)
              xy_dvar_global(:,aa_Lat_index(j,n)) = xy_dvar(:,j)
            end do
            deallocate( xy_dvar )
          case ( 3 )
            allocate( xyz_dvar( a_dimlen(1), a_jmax(n), a_dimlen(3) ) )
            if ( flag_inc_time_dim ) then
              call ni3_get_varss( a_inncid(n), varname, t, xyz_dvar )
            else
              call ni3_get_var  ( a_inncid(n), varname, xyz_dvar )
            end if
            do k = 1, a_dimlen(3)
              do j = 1, a_jmax(n)
                xyz_dvar_global(:,aa_Lat_index(j,n),k) = xyz_dvar(:,j,k)
              end do
            end do
            deallocate( xyz_dvar )
          end select
        end select


      end do loop_proc

      if ( flag_inc_time_dim ) then
        call ni3_put_varss( outncid, 'time', t, a_Time(t) )
      end if

      select case ( xtype )
      case ( NI3_INT )
        select case ( ndims_exc_time )
        case ( 2 )
          if ( flag_inc_time_dim ) then
            call ni3_put_varss( outncid, varname, t, xy_ivar_global )
          else
            call ni3_put_var  ( outncid, varname, xy_ivar_global )
          end if
        case ( 3 )
          if ( flag_inc_time_dim ) then
            call ni3_put_varss( outncid, varname, t, xyz_ivar_global )
          else
            call ni3_put_var  ( outncid, varname, xyz_ivar_global )
          end if
        end select
      case ( NI3_REAL )
        select case ( ndims_exc_time )
        case ( 2 )
          if ( flag_inc_time_dim ) then
            call ni3_put_varss( outncid, varname, t, xy_fvar_global )
          else
            call ni3_put_var  ( outncid, varname, xy_fvar_global )
          end if
        case ( 3 )
          if ( flag_inc_time_dim ) then
            call ni3_put_varss( outncid, varname, t, xyz_fvar_global )
          else
            call ni3_put_var  ( outncid, varname, xyz_fvar_global )
          end if
        end select
      case ( NI3_DOUBLE )
        select case ( ndims_exc_time )
        case ( 2 )
          if ( flag_inc_time_dim ) then
            call ni3_put_varss( outncid, varname, t, xy_dvar_global )
          else
            call ni3_put_var  ( outncid, varname, xy_dvar_global )
          end if
        case ( 3 )
          if ( flag_inc_time_dim ) then
            call ni3_put_varss( outncid, varname, t, xyz_dvar_global )
          else
            call ni3_put_var  ( outncid, varname, xyz_dvar_global )
          end if
        end select
      end select


    end do loop_time


    select case ( xtype )
    case ( NI3_INT )
      select case ( ndims_exc_time )
      case ( 2 )
        deallocate( xy_ivar_global )
      case ( 3 )
        deallocate( xyz_ivar_global )
      end select
    case ( NI3_REAL )
      select case ( ndims_exc_time )
      case ( 2 )
        deallocate( xy_fvar_global )
      case ( 3 )
        deallocate( xyz_fvar_global )
      end select
    case ( NI3_DOUBLE )
      select case ( ndims_exc_time )
      case ( 2 )
        deallocate( xy_dvar_global )
      case ( 3 )
        deallocate( xyz_dvar_global )
      end select
    end select

    deallocate( a_jmax )

    deallocate( aa_Lat )
    deallocate( y_Lat_global )
    deallocate( aa_Lat_index )


    outncfn_b = outncfn
    inncfn_b  = inncfn

    varname = 'ZZZ'
    inncfn  = '-----'
    outncfn = '-----'
    read( ctlfu, nml = item, iostat = ios )
    if ( ios /= 0 ) exit

    write( 6, item )

    if ( inncfn  == '-----' ) inncfn  = trim( varname ) // '.nc'
    if ( outncfn == '-----' ) outncfn = inncfn

    if ( outncfn_b == outncfn ) then
      if ( inncfn_b /= inncfn ) then
        write( 6, * ) 'Input file has to be the same when the output file is the same.'
        write( 6, * ) trim( outncfn_b ), trim( outncfn )
        write( 6, * ) trim( inncfn_b ), trim( inncfn )
        stop
      end if
    else
      call ni3_close( outncid )
      flag_output_file_opened = .false.
    end if

    deallocate( a_Time  )

    deallocate( a_dimnames )
    deallocate( a_dimlen )

    do n = 0, nprocs-1
      call ni3_close( a_inncid(n) )
    end do

  end do loop_namelist

  if ( flag_output_file_opened ) then
    call ni3_close( outncid )
  end if

  close( ctlfu )

  stop
end program merge

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

subroutine setupoutputfile( &
  & inncid, outncfn, &
  & jmax, a_Lat, &
  & outncid &
  & )

  use vtype_module
  use ni3_module

  implicit none

  integer     , intent(in ) :: inncid
  character(*), intent(in ) :: outncfn
  integer     , intent(in ) :: jmax
  real(DP)    , intent(in ) :: a_Lat(jmax)
  integer     , intent(out) :: outncid


  character(extstr)     :: mode
  character(extstr), allocatable :: a_dimnames(:)
  integer               :: ndims, nvars, natts
  integer               :: xtype
  integer , allocatable :: a_idim(:)
  real(sp), allocatable :: a_fdim(:)
  real(dp), allocatable :: a_ddim(:)

  integer :: dimlen
  integer :: l



  call ni3_inq( inncid, ndims, nvars, natts )
  allocate( a_dimnames( ndims ) )
  call ni3_inq_dimnames( inncid, ndims, a_dimnames )


  mode = 'new'
  call ni3_open( outncfn, mode, outncid )

  call ni3_cp_atts( inncid, outncid, 'global' )

  do l = 1, ndims
    call ni3_inq_var( inncid, a_dimnames(l), xtype = xtype )
    call ni3_inq_dimlen( inncid, a_dimnames(l), dimlen )

    if ( a_dimnames(l) == 'lat' ) then

      call ni3_set_dim( outncid, a_dimnames(l), xtype, a_lat )

    else

      select case ( xtype )
      case ( NI3_INT )
        allocate( a_idim( dimlen ) )
        call ni3_get_var( inncid, a_dimnames(l), a_idim )
        call ni3_set_dim( outncid, a_dimnames(l), xtype, a_idim )
        deallocate( a_idim )
      case ( NI3_REAL )
        allocate( a_fdim( dimlen ) )
        call ni3_get_var( inncid, a_dimnames(l), a_fdim )
        call ni3_set_dim( outncid, a_dimnames(l), xtype, a_fdim )
        deallocate( a_fdim )
      case ( NI3_DOUBLE )
        allocate( a_ddim( dimlen ) )
        call ni3_get_var( inncid, a_dimnames(l), a_ddim )
        call ni3_set_dim( outncid, a_dimnames(l), xtype, a_ddim )
        deallocate( a_ddim )
      end select

    end if

    call ni3_cp_atts( inncid, outncid, a_dimnames(l) )
  end do

  deallocate( a_dimnames )


end subroutine setupoutputfile

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

subroutine putweight( &
  & inncid, outncid, &
  & jmax, a_Lat_weight &
  & )

  use vtype_module
  use ni3_module

  implicit none

  integer     , intent(in ) :: inncid
  integer     , intent(in ) :: outncid
  integer     , intent(in ) :: jmax
  real(DP)    , intent(in ) :: a_Lat_weight(jmax)


  character(extstr), allocatable :: a_dimnames(:)
  integer               :: ndims, nvars, natts
  integer               :: xtype
  integer , allocatable :: a_idim(:)
  real(sp), allocatable :: a_fdim(:)
  real(dp), allocatable :: a_ddim(:)
  character(extstr)     :: name

  integer :: dimlen
  integer :: l


  call ni3_inq( inncid, ndims, nvars, natts )


  do l = 1, nvars
    call ni3_inq_varname( inncid, l, name )

    if ( len_trim( name ) <= 6 ) cycle

    if ( name == 'lat_weight' ) then

      call ni3_inq_var( inncid, name, ndims = ndims, xtype = xtype )
      if ( ndims /= 1 ) then
        write( 6, * ) 'ndims has to be 1, but is, ', ndims, '.'
        stop
      end if

      allocate( a_dimnames(ndims) )
      call ni3_inq_vardimnames( inncid, name, ndims, a_dimnames )
      call ni3_def_var( outncid, name, xtype, ndims, a_dimnames )
      call ni3_put_var( outncid, name, a_Lat_weight )
      call ni3_cp_atts( inncid, outncid, name )
      deallocate( a_dimnames )

    else if ( name(len_trim(name)-6:len_trim(name)) == '_weight' ) then

      call ni3_inq_var( inncid, name, ndims = ndims, xtype = xtype )
      if ( ndims /= 1 ) then
        write( 6, * ) 'ndims has to be 1, but is, ', ndims, '.'
        stop
      end if

      allocate( a_dimnames(ndims) )
      call ni3_inq_vardimnames( inncid, name, ndims, a_dimnames )
      call ni3_inq_dimlen( inncid, a_dimnames(1), dimlen )
      call ni3_def_var( outncid, name, xtype, ndims, a_dimnames )

      select case ( xtype )
      case ( NI3_INT )
        allocate( a_idim( dimlen ) )
        call ni3_get_var( inncid, name, a_idim )
        call ni3_put_var( outncid, name, a_idim )
        deallocate( a_idim )
      case ( NI3_REAL )
        allocate( a_fdim( dimlen ) )
        call ni3_get_var( inncid, name, a_fdim )
        call ni3_put_var( outncid, name, a_fdim )
        deallocate( a_fdim )
      case ( NI3_DOUBLE )
        allocate( a_ddim( dimlen ) )
        call ni3_get_var( inncid, name, a_ddim )
        call ni3_put_var( outncid, name, a_ddim )
        deallocate( a_ddim )
      end select

      deallocate( a_dimnames )

      call ni3_cp_atts( inncid, outncid, name )
    end if

  end do


end subroutine putweight
