program spread

  use vtype_module
  use fi_module
  use ni3_module

  implicit none

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

  integer               :: inncid
  integer               :: outncid
  integer               :: inncid_ps

  integer                        :: ndims_all, nvars_all, natts_all
  character(extstr), allocatable :: a_dimnames_all(:)
  integer          , allocatable :: a_dimlen_all(:)
  integer          , allocatable :: a_xtype_all(:)

  integer                        :: imax
  integer                        :: jmax
  integer                        :: kmax_sig
  integer                        :: kmax_sigm
  integer                        :: kmax_in
  integer                        :: kmax_out
  integer                        :: tmax

  integer                        :: dimlen
  integer            :: ndims
  integer            :: xtype

  character(extstr)  :: dimname
  character(extstr)  :: varname
  character(extstr)  :: inncfn
  character(extstr)  :: outncfn
  character(extstr)  :: inncfn_ps
  character(extstr)  :: varname_ps
  character(extstr)  :: ncfn

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

  real(sp), allocatable :: z_Sigma(:)
  real(sp), allocatable :: r_Sigma(:)
  real(sp), allocatable :: a_Sigma(:)

  real(DP), allocatable :: xy_Ps    (:,:)
  real(DP), allocatable :: xya_Press(:,:,:)

  integer , allocatable :: a_iVar(:)
  real(sp), allocatable :: a_fVar(:)
  real(dp), allocatable :: a_dVar(:)

!!$  integer               :: iVar
!!$  integer , allocatable :: xy_iVar (:,:)
!!$  integer , allocatable :: xyz_iVar(:,:,:)
!!$  real(SP), allocatable :: xy_fVar (:,:)
!!$  real(SP), allocatable :: xyz_fVar    (:,:,:)
!!$  real(SP), allocatable :: xyz_fVar_out(:,:,:)
!!$  real(DP), allocatable :: xy_dVar     (:,:)
  real(DP), allocatable :: xyz_dVar    (:,:,:)
  real(DP), allocatable :: xyz_dVar_out(:,:,:)


  logical            :: flag_dim
  logical            :: flag_weight
  logical            :: flag_rst
  logical            :: flag_output

  integer            :: varid
  integer            :: dimid


  integer         , parameter      :: plevmax    = 150
!!$  integer         , parameter      :: d_plevn    =  35
  integer         , parameter      :: d_plevn    =  17
  integer         , parameter      :: d_plevempn = plevmax - d_plevn
  integer                          :: plevn
  real(DP)                         :: plevel( plevmax )
  logical                          :: flaglog

  real(DP)        , parameter      :: missing_value = -999.0d0

  logical           :: FlagIncPlanetLon
  character(extstr) :: inplncfn
  integer           :: inplncid
  character(extstr) :: plvarname
  real(DP)          :: PlanetLon


  integer            :: k, l
  integer            :: t
  integer            :: ts
  integer            :: te


  data plevn     /d_plevn/
