  program main

    use vtype_module
    use ni3_module
    use netcdf

    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)                  :: lon1( 1 ), time1( 1 )
    real(sp)    , allocatable :: &
         buf2d( :, : )   , buf2d_sum( :, : ), &
         buf3d( :, :, : ), buf3d_sum( :, :, : )
    real(sp)    , allocatable :: &
         buf3d2( :, :, : )
    real(sp)                  :: zm


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

    integer(i4b)                       :: ndims
    character(len=extstr), allocatable :: dimname( : )


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


    namelist /input/     ncfn_in, item, ncfn_in2
    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( &
            buf2d( im, jm )    , buf2d_sum( 1, jm )    , &
            buf3d( im, jm, km ), buf3d_sum( 1, jm, km ), &
            buf3d2( im, jm, km ) &
            )


       if( item .eq. 'T2' ) then
          name = 'Temp'
       else if( item .eq. 'VT' ) then
          name  = 'V'
          name2 = 'Temp'

          mode = 'read'
          call ni3_open( ncfn_in2, mode, ncid_in2 )

       else
          name = item
       end if


       call ni3_inq_ndims( ncid_in, name, ndims )

       allocate( dimname( ndims ) )
       if( ndims .eq. 3 ) then
          dimname( 1 ) = 'lon'
          dimname( 2 ) = 'lat'
          dimname( 3 ) = 'time'
       else
          dimname( 1 ) = 'lon'
          dimname( 2 ) = 'lat'
          dimname( 3 ) = 'sig'
          dimname( 4 ) = 'time'
       end if

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

       if( ndims .eq. 3 ) then

          buf2d_sum( :, : )    = 0.0
          do t = ts, te
             call ni3_get_varss( ncid_in, name, t, buf2d )

             buf2d_sum( :, : ) = buf2d_sum( :, : ) + buf2d( :, : )
          end do
          buf2d_sum( :, : ) = buf2d_sum( :, : ) / real( im * ( te - ts + 1 ) )

       else

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

             if( item .eq. 'T2' ) then
                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_sum( 1, j, k ) = buf3d_sum( 1, j, k ) &
                              + ( buf3d( i, j, k ) - zm )**2
                      end do
                   end do
                end do
             else if( item .eq. 'VT' ) then

               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

               call ni3_get_varss( ncid_in2, name2, t, buf3d2 )
               do k = 1, km
                 do j = 1, jm
                   zm = 0.0
                   do i = 1, im
                     zm = zm + buf3d2( i, j, k )
                   end do
                   zm = zm / real( im )
                   do i = 1, im
                     buf3d2( i, j, k ) = buf3d2( i, j, k ) - zm
                   end do
                 end do
               end do

               do k = 1, km
                 do j = 1, jm
                   do i = 1, im
                     buf3d_sum( 1, j, k ) = buf3d_sum( 1, j, k ) &
                       + buf3d( i, j, k ) * buf3d2( i, j, k )
                   end do
                 end do
               end do

             else
!!$                buf3d_sum( :, :, : ) = buf3d_sum( :, :, : ) + 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_sum( 1, j, k ) = buf3d_sum( 1, j, k ) + zm
                      end do
                   end do
                end do
             end if

          end do
          buf3d_sum( :, :, : ) = buf3d_sum( :, :, : ) / real( im * ( te - ts + 1 ) )

       end if



       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, 'sig'  , 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, 'sig'   )
          call ni3_cp_atts( ncid_in, ncid_out, 'time'  )
       end if


       if( item .eq. 'PS' ) then
          varname = 'ps'
          stdname = 'surface_air_pressure'
          unit    = 'Pa'
       else if( item .eq. 'U' ) then
          varname = 'ua'
          stdname = 'eastward_wind'
          unit    = 'm s-1'
       else if( item .eq. 'V' ) then
          varname = 'va'
          stdname = 'northward_wind'
          unit    = 'm s-1'
       else if( item .eq. 'Temp' ) then
          varname = 'ta'
          stdname = 'air_temperature'
          unit    = 'K'
       else if( item .eq. 'VT' ) then
          varname = 'vt'
          stdname = 'northward_wind_air_temperature_eddy_covariance'
          unit    = 'm s-1 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, units = unit )
!!$    call ni3_cp_atts( ncid_in, ncid_out, name )
       t = 1
       if( ndims .eq. 3 ) then
          call ni3_put_varss( ncid_out, varname, t, buf2d_sum )
       else
          call ni3_put_varss( ncid_out, varname, t, buf3d_sum )
       end if


       deallocate( dimname )
       deallocate( lon, lat, sigma )
       deallocate( buf2d, buf2d_sum, buf3d, buf3d_sum )
       deallocate( buf3d2 )

       call ni3_close( ncid_in )
       if ( item == 'VT' ) then
         call ni3_close( ncid_in2 )
       end if

    end do loop


    call ni3_close( ncid_out )


  end program main
