  program main

    use vtype_module
    use ni3_module
    use netcdf
    use fft6_module

    implicit none

    interface
       subroutine findfu( fn, ios, fu, mode )
         use vtype_module
         implicit none
         character(len=*), intent(in )           :: fn
         integer(i4b)    , intent(out)           :: ios, fu
         character(len=*), intent(in ), optional :: mode
       end subroutine findfu
    end interface


    character(len=extstr), parameter :: ctlfn = 'mknc.cntl'
    integer(i4b)                     :: ctlfu, ios
    integer(i4b)                     :: loop_cnt

    integer(i4b)              :: im, jm, km, tm
    real(sp)    , allocatable :: lon( : ), lat( : ), sigma( : )
    real(sp)                  :: time1( 1 )
    real(sp)    , allocatable :: buf3d( :, :, : )
    real(sp)                  :: zm


    character(len=extstr)     :: ncfn_in, ncfn_out
    integer(i4b)              :: ncid_in, ncid_out
    character(len=extstr)     :: item, mode, name, varname, stdname, unit

    integer(i4b)         , parameter :: ndims = 3
!!$    integer(i4b)         , parameter :: ndims = 4
    character(len=extstr)            :: dimname( ndims )


    integer(i4b)              :: mm, nn
    real(dp)    , allocatable :: data ( :, : )
    real(dp)    , allocatable :: datar( :, : ), datai( :, : )
    real(sp)    , allocatable :: power_ave( :, : )
    real(sp)    , allocatable :: zwn( : )
    real(sp)    , allocatable :: hsgm( : ), dsgm( : )


    integer(i4b)              :: ita( 5 )
    real(dp)    , allocatable :: dta( :, : )
    integer(i4b)              :: isign


    integer(i4b)              :: i, j, k, m, n, t
    integer(i4b)              :: ts , te
    integer(i4b)              :: ts0, te0


    namelist /input/      ncfn_in, item
    namelist /output/     ncfn_out
    namelist /time_info/  ts, te


    call findfu( ctlfn, ios, ctlfu )
    if( ios /= 0 ) then
       write( 6, * ) 'STOP: ', ios
       stop
    end if
    open( ctlfu, file = ctlfn, status = 'unknown' )


    ncfn_in  = '../ncfiles/T21L48_hs94m_ls270_tau0.2_zm.nc'
    ncfn_out = '../ncfiles_submit/T21L48_hs94m_ls270_tau0.2_zm.nc'
    name = 'U'


    rewind( ctlfu )
    read( ctlfu, output    , iostat = ios )
    rewind( ctlfu )
    read( ctlfu, time_info , iostat = ios )
    ts0 = ts; te0 = te


    loop_cnt = 0
    rewind( ctlfu )
    loop : do
       loop_cnt = loop_cnt + 1

       item = '###'
       read( ctlfu, input, iostat = ios )
       if( ios .ne. 0 ) exit loop

       write( 6, * ) 'input file: ', trim( ncfn_in )
       write( 6, * ) 'item      : ', trim( item    )

       mode = 'read'
       call ni3_open( ncfn_in , mode, ncid_in  )

       call ni3_inq_dimlen( ncid_in, 'lon'  , im )
       call ni3_inq_dimlen( ncid_in, 'lat'  , jm )
       call ni3_inq_dimlen( ncid_in, 'sig'  , km )
       call ni3_inq_dimlen( ncid_in, 'time' , tm )
       allocate( lon( im ), lat( jm ), sigma( km ) )
       call ni3_get_var( ncid_in, 'lon'  , lon   )
       call ni3_get_var( ncid_in, 'lat'  , lat   )
       call ni3_get_var( ncid_in, 'sig'  , sigma )

       allocate( buf3d( im, jm, km ) )


       mm = jm * km
       nn = im

       allocate( dta( nn, 2 ) )
       allocate( data     ( mm, nn ) )
       allocate( datar    ( mm, nn ), datai    ( mm, nn ) )
       allocate( power_ave( nn, jm ) )
       allocate( zwn( nn ) )

       call fft6_fft_init( nn, ita, dta )


       do n = 1, nn
          zwn( n ) = fft6_i2m( nn, n )
       end do


       allocate( hsgm( km+1 ), dsgm( km ) )
       call ni3_get_var( ncid_in, 'sigm', hsgm )
       do k = 1, km
          dsgm( k ) = hsgm( k ) - hsgm( k+1 )
       end do


       ts = max( 1 , ts0 )
       te = min( tm, te0 )

       power_ave( :, : ) = 0.0
       do t = ts, te
          write( 6, * ) ts, t, te
          call ni3_get_varss( ncid_in, name, t, buf3d )

          do k = 1, km
             do j = 1, jm
                zm = 0.0
                do i = 1, im
                   zm = zm + buf3d( i, j, k )
                end do
                zm = zm / real( im )
                do i = 1, im
                   buf3d( i, j, k ) = buf3d( i, j, k ) - zm
                end do
             end do
          end do