!!$  data plevel &
!!$    & / 1.0d-6,   3.0d-6,   5.0d-6,   7.0d-6,   &
!!$    &   1.0d-5,   3.0d-5,   5.0d-5,   7.0d-5,   &
!!$    &   1.0d-4,   3.0d-4,   5.0d-4,   7.0d-4,   &
!!$    &   1.0d-3,   3.0d-3,   5.0d-3,   7.0d-3,   &
!!$    &   1.0d-2,   3.0d-2,   5.0d-2,   7.0d-2,   &
!!$    &   1.0d-1,   2.0d-1,   3.0d-1,   5.0d-1,   7.0d-1,   &
!!$    &   1.0d0 ,   2.0d0 ,   3.0d0 ,   4.0d0 ,   5.0d0 ,   &
!!$    &   6.0d0 ,   7.0d0 ,   8.0d0 ,   9.0d0 ,   1.0d1 ,   &
!!$    &   d_plevempn*0.0d0 /
  data plevel &
    & / 1000d2, 925d2, 850d2, 700d2, 600d2, 500d2, 400d2, 300d2, 250d2, 200d2, 150d2, 100d2, 70d2, 50d2, 30d2, 20d2, 10d2, &
    &   d_plevempn*0.0d0 /


  namelist /file/        inncfn, outncfn
  namelist /file_ps/     inncfn_ps, varname_ps
  namelist /out_vlev_p/  plevn, plevel, flaglog
  namelist /time/        ts, te
  namelist /PlanetLongitude/   FlagIncPlanetLon, inplncfn, plvarname


  inncfn_ps  = '-----'
  varname_ps = 'Ps'
  ts         = -1
  te         = -1

  FlagIncPlanetLon = .false.
  inplncfn         = '------'
  plvarname        = '------'

  flaglog    = .false.

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

  rewind( ctlfu )
  read( ctlfu, nml = out_vlev_p, iostat = ios )
  write( 6, out_vlev_p )
  !-----
  rewind( ctlfu )
  read( ctlfu, nml = file_ps, iostat = ios )
  if ( ios /= 0 ) stop 'Unable to read namelist file.'
  write( 6, file_ps )
  !-----
  rewind( ctlfu )
  read( ctlfu, nml = time, iostat = ios )
  write( 6, time )
  !-----
  rewind( ctlfu )
  read( ctlfu, nml = PlanetLongitude, iostat = ios )
  write( 6, PlanetLongitude )

  if ( inncfn_ps == '-----' ) stop 'inncfn_ps has to be given.'

  rewind( ctlfu )

  loop_file : do

    inncfn    = '-----'
    outncfn   = '-----'

    read( ctlfu, nml = file, iostat = ios )
    if ( ios /= 0 ) exit
    write( 6, file )

    if ( inncfn    == '-----' ) stop 'inncfn has to be given.'
    if ( outncfn   == '-----' ) stop 'outncfn has to be given.'


    ncfn = inncfn_ps
    mode = "read"
    call ni3_open( ncfn, mode, inncid_ps )

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

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



    call ni3_inq( inncid, ndims_all, nvars_all, natts_all )
    allocate( a_dimnames_all( ndims_all ) )
    call ni3_inq_dimnames( inncid, ndims_all, a_dimnames_all )

    allocate( a_xtype_all ( ndims_all ) )
    allocate( a_dimlen_all( ndims_all ) )

    ! Check length of dimensions
    !
    do l = 1, ndims_all
      call ni3_inq_var( inncid, a_dimnames_all(l), xtype = a_xtype_all(l) )
      call ni3_inq_dimlen( inncid, a_dimnames_all(l), a_dimlen_all(l) )
    end do

    imax      = -1
    jmax      = -1
    kmax_sig  = -1
    kmax_sigm = -1
    tmax      = -1
    do l = 1, ndims_all
      if ( a_dimnames_all(l) == 'lon' ) then
        imax      = a_dimlen_all(l)
      else if ( a_dimnames_all(l) == 'lat' ) then
        jmax      = a_dimlen_all(l)
      else if ( a_dimnames_all(l) == 'sig' ) then
        kmax_sig  = a_dimlen_all(l)
      else if ( a_dimnames_all(l) == 'sigm' ) then
        kmax_sigm = a_dimlen_all(l)
      else if ( a_dimnames_all(l) == 'time' ) then
        tmax      = a_dimlen_all(l)
        if ( ts == -1  ) ts = 1
        if ( ts <=  0  ) ts = 1
        if ( te == -1  ) te = tmax
        if ( te <=  0  ) te = 1
        if ( te > tmax ) te = tmax
      end if
    end do
    if ( imax      == -1 ) then
      stop 'longitude cannot be found'
    end if
    if ( jmax      == -1 ) then
      stop 'latitude cannot be found'
    end if
    if ( kmax_sig  == -1 ) then
      stop 'sig cannot be found'
    end if
    if ( kmax_sigm == -1 ) then
      stop 'sigm cannot be found'
    end if
    if ( tmax      == -1 ) then
      stop 'time cannot be found'
    end if

    kmax_out = plevn


    ! Output global attributes
    !

    call ni3_cp_atts( inncid, outncid, 'global' )


    ! Output dimensions
    !

    do l = 1, ndims_all

      dimlen  = a_dimlen_all(l)
      dimname = a_dimnames_all(l)
      xtype   = a_xtype_all(l)

      if ( ( dimname == 'sig' ) .or. ( dimname == 'sigm' ) ) cycle

      select case ( xtype )
      case ( NI3_INT )
        allocate( a_iVar( dimlen ) )
        call ni3_get_var( inncid, dimname, a_iVar )
        if ( dimname == 'time' ) then
          call ni3_def_dim( outncid, dimname, xtype, 0 )
!!$          do t = 1, dimlen
          do t = ts, te
            call ni3_put_varss( outncid, dimname, t-ts+1, a_iVar( t ) )
          end do
        else
          call ni3_set_dim( outncid, dimname, xtype, a_iVar )
        end if
        deallocate( a_iVar )
      case ( NI3_REAL )
        allocate( a_fVar( dimlen ) )
        call ni3_get_var( inncid, dimname, a_fVar )
        if ( dimname == 'time' ) then
          call ni3_def_dim( outncid, dimname, xtype, 0 )
