program merge

  use vtype_module
  use fi_module
  use ni3_module

  implicit none

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

  integer               :: inncid
  integer, allocatable  :: a_outncid(:)

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

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

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

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

  integer            :: jmax_max
  integer            :: tmax
  real(DP), allocatable :: a_Time (:)
  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_output_file_opened
  logical            :: flag_inc_time_dim

  integer            :: jh, jph, js, je, jch, jc

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

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


  ctlfn = 'split.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 )


  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

    ncfn = inncfn
    mode = "read"
    call ni3_open( ncfn, mode, inncid )

    ! Inquire variable
    !
    call ni3_inq_var( inncid, varname, ndims = ndims, xtype = xtype )
    allocate( a_dimnames( ndims ) )
    call ni3_inq_vardimnames( inncid, varname, ndims, a_dimnames )
    allocate( a_dimlen( ndims ) )
    do l = 1, ndims
      call ni3_inq_dimlen( inncid, 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) ) )
      call ni3_get_var( inncid, a_dimnames(l), a_Time )
    end if


    ! read latitude
    !

    l = 2
    allocate( y_Lat_global( a_dimlen(l) ) )
    call ni3_get_var( inncid, a_dimnames(l), y_Lat_global )

    l = 2
    allocate( y_Lat_weight_global( a_dimlen(l) ) )
    if ( ni3_chk_var( inncid, 'lat_weight' ) ) then
      call ni3_get_var( inncid, 'lat_weight', y_Lat_weight_global )
    else
      y_Lat_weight_global = 0.0d0
    end if

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


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



    jh = jmax_global / 2
    jph = ( jh - 1 ) / nprocs + 1

    jmax_max = jph * 2

    ! store latitude
    !
    allocate( a_jmax( 0:nprocs-1 ) )
    allocate( aa_Lat( jmax_max, 0:nprocs-1 ) )
    allocate( aa_Lat_index ( jmax_max, 0:nprocs-1 ) )
    allocate( aa_Lat_weight( jmax_max, 0:nprocs-1 ) )
    !
    do n = 0, nprocs-1

      js = jph * n + 1
      je = min( jph * ( n + 1 ), jh )
      if ( je >= js ) THEN
        jch = je - js + 1
        jc  = jch * 2
      else
        jc = 0
        js = 1
        je = 1
      end if

      a_jmax(n) = jc

      do j = 1, a_jmax(n) / 2
        aa_Lat       (a_jmax(n)/2+1+(j-1),n) = y_Lat_global(jh+1+(js-1)+(j-1))
        aa_Lat       (a_jmax(n)/2  -(j-1),n) = y_Lat_global(jh  -(js-1)-(j-1))
        aa_Lat_index (a_jmax(n)/2+1+(j-1),n) = jh+1+(js-1)+(j-1)
        aa_Lat_index (a_jmax(n)/2  -(j-1),n) = jh  -(js-1)-(j-1)
        aa_Lat_weight(a_jmax(n)/2+1+(j-1),n) = y_Lat_weight_global(jh+1+(js-1)+(j-1))
        aa_Lat_weight(a_jmax(n)/2  -(j-1),n) = y_Lat_weight_global(jh  -(js-1)-(j-1))
      end do

    end do



    if ( .not. flag_output_file_opened ) then

      allocate( a_outncid(0:nprocs-1) )

      do n = 0, nprocs-1
        write( ncfn, '(a,a,i6.6,a)' ) outncfn(1:len_trim(outncfn)-3), &
          & '_rank', n, '.nc'

        call setupoutputfile( &
          & inncid, ncfn, &
          & a_jmax(n), aa_Lat(1:a_jmax(n), n), &
          & a_outncid(n) &
          & )

        if ( ni3_chk_var( inncid, 'lat_weight' ) ) then
          call putweight( &
            & inncid, a_outncid(n), &
            & a_jmax(n), aa_Lat_weight(1:a_jmax(n),n) &
            & )
        end if
      end do

      flag_output_file_opened = .true.
    end if


    ! define variable
    !
    do n = 0, nprocs-1
      call ni3_def_var( a_outncid(n), varname, xtype, ndims, a_dimnames )
      call ni3_cp_atts( inncid, a_outncid(n), varname )
    end do


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


    loop_time : do t = 1, tmax

      select case ( xtype )
      case ( NI3_INT )
        select case ( ndims_exc_time )
        case ( 2 )
          if ( flag_inc_time_dim ) then
            call ni3_get_varss( inncid, varname, t, xy_ivar_global )
          else
            call ni3_get_var  ( inncid, varname, xy_ivar_global )
          end if
        case ( 3 )
          if ( flag_inc_time_dim ) then
            call ni3_get_varss( inncid, varname, t, xyz_ivar_global )
          else
            call ni3_get_var  ( inncid, 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_get_varss( inncid, varname, t, xy_fvar_global )
          else
            call ni3_get_var  ( inncid, varname, xy_fvar_global )
          end if
        case ( 3 )
          if ( flag_inc_time_dim ) then
            call ni3_get_varss( inncid, varname, t, xyz_fvar_global )
          else
            call ni3_get_var  ( inncid, 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_get_varss( inncid, varname, t, xy_dvar_global )
          else
            call ni3_get_var  ( inncid, varname, xy_dvar_global )
          end if
        case ( 3 )
          if ( flag_inc_time_dim ) then
            call ni3_get_varss( inncid, varname, t, xyz_dvar_global )
          else
            call ni3_get_var  ( inncid, varname, xyz_dvar_global )
          end if
        end select
      end select




      loop_proc : do n = 0, nprocs-1

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


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


      end do loop_proc


    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( aa_Lat_index )
    deallocate( aa_Lat_weight )


    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
      do n = 0, nprocs-1
        call ni3_close( a_outncid(n) )
      end do
      deallocate( a_outncid )

      flag_output_file_opened = .false.
    end if

    deallocate( y_Lat_global )
    deallocate( y_Lat_weight_global )
    if ( flag_inc_time_dim ) then
      deallocate( a_Time  )
    end if
    deallocate( a_dimnames )
    deallocate( a_dimlen )


    call ni3_close( inncid  )

  end do loop_namelist

  if ( flag_output_file_opened ) then
    do n = 0, nprocs-1
      call ni3_close( a_outncid(n) )
    end do
    deallocate( a_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
  character(extstr)     :: varname
  integer :: flag_rst
  integer :: l
  integer :: t


  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 )
        if ( a_dimnames(l) == 'time' ) then
          call ni3_def_dim( outncid, a_dimnames(l), xtype, 0 )
          do t = 1, dimlen
            call ni3_put_varss( outncid, a_dimnames(l), t, a_idim( t ) )
          end do
        else
          call ni3_set_dim( outncid, a_dimnames(l), xtype, a_idim )
        end if
        deallocate( a_idim )
      case ( NI3_REAL )
        allocate( a_fdim( dimlen ) )
        call ni3_get_var( inncid, a_dimnames(l), a_fdim )
        if ( a_dimnames(l) == 'time' ) then
          call ni3_def_dim( outncid, a_dimnames(l), xtype, 0 )
          do t = 1, dimlen
            call ni3_put_varss( outncid, a_dimnames(l), t, a_fdim( t ) )
          end do
        else
          call ni3_set_dim( outncid, a_dimnames(l), xtype, a_fdim )
        end if
        deallocate( a_fdim )
      case ( NI3_DOUBLE )
        allocate( a_ddim( dimlen ) )
        call ni3_get_var( inncid, a_dimnames(l), a_ddim )
        if ( a_dimnames(l) == 'time' ) then
          call ni3_def_dim( outncid, a_dimnames(l), xtype, 0 )
          do t = 1, dimlen
            call ni3_put_varss( outncid, a_dimnames(l), t, a_ddim( t ) )
          end do
        else
          call ni3_set_dim( outncid, a_dimnames(l), xtype, a_ddim )
        end if
        deallocate( a_ddim )
      end select

    end if

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


  deallocate( a_dimnames )

  ! output flag_rst
  !
  do l = 1, nvars
    call ni3_inq_varname( inncid, l, varname )
    if ( varname == 'flag_rst' ) then
      call ni3_inq_var( inncid, varname, ndims = ndims, xtype = xtype )
      allocate( a_dimnames(0) )
      call ni3_def_var( outncid, varname, xtype, ndims, a_dimnames )
      call ni3_cp_atts( inncid, outncid, varname )
      call ni3_get_var( inncid , varname, flag_rst )
      call ni3_put_var( outncid, varname, flag_rst )
      deallocate( a_dimnames )
    end if
  end do


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(:)
  character(extstr)     :: name
  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 )

  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 )
      deallocate( a_dimnames )

      call ni3_cp_atts( inncid, outncid, name )

    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
