; Huisheng Bian (6/17/2020)
; read in ATom obs data
; read in GEOS model data that has been sampled along flighttrack
; The GEOS DATA ia a curtain data along flight track


  cobstime = '20180104' ; 20170402, 20170602
  cobstimef2 = '2018-01-04'
  watom = 'atom2'
  caseid = 'atom1to4'

 pickvalue = ['SS']
 ;pickvalue = ['BC','OA','SO4','NO3','NH4','MSA','DU','SS','rMSASO4','HNO3','CO_NOAA','CO_QCLS','CO_UCATS']
 ttrc = n_elements(pickvalue)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Read in the ATom Observational data
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  fdatain = 'Mor.all.123.'+cobstimef2+'.tbl.'+watom
  fileinp = 'atomobs.'+fdatain
  openr, lun, fileinp, /get
; count line
  j = 0L
  str = 'a'
  while(not(eof(lun))) do begin
   readf, lun, str
   j = j+1
  endwhile
  print,'j =',j
; reopen file and extract data
  free_lun, lun
  openr, lun, fileinp, /get
; get header line
  readf, lun, str
  data = fltarr(j-1,24)   ; 5 columns
  nymd = lonarr(j-1)   ; 5 columns
  for i = 0L, j-2 do begin
   readf, lun, str
   tmpstr = strsplit(str,' ',/extract)
   tmptime = strsplit(tmpstr(0),'-T:',/extract)
   nymd[i] = long(tmptime(0))*10000+long(tmptime(1))*100+long(tmptime(2))
   data[i,*] = float(tmpstr(1:24))
  endfor
  free_lun, lun
  datestr = strpad(nymd,10000000)
  fltlon = data[*,0]
  fltlat = data[*,1]
  fltpress = data[*,2]                       ; mb, hPa
  flttemp = data[*,3]                        ; k
  fltalt = data[*,5]/1000.d0                 ; km, G_ALT
  fltstdtovol = data[*,12]                   ; m3/sm3
  vstd = fltpress/1013.25*273.15/flttemp
  fltoaocr = data[*,13]                      ; oa/oc
  fltbc = data[*,6]/fltstdtovol                         ; from ng/sm3 to ng/m3
  ;;;;;;;;;; fltocams is OA data you are interested in ;;;;;;;;;;;;;;;
  fltocams = data[*,7]/fltstdtovol*1000.d0              ; from ug/sm3 to ng/m3

  fltso4ams = data[*,8]/fltstdtovol*1000.d0             ; from ug/sm3 to ng/m3
  fltno3ams = data[*,9]/fltstdtovol*1000.d0             ; from ug/sm3 to ng/m3
  fltnh4ams = data[*,10]/fltstdtovol*1000.d0            ; from ug/sm3 to ng/m3
  fltssams = data[*,11]/fltstdtovol**1000.d0            ; from ug/sm3 to ng/m3
  fltmsaams = data[*,14]/fltstdtovol*1000.d0            ; from ug/sm3 to ng/m3
  fltdu = data[*,15]/fltstdtovol*1000.d0                ; from ug/sm3 to ng/m3
  fltss = data[*,16]/fltstdtovol*1000.0                 ; from ug/sm3 to ng/m3
  fltrmsaso4 = data[*,17]             ; ratio from 0 to 1
  flthno3 = data[*,23]                ; in pptv
  fltconoaa = data[*,19]              ; in ppbv
  fltcoqcls = data[*,20]              ; in ppbv
  fltcoucats = data[*,21]             ; in ppbv

  sizeobs = size(fltocams,/dimensions)
  print,'sizeobs =', sizeobs

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Read in GEOS Model results 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  cdfid = ncdf_open('GEOS5_'+caseid+'_'+cobstime+'_'+watom+'.nc')
  id = ncdf_varid(cdfid,'time')
  ncdf_varget, cdfid, id, time
  id = ncdf_varid(cdfid,'trjLon')
  ncdf_varget, cdfid, id, lon
  id = ncdf_varid(cdfid,'trjLat')
  ncdf_varget, cdfid, id, lat
  id = ncdf_varid(cdfid,'delp')
  ncdf_varget, cdfid, id, delp

  id = ncdf_varid(cdfid,'BCphobic')
  ncdf_varget, cdfid, id,bc1 
  id = ncdf_varid(cdfid,'BCphilic')
  ncdf_varget, cdfid, id,bc2 
  id = ncdf_varid(cdfid,'OCphobic')
  ncdf_varget, cdfid, id,oc1 
  id = ncdf_varid(cdfid,'OCphilic')
  ncdf_varget, cdfid, id,oc2 
  id = ncdf_varid(cdfid,'SO4')
  ncdf_varget, cdfid, id, so4
  id = ncdf_varid(cdfid,'NO3an1')
  ncdf_varget, cdfid, id,no3an1 
  id = ncdf_varid(cdfid,'NH4a')
  ncdf_varget, cdfid, id,nh4 
  id = ncdf_varid(cdfid,'MSA')
  ncdf_varget, cdfid, id,msa 
  id = ncdf_varid(cdfid,'du001')  ; radius: 0.1 - 1.0 um
  ncdf_varget, cdfid, id,du001
  id = ncdf_varid(cdfid,'du002')  ; radious: 1.0 - 1.8 um
  ncdf_varget, cdfid, id,du002
  id = ncdf_varid(cdfid,'ss001')  ; radius: 0.01 - 0.1 um
  ncdf_varget, cdfid, id,ss001
  id = ncdf_varid(cdfid,'ss002')  ; radious: 0.1 - 0.5 um
  ncdf_varget, cdfid, id,ss002
  id = ncdf_varid(cdfid,'ss003')  ; radious: 0.5 - 1.5 um
  ncdf_varget, cdfid, id,ss003
  id = ncdf_varid(cdfid,'HNO3CONC')
  ncdf_varget, cdfid, id,hno3concin
  ;id = ncdf_varid(cdfid,'CO')
  ;ncdf_varget, cdfid, id,coconcin
  id = ncdf_varid(cdfid,'AIRDENS')
  ncdf_varget, cdfid, id, rhoa
  ncdf_close, cdfid

  h = delp/9.8/rhoa/1000.d0   ; km
  print,'h =',h(*,0)

  bcconc = (bc1+bc2)*rhoa*1e12   ; ng/m3
  ;;;;;;;;;; occonc is OA data you are interested in ;;;;;;;;;;;;;;;
  occonc = (oc1+oc2)*rhoa*1e12   ; ng/m3 ; GEOS5 result in OA directly 

  so4conc = so4*rhoa*1e12   ; ng/m3
  no3conc = no3an1*rhoa*1e12   ; ng/m3
  nh4conc = nh4*rhoa*1e12   ; ng/m3
  msaconc = msa*rhoa*1e12   ; ng/m3
  duconc = (du001+du002*0.494)*rhoa*1e12   ; ng/m3
  ssconc = (ss001+ss002+ss003)*rhoa*1e12   ; ng/m3
  rmsaso4 = msaconc/so4conc   ; 
  hno3conc = hno3concin/rhoa*28.96/63.0*1e12   ; from kg/m3 to pptv
  ;coconc = coconcin*1e9   ; mol/mol to ppbv

  sizemod = size(nh4conc,/dimensions)
  print,'sizemod =', sizemod(1)

; make sure obs and mod have the same size
  if sizeobs ne sizemod(1) then begin
     print,'sizeobs, sizemod =', sizeobs, sizemod(1)
     stop
  endif

  end