!!$          do t = 1, dimlen
          do t = ts, te
            call ni3_put_varss( outncid, dimname, t-ts+1, a_fVar( t ) )
          end do
        else
          call ni3_set_dim( outncid, dimname, xtype, a_fVar )
        end if
        deallocate( a_fVar )
      case ( NI3_DOUBLE )
        allocate( a_dVar( dimlen ) )
        call ni3_get_var( inncid, dimname, a_dVar )
        if ( dimname == 'time' ) then
          call ni3_def_dim( outncid, dimname, xtype, 0 )
!!$          do t = 1, dimlen
          do t = ts, te
            call ni3_put_varss( outncid, dimname, t-ts+1, a_dVar( t ) )
          end do
        else
          call ni3_set_dim( outncid, dimname, xtype, a_dVar )
        end if
        deallocate( a_dVar )
      end select

      call ni3_cp_atts( inncid, outncid, dimname )
    end do


    ! Output weights
    !

    do l = 1, ndims_all

      dimlen  = a_dimlen_all(l)
      dimname = a_dimnames_all(l)

      if ( ( dimname == 'lon' ) .or. ( dimname == 'lat' ) ) then

        varname = trim( dimname ) // '_weight'

        call ni3_inq_var( inncid, varname, ndims = ndims, xtype = xtype )
        allocate( a_dimnames( ndims ) )
        call ni3_inq_vardimnames( inncid, varname, ndims, a_dimnames )
        call ni3_def_var( outncid, varname, xtype, ndims, a_dimnames )
        deallocate( a_dimnames )

        select case ( xtype )
        case ( NI3_INT )
          allocate( a_iVar( dimlen ) )
          call ni3_get_var( inncid , varname, a_iVar )
          call ni3_put_var( outncid, varname, a_iVar )
          deallocate( a_iVar )
        case ( NI3_REAL )
          allocate( a_fVar( dimlen ) )
          call ni3_get_var( inncid , varname, a_fVar )
          call ni3_put_var( outncid, varname, a_fVar )
          deallocate( a_fVar )
        case ( NI3_DOUBLE )
          allocate( a_dVar( dimlen ) )
          call ni3_get_var( inncid , varname, a_dVar )
          call ni3_put_var( outncid, varname, a_dVar )
          deallocate( a_dVar )
        end select

        call ni3_cp_atts( inncid, outncid, varname )

      end if
    end do


    ! Output pressure level
    !

    xtype = NI3_REAL
    dimname = 'level'
    if ( flaglog ) then
      call ni3_set_dim( outncid, dimname, xtype, log10( plevel(1:plevn) ) )
    else
      call ni3_set_dim( outncid, dimname, xtype, plevel(1:plevn) )
    end if
    call ni3_put_att( outncid, dimname, 'units', 'Pa' )
    call ni3_put_att( outncid, dimname, 'long_name', 'Level' )
    call ni3_put_att( outncid, dimname, 'positive', 'down' )


    ! Obtain sigma levels
    !

    allocate( z_Sigma( kmax_sig  ) )
    allocate( r_Sigma( kmax_sigm ) )
    call ni3_get_var( inncid, 'sig' , z_Sigma )
    call ni3_get_var( inncid, 'sigm', r_Sigma )



    ! Inquire variable
    !

    loop_variables : do varid = 1, nvars_all

      call ni3_inq_varname( inncid, varid, varname )


