; Huisheng Bian, Dec., 2012
; read in ARCTAS aircraft measurement and plot timeseries of the aerosol fields
; User's choice :

;############# Choose time ##############
yr='2008'
mn1=7 & mn2=mn1
mm = '0'+strcompress(string(mn1),/rem)
if mm eq '04' then dd =['01','04','05','08','09','12','16','17','19']
if mm eq '06' then dd =['18','20','22','24','26','29']
if mm eq '07' then dd =['01','04','05','08','09','10','13']
if mm eq '04' then mon='April'
if mm eq '06' then mon='June'
if mm eq '07' then mon='July'

;############# Choose resolution ##############
res='p5xp6'  ; 1x125 or 2x25 or p5xp6
dt='3hr'
zonal=2  ; plot timeseries
tag = 1
ifall = 0

;########################################


; Color index: 0 --- black          1 --- Magenta       2 --- Cyan
;              3 --- Yellow         4 --- Green         5 --- Red
;              6 --- Blue           7 --- Navy          8 --- Gold
;              9 --- Pink          10 --- Aqua         11 --- Orchid
;             12 --- Gray          13 --- Sky          14 --- Beige
;             15 --- White
;Index 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
RED = [   0, 255,   0, 255,   0, 255,   0,   0, 255, 255, 112, 219, 127,   0, 255, 255]
GREEN=[   0,   0, 255, 255, 255,   0,   0,   0, 187, 127, 219, 112, 127, 163, 171, 255]
BLUE =[   0, 255, 255,   0,   0,   0, 255, 115,   0, 127, 147, 219, 127, 255, 127, 255]
TVLCT, red, green, blue

x1=fltarr(3)
x2=fltarr(3)
x1(0)=0.04 & x2(0)=0.32
x1(1)=0.36 & x2(1)=0.64
x1(2)=0.68 & x2(2)=0.96
y1=fltarr(3)
y2=fltarr(3)
y1(2)=0.05 & y2(2)=0.30
y1(1)=0.35 & y2(1)=0.60
y1(0)=0.65 & y2(0)=0.90


plt='ps'  ; postscript
set_plot,plt
if plt eq 'ps' then begin
 !psym=0
 !p.font=0
 !p.multi=[0,3,3,0,0]
 device,/symbol,font_index=4
device,/landscape,xoff=0.5,yoff=26,xsize=25,ysize=20
endif

; --- Read ARCTAS measurements:  ---

ndays=n_elements(dd(*))

dir = '/misc/mcc20/common/arctas/MERGE/1_MINUTE/DC8/'

; plot which ARCTAS measured aerosol quantities
  trc =['bc','ext']  ; tot and ext, or abs
  nns = n_elements(trc(*))
  tunits=['ppbv','Mm-1']
  if trc(1) eq 'ext' and ifall eq 0 then fign='multi-time-'+trc(0)+'.'+trc(1)+'.'+mm+'.R13.ps'
  device,file=fign,/color


for nd=0,ndays-1 do begin
    openr,1,dir+'mrg60_dc8_2008'+mm+dd(nd)+'_R13.ict'
    head=strarr(1)
    s1=''
    s2=''
    for i=0,385 do begin    ; R13 spring & summer
        readf,1,head
    endfor
    readf,1,s1
    j=0
    var=fltarr(150)
