pro plot_spcza,species,scale=scale,_EXTRA=e


   if (n_elements(scale) eq 0) then $
      scale = 1.e12


   ; read observations

   gte_readbin,'c:/data/gte/pem-tb/dc8/merges/mrg_bwhoxd10.pmb.bdt',  $
        odata,ovardesc,time=otime
   gte_readbin,'c:/data/gte/pem-tb/dc8/merges/mrg_thatgd10.pmb.bdt',  $
        hdata,hvardesc,time=htime


   oheader = ovardesc.name
   ospec = [ 'ALTP', 'Ozone',  $
            'Carbon Monoxide mixing ratio', $
            'Carbon Dioxide mixing ratio'  ]
   osel = make_selection(oheader,ospec)
   if (min(osel) lt 0) then stop

   hheader = hvardesc.name
   hspec = [ 'ALTP', 'HNO3' ]
   hsel = make_selection(hheader,hspec)
   if (min(hsel) lt 0) then stop


   ; identify time intervals for the two SPCZ profiles
   timerange0 = [ 94888., 97365. ]
   timerange1 = [ 100768, 103148. ]

   oxind0 = where(otime ge timerange0[0] AND otime le timerange0[1])
   odata0 = odata[oxind0,*]
   oxind1 = where(otime ge timerange1[0] AND otime le timerange1[1])
   odata1 = odata[oxind1,*]

   hxind0 = where(htime ge timerange0[0] AND htime le timerange0[1])
   hdata0 = hdata[hxind0,*]
   hxind1 = where(htime ge timerange1[0] AND htime le timerange1[1])
   hdata1 = hdata[hxind1,*]


   oalt0 = odata0[*,osel[0]]
   o30 = odata0[*,osel[1]]
   co0 = odata0[*,osel[2]]
   co20 = odata0[*,osel[3]]
   halt0 = hdata0[*,hsel[0]]
   hno30 = hdata0[*,hsel[1]]

   oalt1 = odata1[*,osel[0]]
   o31 = odata1[*,osel[1]]
   co1 = odata1[*,osel[2]]
   co21 = odata1[*,osel[3]]
   halt1 = hdata1[*,hsel[0]]
   hno31 = hdata1[*,hsel[1]]


   ; do 4 panel plot with vertical profiles
   xyouts,0,0,'!6',/norm
   !p.charsize=1.3
   if (!d.name ne 'WIN') then !p.charthick = 2 else !p.charthick=1.2
   !x.charsize=1.1
   !y.charsize=1.1
   !x.minor = 1
   !y.minor = 2

   multipanel,4,margin=[0.06,0.08,0.,0.],omargin=[0.12,0.04,0.02,0.01]


   empty = replicate(' ',30)



   multipanel,position=p
   plot,o31,oalt1,/nodata,color=1,  $
        xtitle='O!L3!N [ppb]',xrange=[0.,60.], $
        ytitle='Altitude [km]',yrange=[0.,12.],xstyle=1,position=p,/noerase

   oplot,o30,oalt0,color=3,psym=sym(1),symsize=1.
   oplot,o31,oalt1,color=2,psym=sym(4),symsize=1.

   multipanel,/advance,position=p
   plot,hno31,halt1,/nodata,color=1,  $
        ytickname=empty,xtitle='HNO!L3!N [ppt]',xrange=[0.,300.], $
        yrange=[0.,12.],xstyle=1,position=p,/noerase

   oplot,hno30,halt0,color=3,psym=sym(1),symsize=1.6
   oplot,hno31,halt1,color=2,psym=sym(4),symsize=1.6

   multipanel,/advance,position=p
   plot,co1,oalt1,/nodata,color=1,  $
        xtitle='CO [ppb]',xrange=[30.,60.], $
        ytitle='Altitude [km]',yrange=[0.,12.],xstyle=1,position=p,/noerase

   oplot,co0,oalt0,color=3,psym=sym(1),symsize=1.
   oplot,co1,oalt1,color=2,psym=sym(4),symsize=1.

   multipanel,/advance,position=p
   plot,co21,oalt1,/nodata,color=1,  $
        ytickname=empty,xtitle='CO!L2!N [ppb]',xrange=[364., 368.], $
        yrange=[0.,12.],xstyle=1,position=p,/noerase

   oplot,co20,oalt0,color=3,psym=sym(1),symsize=1.
   oplot,co21,oalt1,color=2,psym=sym(4),symsize=1.




   xyouts,0.5,0.97,'DC-8 Fiji local 1', $
      /norm,align=0.5,charsize=1.8,color=1
   xyouts,0.5,0.94,'SPCZ characterization : 15!Uo!NS and 21!Uo!NS', $
      /norm,align=0.5,charsize=1.5,color=1



   return
end