;********************************************************************
;
; idl_plot
;
; IDL program to read in line-by-line fluxes and cooling ratesand
; plot selected data interactively.

; Formatted data (fm93) can be converted to IDL binary interactively.

; Data are available for AER, GLA, GFDL, AND KIAE line-by-line models
; from simulations of various ICRCCM cases.  Fluxes and cooling rates
; can be plotted as spectral sums over any selected wavenumber band, or
; as spectral cross-sections at single pressure levels, and multiple
; plots can be overlaid for comparison.
;
; Mike Iacono,   Atmospheric and Environmental Research, Inc.
;      V1.0  Initial release; integrated fluxes and cooling rates.
;            March 25, 1993
;
;      V2.0  Updated to include multiple plots, adjustible X and Y
;            axis scaling, and colorization.
;            April 2, 1993
;
;      V3.0  Updated to include spectral flux and cooling rate
;            cross-section plots, input of parameters from datafile,
;            multiple windows, and an option to perform formatted
;            to binary data conversion.
;            April 21, 1993
;
;            Future changes:
;            1) Plots of flux and cooling rate difference.
;            2) Postscript output for any generated plot.
;            3) Write p-level integrated total for cross-section plots.
;
; Bill Ridgway,  Applied Research Corporation &
;                NASA Goddard Laboratory for Atmospheres
;                Email: ridgway@climate.gsfc.nasa.gov
;                Phone: (301) 286-9138
;
;      V3.1  Incorporated lists for all ICRCCM cases, all groups
;            Modified design of IDL binary file
;            Created directory structure convention for file names
;            May 26, 1993
;
;********************************************************************

; Initialize control parameters and strings

   fnmstr = strarr(10) & titlenm = strarr(10)
   ovrin = '' & sclin = 1
   wn1s = '' & wn2s = ''
   pltflg = 0 & crfac = 1
   winflg = intarr(4) & for i=0,3 do winflg(i) = 0

; Set up output display for multiple plots
;   !p.multi=[0,1,2]
; Send plot to X window and define its size
   set_plot, 'X'
   window, xsize=700 & winflg(0) = 1
;   window, xsize=640, ysize=800
;   window, xsize=512, ysize=640

; Top of program loop
; Select LBL case for input
lbltop:
   print, ' '
   print, ' Select LBL data source:'
   print, ' 1) AER '
   print, ' 2) GLA '
   print, ' 3) GFDL'
   print, ' 4) KIAE'
   print, ' 5) Data conversion (fmt to bin)'
   print, ' 6) Exit program'
   read, gname
   if gname lt 1 or gname gt 6 then goto, lbltop
   if gname eq 6 then goto, pltend
;
; Data conversion; formatted to binary
;
   if gname eq 5 then begin

      filename = ' '
      print, ' '
      print, 'Enter input data filename, or, x to exit:'
      read, filename
      if filename eq 'x' then goto, lbltop
      inname = '../data/'+filename
      outname = '../idl_data/'+filename

; Read in line-by-line data in GLA 'fm93' format and
; write as IDL binary
; reads nlyrs, nlvls, and nbnds from datafile

; Open input file
    openr, 11, inname
; Open output file
    openw, 12, outname

