function read_L2_no2,filename

data_object={scd:-1.e10,vcd:-1.e10,cp_ll:fltarr(2),cp_lr:fltarr(2),cp_ul:fltarr(2),cp_ur:fltarr(2),amf:-1.e10,$
  sza:-1.e10,scdprec:-1.e10,scdstrat:-1.e10,fltrop:0,rms:-1.e10,chi:-1.e10,fqf:0,vza:-1.e10,cld_frac:-1.e10,cld_pres:-1.e10,amfclear:-1.e10,amfcloud:-1.e10,$
  gpxQAflags:-32767L,vcdtrop:-1.e10,vcdstrat:-1.e10,lat:-1.e10,lon:-1.e10,amfgeo:-1.e10,amftrop:-1.e10,albedo:-1.e10,sigvcdtrop:-1.e10,$
  cld_rfrac:-1.e10,data:-1.e10,scd_orig:-1.e10,ghostcolumn:-1.e10,vcdtropm:-1.e10,terrainheight:0,pixelnr:0,time:0.D,psurf:-1.e10,$
  tm4terrainheight:0,kernel:ptr_new()}

res=h5_parse(filename,/read_data)

tagnames=tag_names(res.hdfeos.swaths)
ii=-1
dum=where(strlowcase(tagnames) eq 'dominono2' ,count)
if (count eq 1) then begin
  ii=dum[0]
endif

if ii le 0 then begin
  print,'no valid swaths present'
  stop
endif

dum=tag_names(res.hdfeos.swaths.(ii).data_fields)
dres=res.hdfeos.swaths.(ii) & res=dres & dres=0

amf=res.data_fields.airmassfactor._data * res.data_fields.airmassfactor.scalefactor._data[0]
amfgeo=res.data_fields.airmassfactorgeometric._data * res.data_fields.airmassfactorgeometric.scalefactor._data[0]
amftrop=res.data_fields.airmassfactortropospheric._data * res.data_fields.airmassfactortropospheric.scalefactor._data[0]
albedo=res.data_fields.surfacealbedo._data * res.data_fields.surfacealbedo.scalefactor._data[0]
vcd=res.data_fields.totalverticalcolumn._data * res.data_fields.totalverticalcolumn.scalefactor._data[0]
vcdtrop=res.data_fields.troposphericverticalcolumn._data * res.data_fields.troposphericverticalcolumn.scalefactor._data[0]
vcdstrat=res.data_fields.assimilatedstratosphericverticalcolumn._data * res.data_fields.assimilatedstratosphericverticalcolumn.scalefactor._data[0]
vcdtropstd=res.data_fields.troposphericverticalcolumnerror._data * res.data_fields.troposphericverticalcolumnerror.scalefactor._data[0]
ghostcolumn=res.data_fields.ghostcolumn._data * res.data_fields.ghostcolumn.scalefactor._data[0]
vcdtropm=res.data_fields.troposphericverticalcolumnmodel._data * res.data_fields.troposphericverticalcolumnmodel.scalefactor._data[0]
fltrop=res.data_fields.troposphericcolumnflag._data
kernel=res.data_fields.averagingkernel._data * res.data_fields.averagingkernel.scalefactor._data[0]
;  scd_orig=res.data_fields.originalslantcolumnamountno2._data * res.data_fields.originalslantcolumnamountno2.scalefactor._data[0]
psurf=res.data_fields.tm4surfacepressure._data * res.data_fields.tm4surfacepressure.scalefactor._data[0]
tm4terrainheight=res.data_fields.tm4terrainheight._data
rms=0. & chi=0. & fqf=0 & amfclear=0. & amfcloud=0.
dum=reform(res.geolocation_fields.longitudecornerpoints._data[*,*,0])
cp_ll=fltarr(2,n_elements(dum[*,0]),n_elements(dum[0,*])) & cp_ll[*]=!values.f_nan
cp_ul=cp_ll & cp_lr=cp_ll & cp_ur=cp_ll
cp_ll[0,*,*]= res.geolocation_fields.latitudecornerpoints._data[*,*,3]
cp_ll[1,*,*]=res.geolocation_fields.longitudecornerpoints._data[*,*,3]
cp_ul[0,*,*]= res.geolocation_fields.latitudecornerpoints._data[*,*,2]
cp_ul[1,*,*]=res.geolocation_fields.longitudecornerpoints._data[*,*,2]
cp_lr[0,*,*]= res.geolocation_fields.latitudecornerpoints._data[*,*,1]
cp_lr[1,*,*]=res.geolocation_fields.longitudecornerpoints._data[*,*,1]
cp_ur[0,*,*]= res.geolocation_fields.latitudecornerpoints._data[*,*,0]
cp_ur[1,*,*]=res.geolocation_fields.longitudecornerpoints._data[*,*,0]

