; This procedure plots the spatial distribution of tropospheric 
; column amount of desired species over CO

PRO spl_trop_col_diff, outfile, trop_column_all, true_trop_column_all, $
	lat_all, lon_all, dd, hh, scenario, spc_ttl, sp, choice

 shapefile = 'inputs/Colorado_County_Boundaries.shp'

 psopen, file = outfile, ticklen = -200, font = 6, tfont = 6, $
 tcharsize = 120, charsize = 120, space1 = 150, space2 = 150, $
 space3 = 150, /PORTRAIT 

  cs, scale = 1, ncols = 22

  xsize = 16000
  ysize = 11000
  ;rlon = [-109.5, -101.5]
  ;rlat = [36.6, 41.4]
  rlon  = [-106.0, -103.5]
  rlat  = [39.2, 40.8]

  ; --------------------------------------------------------------------------
  ; Section-1 for TEMPO
  ; --------------------------------------------------------------------------
  if choice eq 'TEMPO' then begin 
  ; --------------------------------------------------------------------------
  ; variables for GLYX

  if sp eq 0 then begin ; this will be replaces by CASES for sp
   levs  = findgen(21)-10
   trop_column_all1 = trop_column_all/1.0e+14
   true_trop_column_all1 = true_trop_column_all/1.0e+14
   ;print, minmax(trop_column_all, /NAN)
   ;print, minmax(true_trop_column_all, /NAN)
   lev_min = -10
   lev_max = 10
   lev_stp = 1
   lev_dec = 0
   spc_title = 'Tropospheric Column '+spc_ttl+' ( x 10!E14!N mol-cm!E-2!N)'
  endif


  ; variables for H2O
  if sp eq 1 then begin ; this will be replaces by CASES for sp
   levs  = findgen(21)-10
   trop_column_all1 = trop_column_all/1.0e+22
   true_trop_column_all1 = true_trop_column_all/1.0e+22
   ;print, minmax(trop_column_all, /NAN)
   ;print, minmax(true_trop_column_all, /NAN)
   lev_min = -10
   lev_max = 10
   lev_stp =  1
   spc_title = 'Tropospheric Column '+spc_ttl+' ( x 10!E22!N mol-cm!E-2!N)'
  endif


  ; variables for HCHO
  if sp eq 2 then begin ; this will be replaces by CASES for sp
   levs  = findgen(27)-10
   trop_column_all1 = trop_column_all/1.0e+15
   true_trop_column_all1 = true_trop_column_all/1.0e+15
   ;print, minmax(trop_column_all, /NAN)
   ;print, minmax(true_trop_column_all, /NAN)
   lev_min = -10
   lev_max = 10
   lev_stp =  1
   spc_title = 'Tropospheric Column '+spc_ttl+' ( x 10!E15!N mol-cm!E-2!N)'
  endif


  ; variables for NO2
  if sp eq 3 then begin	; this will be replaces by CASES for sp
   levs  = findgen(21)-10
   trop_column_all1 = trop_column_all/1.0e+15
   true_trop_column_all1 = true_trop_column_all/1.0e+15
   ;print, minmax(trop_column_all1, /NAN)
   ;print, minmax(true_trop_column_all1, /NAN)
   lev_min = -10
   lev_max = 10
   lev_stp = 1
   spc_title = 'Tropospheric Column '+spc_ttl+' ( x 10!E15!N mol-cm!E-2!N)'
  endif


  ; variables for ozone
  if sp eq 4 then begin ; this will be replaces by CASES for sp
   levs  = findgen(21)*2 - 20
   trop_column_all1 = trop_column_all/(2.6867*1.0e+16)
   true_trop_column_all1 = true_trop_column_all/(2.6867*1.0e+16)
   ;print, minmax(trop_column_all, /NAN)
   ;print, minmax(true_trop_column_all, /NAN)
   lev_min = -20
   lev_max = 20
   lev_stp = 2
   spc_title = 'Tropospheric Column '+spc_ttl+' (DU)'  
  endif


  ; variables for SO2
  if sp eq 5 then begin ; this will be replaces by CASES for sp
   levs  = findgen(27)-12
   trop_column_all1 = trop_column_all/1.0e+14
   true_trop_column_all1 = true_trop_column_all/1.0e+14
   ;print, minmax(trop_column_all, /NAN)
   ;print, minmax(true_trop_column_all, /NAN)
   lev_min = -12
   lev_max = 12
   lev_stp = 1
   spc_title = 'Tropospheric Column '+spc_ttl+' ( x 10!E14!N mol-cm!E-2!N)'
  endif
  ; --------------------------------------------------------------------------

  ; ----------- TRUE Column ----------
  pos, xsize = xsize, ysize = ysize, xoffset = 2500, yoffset = 16000
  levs, min = lev_min, max = lev_max, step = lev_stp;, ndecs=lev_dec
  con_lc, reform(true_trop_column_all1), lon_all, lat_all, levs, rlon, rlat
  cgDrawShapes, shapefile ,/AutoDraw
  gplot, x = -103.3, y = 40.0, text = 'TRUE Column', col = 1, orientation=90
  axes, xstep = 1, ystep = 0.5, ndecs=1

  ; ----------- Derived Column ----------
  pos, xsize = xsize, ysize = ysize, xoffset = 2500, yoffset = 4000
  levs, min = lev_min, max = lev_max, step = lev_stp;, ndecs=lev_dec
  con_lc, reform(trop_column_all1), lon_all, lat_all, levs, rlon, rlat
  cgDrawShapes, shapefile ,/AutoDraw
  gplot, x = -103.3, y = 40.0, text = 'DERIVED Column', col = 1, orientation=90
  axes, xstep = 1, ystep = 0.5, ndecs=1

  colbar, coords = [2500, 2400, 18500, 3000], title = '', nth=2,/nolines
  gplot, x=10000, y=900, text = spc_title, col=1,/DEVICE
  gplot, x = 10000, y = 28000, text = string(scenario+' ['+dd+'-Aug-2014, '+hh+' UTC]'), col = 1, /device


 psclose, /noview
 ;spawn, 'convert -density 400 -trim idl.eps ' + outfile
 print, 'Plot is located here', outfile

 endif

 ; --------------------------------------------------------------------------
 ; Section-2 : For WRF variables
 ; --------------------------------------------------------------------------
 if choice eq 'WRF' then begin
 ; --------------------------------------------------------------------------

  ; ------------ Surface NO2 VMR ---------- 
   levs  = findgen(21)*2 - 20
   lev_min = -20
   lev_max = 20
   lev_stp = 2
   vmr_title = spc_ttl+' (ppb)'
   emi_title = spc_ttl+' (mol-km!E-2!N-hr!E-1!N)'

   ;- Mixing Ratio -
   pos, xsize = xsize, ysize = ysize, xoffset = 2500, yoffset = 16000
   levs, min = lev_min, max = lev_max, step = lev_stp;, ndecs=lev_dec, /UPPER;,/V
   con_lc, reform(trop_column_all), lon_all, lat_all, levs, rlon, rlat
   cgDrawShapes, shapefile ,/AutoDraw
   gplot, x = -103.3, y = 40.0, text = 'Mixing Ratio (ppb)', col = 1, orientation=90
   axes, xstep = 1, ystep = 0.5, ndecs=1

   ;- Emissions -
   pos, xsize = xsize, ysize = ysize, xoffset = 2500, yoffset = 4000
   levs, min = lev_min, max = lev_max, step = lev_stp;, ndecs=lev_dec, /UPPER;,/V
   con_lc, reform(true_trop_column_all), lon_all, lat_all, levs, rlon, rlat
   cgDrawShapes, shapefile ,/AutoDraw
;x 10!E-1!N
   gplot, x = -103.3, y = 40.0, text = 'Emissions ( mol-km!E-2!N-hr!E-1!N)',$
         col = 1, orientation=90
   axes, xstep = 1, ystep = 0.5, ndecs=1

  colbar, coords = [2500, 2400, 18500, 3000], title = '', nth=2,/nolines
  gplot, x=10000, y=900, text = 'Surface '+spc_ttl, col=1,/DEVICE
  gplot, x = 10000, y = 28000, text = string(scenario+' ['+dd+'-Aug-2014, '+hh+' UTC]'), col = 1, /device

  psclose, /noview
  ;spawn, 'convert -density 400 -trim idl.eps ' + outfile
  print, 'Plot is located here', outfile

 endif

END