  program main

    use vtype
    use gi_module
    use netcdf
    use ni3_module

    implicit none


    integer(i4b)                          :: im, jm
    integer(i4b)                          :: kup, km
    real(dp)                              :: pint, pt
    real(dp)                , allocatable :: sigma( : )
    real(dp)                , allocatable :: gps( :, : ), gp( :, :, : )
    real(dp)                , allocatable :: vlev8( : )
    real(dp)                , allocatable :: &
         buf_in3d8( :, :, : ), buf_out3d8( :, :, : )

    integer(i4b)                          :: item_index, ps_index

    real(sp)                , allocatable :: lon( : ), lat( : ), sigmam( : )
    integer(i4b)                          :: vlevn
    real(sp)                , allocatable :: vlev( : ), vlev_inv( : )
    real(sp)                              :: time

    real(sp)                , allocatable :: buf_2d( :, : ), &
         buf_in3d ( :, :, : ), buf_out3d( :, :, : ), buf_out3d_inv( :, :, : )


    character( len=extstr )               :: fn

    integer( i4b )                        :: i, k, l, t
    integer( i4b )                        :: ts, te

    logical                               :: ex
    integer(i4b)                          :: ios

    type( gi_file_att_type )              :: gfa_ps, gfa_in

    character(len=extstr)                 :: ctlfn = "gcm_gr2nc.cntl"
    integer(i4b)                          :: ctlfu

    integer(i4b)                          :: ncid
    character(len=extstr)                 :: mode, title, com, &
         name, stdname, longname, units, &
         conv, title, hist, inst, src
    integer(i4b)            , parameter   :: ndims2dvar = 3, ndims3dvar = 4
    character(len=extstr)                 :: &
         dimname2dvar( ndims2dvar ), dimname3dvar( ndims3dvar )
    real(sp)                              :: fillvalue

    character(len=extstr)                 :: dir_ps , gfn_ps, item_ps
    character(len=extstr)                 :: dir_in , gfn_in, item
    character(len=extstr)                 :: ncfn_out
    logical                               :: sw_pout
    character(len=extstr)                 :: vdim_type
    integer(i4b)         , parameter      :: plevmax    = 50
    integer(i4b)         , parameter      :: d_plevn    = 35
    integer(i4b)         , parameter      :: d_plevempn = plevmax - d_plevn
    integer(i4b)                          :: plevn
    real(dp)                              :: plevel( plevmax )

    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 /


    namelist /p_out/     sw_pout, dir_ps, gfn_ps, item_ps, vdim_type, &
         plevn, plevel
    namelist /input/     dir_in , gfn_in, item
    namelist /output/    ncfn_out
    namelist /timeinfo/  ts, te

!    namelist /varinfo/   var_varname, var_stdname, var_longname, var_units
    namelist /varinfo/   name, stdname, longname, units

    fn = ctlfn
    inquire( file = fn, exist = ex )
    if( .not. ex ) then
       write( 6, * ) "Control file ", trim( ctlfn ), " does not exit."
       write( 6, * ) "Control file should contain namelists listed below."
       write( 6, * ) "    &p_out"
       write( 6, * ) "        sw_pout   = .true."
       write( 6, * ) "        dir_ps    = 'psdir'"
       write( 6, * ) "        gfn_ps    = 'ps'"
       write( 6, * ) "        item_ps   = 'ps'"
       write( 6, * ) "        vdim_type = 'mymodel'"
       write( 6, * ) "        plevn     = 4"
       write( 6, * ) "        plevel    = 0.18d0   0.5d0"
       write( 6, * ) "                    1.36d0   3.7d0"
       write( 6, * ) "    &end"
       write( 6, * ) "    &input"
       write( 6, * ) "        dir_in    = 'indir'"
       write( 6, * ) "        gfn_in    = 'ingrd'"
       write( 6, * ) "        item      = 't'"
       write( 6, * ) "    &end"
       write( 6, * ) "    &output"
       write( 6, * ) "        ncfn_out  = 'out.nc'"
       write( 6, * ) "    &end"
       write( 6, * ) "plev is optional."
       stop
    end if
    call findfu( fn, ios, ctlfu )
    if( ios /= 0 ) then
       write( 6, * ) 'STOP: ', ios
       stop
    end if
    open( ctlfu, file = fn, status = 'unknown' )


    !
    ! check size of arrays
    !
    rewind( ctlfu )
    read( ctlfu, input )
    write( 6, * ) "Checking size of array using GrADs file ", &
         trim( dir_in ) // "/" // trim( gfn_in ), ".ctl ..."
    write( fn, '(3a)' ) trim( dir_in ), "/", trim( gfn_in )
    call gi_parse_ctl( fn, gfa_in )
    write( 6, * ) "xdef : ", gfa_in % xdef
    write( 6, * ) "ydef : ", gfa_in % ydef
    write( 6, * ) "zdef : ", gfa_in % zdef
    write( 6, * ) "tdef : ", gfa_in % tdef