!!$      write( 6, * ) varid, ':', trim( varname )

      flag_dim = .false.
      do dimid = 1, ndims_all
        dimname = a_dimnames_all(dimid)
        if ( dimname == varname ) then
          flag_dim = .true.
        end if
      end do

      flag_weight = .false.
      do dimid = 1, ndims_all
        dimname = a_dimnames_all(dimid)
        if ( trim( dimname ) // '_weight' == varname ) then
          flag_weight = .true.
        end if
      end do


      flag_output = .true.
      flag_rst    = .false.


      if ( ( flag_dim ) .or. ( flag_weight ) ) then

        flag_output = .false.

      else

        call ni3_inq_var( inncid, varname, ndims = ndims, xtype = xtype )

        if ( ndims >= 4 ) then

          allocate( a_dimnames( ndims ) )
          call ni3_inq_vardimnames( inncid, varname, ndims, a_dimnames )


          ! check dimensions
          !
          if ( a_dimnames(1) /= 'lon' ) then
            write( 6, * ) '      1st dimension is not lon, but ', trim( a_dimnames(1) )
            flag_output = .false.
          end if
          if ( a_dimnames(2) /= 'lat' ) then
            write( 6, * ) '      2nd dimension is not lat, but ', trim( a_dimnames(2) )
            flag_output = .false.
          end if
          if ( ( a_dimnames(3) /= 'sig' ) .and. ( a_dimnames(3) /= 'sigm' ) ) then
            write( 6, * ) '      3rd dimension is not sig or sigm, but ', trim( a_dimnames(3) )
            flag_output = .false.
          end if
          if ( a_dimnames(4) /= 'time' ) then
            write( 6, * ) '      4th dimension is not time, but ', trim( a_dimnames(4) )
            flag_output = .false.
          end if

          if ( .not. flag_output ) then
            deallocate( a_dimnames )
          end if

          ! Lines below would be unnecessary.
        else

!!$          if ( varname == 'flag_rst' ) then
!!$            ! define variable
!!$            !
!!$            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, ivar )
!!$            call ni3_put_var( outncid, varname, ivar )
!!$            deallocate( a_dimnames )
!!$
!!$            flag_rst = .true.
!!$          end if

          flag_output = .false.

        end if

      end if

      if ( flag_output ) then

        if ( a_dimnames(3) == 'sig' ) then
          kmax_in = kmax_sig
        else if ( a_dimnames(3) == 'sigm' ) then
          kmax_in = kmax_sigm
        else
          write( 6, * ) 'Unexpected vertical dimension, ', trim( a_dimnames(3) )
          stop
        end if
        allocate( a_Sigma     ( kmax_in ) )
        allocate( xyz_dVar    ( imax, jmax, kmax_in ) )
        allocate( xy_Ps       ( imax, jmax ) )
        allocate( xya_Press   ( imax, jmax, kmax_in ) )
        allocate( xyz_dVar_out( imax, jmax, kmax_out ) )
        if ( a_dimnames(3) == 'sig' ) then
          a_Sigma = z_Sigma
        else if ( a_dimnames(3) == 'sigm' ) then
          a_Sigma = r_Sigma
        else
          write( 6, * ) 'Unexpected vertical dimension, ', trim( a_dimnames(3) )
          stop
        end if


        ! define variable
        !
        a_dimnames(3) = 'level'
        call ni3_def_var( outncid, varname, xtype, ndims, a_dimnames )
        call ni3_cp_atts( inncid, outncid, varname )
        call ni3_put_att( outncid, varname, 'missing_value', missing_value )

        deallocate( a_dimnames )


!!$        loop_time : do t = 1, tmax
        loop_time : do t = ts, te

          call ni3_get_varss( inncid   , varname   , t, xyz_dVar )
          call ni3_get_varss( inncid_ps, varname_ps, t, xy_Ps    )

          do k = 1, kmax_in
            xya_Press(:,:,k) = xy_Ps * a_Sigma(k)
          end do

          call logintp2p( &
            & imax, jmax, kmax_in, xya_Press, xyz_dVar,   &
            & kmax_out, plevel(1:kmax_out), xyz_dVar_out, &
            & missing_value &
            & )

          call ni3_put_varss( outncid, varname, t-ts+1, xyz_dVar_out )

        end do loop_time

        deallocate( a_Sigma      )
        deallocate( xyz_dVar     )
        deallocate( xy_Ps        )
        deallocate( xya_Press    )
        deallocate( xyz_dVar_out )


      else
        if ( ( .not. flag_dim ) .and. ( .not. flag_weight ) .and. ( .not. flag_rst ) ) then
!!$          write( 6, * ) '      Skip'
          write( 6, * ) trim( varname ), ' is skipped.'
        end if
      end if


    end do loop_variables



    if ( FlagIncPlanetLon ) then
      ncfn = inplncfn
      mode = "read"
      call ni3_open( ncfn, mode, inplncid )

      xtype = NI3_REAL
      ndims = 1
      allocate( a_dimnames(ndims) )
      a_dimnames(1) = 'time'
      varname = plvarname
      call ni3_def_var( outncid, varname, xtype, ndims, a_dimnames )
      deallocate( a_dimnames )

      call ni3_cp_atts( inplncid, outncid, varname )

      loop_time_pl : do t = ts, te
        call ni3_get_varss( inplncid, varname, t     , PlanetLon )
        call ni3_put_varss( outncid , varname, t-ts+1, PlanetLon )
      end do loop_time_pl

      call ni3_close( inplncid )
    end if



    deallocate( z_Sigma )
    deallocate( r_Sigma )

    deallocate( a_dimnames_all )

    deallocate( a_xtype_all  )
    deallocate( a_dimlen_all )


    call ni3_close( outncid   )
    call ni3_close( inncid    )
    call ni3_close( inncid_ps )

  end do loop_file

  close( ctlfu )


end program spread