; var in locar are utc=0, JDAY =1, localtime=4, lat=5, lon=6, altp=7, pressure=8, temp=9, GPS_Alt=15, CO (ppbv)=49, SAGA SO4 (pptv) (Fine-Aerosol_Sulfate, ppt)=64, CIT-CIMS SO2 (pptv)=65, GT-CIMS SO2 (pptv)=66, H2SO4_CIMS=72 (molecule/cm3), Na-=117 (pptv), NH4= 118 (pptv), K=119 (pptv), Mg= 120 (pptv), Ca= 121 (pptv), Cl=122 (pptv), NO3=124 (pptv), SAGA SO4 = 125 (pptv), Chloride=130 (ug/m3), Ammonium= 131 (ug/m3), Nitrate= 132 (ug/m3), Organics<213mz= 133 (ug/m3), AMS SO4(ug/m3)=134, Chloride_PILS-IONS=138 (ug/m3), Nitrate_PILS-IONS=139 (ug/m3), Nitrite_PILS-IONS=140 (ug/m3), Sulfate_PILS-IONS=141, (ug/m3) Ammonium_PILS-IONS=142 (ug/m3), Calcium_PILS-IONS=143 (ug/m3), Sodium_PILS-IONS=144 (ug/m3), Potassium_PILS-IONS= 145 (ug/m3), SP2 BC=150 (ng/m3), Total_Scatter@550_nm = 176 (Mm-1), Backscatter@550_nm= 179 (Mm-1), Submicron_Scattering@550_nm= 182 (Mm-1), Angstron_Exponent_(Scattering)_450/700_nm= 184, Total_Absorption@532_nm= 186 (Mm-1), Submicron_Absorption@532_nm= 189 (Mm-1), Single_Scatter_Albedo_(Blue)= 198, Single_Scatter_Albedo_(Green)= 199, Single_Scatter_Albedo_(Red)= 200, Acetonitril_PTRMS = 215, Acetonitrile_TOGA = 236, OCS= 247 (pptv), DMS= 248 (pptv), Relative_Humidity_Water_(RHw)_LARGE-fRH= 326 (%), Relative_Humidity_Ice_(RHi)_LARGE-fRH= 327 (%), Wet-Scattering@532nm_LARGE-fRH= 328 (mm-1), Dry_Scattering@532nm_LARGE-fRH= 329 (Mm-1)   
    namevar =['utc','jday','lat','lon','altp','pres','temp','gpsalt','co','nh4s','cas','no3s','so4s','nh4a','no3a','orga','so4a','no3p','nitritep','so4p','nh4p','cap','nap','bc','tsca','subsca','tabs','subabs','ch3cn','rhw','rhi']