scd=res.data_fields.slantcolumnamountno2._data*res.data_fields.slantcolumnamountno2.scalefactor._data[0]
scdstrat=res.data_fields.assimilatedstratosphericslantcolumn._data*res.data_fields.assimilatedstratosphericslantcolumn.scalefactor._data[0]
scdprec=res.data_fields.slantcolumnamountno2std._data* res.data_fields.slantcolumnamountno2std.scalefactor._data[0]
lat=res.geolocation_fields.latitude._data
lon=res.geolocation_fields.longitude._data
sza=res.geolocation_fields.solarzenithangle._data
vza=res.geolocation_fields.viewingzenithangle._data
meas_time=res.geolocation_fields.time._data
cld_frac=res.data_fields.CloudFraction._data * res.data_fields.CloudFraction.scalefactor._data[0]
cld_pres=res.data_fields.CloudPressure._data * res.data_fields.CloudPressure.scalefactor._data[0]
cld_rfrac=res.data_fields.Cloudradiancefraction._data * res.data_fields.Cloudradiancefraction.scalefactor._data[0]

;stop
dum=where(res.data_fields.CloudFraction._data eq res.data_fields.CloudFraction._fillvalue._data[0],count) 
if count gt 0 then cld_frac[dum]=res.data_fields.CloudradianceFraction._fillvalue._data[0]
dum=where(res.data_fields.Cloudpressure._data eq res.data_fields.Cloudpressure._fillvalue._data[0],count) 
if count gt 0 then cld_pres[dum]=res.data_fields.Cloudradiancefraction._fillvalue._data[0]
dum=where(res.data_fields.CloudradianceFraction._data eq res.data_fields.CloudradianceFraction._fillvalue._data[0],count) 
if count gt 0 then cld_rfrac[dum]=res.data_fields.CloudRadianceFraction._fillvalue._data[0]

dum=where(strlowcase(tag_names(res.data_fields)) eq 'terrainheight',count)
if (count gt 0) then begin
  terrainheight=res.data_fields.TerrainHeight._data
endif else begin
  terrainheight=res.geolocation_fields.TerrainHeight._data
endelse

QAflags=0. & dum=tag_names(res.geolocation_fields)
dum=where(strlowcase(dum) eq 'groundpixelqualityflags',count)
if (count gt 0) then QAflags=res.geolocation_fields.GroundPixelQualityFlags._data
pixel=replicate(data_object,n_elements(sza[*,0]),n_elements(sza[0,*]))

if n_elements(albedo) gt 100 then pixel.albedo     = albedo
if n_elements(amf) gt 100 then pixel.amf           = amf
if n_elements(amfgeo) gt 100 then pixel.amfgeo     = amfgeo
if n_elements(amftrop) gt 100 then pixel.amftrop   = amftrop
if n_elements(amfclear) gt 100 then pixel.amfclear = amfclear
if n_elements(amfcloud) gt 100 then pixel.amfcloud = amfcloud
if n_elements(sza) gt 100 then pixel.sza           = sza
if n_elements(scd) gt 100 then pixel.scd           = scd
if n_elements(scdstrat) gt 100 then pixel.scdstrat = scdstrat
if n_elements(vcd) gt 100 then pixel.vcd           = vcd
if n_elements(rms) gt 100 then pixel.rms           = rms
if n_elements(chi) gt 100 then pixel.chi           = chi
if n_elements(fqf) gt 100 then pixel.fqf           = fqf
 if n_elements(scdprec) gt 100 then pixel.scdprec   = scdprec