!!$          do k = 1, km
!!$             do j = 1, jm
!!$                do i = 1, im
!!$                   n = i
!!$                   m = ( k - 1 ) * jm + j
!!$                   datar( m, n ) = buf3d( i, j, k )
!!$                   datai( m, n ) = 0.0d0
!!$                end do
!!$             end do
!!$          end do
!!$
!!$          isign = 1
!!$          call fft6_fft( mm, nn, ita, dta, datar, datai, isign )

          do k = 1, km
             do j = 1, jm
                do i = 1, im
                   n = i
                   m = ( k - 1 ) * jm + j
                   data( m, n ) = buf3d( i, j, k )
                end do
             end do
          end do
          isign = 1
          call fft6_fftreal( mm, nn, ita, dta, data, datar, datai, isign )

          do k = 1, km
             do j = 1, jm
                do i = 1, im
                   n = i
                   m = ( k - 1 ) * jm + j
                   power_ave( n, j ) = power_ave( n, j ) &
                        + ( datar( m, n )**2 + datai( m, n )**2 ) * dsgm( k )
!!$                   power_ave( n, j ) = power_ave( n, j ) &
!!$                        + ( datar( m, n )**2 + datai( m, n )**2 )
                end do
             end do
          end do

       end do

       power_ave( :, : ) = power_ave( :, : ) / real( te - ts + 1 )


!!$       power_ave( :, : ) = power_ave( :, : ) * 3389.0e3
!!$       power_ave( :, : ) = power_ave( :, : ) * 6371.0e3
       power_ave( :, : ) = power_ave( :, : )

       power_ave( :, : ) = power_ave( :, : ) * 2.0




       if( loop_cnt .eq. 1 ) then
          mode = 'new'
          call ni3_open( ncfn_out, mode, ncid_out )

!!$          lon1( 1 )  = 0.0
!!$          call ni3_set_dim( ncid_out, 'lon'  , NF90_REAL, lon1  )
!!$          call ni3_set_dim( ncid_out, 'lat'  , NF90_REAL, lat   )
!!$          call ni3_set_dim( ncid_out, 'sigma', NF90_REAL, sigma )
!!$          time1( 1 ) = 0.0
!!$          call ni3_set_dim( ncid_out, 'time' , NF90_REAL, time1 )
!!$
!!$          call ni3_cp_atts( ncid_in, ncid_out, 'lon'   )
!!$          call ni3_cp_atts( ncid_in, ncid_out, 'lat'   )
!!$          call ni3_cp_atts( ncid_in, ncid_out, 'sigma' )
!!$          call ni3_cp_atts( ncid_in, ncid_out, 'time'  )

!!$          call ni3_set_dim( ncid_out, 'zwn'  , NF90_REAL, zwn   )
          call ni3_set_dim( ncid_out, 'zwn'  , NF90_REAL, zwn( 1:(nn-1)/3+1 ) )
          call ni3_set_dim( ncid_out, 'lat'  , NF90_REAL, lat   )
          time1( 1 ) = 0.0
          call ni3_set_dim( ncid_out, 'time' , NF90_REAL, time1 )

          call ni3_cp_atts( ncid_in, ncid_out, 'lat'   )
          call ni3_cp_atts( ncid_in, ncid_out, 'time'  )
       end if


!!$       dimname( 1 ) = 'lon'
!!$       dimname( 2 ) = 'lat'
!!$       dimname( 3 ) = 'sigma'
!!$       dimname( 4 ) = 'time'
       dimname( 1 ) = 'zwn'
       dimname( 2 ) = 'lat'
       dimname( 3 ) = 'time'
       if( item .eq. 'U' ) then
          varname = 'uaspe'
          stdname = 'eastward_wind_eddy_variance_spectra'
          unit    = 'm2 s-2'
       else if( item .eq. 'V' ) then
          varname = 'vaspe'
          stdname = 'northward_wind'
          unit    = 'm s-1'
       else if( item .eq. 'T' ) then
          varname = 'ta'
          stdname = 'air_temperature'
          unit    = 'K'
       else if( item .eq. 'T2' ) then
          varname = 'ta2'
          stdname = 'air_temperature_eddy_variance'
          unit    = 'K2'
       end if
       call ni3_def_var( ncid_out, varname, NF90_REAL, ndims, dimname, &
            stdname = stdname, longname = stdname, units = unit )
       t = 1
!!$       call ni3_put_varss( ncid_out, varname, t, buf3d_sum )
!!$       call ni3_put_varss( ncid_out, varname, t, power_ave )
       call ni3_put_varss( ncid_out, varname, t, power_ave( 1:(nn-1)/3+1, 1:jm ) )


!!$       do j = 1, jm
!!$          do n = 1, (nn-1)/3+1
!!$             write( 60, * ) zwn( n ), lat( j ), power_ave( n, j )
!!$          end do
!!$          write( 60, * )
!!$       end do
!!$
!!$
!!$       do j = jm/2+5, jm/2+5
!!$          do n = 1, nn
!!$             write( 61, * ) zwn( n ), lat( j ), power_ave( n, j )
!!$          end do
!!$          write( 61, * )
!!$          write( 61, * )
!!$       end do



       deallocate( lon, lat, sigma )
       deallocate( buf3d )

       deallocate( hsgm, dsgm )

       deallocate( dta )
       deallocate( data )
       deallocate( datar, datai )
       deallocate( power_ave )
       deallocate( zwn )


    end do loop



    call ni3_close( ncid_in  )

    call ni3_close( ncid_out )


  end program main