locvar=[0,1,5,6,7,8,9,15,49,118,121,124,125,131,132,133,134,137,138,139,140,141,142,150,176,182,186,189,215,326,327] ; co. bc, so2, so2, so4, so4, note: start from 0
    fct=[1,1,1,1,1,1,1,1,1,0.018/22.4,0.040078/22.4/0.06,0.0620049/22.4,0.09606/22.4,1,1,1,1,1,1,1,1,1./0.06,1,.001,1,1,1,1,1000.,1,1]
    nvar=n_elements(locvar(*))
    fvar1=fltarr(3600,nvar)
    var = dblarr(333)
    j=0

    WHILE NOT EOF(1) DO BEGIN
     readf,1,s2
     reads,s2,var

     for i=0,n_elements(locvar(*))-1 do begin
         fvar1(j,i)=var(locvar(i))
     endfor

     j=j+1
    ENDWHILE
    print,'j-1 =',j-1
    close,1

    fvara=fltarr(j,nvar)
    fvara(*,0)=fvar1(0:j-1,0)/3600.0
    for i=1,n_elements(locvar(*))-1 do begin
        fvara(*,i)=fvar1(0:j-1,i)*fct(i)
    endfor


    if zonal eq 2 then begin
       ; plot timeseries, set x-axis every 1800s (half hour)
       xv=fltarr(j-1)
       xv = indgen(j-1)+1
       xtkn = (fvara(j-1,0)-fvara(0,0))*2+1
       xtk = fltarr(xtkn)
       pxtk = fltarr(xtkn+1)
       xtk(0) = fvara(0,0)
       pxtk(0) = 1
       for i = 1,xtkn-1 do begin
           xtk(i) = xtk(i-1)+0.5
           pxtk(i) = pxtk(i-1)+30
       endfor
       pxtk(xtkn) = pxtk(xtkn-1)+30
       xlonnm = strcompress(string(xtk,format='(f4.1)'),/remove_all)
       xnon = replicate(' ',xtkn+1)
    endif

    ; cycle for plotting quantities
    for ns = 0,nns-1 do begin
    ; cal total aerosol mass in ug/m3
        ttot = fltarr(j) & ttot(*) = 0.0

        if ifall eq 0 then begin

        if trc(ns) ne 'ext' and trc(ns) eq 'abs' then begin
           vnh4a = fltarr(j) & vnh4a(*) = 0.0
           vno3a = fltarr(j) & vno3a(*) = 0.0
           vorga= fltarr(j) & vorga(*) = 0.0
           vso4a= fltarr(j) & vso4a(*) = 0.0
           vcas = fltarr(j) & vcas(*) = 0.0
           vbc  = fltarr(j) & vbc(*) = 0.0
           ccn  = fltarr(j) & ccn(*) = 0.0
           for i=0,nvar-1 do begin
              if namevar(i) eq 'nh4a' then begin
                 vnh4a(*) = fvara(*,i)
              endif
              if namevar(i) eq 'no3a' then begin
                 vno3a(*) = fvara(*,i)
              endif
              if namevar(i) eq 'orga' then begin
                 vorga(*) = fvara(*,i)
              endif
              if namevar(i) eq 'so4a' then begin
                 vso4a(*) = fvara(*,i)
              endif
              if namevar(i) eq 'cas' then begin
                 vcas(*) = fvara(*,i)
              endif
              if namevar(i) eq 'bc' then begin
                 vbc(*) = fvara(*,i)
              endif
              if namevar(i) eq 'ch3cn' then begin
                 ccn(*) = fvara(*,i)
              endif
              loc = where(vnh4a(*) > 0.0 and vno3a(*) > 0.0 and vorga(*) > 0.0 and vso4a(*) > 0.0 and vcas(*) > 0.0 and vbc(*) > 0.0)
              ttot(loc) = vnh4a(loc)+vno3a(loc)+vorga(loc)+vso4a(loc)+vcas(loc)+vbc(loc)
           endfor
           loc = where(ttot(*) eq 0.0)
           ttot(loc) = -1.0e+9
           ;print,'ttot=',ttot(*)
	endif  


        if trc(ns) eq 'ext' or trc(ns) eq 'abs' then begin
           ttot = fltarr(j) & ttot(*) = 0.0  
           ccn = fltarr(j) & ccn(*) = 0.0

           vtsca = fltarr(j) & vtsca(*) = 0.0
           vtabs = fltarr(j) & vtabs(*) = 0.0
           for i=0,nvar-1 do begin
               if namevar(i) eq 'tsca' then begin
                  vtsca(*) = fvara(*,i)
               endif 
               if namevar(i) eq 'tabs' then begin
                  vtabs(*) = fvara(*,i)
	          ;print,'vtabs =',vtabs(*)
               endif 
	       if namevar(i) eq 'ch3cn' then begin
                  ccn(*) = fvara(*,i)
               endif
           endfor
           if trc(ns) eq 'ext' then loc = where(vtsca(*) > 0.0 and vtabs(*) > 0.0)
           if trc(ns) eq 'ext' then ttot(loc) = vtsca(loc) + vtabs(loc) 
           if trc(ns) eq 'abs' then begin
	      loc = where(vtabs(*) > 0.0)
              ttot(loc) = vtabs(loc) 
	   endif
           ;;;;; bianfix needed ;;; sum sca and abs when both measurements are available !!!!!!
           loc = where(ttot(*) eq 0.0)
           ttot(loc) = -1.e9
        endif  ; ext
     endif ; ifall = 0
	print,'ttot =',ttot(*)

     ; check the max/min value of plotting quantities
     for i = 0,nvar-1 do begin
         if namevar(i) eq 'co' or namevar(i) eq 'nh4s' or namevar(i) eq 'no3s' or namevar(i) eq 'so4s' or namevar(i) eq 'nh4a' or namevar(i) eq 'no3a' or namevar(i) eq 'orga' or namevar(i) eq 'so4a' or namevar(i) eq 'no3p' or namevar(i) eq 'nitritep' or namevar(i) eq 'so4p' or namevar(i) eq 'nh4p' or namevar(i) eq 'bc' or namevar(i) eq 'tsca' or namevar(i) eq 'tabs' or namevar(i) eq 'rhw' or namevar(i) eq 'rhi' then begin
            print, 'max and min of '+namevar(i)+' =',max(fvara(*,i)),min(where(fvara(*,i) gt 0.0))
         endif
         if namevar(i) eq 'cas' or namevar(i) eq 'cap' then $
            print, 'max and min of '+namevar(i)+' =',max(fvara(*,i))/0.06,min(where(fvara(*,i)/0.06 gt 0.0))
         endfor