if n_elements(cld_frac) gt 100 then pixel.cld_frac = cld_frac
if n_elements(cld_pres) gt 100 then pixel.cld_pres = cld_pres
if n_elements(cld_rfrac) gt 100 then pixel.cld_rfrac= cld_rfrac
if n_elements(vza) gt 100 then pixel.vza           = vza
if n_elements(lat) gt 100 then pixel.lat           = lat
if n_elements(lon) gt 100 then pixel.lon           = lon
if n_elements(qaflags) gt 100 then pixel.gpxQAflags= QAflags
if n_elements(vcdstrat) gt 100 then pixel.vcdstrat = vcdstrat
; repair idl failure to handle byte datatypes
if n_elements(vcdtrop) gt 100 then begin
  fltrop=fix(fltrop)
  dum=where(fltrop ne 0,count)
  if count gt 0 then fltrop[dum]=fltrop[dum]-256
  pixel.vcdtrop=vcdtrop
endif
if n_elements(fltrop) gt 100 then pixel.fltrop               = fltrop
if n_elements(scd_orig) gt 100 then pixel.scd_orig           = scd_orig
if n_elements(vcdtropstd) gt 100 then pixel.sigvcdtrop       = vcdtropstd
if n_elements(terrainheight) gt 100 then pixel.terrainheight = terrainheight
if n_elements(ghostcolumn) gt 100 then pixel.ghostcolumn     = ghostcolumn
if n_elements(vcdtropm) gt 100 then pixel.vcdtropm           = vcdtropm
if n_elements(psurf) gt 100 then pixel.psurf                 = psurf
if n_elements(tm4terrainheight) gt 100 then pixel.tm4terrainheight=tm4terrainheight
if n_elements(meas_time) gt 100 then begin
  nrows=n_elements(pixel[*,0]) & nims=n_elements(pixel[0,*])
  time=dblarr(nrows,nims)
  for row=0,nrows-1 do time[row,*]=meas_time
  pixel.time=time
endif
if n_elements(kernel) gt 100 then begin  
  dum=where(kernel lt -1e10,count)
  if count gt 0 then kernel[dum]=!values.f_nan
  n_layers=n_elements(kernel[0,0,*])
  for i=0L,n_elements(pixel)-1 do pixel[i].kernel=ptr_new(fltarr(n_layers))
  for i=0,n_elements(pixel[*,0])-1 do begin
    for j=0,n_elements(pixel[0,*])-1 do *pixel[i,j].kernel=kernel[i,j,*]
  endfor
  ;stop
endif

for row=0,n_elements(pixel[*,0])-1 do pixel[row,*].pixelnr=row

if (n_elements(cp_ll) gt 1) then begin
  pixel.cp_ll=cp_ll & pixel.cp_ul=cp_ul & pixel.cp_lr=cp_lr & pixel.cp_ur=cp_ur
endif else begin
  calculate_pixel_cornerpoints,pixel
endelse

tagnames=strlowcase(tag_names(pixel))
for i=0,n_tags(pixel)-1 do begin
  if (tagnames[i] eq 'kernel') then begin
    continue
  endif else begin
    dum=where(pixel.(i) le -1.e10,count)
    if count gt 0 then pixel[dum].(i)=!values.f_nan
  endelse
  
endfor

return,pixel

end ; function read_L2_no2

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

pro repair_nans,data,aux_data,admin,dimension,name

  nxtrack=n_elements(data[*,0])
  n=n_elements(data[0,*])
  dum=where(finite(data) eq 0,count)
  if count le 0 then return
    for i=0L,count-1 do begin
    ii=dum[i] mod nxtrack
    jj=dum[i]/nxtrack
    if dimension eq 1 then begin
      if (ii eq 0) then continue
      data[ii,jj]=max(aux_data[ii-1:ii,jj],/nan)
      admin[ii,jj]=admin[ii,jj]-1.
    endif else begin
      if (jj eq 0) then continue
      data[ii,jj]=max(aux_data[ii,jj-1:jj],/nan)
      admin[ii,jj]=admin[ii,jj]-2.
    endelse
  endfor
  dum=where(finite(data) eq 0,count1)
;  print,name+': '+st(count-count1)+' of '+st(count) +' NaN pixels removed'

end

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

pro calculate_pixel_cornerpoints,pixel

lat=pixel.lat
lon=pixel.lon