!    if( gfa_in % nvars .ne. 1 ) &
!         stop "STOP: Number of variables included in GrADs file is not 1. "

    im = gfa_in % xdef
    jm = gfa_in % ydef
    km = gfa_in % zdef

    allocate( lon( im ) )
    allocate( lat( jm ) )
    do i = 1, im
       lon( i ) = gfa_in % x1 + gfa_in % x2 * real( i - 1 )
    end do
    lat( : ) = gfa_in % ylevs( : )
    allocate( sigmam( km ) )
    do k = 1, km
       sigmam( k ) = gfa_in % zlevs( km - k + 1 )
    end do

    call gi_close_data( gfa_in )


    !
    ! set vertical levels
    !
    sw_pout   = .false.
    item_ps   = "ps"
    vdim_type = "mymodel"
    rewind( ctlfu )
    read( ctlfu, p_out, iostat = ios )
    if( sw_pout ) then

       if( ( vdim_type .eq. "mymodel" ) &
            .or. ( vdim_type .eq. "MYMODEL" ) &
            .or. ( vdim_type .eq. "Mymodel" ) ) then
          allocate( sigma( km+1 ) )
          k = 1
          sigma( k ) = -1.0d0

!          sigmam( k ) = ( sigma( k ) + sigma( k+1 ) ) / 2.0d0

          do k = 1+1, km+1-1
             sigma( k ) = 2.0d0 * sigmam( k-1 ) - sigma( k-1 )
          end do
          k = km+1
          sigma( k ) = 1.0d0
       end if


       vlevn = plevn
       allocate( vlev( vlevn ) )
       do k = 1, vlevn
          vlev( k ) = plevel( k )
       end do
    else
       vlevn = km
       allocate( vlev( vlevn ) )
       vlev( : ) = sigmam( : )
    end if
    do k = 1, vlevn
       write( 6, * ) k, vlev( k )
    end do
    allocate( vlev8( vlevn ) )
    vlev8( : ) = vlev( : )


    !
    ! output file
    !
    rewind( ctlfu )
    read( ctlfu, output )
    write( 6, * ) "Output NetCDF file  : ", trim( ncfn_out )

    mode = "new"
    call ni3_open( ncfn_out, mode, ncid )
    !
    conv  = "CF-1.0"
    title = "Aqua Planet: Transients from 'control' Experiment" 
    hist  = "Original data produced: 2005/07/01 on ES; translated form at to conform with APE TR standard output" 
    inst = "AGU for APE; AFES Working Team, GFD Dennou Club, University of Tokyo" 
    src  = "AFES v1.15" 
    call ni3_set_ga( ncid, conv, title, hist, inst, src )
    !
    name    = "lon"
    stdname = "longitude"
    units   = "degrees_east"
    call ni3_set_dim( ncid, name, NF90_REAL, lon, &
         stdname = stdname, units = units )