; read initial header, extract key variables
      text = strarr(1)
      for i= 1,13 do readf, 11, text ; 13 lines of header
      readf, 11, text & nsav = strmid(text,30,10) & nlyrs = fix(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & nlvls = fix(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & nbnds = fix(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & psfc = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & tsfc = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & ptrp = float(nsav(0))
      for i=20,29 do readf, 11, text ; 12 lines of unneed variables
      readf, 11, text & nsav = strmid(text,30,10) & vfrst = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & vlast = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & fusfc = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & fdsfc = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & futrp = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & fdtrp = float(nsav(0))
      readf, 11, text & nsav = strmid(text,30,10) & futoa = float(nsav(0))

      data1 = fltarr(6,nlyrs)
      data2 = fltarr(5,nlvls)

      for i= 1, 2 do readf, 11, text ; skip profile data header
      readf, 11, data1
      for i= 1, 2 do readf, 11, text ; skip profile data header
      readf, 11, data2

; print key variables
      print,'...................................'
      print,'..... Data conversion for file : ',inname
      print,'...................................'
      print,'..... nlyrs =',nlyrs
      print,'..... nlvls =',nlvls
      print,'..... nbnds =',nbnds
      print,'..... psfc  =',psfc
      print,'..... tsfc  =',tsfc
      print,'..... ptrp  =',ptrp
      print,'...................................'
      print,'..... vfrst =',vfrst
      print,'..... vlast =',vlast
      print,'..... fusfc =',fusfc
      print,'..... fdsfc =',fdsfc
      print,'..... futrp =',futrp
      print,'..... fdtrp =',fdtrp
      print,'..... futoa =',futoa
      print,'...................................'

; write header for binary IDL file
      writeu, 12, nlyrs, nlvls, nbnds, psfc, tsfc, ptrp
      writeu, 12, vfrst,vlast,fusfc,fdsfc,futrp,fdtrp,futoa
      writeu, 12, data1
      writeu, 12, data2

; read and output spectral band data
      fup = fltarr(nlvls)
      fdn = fltarr(nlvls)
      fnt = fltarr(nlvls)
      htr = fltarr(nlyrs)
      for i=0,nbnds-1 do begin
         readf, 11, vfrst, vlast
         readf, 11, fusfc,fdsfc,futrp,fdtrp,futoa
         readf, 11, format='(5e15.6)', fup
         readf, 11, format='(5e15.6)', fdn
         readf, 11, format='(5e15.6)', fnt
         readf, 11, format='(5e15.6)', htr
         writeu, 12, vfrst,vlast,fusfc,fdsfc,futrp,fdtrp,futoa
         writeu, 12, fup
         writeu, 12, fdn
         writeu, 12, fnt
         writeu, 12, htr
      endfor

      print,'..... Data conversion complete : ',outname
      print,'...................................'
      print,' '
      print,' '
      print,' '

   close,11,12
   goto, lbltop
   endif
;
; End of data conversion
;

; Binary data input

; Set data ID for file name, plot title, and plotting color

; *****AER*****
aertop:
   if gname eq 1 then begin
      fnmstr = ['aer.ph0.c0a', $
                'aer.ph0.c0b', $
                'aer.ph1.c19', $
                'aer.ph1.c20', $
                'aer.ph1.c25a', $
                'aer.ph1.c27a', $
                'aer.ph1.c28a', $
                'aer.ph1.c29a']

      titlenm = [ $
        '(1) AER (0,0a),  TRP, H2O + Cont   ', $
        '(2) AER (0,0b),  SAW, H2O + Cont   ', $
        '(3) AER (1,19),  MLS, H2O + Cont   ',$
        '(4) AER (1,20),  MLS, H2O Lines    ', $
        '(5) AER (1,25a),  TRP, H2O Cont + CO2 (355ppm) + O3   ', $
        '(6) AER (1,27a),  MLS, H2O Cont + CO2 (355ppm) + O3   ', $
        '(7) AER (1,28a),  MLS, H2O Cont + CO2 (710ppm) + O3   ', $
        '(8) AER (1,29a),  MLW, H2O Cont + CO2 (355ppm) + O3   ']

      print, ' '
      print, ' Select AER LBL case:'
      print, ' '
      for i=1,8 do print, '     ',titlenm(i-1)
      print, ' '
      read, fname
      if fname lt 1 or fname gt 8 then goto, aertop
      clrid = 2
   endif

; *****GLA*****
glatop:
   if gname eq 2 then begin
      fnmstr = ['gla.ph0.c00', $
                'gla.ph1.c07', $
                'gla.ph1.c08', $
                'gla.ph1.c09', $
                'gla.ph1.c10', $
                'gla.ph1.c11', $
                'gla.ph1.c12', $
                'gla.ph1.c13', $
                'gla.ph1.c14', $
                'gla.ph1.c15', $
                'gla.ph1.c16', $
                'gla.ph1.c17', $
                'gla.ph1.c18', $
                'gla.ph1.c19', $
                'gla.ph1.c20', $
                'gla.ph1.c21', $
                'gla.ph1.c22', $
                'gla.ph1.c23', $
                'gla.ph1.c24a', $
                'gla.ph1.c24b', $
                'gla.ph1.c27', $
                'gla.ph1.c27x', $
                'gla.ph1.c28', $
                'gla.ph1.c29', $
                'gla.ph1.c30', $
                'gla.ph1.c31', $
                'gla.ph1.c32', $
                'gla.ph1.c33', $
                'gla.ph1.c34', $
                'gla.ph1.c35', $
                'gla.ph1.c36']

      titlenm = [ $
        '(1)  GLA (0, 0),  TRP, H2O + Continuum   ', $
        '(2)  GLA (1, 7),  TRP, CO2 (300ppm)   ', $
        '(3)  GLA (1, 8),  TRP, CO2 (600ppm)   ', $
        '(4)  GLA (1, 9),  MLS, CO2 (300ppm)   ', $
        '(5)  GLA (1,10),  MLS, CO2 (600ppm)   ', $
        '(6)  GLA (1,11),  MLW, CO2 (300ppm)   ', $
        '(7)  GLA (1,12),  MLW, CO2 (600ppm)   ', $
        '(8)  GLA (1,13),  SAS, CO2 (300ppm)   ', $
        '(9)  GLA (1,14),  SAS, CO2 (600ppm)   ', $
        '(10) GLA (1,15),  SAW, CO2 (300ppm)   ', $
        '(11) GLA (1,16),  SAW, CO2 (600ppm)   ', $
        '(12) GLA (1,17),  MLS, H2O + Cont (75% concentration)   ', $
        '(13) GLA (1,18),  MLS, H2O Lines (75% concentration)   ', $
        '(14) GLA (1,19),  MLS, H2O + Cont   ', $
        '(15) GLA (1,20),  MLS, H2O Lines only   ', $
        '(16) GLA (1,21),  MLS, H2O + Cont (125% concentration)   ', $
        '(17) GLA (1,22),  MLS, H2O Lines (125% concentration)   ', $
        '(18) GLA (1,23),  MLS, Ozone   ', $
        '(19) GLA (1,24a), MLS, Ozone (75% concentration)   ', $
        '(20) GLA (1,24b), MLS, Ozone (125% concentration)   ', $
        '(21) GLA (1,27),  MLS, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(22) GLA (1,27x), MLS, H2O + Cont + CO2 + O3 + CH4 + N2O   ', $
        '(23) GLA (1,28),  MLS, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(24) GLA (1,29),  MLW, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(25) GLA (1,30),  MLW, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(26) GLA (1,31),  SAS, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(27) GLA (1,32),  SAS, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(28) GLA (1,33),  SAW, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(29) GLA (1,34),  SAW, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(30) GLA (1,35),  MLS, H2O Lines + CO2 (300ppm) + O3   ', $
        '(31) GLA (1,36),  MLS, H2O Lines + CO2 (600ppm) + O3   ']

      print, ' '
      print, ' Select GLA LBL case:'
      print, ' '
      for i=1,31 do print, '     ',titlenm(i-1)
      print, ' '
      read, fname
      if fname lt 1 or fname gt 31 then goto, glatop
      clrid = 4
   endif

; *****GFDL*****
gfdltop:
   if gname eq 3 then begin
      fnmstr = ['gfdl.ph1.c01', $
                'gfdl.ph1.c02', $
                'gfdl.ph1.c03', $
                'gfdl.ph1.c04', $
                'gfdl.ph1.c05', $
                'gfdl.ph1.c06', $
                'gfdl.ph1.c07', $
                'gfdl.ph1.c08', $
                'gfdl.ph1.c09', $
                'gfdl.ph1.c10', $
                'gfdl.ph1.c11', $
                'gfdl.ph1.c12', $
                'gfdl.ph1.c13', $
                'gfdl.ph1.c14', $
                'gfdl.ph1.c15', $
                'gfdl.ph1.c16', $
                'gfdl.ph1.c17', $
                'gfdl.ph1.c18', $
                'gfdl.ph1.c19', $
                'gfdl.ph1.c20a', $
                'gfdl.ph1.c20b', $
                'gfdl.ph1.c21', $
                'gfdl.ph1.c22', $
                'gfdl.ph1.c23', $
                'gfdl.ph1.c24a', $
                'gfdl.ph1.c24b', $
                'gfdl.ph1.c25', $
                'gfdl.ph1.c26', $
                'gfdl.ph1.c27', $
                'gfdl.ph1.c28', $
                'gfdl.ph1.c29', $
                'gfdl.ph1.c30', $
                'gfdl.ph1.c31', $
                'gfdl.ph1.c32', $
                'gfdl.ph1.c33', $
                'gfdl.ph1.c34', $
                'gfdl.ph1.c35', $
                'gfdl.ph1.c36', $
                'gfdl.ph1.c37', $
                'gfdl.ph2.c01', $
                'gfdl.ph2.c02', $
                'gfdl.ph2.c03', $
                'gfdl.ph2.c04', $
                'gfdl.ph2.c05', $
                'gfdl.ph2.c06', $
                'gfdl.ph2.c07', $
                'gfdl.ph2.c15', $
                'gfdl.ph2.c16', $
                'gfdl.ph2.c17', $
                'gfdl.ph2.c18']

      titlenm = [ $
        '(1)  GFDL (1, 1),  T=200K Isothermal, CO2 (300ppm)   ', $
        '(2)  GFDL (1, 2),  T=200K Isothermal, CO2 (600ppm)   ', $
        '(3)  GFDL (1, 3),  T=250K Isothermal, CO2 (300ppm)   ', $
        '(4)  GFDL (1, 4),  T=250K Isothermal, CO2 (600ppm)   ', $
        '(5)  GFDL (1, 5),  T=300K Isothermal, CO2 (300ppm)   ', $
        '(6)  GFDL (1, 6),  T=300K Isothermal, CO2 (600ppm)   ', $
        '(7)  GFDL (1, 7),  TRP, CO2 (300ppm)   ', $
        '(8)  GFDL (1, 8),  TRP, CO2 (600ppm)   ', $
        '(9)  GFDL (1, 9),  MLS, CO2 (300ppm)   ', $
        '(10) GFDL (1,10),  MLS, CO2 (600ppm)   ', $
        '(11) GFDL (1,11),  MLW, CO2 (300ppm)   ', $
        '(12) GFDL (1,12),  MLW, CO2 (600ppm)   ', $
        '(13) GFDL (1,13),  SAS, CO2 (300ppm)   ', $
        '(14) GFDL (1,14),  SAS, CO2 (600ppm)   ', $
        '(15) GFDL (1,15),  SAW, CO2 (300ppm)   ', $
        '(16) GFDL (1,16),  SAW, CO2 (600ppm)   ', $
        '(17) GFDL (1,17),  MLS, H2O + Cont (75% concentration)   ', $
        '(18) GFDL (1,18),  MLS, H2O Lines (75% concentration)   ', $
        '(19) GFDL (1,19),  MLS, H2O + Cont (100% concentration)  ', $
        '(20) GFDL (1,20a), MLS, H2O Lines (100% concentration)  ', $
        '(21) GFDL (1,20b), MLS, H2O Lines (100% concentration)  ', $
        '(22) GFDL (1,21),  MLS, H2O + Cont (125% concentration)   ', $
        '(23) GFDL (1,22),  MLS, H2O Lines (125% concentration)   ', $
        '(24) GFDL (1,23),  MLS, Ozone   ', $
        '(25) GFDL (1,24a), MLS, Ozone (75% concentration)   ', $
        '(26) GFDL (1,24b), MLS, Ozone (125% concentration)   ', $
        '(27) GFDL (1,25),  TRP, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(28) GFDL (1,26),  TRP, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(29) GFDL (1,27),  MLS, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(30) GFDL (1,28),  MLS, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(31) GFDL (1,29),  MLW, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(32) GFDL (1,30),  MLW, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(33) GFDL (1,31),  SAS, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(34) GFDL (1,32),  SAS, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(35) GFDL (1,33),  SAW, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(36) GFDL (1,34),  SAW, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(37) GFDL (1,35),  MLS, H2O Lines + CO2 (300ppm) + O3   ', $
        '(38) GFDL (1,36),  MLS, H2O Lines + CO2 (600ppm) + O3   ', $
        '(39) GFDL (1,37),  MLS, N2O (0.28ppm)   ', $
        '(40) GFDL (2, 1),  T=200K Isothermal, H2O + Cont   ', $
        '(41) GFDL (2, 2),  T=200K Isothermal, H2O Lines   ', $
        '(42) GFDL (2, 3),  T=250K Isothermal, H2O + Cont   ', $
        '(43) GFDL (2, 4),  T=250K Isothermal, H2O Lines   ', $
        '(44) GFDL (2, 5),  T=300K Isothermal, H2O + Cont   ', $
        '(45) GFDL (2, 6),  T=300K Isothermal, H2O Lines   ', $
        '(46) GFDL (2, 7),  MLS, H2O + e-type, 400-1200 cm-1   ', $
        '(47) GFDL (2,15),  MLS, H2O (75%) + Cont + CO2 (300ppm) + O3  ', $
        '(48) GFDL (2,16),  MLS, H2O Lines (75%) + CO2 (300ppm) + O3  ', $
        '(49) GFDL (2,17),  MLS, H2O (125%) + Cont + CO2 (300ppm) + O3  ', $
        '(50) GFDL (2,18),  MLS, H2O Lines (125%) + CO2 (300ppm) + O3  ']

      print, ' '
      print, ' Select GFDL LBL case:'
      print, ' '
      for i=1,50 do print, '     ',titlenm(i-1)
      print, ' '
      read, fname
      if fname lt 1 or fname gt 50 then goto, gfdltop
      clrid = 3
   endif

; *****KIAE*****
kiaetop:
   if gname eq 4 then begin
      fnmstr = [ $
                'kiae.ph1.c09', $
                'kiae.ph1.c18', $
                'kiae.ph1.c20', $
                'kiae.ph1.c22', $
                'kiae.ph1.c23', $
                'kiae.ph1.c25', $
                'kiae.ph1.c27', $
                'kiae.ph1.c27x', $
                'kiae.ph1.c29', $
                'kiae.ph1.c31', $
                'kiae.ph1.c33', $
                'kiae.ph1.c35'  ]

      titlenm = [ $
        '(1)  KIAE (1, 9),  MLS, CO2 (300ppm)   ', $
        '(2)  KIAE (1,18),  MLS, H2O Lines (75% concentration)   ', $
        '(3)  KIAE (1,20),  MLS, H2O Lines (100% concentration)  ', $
        '(4)  KIAE (1,22),  MLS, H2O Lines (125% concentration)   ', $
        '(5)  KIAE (1,23),  MLS, Ozone   ', $
        '(6)  KIAE (1,25),  TRP, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(7)  KIAE (1,26),  TRP, H2O + Cont + CO2 (600ppm) + O3   ', $
        '(8)  KIAE (1,27),  MLS, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(9)  KIAE (1,27x), MLS, H2O + Cont + CO2 + O3 + CH4 + N2O   ', $
        '(10) KIAE (1,29),  MLW, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(11) KIAE (1,31),  SAS, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(12) KIAE (1,33),  SAW, H2O + Cont + CO2 (300ppm) + O3   ', $
        '(13) KIAE (1,35),  MLS, H2O Lines + CO2 (300ppm) + O3   ' ]

      print, ' '
      print, ' Select KIAE LBL case:'
      print, ' '
      for i=1,13 do print, '     ',titlenm(i-1)
      print, ' '
      read, fname
      if fname lt 1 or fname gt 13 then goto, kiaetop
      clrid = 1
   endif

   print, 'Please wait for data input ...'

; Open data file
openr, 1, '../idl_data/' + fnmstr(fname-1)

; read initial header, including control parameters
   text = strarr(1)

   nlyrs=0 & nlvls=0 & nbnds=0 & tsfc=0. & ptrp=0.
   readu, 1, nlyrs, nlvls, nbnds, psfc, tsfc, ptrp

   wnl=0. & wnu=0. & fus=0. & fds=0. & fur=0. & fdr=0. & fut=0.
   readu, 1, wnl, wnu, fus, fds, fur, fdr, fut

   if gname eq 2 then wnu = wnu+0.01
   wni = (wnu-wnl)/float(nbnds)

   wnls = fix(wnl)
   wnus = fix(wnu)
   wnis = fix(wni)

; read layer data
   data1 = fltarr(6,nlyrs)
   readu, 1, data1
   plyr = data1(0,*)
   tlyr = data1(1,*)
   air  = data1(2,*)
   h2o  = data1(3,*)
   o3   = data1(4,*)
   htrt = data1(5,*)

; read level data
   data2 = fltarr(5,nlvls)
   readu, 1, data2
   plvl = data2(0,*)
   tlvl = data2(1,*)
   flxu = data2(2,*)
   flxd = data2(3,*)
   flxn = data2(4,*)

; read spectral band data
   band1 = fltarr(nbnds)
   band2 = fltarr(nbnds)
   fusfc = fltarr(nbnds)
   fdsfc = fltarr(nbnds)
   futrp = fltarr(nbnds)
   fdtrp = fltarr(nbnds)
   futoa = fltarr(nbnds)
   fup = fltarr(nbnds,nlvls)
   fdn = fltarr(nbnds,nlvls)
   fnt = fltarr(nbnds,nlvls)
   htr = fltarr(nbnds,nlyrs)
   data3 = fltarr(nlvls)
   data4 = fltarr(nlyrs)
   for i=0,nbnds-1 do begin
      readu, 1, bnd1, bnd2, fus, fds, fur, fdr, fut
        band1(i) = bnd1
        band2(i) = bnd2
        fusfc(i) = fus
        fdsfc(i) = fds
        futrp(i) = fur
        fdtrp(i) = fdr
        futoa(i) = fut
      readu, 1, data3
        fup(i,*) = data3
      readu, 1, data3
        fdn(i,*) = data3
      readu, 1, data3
        fnt(i,*) = data3
      readu, 1, data4
        htr(i,*) = data4
   endfor
   close,1

; top of plotting loop

; set up color table (b,w,r,g,b)
  r = [0,1,1,0,0]
  g = [0,1,0,1,0.6]
  b = [0,1,0,0,1]
  tvlct, r*255, g*255, b*255

plttop:

!p.color = 1
print, ' '
print, 'Select plotting option:'
print, ' 1) Atmosphere profiles'
print, ' 2) Spectrally integrated flux profiles'
print, ' 3) Spectrally integrated cooling rate profiles'
print, ' 4) Spectral flux cross-sections'
print, ' 5) Spectral cooling rate cross-sections'
print, ' 6) Plot another LBL case'
print, ' 7) Set X/Y axis scaling'
print, ' 8) Change plotting window'
print, ' 9) Exit program'
read, typin
if typin eq 6 then goto, lbltop
if typin eq 9 then goto, pltend
if typin lt 1 or typin gt 9 then goto, plttop
scltop:
if typin eq 7 then begin
   print, ' '
   print, 'Select scaling for X/Y axes:'
   print, ' 1) Linear/Linear'
   print, ' 2) Linear/Log'
;   print, ' 3) Log/Linear'
;   print, ' 4) Log/Log'
   read, sclin
   if sclin eq 1 then  print, 'Plot axis scaling set to: X - linear; Y - linear'
   if sclin eq 2 then  print, 'Plot axis scaling set to: X - linear; Y - log'
   if sclin lt 1 or sclin gt 2 then goto, scltop
;   if sclin eq 1 then plotstr = 'plot'
;   if sclin eq 2 then plotstr = 'plot_io'
;   if sclin eq 3 then plotstr = 'plot_oi'
;   if sclin eq 4 then plotstr = 'plot_oo'
   goto, plttop
endif
if typin eq 8 then begin
windin:
   print, ' '
   print, 'Select output window number:'
   print, ' 0) Window 0'
   print, ' 1) Window 1'
   print, ' 2) Window 2'
   print, ' 3) Window 3'
   read, wndin
   if wndin lt 0 or wndin gt 3 then goto, windin
   if winflg(wndin) eq 0 then begin
      window, wndin, xsize=420, ysize=300
      winflg(wndin) = 1
      crfac=0.7
   endif
   if winflg(wndin) eq 1 then begin
      wset, wndin
      if wndin eq 0 then crfac=1.0
      if wndin ge 1 then crfac=0.7
   endif
   pltflg = 0
   goto, plttop
endif

; Atmosphere profiles
if typin eq 1 then begin
dtin:
   print, ' '
   print, 'Choose one:'
   print, ' 1) Temperature'
   print, ' 2) Air column amount'
   print, ' 3) H2O mixing ratio'
   print, ' 4) O3 mixing ratio'
   read, datin
   if datin lt 1 or datin gt 4 then goto, dtin
   if datin eq 1 then begin
      datax = tlvl
      xt = 'Temperature (K)'
   endif
   if datin eq 2 then begin
      datax = air
      xt = 'Air Column Amount (molecules/cm!e2!n)'
   endif
   if datin eq 3 then begin
      datax = h2o
      xt = 'H2O mixing ratio (ppm)'
   endif
   if datin eq 4 then begin
      datax = o3
      xt = 'O3 mixing ratio (ppm)'
   endif

   if pltflg eq 1 then begin
      print, ' '
      print, 'Overlay new plot on current plot? (Y/N)'
      read, ovrin
      if ovrin ne 'y' then pltflg = 0
      if ovrin eq 'y' and sclin eq 1 then goto, ovrly1a
      if ovrin eq 'y' and sclin eq 2 then goto, ovrly2a
      if sclin eq 2 then goto, pltloga
   endif

;!p.region = [0,0,1,0.5]
   plot, datax, plvl, $
     charsize = 1.3*crfac, $
     title = strmid(titlenm(fname-1),0,9)+' Atmosphere', $
     xtitle = xt, $
     ytitle = 'Pressure (mb)', $
     yrange = [1013,0], $
     ytickv = [1000,800,600,400,200,0], $
     yticks = 5, $
     yminor = 4
   pltflg = 1

ovrly1a:
!p.color = clrid
   oplot, datax, plvl
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

pltloga:
;!p.region = [0,0.5,1,1]
   plot_io, datax, plvl, $
     charsize = 1.3*crfac, $
     title = strmid(titlenm(fname-1),0,9)+' Atmosphere', $
     xtitle = xt, $
     ytitle = 'Pressure (mb)', $
     yrange = [1013,0.1], $
     ytickv = [1000,100,10,1,0.1], $
     yticks = 4, $
     yminor = 9
   pltflg = 1

ovrly2a:
!p.color = clrid
   oplot, datax, plvl
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

endif

; Spectrally integrated flux profiles
if typin eq 2 then begin
   print, ' '
   print, 'Selected flux type:'
   print, ' 1) Up'
   print, ' 2) Down'
   print, ' 3) Net'
   read, flxin
wnin:
   print, ' '
   print, 'Input wavenumber range to plot, (RETURN for full spectral range):'
   read, wn1s, wn2s
   wn1 = float(wn1s) & wn2 = float(wn2s)
   if wn1s eq '' and wn2s eq '' then begin
      wn1 = wnl & wn2 = wnu
   endif
   if wn2 le wn1 then begin
      print, ' '
      print, 'Upper wavenumber must be greater than lower wavenumber.'
      goto, wnin
   endif
   if wn1 lt wnl or wn2 gt wnu then begin
      print, ' '
      print, 'Available spectral range is ',wnls,' - ',wnus
      goto, wnin
   endif
   if wn1/wni ne fix(wn1/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin
   endif
   if wn2/wni ne fix(wn2/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin
   endif
; Compute flux sum for desired region
   if wn1 eq wnl and wn2 eq wnu then begin
      if flxin eq 1 then sumx = flxu
      if flxin eq 2 then sumx = flxd
      if flxin eq 3 then sumx = flxn
      goto, flxlbl
   endif
   idx1 = (wn1-wnl)/wni
   idx2 = (wn2-wnl)/wni
   sumx = fltarr(nlvls)
   for i=0,nlvls-1 do sumx(i) = 0.
   if flxin eq 1 then arrx = fup
   if flxin eq 2 then arrx = fdn
   if flxin eq 3 then arrx = fnt
   for i = 0,nlvls-1 do begin
      for j=idx1,idx2-1 do begin
         sumx(i) = sumx(i) + arrx(j,i)
      endfor
   endfor
flxlbl:
   if flxin eq 1 then fstr = 'Up '
   if flxin eq 2 then fstr = 'Down '
   if flxin eq 3 then fstr = 'Net '
   ttl = strarr(1)
   ttl = titlenm(fname-1) + string(fix(wn1),format='(i4)') + ' - ' + string(fix(wn2),format='(i4)') + ' cm!e-1!n'

   if pltflg eq 1 then begin
      print, ' '
      print, 'Overlay new plot on current plot? (Y/N)'
      read, ovrin
      if ovrin ne 'y' then pltflg = 0
      if ovrin eq 'y' and sclin eq 1 then goto, ovrly1b
      if ovrin eq 'y' and sclin eq 2 then goto, ovrly2b
      if sclin eq 2 then goto, pltlogb
   endif

;!p.region = [0,0,1,0.5]
   plot, sumx, plvl, $
     charsize = 1.3*crfac, $
     title = ttl, $
     xtitle = fstr+'Flux (W/m2)', $
     ytitle = 'Pressure (mb)', $
     yrange = [1013,0], $
     ytickv = [1000,800,600,400,200,0], $
     yticks = 5, $
     yminor = 4
   pltflg = 1

ovrly1b:
!p.color = clrid
   oplot, sumx, plvl
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

pltlogb:
;!p.region = [0,0.5,1,1]
   plot_io, sumx, plvl, $
     charsize = 1.3*crfac, $
     title = ttl, $
     xtitle = fstr+'Flux (W/m2)', $
     ytitle = 'Pressure (mb)', $
     yrange = [1013,0.1], $
     ytickv = [1000,100,10,1,0.1], $
     yticks = 4, $
     yminor = 9
   pltflg = 1

ovrly2b:
!p.color = clrid
   oplot, sumx, plvl
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

endif

; Spectrally integrated cooling rate profiles
if typin eq 3 then begin
wnin1:
   print, ' '
   print, 'Input wavenumber range to plot, (RET for full spectral range):'
   read, wn1s, wn2s
   wn1 = float(wn1s) & wn2 = float(wn2s)
   if wn1s eq '' and wn2s eq '' then begin
      wn1 = wnl & wn2 = wnu
   endif
   if wn2 le wn1 then begin
      print, ' '
      print, 'Upper wavenumber must be greater than lower wavenumber.'
      goto, wnin1
   endif
   if wn1 lt wnl or wn2 gt wnu then begin
      print, ' '
      print, 'Available spectral range is ',wnls,' - ',wnus
      goto, wnin1
   endif
   if wn1 eq '' then wn1 = wnl
   if wn2 eq '' then wn2 = wnu
   if wn1/wni ne fix(wn1/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin1
   endif
   if wn2/wni ne fix(wn2/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin1
   endif
; Compute flux sum for desired region
   if wn1 eq wnl and wn2 eq wnu then begin
      sumx1 = -htrt
      goto, htrlbl
   endif
   idx1 = (wn1-wnl)/wni
   idx2 = (wn2-wnl)/wni
   sumx1 = fltarr(nlyrs)
   for i=0,nlyrs-1 do sumx1(i) = 0.
; Change sign for cooling rate
   arrx1 = -htr
   for i=0,nlyrs-1 do begin
      for j=idx1,idx2-1 do begin
         sumx1(i) = sumx1(i) + arrx1(j,i)
      endfor
   endfor
htrlbl:
   ttl = titlenm(fname-1) + string(fix(wn1),format='(i4)') + ' - ' + string(fix(wn2),format='(i4)') + ' cm!e-1!n'

   if pltflg eq 1 then begin
      print, ' '
      print, 'Overlay new plot on current plot? (Y/N)'
      read, ovrin
      if ovrin ne 'y' then pltflg = 0
      if ovrin eq 'y' and sclin eq 1 then goto, ovrly1c
      if ovrin eq 'y' and sclin eq 2 then goto, ovrly2c
      if sclin eq 2 then goto, pltlogc
   endif

;!p.region = [0,0,1,0.5]
   plot, sumx1, plyr, $
     charsize = 1.3*crfac, $
     title = ttl, $
     xtitle = 'Cooling Rate (K/day)', $
     ytitle = 'Pressure (mb)', $
     yrange = [1013,0], $
     ytickv = [1000,800,600,400,200,0], $
     yticks = 5, $
     yminor = 4
   pltflg = 1
   oplot, [0,0], [1013,0]

ovrly1c:
!p.color = clrid
   oplot, sumx1, plyr
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

pltlogc:
;!p.region = [0,0.5,1,1]
   plot_io, sumx1, plyr, $
     charsize = 1.3*crfac, $
     title = ttl, $
     xtitle = 'Cooling Rate (K/day)', $
     ytitle = 'Pressure (mb)', $
     yrange = [1013,0.1], $
     ytickv = [1000,100,10,1,0.1], $
     yticks = 4, $
     yminor = 9
   pltflg = 1
   oplot, [0,0], [1013,0.1]

ovrly2c:
!p.color = clrid
   oplot, sumx1, plyr
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

endif

; Spectral flux cross-sections
if typin eq 4 then begin
wnin2:
   print, ' '
   print, 'Input wavenumber range to plot, (RET for full spectral range):'
   read, wn1s, wn2s
   wn1 = float(wn1s) & wn2 = float(wn2s)
   if wn1s eq '' and wn2s eq '' then begin
      wn1 = wnl & wn2 = wnu
   endif
   if wn2 le wn1 then begin
      print, ' '
      print, 'Upper wavenumber must be greater than lower wavenumber.'
      goto, wnin2
   endif
   if wn1 lt wnl or wn2 gt wnu then begin
      print, ' '
      print, 'Available spectral range is ',wnls,' - ',wnus
      goto, wnin2
   endif
   if wn1 eq '' then wn1 = wnl
   if wn2 eq '' then wn2 = wnu
   if wn1/wni ne fix(wn1/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin2
   endif
   if wn2/wni ne fix(wn2/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin2
   endif

   print, ' '
   print, 'Selected flux type:'
   print, ' 1) Up'
   print, ' 2) Down'
   print, ' 3) Net'
   read, flxin

prin:
   print, ' '
   print, 'Input pressure level (in mb) for cross-section plot;'
   print, '  (i.e. 500; the nearest available level will be plotted), or:'
   print, ' 1) Top of atmosphere'
   print, ' 2) Tropopause'
   print, ' 3) Surface'
   read, prsin
   if prsin lt 0 or prsin gt 1050 then goto, prin
   if prsin eq 1 then prsin = 0.
   if prsin eq 2 then prsin = ptrp
   if prsin eq 3 then prsin = psfc


; Find nearest pressure level to requested level
   psxn = min(abs(plvl-prsin),pmnind)
   psxn = plvl(pmnind)
   if flxin eq 1 then arry = fup(*,pmnind)/wni
   if flxin eq 1 then fstr = 'Up '
   if flxin eq 2 then arry = fdn(*,pmnind)/wni
   if flxin eq 2 then fstr = 'Down '
   if flxin eq 3 then arry = fnt(*,pmnind)/wni
   if flxin eq 3 then fstr = 'Net '
   ttl = strarr(1)
   ttl = titlenm(fname-1) + string(psxn,format='(f8.2)') + ' mb'

   idx1 = (wn1-wnl)/wni
   idx2 = (wn2-wnl)/wni
   xwn = fltarr((wn2-wn1)/wni)
   for i=0,((wn2-wn1)/wni)-1 do xwn(i)=(wn1+i*wni+wni/2.)

   if pltflg eq 1 then begin
      print, ' '
      print, 'Overlay new plot on current plot? (Y/N)'
      read, ovrin
      if ovrin ne 'y' then pltflg = 0
      if ovrin eq 'y' then goto, ovrly1d
   endif

;!p.region = [0,0,1,0.5]
   plot, xwn, arry(idx1:idx2-1), $
     charsize = 1.3*crfac, $
     title = ttl, $
     xtitle = 'Wavenumber (cm!e-1!n)', $
     xrange = [wn1,wn2], $
     ytitle = fstr+'Flux (W m!e-2!n/cm!e-1!n)'
   pltflg = 1

ovrly1d:
!p.color = clrid
   oplot, xwn, arry(idx1:idx2-1)
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

endif

; Spectral cooling rate cross-sections
if typin eq 5 then begin
wnin3:
   print, ' '
   print, 'Input wavenumber range to plot, (RET for full spectral range):'
   read, wn1s, wn2s
   wn1 = float(wn1s) & wn2 = float(wn2s)
   if wn1s eq '' and wn2s eq '' then begin
      wn1 = wnl & wn2 = wnu
   endif
   if wn2 le wn1 then begin
      print, ' '
      print, 'Upper wavenumber must be greater than lower wavenumber.'
      goto, wnin3
   endif
   if wn1 lt wnl or wn2 gt wnu then begin
      print, ' '
      print, 'Available spectral range is ',wnls,' - ',wnus
      goto, wnin3
   endif
   if wn1 eq '' then wn1 = wnl
   if wn2 eq '' then wn2 = wnu
   if wn1/wni ne fix(wn1/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin3
   endif
   if wn2/wni ne fix(wn2/wni) then begin
      print, 'Wavenumber value must be a multiple of ',wnis
      goto, wnin3
   endif

prin2:
   print, ' '
   print, 'Input pressure layer (in mb) for cross-section plot;'
   print, '  (i.e. 500; the nearest available layer will be plotted), or:'
   print, ' 1) Top of atmosphere'
   print, ' 2) Tropopause'
   print, ' 3) Surface'
   read, prsin
   if prsin lt 0 or prsin gt 1050 then goto, prin2
   if prsin eq 1 then prsin = 0.
   if prsin eq 2 then prsin = ptrp
   if prsin eq 3 then prsin = psfc


; Find nearest pressure level to requested level
   psxn = min(abs(plyr-prsin),pmnind)
   psxn = plyr(pmnind)
   arry = -htr(*,pmnind)/wni
   ttl = strarr(1)
   ttl = titlenm(fname-1) + string(psxn,format='(f8.2)') + ' mb'

   idx1 = (wn1-wnl)/wni
   idx2 = (wn2-wnl)/wni
   xwn = fltarr((wn2-wn1)/wni)
   for i=0,((wn2-wn1)/wni)-1 do xwn(i)=(wn1+i*wni+wni/2.)

   if pltflg eq 1 then begin
      print, ' '
      print, 'Overlay new plot on current plot? (Y/N)'
      read, ovrin
      if ovrin ne 'y' then pltflg = 0
      if ovrin eq 'y' then goto, ovrly1e
   endif

;!p.region = [0,0,1,0.5]
   plot, xwn, arry(idx1:idx2-1), $
     charsize = 1.3*crfac, $
     title = ttl, $
     xtitle = 'Wavenumber (cm!e-1!n)', $
     xrange = [wn1,wn2], $
     ytitle = 'Cooling rate (K/day/cm!e-1!n)'
   pltflg = 1
   oplot, [wn1,wn2], [0,0]

ovrly1e:
!p.color = clrid
   oplot, xwn, arry(idx1:idx2-1)
   if clrid eq 2 then xyouts, 0.70,0.03, 'AER', /normal, size=1.2
   if clrid eq 4 then xyouts, 0.77,0.03, 'GLA', /normal, size=1.2
   if clrid eq 3 then xyouts, 0.84,0.03, 'GFDL', /normal, size=1.2
   if clrid eq 1 then xyouts, 0.91,0.03, 'KIAE', /normal, size=1.2
   goto, plttop

endif

pltend:

end