n=n_elements(lat[0,*])
nxtrack=n_elements(lat[*,0])

dum=where(abs(lat) gt 200.,count) & if count gt 0 then lat[dum]=!values.f_nan
dum=where(abs(lon) gt 200.,count) & if count gt 0 then lon[dum]=!values.f_nan

;***************Repair missing geolocation coordinates*********************
;
; Find the datapoints for which the geolocation coordinates are missing
; and try to fill in the gaps by linear interpolation
;
; Do this for the latitude field
;
dum=where(finite(lat) eq 0,count)
if count gt 0 then begin
  for i=0L,count-1 do begin
    jj=dum[i]/nxtrack
    ii=dum[i] mod nxtrack
    if ((ii lt 1) OR (jj lt 1) OR (ii gt nxtrack-2) OR (jj gt n-2)) then continue
    lat[ii,jj]=mean([reform(lat[ii-1:ii+1,jj]),reform(lat[ii,jj-1:jj+1])],/nan)
    if finite(lon[ii,jj] eq 0) then lon[ii,jj]=geomean([reform(lon[ii-1:ii+1,jj]),reform(lon[ii,jj-1:jj+1])])
  endfor
endif
;
; and for the longitude field
;
dum=where(finite(lon) eq 0,count)
if count gt 0 then begin
  for i=0L,count-1 do begin
    jj=dum[i]/nxtrack
    ii=dum[i] mod nxtrack
    if ((ii lt 1) OR (jj lt 1) OR (ii gt nxtrack-2) OR (jj gt n-2)) then continue
    lon[ii,jj]=geomean([reform(lon[ii-1:ii+1,jj]),reform(lon[ii,jj-1:jj+1])])
  endfor
endif
pixel.lat=lat
pixel.lon=lon

;
; *********************Finished repairing missing geolocation coordinates********************
;

; ##############Here begins the time critical part of the code##############################################
;
; What needs to be done is calculating the lat/lon value for each pixel by taking the mean of its four neighbours
; The old (and slow) method is to do this in a double for-loop and call the mean and geomean function for each row and image number.
;
; The fast approach is to do this vector wise, it makes the code a bit more complicated to read, but at the same time about
; ten times faster!
; Basically the method boils down to calculate the average for an entire column by adding two columns and dividing by 2.
; However we need to take care of NaN ourselves, something which in the old method was taken care of by the mean function.
;

; create and initialise temporary bookkeeping arrays.
;
  tlat=fltarr(nxtrack,n) & tlat[*,0]=!values.f_nan & tlat1=tlat & tlon=tlat & tlon1=tlon & deltalon=lon
  admin_lat=tlat & admin_lat[*]=4. & admin_lon=admin_lat
;
; First step of calculating the average of four points: add two neighboring columns for all columns in the swath
; Note: an array column corresponds to an OMI row.
;
  for i=1,nxtrack-1 do begin
    tlat[i,1:*]     = lat[i,1:*] + lat[i-1,1:*]
    tlon[i,1:*]     = lon[i,1:*] + lon[i-1,1:*]
    deltalon[i,1:*] = lon[i,1:*] - lon[i-1,1:*]
  endfor
  deltalon=abs(deltalon)

;
; Fix NaN occurrences
;
  repair_nans,tlat,lat,admin_lat,1,'lat1'
  repair_nans,tlon,lon,admin_lon,1,'lon1'

;
;Repair dateline issues for the longitude parameter
;
  dum=where(deltalon gt 180.,count)
  if count gt 0 then begin
    for i=0,count -1 do begin
      ii=dum[i] mod nxtrack
      jj=dum[i] / nxtrack
      tlon[ii,jj] = (tlon[ii,jj] gt 0.) ? tlon[ii,jj]-360. : tlon[ii,jj]+360.
    endfor
  endif

;
; Second step of calculating the average: add two adjacent rows for all rows (images) in the orbit
; Note: an array row corresponds to an OMI image
;
  for i=1,n-1 do begin
    tlat1[1:*,i]=tlat[1:*,i-1]+tlat[1:*,i]
    tlon1[1:*,i]=tlon[1:*,i-1]+tlon[1:*,i]
    deltalon[1:*,i]=tlon[1:*,i-1]-tlon[1:*,i]
  endfor
  deltalon=0.5*abs(deltalon)