!    call ni3_put_att( ncid, "lon", "topology", "circular" )
!    call ni3_put_att( ncid, "lon", "modulo", 360.0 )

    !
    name    = "lat"
    stdname = "latitude"
    units   = "degrees_north"
    call ni3_set_dim( ncid, name, NF90_REAL, lat, &
         stdname = stdname, units = units )
    !
    allocate( vlev_inv( vlevn ) )
    do k = 1, vlevn
       vlev_inv( k ) = vlev( vlevn - k + 1 )
    end do
    if( sw_pout ) then
       name    = "plev"
       stdname = "air_pressure"
       units   = "Pa"
       call ni3_set_dim( ncid, name, NF90_REAL, vlev_inv, &
            stdname = stdname, units = units )
    else
       name    = "sigma"
       stdname = "sigma"
       units   = "sigma_level"
       call ni3_set_dim( ncid, name, NF90_REAL, vlev_inv, &
            stdname = stdname, units = units )
    end if
    !
    name    = "time"
    stdname = "time"
    units   = "days since 0000-01-01"
    call ni3_def_dim( ncid, name, NF90_REAL, NF90_UNLIMITED, &
         stdname = stdname, units = units )


    dimname2dvar( 1 ) = "lon"
    dimname2dvar( 2 ) = "lat"
    dimname2dvar( 3 ) = "time"
    dimname3dvar( 1 ) = "lon"
    dimname3dvar( 2 ) = "lat"
    if( sw_pout ) then
       dimname3dvar( 3 ) = "plev"
    else
       dimname3dvar( 3 ) = "sigma"
    end if
    dimname3dvar( 4 ) = "time"


    allocate( gps          ( im, jm ) )
    allocate( gp           ( im, jm, km    ) ) 
    allocate( buf_2d       ( im, jm ) ) 
    allocate( buf_in3d     ( im, jm, km    ) ) 
    allocate( buf_out3d    ( im, jm, vlevn ) ) 
    allocate( buf_out3d_inv( im, jm, vlevn ) ) 

    allocate( buf_in3d8    ( im, jm, km    ) ) 
    allocate( buf_out3d8   ( im, jm, vlevn ) ) 


    !
    ! read time information
    !
    ts = 1
    te = gfa_in % tdef
    rewind( ctlfu )
    read( ctlfu, timeinfo, iostat = ios )

    write( 6, * ) 'Converting time index from ', ts, ' to ', te, '.'
    if( te .gt. gfa_in % tdef ) then
       te = min( te, gfa_in % tdef )
       write( 6, * ) 'End time index is modified to ', te, '.'
    end if


    rewind( ctlfu )

    loop : do
       item = "###"
       read( ctlfu, input, iostat = ios )
       if( ios .ne. 0 ) exit loop
       if( item .eq. "###" ) item = gfn_in

       write( fn, '(3a)' ) trim( dir_in ), "/", trim( gfn_in )
       write( 6, * ) "Reading GrADs control file : ", trim( fn ), ".ctl"
       call gi_parse_ctl( fn, gfa_in )


       search_item : do item_index = 1, gfa_in % nvars
          if( gfa_in % vars( item_index ) % var_name .eq. item ) &
               exit search_item
       end do search_item
       if( item_index .gt. gfa_in % nvars ) then
          write( 6, * ) "Unable to find ", trim( item ), &
               " in ", trim( fn ), "."
          stop
       end if


       l = item_index

!       name      = gfa_in % vars( l ) % var_name
!       stdname   = gfa_in % vars( l ) % var_name
!       longname  = gfa_in % vars( l ) % var_title
!       units     = gfa_in % vars( l ) % var_unit

       read( ctlfu, varinfo )
       fillvalue = gfa_in % undef

       write( 6, * ) "Variable : ", trim( name )


       !
       ! Check array size
       !
       if( ( im .ne. gfa_in % xdef ) .or. ( jm .ne. gfa_in % ydef ) &
            .or. &
            ( ( km .ne. gfa_in % vars( item_index ) % nz ) &
            .and. ( 1 .gt. gfa_in % vars( item_index ) % nz ) ) ) then
          write( 6, * ) "Array size is not appropriate."
          write( 6, * ) im, jm, km
          write( 6, * ) gfa_in % xdef, gfa_in % ydef, &
               gfa_in % vars( item_index ) % nz
          stop
       end if

       write( 6, * ) "Opening GrADs data file : ", trim( fn ), ".grd"
       call gi_open_data( dir_in , gfa_in , "read"  )
       write( 6, * ) "File unit of input  GrADs data file : ", gfa_in % un



       if( gfa_in % vars( l ) % nz .le. 1 ) then

          call ni3_def_var( ncid, name, NF90_REAL, &
               ndims2dvar, dimname2dvar, &
               stdname, longname, units, fillvalue )