; assign points along flight track here
       nc = nd/3
       nr = nd-nd/3*3
       print,'nc=',nc,nr
       set_viewport,x1(nr),x2(nr),y1(nc),y2(nc)

       ;yr1 = max([0.0, min(tvar(*))*0.9])  & yr2=max(tvar(*))*1.1

       yr11 = 0.0  & yr21=max(vbc(*))*1.1
       if trc(ns) eq 'tot' or trc(ns) eq 'ext' or trc(ns) eq 'abs' then begin
          yr12 = 0.0  & yr22=max(ttot(*))*1.1
	  if yr22 le 0. then yr22 = 70.
       endif
       if nr eq 0 then begin
          yhd = trc(ns)+' ('+tunits(ns)+')'
       endif else begin
          yhd = ' '
       endelse
       print,'yr1=',yr1,yr2
       if ns eq 0 then begin
          plot,xv,vbc(*),/nodata,/noerase,xthick=3,xrange=[1,pxtk(xtkn)],xstyle=1, $
           yrange=[yr11,yr21],ystyle=4, $
	   xtit = '', ytit = '', title = '', $
	    charsize=1.5, $
           xtickname=xnon,xticks=xtkn, min_value = 0.0
          if ns eq 0 and nc eq 1 then xyouts, 0.15*((pxtk(xtkn)-1)),-0.25*yr21, 'hours within the day',charsize = 1.3
       endif
       for i = 0,nvar-1 do begin
           if namevar(i) eq trc(ns) then begin
              plot,xv,fvara(*,i),xrange=[1,pxtk(xtkn)],yrange=[yr11,yr21],xstyle=5,/noe,ystyle=9,COLOR=6,th=5,min_value=0,max_value=10000.,charsize=1.5
           endif
       endfor

       if ns eq 0 then xyouts,xv(j-2)*0.7,(yr21-yr11)*0.7,mm+dd(nd),charsize=1.0
       if ns eq 1 and nc eq 0 and nr eq 0 then xyouts,xv(j-2)*0.03,(yr22-yr12)*0.9,'BC (ng/m3)',color=6,charsize=1.0
       if ns eq 1 and nc eq 0 and nr eq 0 and trc(ns) eq 'ext' then xyouts,xv(j-2)*0.03,(yr22-yr12)*0.8,'Ext 532-550nm (Mm-1)',color=4,charsize=1.0
       if ns eq 1 and nc eq 0 and nr eq 0 and trc(ns) eq 'abs' then xyouts,xv(j-2)*0.03,(yr22-yr12)*0.8,'Abs 532-550nm (Mm-1)',color=4,charsize=1.0

       if ns eq 0 then begin
          for i=0,xtkn-1 do begin
              xyouts,pxtk(i),-(yr21-yr11)*0.08,xlonnm(i),charsize=0.7,orientation=40,alignment=0.5
          endfor
       endif

print,'Done ',fign
endfor  ; trc
 if ifall eq 0 then openw,lun,'ARCTAS-co-'+trc(1)+'-'+mm+dd(nd)+'.dat',/get_lun
 if ifall eq 1 then openw,lun,'ARCTAS-co-'+trc(1)+'-'+mm+dd(nd)+'.all.dat',/get_lun
 print,'j=',j
    for i = 0,j-1 do begin
        if fvara(i,8) gt 0.0 and ttot(i) gt 0.0 and ccn(i) gt 0.0 then begin
           printf,lun,format='(i5,3f10.3)',i,fvara(i,8),ttot(i),ccn(i)
	endif
    endfor
    free_lun,lun
endfor ; nday
device,/close
end