;
;Fix NaN occurrences
;
  repair_nans,tlat1,tlat,admin_lat,2,'lat2'
  repair_nans,tlon1,tlon,admin_lon,2,'lon2'
;
; Repair dateline issues for the longitude parameter
;
  dum=where(deltalon gt 180.,count)
  if count gt 0 then begin
    for i=0,count -1 do begin
      ii=dum[i] mod nxtrack
      jj=dum[i] / nxtrack
      tlon1[ii,jj] = (tlon1[ii,jj] gt 0.) ? tlon1[ii,jj]-720. : tlon1[ii,jj]+720.
    endfor
  endif
;
; the admin_ array contains the bookkeeping parameter N by which we are going to divide.
; So prevent division by zero.
;
  dum=where(admin_lat le 0.,count) & if dum gt 0 then admin_lat[dum]=!values.f_nan
  dum=where(admin_lon le 0.,count) & if dum gt 0 then admin_lon[dum]=!values.f_nan

  tlat=tlat1/admin_lat & tlat1=0.
  tlon=tlon1/admin_lon & tlon1=0.
;
; the cornerpoints have been calculated. Assign them to the proper array elements in the pixel struct.
;
  for j=1,nxtrack-1 do begin
    for i=1,n-1 do begin
;      cp_ll=[mean(lat[j-1:j,i-1:i],/nan),geomean(lon[j-1:j,i-1:i])]
      cp_ll=[tlat[j,i],tlon[j,i]]
      pixel[j,i].cp_ll=cp_ll
      pixel[j,i-1].cp_ul=cp_ll
      pixel[j-1,i].cp_lr=cp_ll
      pixel[j-1,i-1].cp_ur=cp_ll
    endfor
  endfor
;
; #################################end of the time critical part##############################
;

; process extreme rows of the swath
  cp_el=fltarr(2) & cp_er=cp_el
  for i=1,n-1 do begin
    cp_el[*]=!values.f_nan & cp_er[*]=!values.f_nan
    cp_el[0]=2.*pixel[1,i].cp_ll[0]-pixel[1,i].cp_lr[0]
    cp_er[0]=2.*pixel[nxtrack-2,i].cp_lr[0]-pixel[nxtrack-2,i].cp_ll[0]
    cp_el[1]=longitude_extrapolate(pixel[1,i].cp_ll[1],pixel[1,i].cp_lr[1])
    cp_er[1]=longitude_extrapolate(pixel[nxtrack-2,i].cp_lr[1],pixel[nxtrack-2,i].cp_ll[1])
    pixel[0,i].cp_ll=cp_el
    pixel[0,i-1].cp_ul=cp_el
    pixel[nxtrack-1,i].cp_lr=cp_er
    pixel[nxtrack-1,i-1].cp_ur=cp_er
  endfor

;
; Process first and last image of the orbit
;
  for j=1,nxtrack-1 do begin
    pixel[j,0].cp_ll[0]=2.*pixel[j,1].cp_ll[0] - pixel[j,2].cp_ll[0]
    pixel[j,0].cp_ll[1]=longitude_extrapolate(pixel[j,1].cp_ll[1],pixel[j,2].cp_ll[1])
    pixel[j-1,0].cp_lr=pixel[j,0].cp_ll
    pixel[j,n-1].cp_ul[0]=2.*pixel[j,n-2].cp_ul[0] - pixel[j,n-3].cp_ul[0]
    pixel[j,n-1].cp_ul[1]=longitude_extrapolate(pixel[j,n-2].cp_ul[1],pixel[j,n-3].cp_ul[1])
    pixel[j-1,n-1].cp_ur=pixel[j,n-1].cp_ul
  endfor