!          do t = 1, gfa_in % tdef
          do t = ts, te
             time = real( t )

             call gi_read_data( gfa_in, l, t, 1, buf_2d )
!!$             call ni3_put_varss( ncid, "time", t, time )
!!$             call ni3_put_varss( ncid, name  , t, buf_2d )
             call ni3_put_varss( ncid, "time", t-ts+1, time )
             call ni3_put_varss( ncid, name  , t-ts+1, buf_2d )
          end do

       else

          if( sw_pout ) then
             write( fn, '(3a)' ) trim( dir_ps ), "/", trim( gfn_ps )
             call gi_parse_ctl( fn, gfa_ps )

             call gi_open_data( dir_ps, gfa_ps, "read" )

             !
             ! Check array size
             !
             if( ( im .ne. gfa_ps % xdef ) &
                  .or. ( jm .ne. gfa_ps % ydef ) ) then
                write( 6, * ) &
                     "Number of vertical levels is not appropriate."
                write( 6, * ) im, jm
                write( 6, * ) gfa_ps % xdef, gfa_ps % ydef
                stop
             end if

             write( 6, * ) "File unit of ps     GrADs data file : ", &
                  gfa_ps  % un

             search_ps : do ps_index = 1, gfa_ps % nvars
                if( gfa_ps % vars( ps_index ) % var_name .eq. item_ps ) &
                     exit search_ps
             end do search_ps
             if( ps_index .gt. gfa_ps % nvars ) then
                write( 6, * ) "Unable to find ", trim( item_ps ), &
                     " in ", trim( fn ), "."
                stop
             end if

          end if


          call ni3_def_var( ncid, name, NF90_REAL, &
               ndims3dvar, dimname3dvar, &
               stdname, longname, units, fillvalue )

!          do t = 1, gfa_in % tdef
          do t = ts, te
             time = real( t )

             call gi_read_3ddata( gfa_in, buf_in3d, l, t, .false. )

             !
             ! convert vertical grid to pressure
             !
             if( sw_pout ) then
                call gi_read_data( gfa_ps, ps_index, t, 1, gps )
!                gps( :, : ) = gps( :, : ) * 1.0d2

                if( ( vdim_type .eq. "mymodel" ) &
                     .or. ( vdim_type .eq. "MYMODEL" ) &
                     .or. ( vdim_type .eq. "Mymodel" ) ) then
                   call calcgp( im, jm, kup, km, pint, pt, sigma, gps, gp )
                else if( ( vdim_type .eq. "traditional" ) &
                     .or. ( vdim_type .eq. "TRADITIONAL" ) &
                     .or. ( vdim_type .eq. "Traditional" ) ) then
                   do k = 1, km
                      gp( :, :, k ) = gps( :, : ) * sigmam( k )
                   end do
                else
                   write( 6, * ) 'Undefined vdim_type, ', &
                        trim( vdim_type ), '.'
                   stop
                end if
!                gp( :, :, : ) = gp( :, :, : ) * 1.0d-2
                buf_in3d8( :, :, : ) = buf_in3d  ( :, :, : )
                call logintp2p( im, jm, km, gp, buf_in3d8, &
                     vlevn, vlev8, buf_out3d8, dble( fillvalue ) )
                buf_out3d( :, :, : ) = buf_out3d8( :, :, : )
             else
                buf_out3d( :, :, : ) = buf_in3d  ( :, :, : )
             end if

!!$             call ni3_put_varss( ncid, "time", t, time )
             call ni3_put_varss( ncid, "time", t-ts+1, (time-ts+1)/4+1160 )
             do k = 1, vlevn
                buf_out3d_inv( :, :, k ) = buf_out3d( :, :, vlevn - k + 1 )
             end do
!!$             call ni3_put_varss( ncid, name  , t, buf_out3d_inv )
             call ni3_put_varss( ncid, name  , t-ts+1, buf_out3d_inv )
          end do

          if( sw_pout ) call gi_close_data( gfa_ps )

       end if

       call gi_close_data( gfa_in )

    end do loop

    close( ctlfu )

    deallocate( buf_2d    )
    deallocate( buf_in3d  )
    deallocate( buf_out3d )

    call ni3_close( ncid )


  end program main
