program main

  use vtype_module
  use ni3_module
  use netcdf

  implicit none

  integer, parameter :: NumCases = 21

  integer               :: a_CaseNum     (NumCases)
  character(len=extstr) :: a_CaseChar    (NumCases)
  character(len=extstr) :: a_CaseTailChar(NumCases)

  integer               :: CaseNum
  character(len=extstr) :: CaseChar
  character(len=extstr) :: CaseTailChar

  integer         , parameter :: kmax = 34-2
  real(sp)                    :: z_PressTop(1:kmax)
  real(sp)                    :: z_PressBot(1:kmax)
  real(sp)                    :: z_Press   (1:kmax)
  real(sp)                    :: z_DTempDt (1:kmax)


  character(len=extstr)            :: line

  character(len=extstr)            :: ncfn
  character(len=extstr)            :: mode
  character(len=4*extstr)          :: title, com, ref, inst, conv, src, hist
  character(len=extstr)            :: name, stdname, longname, units

  integer(i4b)                     :: ncid

  integer(i4b)                       :: NDims
  character(len=extstr), allocatable :: a_Dimname( : )

  real(sp)             , parameter :: FillValue = -99999.0

  character(len=extstr)            :: dir0
  character(len=extstr)            :: dir
  character(len=extstr)            :: fn
  character(len=extstr)            :: path
  integer                 :: fu
  integer                 :: ios
  logical                 :: ex
  logical                 :: od

  integer(i4b)          :: k
  integer(i4b)          :: l
  integer(i4b)          :: n

  data a_CaseNum      / 9,     10,    19,    20,    23,    23,    23,    25,    27,    28,    29,    31,    33,    35,    37,    37,    37,    52,    53,    54,    55 /
  data a_CaseChar     / '',    '',    '',    '',   'a',   'b',   'c',    '',    '',    '',    '',    '',    '',    '',   'a',   'b',   'g',    '',    '',    '',    '' /
  data a_CaseTailChar / 'mls', 'mls', 'mls', 'mls', 'mls', 'mls', 'mls', 'trp', 'mls', 'mls', 'mlw', 'sas', 'saw', 'mls', 'mls', 'mls', 'mls', 'mls', 'mls', 'mls', 'mls' /


  do n = 1, NumCases
    write( 6, * ) n

    CaseNum      = a_CaseNum(n)
    CaseChar     = a_CaseChar(n)
    CaseTailChar = a_CaseTailChar(n)

    !
    ! Output to a NetCDF file
    !
    mode = "new"
    title = "Radiative heating rate calculated by LBLRTM for ICRCCM"
    inst  = "GFD Dennou Club, DCPAM project"
    src   = "http://www.met.fsu.edu/people/ellingson/LBLWeb/"
    conv  = "CF-1.5"
    hist  = "An electronic (ASCII) data file was downloaded from http://www.met.fsu.edu/people/ellingson/LBLWeb/, and was converted into netcdf file by GFD Dennou Club."
    ref   = "Ellingson, R.G., J. Ellis, and S. Fels, 1991: "
    ref   = trim(ref) // "The intercomparison of radiation codes used in climate models: Long wave results, J. Geophys. Res., 96, D5, 8929."
    com   = ""!Please note that all data are not included."
    !
    write( ncfn, '(a,i2.2,a)' ) "out/Case", CaseNum, trim( CaseChar ) // "hr.nc"
    write( 6, * ) 'Creating ', trim( ncfn )
    call ni3_open( ncfn, mode, ncid )
    call ni3_set_ga( ncid, title = title, conv = conv, inst = inst, src = src, hist = hist, ref = ref, com = com )
    !
    !============================================================================
    name    = "level"
    stdname = "level"
    units   = "Pa"
!!$  call ni3_set_dim( ncid, name, NF90_REAL, z_LevelNum, &
!!$    & stdname = stdname, units  = units )
    call ni3_def_dim( ncid, name, NF90_REAL, kmax, &
      & stdname = stdname, units  = units )
    call ni3_put_att( ncid, name, 'positive', 'down' )
    !
    NDims = 1
    allocate( a_Dimname( NDims ) )
    a_DimName( 1 ) = "level"
    !
    name         = "plevbot"
    stdname      = "pressure of layer bottom"
    units        = "Pa"
    call ni3_def_var( ncid, name, NF90_REAL, ndims, a_DimName, &
      & stdname = stdname, units = units, fillvalue = FillValue )
    !
    name         = "plevtop"
    stdname      = "pressure of layer top"
    units        = "Pa"
    call ni3_def_var( ncid, name, NF90_REAL, ndims, a_DimName, &
      & stdname = stdname, units = units, fillvalue = FillValue )
    !
    name         = "hr"
    stdname      = "radiative_heating_rate"
    units        = "K s-1"
    call ni3_def_var( ncid, name, NF90_REAL, ndims, a_DimName, &
      & stdname = stdname, units = units, fillvalue = FillValue )
    !
    deallocate( a_Dimname )
    !============================================================================


    dir0 = '../data/www.met.fsu.edu/people/ellingson/LBLWeb/results'
    write( dir, '(a,i2.2,a)' ) trim( dir0 ) // '/' // 'Case', CaseNum, trim( CaseChar ) // '-' // trim( CaseTailChar )
    write( fn, '(a,i2.2,a)' ) 'Case', CaseNum, trim( CaseChar ) // 'coolrt_' // trim( CaseTailChar ) // '.txt'
    path = trim( dir ) // '/' // trim( fn )

    write( 6, * ) 'Reading ', trim( path )

    inquire( file = path, opened = od, exist = ex )

    fu = 50
    open( fu, file=path, status='unknown' )

    ! header
    do l = 1, 2
      read( fu, * ) line
    end do

    do k = 1, kmax
!    mb          mb        deg K/day
      read( fu, *, iostat = ios ) z_PressBot(k), z_PressTop(k), z_DTempDt(k)
      if ( ios == -1 ) stop 'Error during reading a file'
    end do

    z_PressBot = z_PressBot * 1.0d2
    z_PressTop = z_PressTop * 1.0d2
    z_DTempDt  = - z_DTempDt  / 86400.0d0

    z_Press = ( z_PressBot + z_PressTop ) / 2.0d0

    name = "level"
    call ni3_put_var( ncid, name, z_Press )
    name = "plevbot"
    call ni3_put_var( ncid, name, z_PressBot )
    name = "plevtop"
    call ni3_put_var( ncid, name, z_PressTop )
    name = "hr"
    call ni3_put_var( ncid, name, z_DTempDt )

    close( fu )

    call ni3_close( ncid )

  end do


end program main