;
; process extreme cornerpoints
;
  pixel[0,0].cp_ll[0]=2.*pixel[1,0].cp_ll[0] - pixel[2,0].cp_ll[0]
  pixel[0,n-1].cp_ul[0]=2.*pixel[1,n-1].cp_ul[0] - pixel[2,n-1].cp_ul[0]
  pixel[nxtrack-1,0].cp_lr[0] = 2.*pixel[nxtrack-2,0].cp_lr[0] - pixel[nxtrack-3,0].cp_lr[0]
  pixel[nxtrack-1,n-1].cp_ur[0] = 2.*pixel[nxtrack-2,n-1].cp_ur[0] - pixel[nxtrack-3,n-1].cp_ur[0]

  pixel[0,0].cp_ll[1]=longitude_extrapolate(pixel[1,0].cp_ll[1],pixel[2,0].cp_ll[1])
  pixel[0,n-1].cp_ul[1]=longitude_extrapolate(pixel[1,n-1].cp_ul[1],pixel[2,n-1].cp_ul[1])
  pixel[nxtrack-1,0].cp_lr[1] =longitude_extrapolate( pixel[nxtrack-2,0].cp_lr[1],pixel[nxtrack-3,0].cp_lr[1])
  pixel[nxtrack-1,n-1].cp_ur[1] = longitude_extrapolate(pixel[nxtrack-2,n-1].cp_ur[1],pixel[nxtrack-3,n-1].cp_ur[1])

;
; is there a large gap in the orbit?
;
  dlat=lat[30,*] & n=n_elements(dlat)-1 & dlat=dlat[0:n-1] & dlat[*]=!values.f_nan
  for i=1,n-1 do dlat[i]=lat[30,i]-lat[30,i-1]
  dum=where(abs(dlat) gt 2,count)
  if count gt 0 then begin
    for index=0,count-1 do begin
      image=dum[index]-1
; Process first and last image of the orbitparts
      for j=1,nxtrack-1 do begin
        pixel[j,image+1].cp_ll[0] = 2.*pixel[j,image+2].cp_ll[0] - pixel[j,image+3].cp_ll[0]
        pixel[j,image+1].cp_ll[1] = longitude_extrapolate(pixel[j,image+2].cp_ll[1],pixel[j,image+3].cp_ll[1])
        pixel[j-1,image+1].cp_lr  = pixel[j,image+1].cp_ll
        pixel[j,image].cp_ul[0]   = 2.*pixel[j,image-1].cp_ul[0] - pixel[j,image-2].cp_ul[0]
        pixel[j,image].cp_ul[1]   = longitude_extrapolate(pixel[j,image-1].cp_ul[1],pixel[j,image-2].cp_ul[1])
        pixel[j-1,image].cp_ur    = pixel[j,image].cp_ul
      endfor
; process extreme cornerpoints
      pixel[0,image+1].cp_ll[0]         = 2.*pixel[1,image+1].cp_ll[0]         - pixel[2,image+1].cp_ll[0]
      pixel[0,image].cp_ul[0]           = 2.*pixel[1,image].cp_ul[0]           - pixel[2,image].cp_ul[0]
      pixel[nxtrack-1,image+1].cp_lr[0] = 2.*pixel[nxtrack-2,image+1].cp_lr[0] - pixel[nxtrack-3,image+1].cp_lr[0]
      pixel[nxtrack-1,image].cp_ur[0]   = 2.*pixel[nxtrack-2,image].cp_ur[0]   - pixel[nxtrack-3,image].cp_ur[0]
      pixel[0,image+1].cp_ll[1]         = longitude_extrapolate(pixel[1,image+1].cp_ll[1],pixel[2,image+1].cp_ll[1])
      pixel[0,image].cp_ul[1]           = longitude_extrapolate(pixel[1,image].cp_ul[1],pixel[2,image].cp_ul[1])
      pixel[nxtrack-1,image+1].cp_lr[1] = longitude_extrapolate(pixel[nxtrack-2,image+1].cp_lr[1],pixel[nxtrack-3,image+1].cp_lr[1])
      pixel[nxtrack-1,image].cp_ur[1]   = longitude_extrapolate(pixel[nxtrack-2,image].cp_ur[1],pixel[nxtrack-3,image].cp_ur[1])
    endfor
  endif
end

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

function st,dumin
  dums=''
  if n_elements(dumin) eq 0 then return,dums
  if n_elements(dumin) eq 1 then dums=strtrim(string(dumin),2)
  if n_elements(dumin) gt 1 then begin
    for i=0,n_elements(dumin)-1 do dums=dums+strtrim(string(dumin[i]),2)+' '
  endif
return,dums
end