; Program: convert_hdfeos_hdf
; Author: Ruud Dirksen
; Date:   24 April 2008
; Goal: convert DOMINO hdf-eos5 files back to old hdf4 format
;
; To start the program type at the idl prompt:
;
; .com convert_hdfeos_hdf
; convert_hdfeos_hdf
;

function user_inputs

;###########User input section########################
;
; The program looks for the input files in the path below,
; the output files will be stored in this directory as well.

path='/nobackup/users/boersma/omi/coll3/2004/11/'

;##########End of user input section#######################

  return,path
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

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

pro sort_attributes,attr

dum=where(strpos(strlowcase(attr.name),'input_pointer') ge 0,count )
if count lt 1 then return
first_inputpointer=dum[0]

; find the Input_pointer entries that contain a L2 NO2 file, sort these and store them 
; directly below the generic attribute list
i1=0 & index=intarr(count) & index[*]=-1 
for ii=0,count-1 do begin
  i=dum[ii]
  if (ptr_valid(attr[i].data) eq 0) then continue
  dumstring=*attr[i].data
  if (strpos(strlowcase(dumstring),'omno2') ge 0) then begin & index[i1]=i & i1=i1+1 & endif
endfor
dum=where(index ge 0,count) 
index = (count gt 0) ? index[dum] : -1
Nno2=0 
if (index[0] ge 0) then begin
  Nno2=n_elements(index)
  for i=0,n_elements(index)-1 do begin
    si=index[i] & ti=first_inputpointer+i
    swap=attr[si] & attr[si]=attr[ti] & attr[ti]=swap
  endfor
endif

dum=where(strpos(strlowcase(attr.name),'input_pointer') ge 0,count )
i1=0 & index=intarr(count) & index[*]=-1 
for ii=0,count-1 do begin
  i=dum[ii]
  if (ptr_valid(attr[i].data) eq 0) then continue
  dumstring=*attr[i].data
  if (strpos(strlowcase(dumstring),'amf_lut') ge 0) then begin & index[i1]=i & i1=i1+1 & endif
endfor
dum=where(index ge 0,count) 
index = (count gt 0) ? index[dum] : -1
Namf=0 
if (index[0] ge 0) then begin
  Namf=n_elements(index)
  for i=0,n_elements(index)-1 do begin
    si=index[i] & ti=first_inputpointer+i+Nno2
    swap=attr[si] & attr[si]=attr[ti] & attr[ti]=swap
  endfor
endif

dum=where(strpos(strlowcase(attr.name),'input_pointer') ge 0,count )
for i=0,count-1 do attr[dum[i]].name='Input_pointer_'+st(i+1)

end

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

function fill_datacontainer

  datacontainer_struct={hdf_name:'',hdfeos_name:'',target:-1,origin:-1,data:ptrarr(20),fixedvalue:0.,$
                        attributes:ptrarr(20),global_attributes:ptr_new(),use_geo:-1,datatype:-1}
  datacontainer=replicate(datacontainer_struct,100)
  datacontainer[*].fixedvalue=!values.f_nan

  i=0
  datacontainer[i].hdf_name='date'
  datacontainer[i].hdfeos_name='date'
  datacontainer[i].target=0
  datacontainer[i].origin=-99

  i=i+1
  datacontainer[i].hdf_name='time'
  datacontainer[i].hdfeos_name='Time'
  datacontainer[i].target=0
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='lon'
  datacontainer[i].hdfeos_name='Longitude'
  datacontainer[i].target=0
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='lat'
  datacontainer[i].hdfeos_name='Latitude'
  datacontainer[i].target=0
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='vcd'
  datacontainer[i].hdfeos_name='TotalVerticalColumn'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='sigvcd'
  datacontainer[i].hdfeos_name='TotalVerticalColumnError'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='vcdtrop'
  datacontainer[i].hdfeos_name='TroposphericVerticalColumn'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='sigvcdt'
  datacontainer[i].hdfeos_name='TroposphericVerticalColumnError'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='vcdstrat'
  datacontainer[i].hdfeos_name='AssimilatedStratosphericVerticalColumn'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='sigvcds'
  datacontainer[i].hdfeos_name='AssimilatedStratosphericVerticalColumnError'
  datacontainer[i].target=0
  datacontainer[i].origin=-1
  datacontainer[i].fixedvalue=0.25

  i=i+1
  datacontainer[i].hdf_name='fltrop'
  datacontainer[i].hdfeos_name='TroposphericColumnFlag'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='psurf'
  datacontainer[i].hdfeos_name='TM4SurfacePressure'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='sigvcdak'
  datacontainer[i].hdfeos_name='VCDErrorUsingAvKernel'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='sigvcdtak'
  datacontainer[i].hdfeos_name='VCDTropErrorUsingAvKernel'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='kernel'
  datacontainer[i].hdfeos_name='AveragingKernel'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='ghostcol'
  datacontainer[i].hdfeos_name='GhostColumn'
  datacontainer[i].target=0
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='sza'
  datacontainer[i].hdfeos_name='SolarZenithAngle'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='vza'
  datacontainer[i].hdfeos_name='ViewingZenithAngle'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='saa'
  datacontainer[i].hdfeos_name='SolarAzimuthAngle'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='vaa'
  datacontainer[i].hdfeos_name='ViewingAzimuthAngle'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='ssc'
  datacontainer[i].hdfeos_name='InstrumentConfigurationID'
  datacontainer[i].target=1
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='loncorn'
  datacontainer[i].hdfeos_name='LongitudeCornerPoints'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='latcorn'
  datacontainer[i].hdfeos_name='LatitudeCornerPoints'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='pixelnr'
  datacontainer[i].hdfeos_name='pixelnr'
  datacontainer[i].target=1
  datacontainer[i].origin=-1

  i=i+1
  datacontainer[i].hdf_name='imagenr'
  datacontainer[i].hdfeos_name='imagenr'
  datacontainer[i].target=1
  datacontainer[i].origin=-1

  i=i+1
  datacontainer[i].hdf_name='gpxflag'
  datacontainer[i].hdfeos_name='GroundPixelQualityFlags'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='theight'
  datacontainer[i].hdfeos_name='TerrainHeight'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='mheight'
  datacontainer[i].hdfeos_name='TM4TerrainHeight'
  datacontainer[i].target=1
  datacontainer[i].origin=1

  i=i+1
  datacontainer[i].hdf_name='scd'
  datacontainer[i].hdfeos_name='SlantColumnAmountNO2'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='sigscd'
  datacontainer[i].hdfeos_name='SlantColumnAmountNO2Std'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='amf'
  datacontainer[i].hdfeos_name='AirmassFactor'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='amftrop'
  datacontainer[i].hdfeos_name='AirmassFactorTropospheric'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='amfgeo'
  datacontainer[i].hdfeos_name='AirmassFactorGeometric'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='scdstr'
  datacontainer[i].hdfeos_name='AssimilatedStratosphericSlantcolumn'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='clfrac'
  datacontainer[i].hdfeos_name='CloudFraction'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='cltpres'
  datacontainer[i].hdfeos_name='CloudPressure'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='albclr'
  datacontainer[i].hdfeos_name='SurfaceAlbedo'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='crfrac'
  datacontainer[i].hdfeos_name='CloudRadianceFraction'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='ltropo'
  datacontainer[i].hdfeos_name='TM4TropopauseLevel'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  i=i+1
  datacontainer[i].hdf_name='fcvcdtr'
  datacontainer[i].hdfeos_name='TroposphericVerticalColumnModel'
  datacontainer[i].target=2
  datacontainer[i].origin=0

  dum=where(strlen(datacontainer.hdfeos_name) gt 0,count)
  if count le 0 then return,0
  datacontainer=datacontainer[dum]

  return,datacontainer
end

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

pro read_global_attributes,data,attr

;attributes_struct={name:'',data:ptr_new(),count:0,type:'dfnt_char88'}
;attr=replicate(attributes_struct,150)
months=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']

i=0 & eol=0 & ip=1
if (ptr_valid(attr[0].data) eq 0) then begin
; if this is the first call the list of attributes still is completely empty
; Fill the top of the list with general information
  attr[i].name='Authors'            & attr[i].data=ptr_new(data.authors._data)   & i=i+1
  attr[i].name='Affiliation'        & attr[i].data=ptr_new(data.affiliation._data) & i=i+1
  attr[i].name='Email'              & attr[i].data=ptr_new(data.email._data)     & i=i+1
  text_string=data.pge_name._data+': version '+data.pge_version._data
  attr[i].name='Data_created_by'    & attr[i].data=ptr_new(text_string)          & i=i+1
  dum_string=strmid(data._file,19,9,/reverse)
  text_string=strmid(dum_string,7,2)+' '+months[fix(strmid(dum_string,5,2))-1]+' '+strmid(dum_string,0,4)
  attr[i].name='Creation_date'      & attr[i].data=ptr_new(text_string)          & i=i+1
  attr[i].name='Unit_of_NO2_column' & attr[i].data=ptr_new('1e15 molecules/cm2') & i=i+1
  eol=i
endif else begin
; this is not the first call. So the attributes list is not empty
; count the number of entries that carry the name Input_pointer
  for i=0,n_elements(attr)-1 do if ptr_valid(attr[i].data) eq 0 then break
  eol=i
  dum=where(strmid(attr.name,0,13) eq 'Input_pointer',count)
  ip=count+1
endelse

; now we know the number for the next (to be defined) inputpointer on the list
attr[eol].name='Input_pointer_'+st(ip) & attr[eol].data=ptr_new(data.no2_l2_file._data ) & eol=eol+1 & ip=ip+1

; prevent duplication of attribute list entries. 
; If there already exists an identical entry for amf_lut then don't add a new one
i=0 & found=0
while ((i lt eol) AND (found eq 0)) do begin
  dumstring=*attr[i].data
  if (dumstring eq data.amf_lut._data) then found=1
  i=i+1
endwhile
if (found eq 0) then begin
  attr[eol].name='Input_pointer_'+st(ip) & attr[eol].data=ptr_new(data.amf_lut._data) & eol=eol+1 & ip=ip+1
endif

; prevent duplicate entries for attributes dealing with meteodata 
; If there already exists an identical entry for amf_lut then don't add a new one
dum=where(strmid(tag_names(data),0,5) eq 'METEO',count)
for index=0,count-1 do begin
  i=0 & found=0
  while ((i lt eol) AND (found eq 0)) do begin
    dumstring=*attr[i].data
    if (dumstring eq data.(dum[index])._data) then found=1
    i=i+1
  endwhile
  if found eq 0 then begin
    attr[eol].name='Input_pointer_'+st(ip)
    attr[eol].data=ptr_new(data.(dum[index])._data)
    eol=eol+1 & ip=ip+1
  endif
endfor

end ; read_global_attributes

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

pro write_global_attributes,fn,attributes

  if (file_test(fn) eq 1) then begin
    fid=hdf_open(fn,/rdwr)
  endif else begin
    print,'The output file doesnt exist'
    stop
  endelse

  sd_id=hdf_sd_start(fn,/rdwr)

  for i=0,n_elements(attributes)-1 do begin
    name=attributes[i].name
    count=attributes[i].count
    data=*attributes[i].data
    hdf_sd_attrset,sd_id,name,data,count
    ptr_free,attributes[i].data
  endfor

  hdf_sd_end,sd_id
  hdf_close,fid

end ; write_global_attributes

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

pro write_swath,fn,data,datestring,no2=no2,anc=anc,geo=geo

type_lut=intarr(20) & type_lut[*]=-1
type_lut[1]= 21   ;byte
type_lut[2]= 22   ;integer
type_lut[3]= 24   ;long
type_lut[4]= 5    ;float
type_lut[5]= 6    ;double
type_lut[12]= 23  ;
type_lut[13]= 25  ;

;
; collect fieldnames that are present in the HDFEOS file
;
datafieldnames  =strlowcase(tag_names(data.data_fields))
geolocationnames=strlowcase(tag_names(data.geolocation_fields))
geolocationnames=[geolocationnames,'date']
;
; make index array for valid datapoints.
;
ii=where(datafieldnames eq 'totalverticalcolumn')
dumdata=data.data_fields.(ii)._data
data_index_array=where(dumdata gt -1.e10,tot_count)
if tot_count le 0 then begin
  print,'no valid data present'
  return
endif

imagenr=fix(dumdata) & imagenr[*]=0 & pixelnr=imagenr
n1=n_elements(imagenr[*,0])
n2=n_elements(imagenr[0,*])
for i=0,n1-1 do begin
  imagenr[i,*]=indgen(n2)+1
  pixelnr[i,*]=replicate(i+1,n2)
endfor

;
; special treatment for the date and time field. In HDFEOS the time array is one dimensional,
; one element per image.
; In the HDF file the date and time is provided for each datapoint/pixel
;
; timedata_orig is 1644x60 of float containing TAI time for each pixel in the swath
;
t=data.geolocation_fields.time._data
timedata_orig=double(dumdata) & timedata_orig[*]=0.
for i=0,n_elements(timedata_orig[*,0])-1 do timedata_orig[i,*]=t

;
; time_mask is 1D integer array. It contains the image numbers for which a valid measurement
; is present.
;
time_mask=fix(dumdata) & time_mask[*]=0 & time_mask[data_index_array]=1
for i=0,n_elements(time_mask[0,*])-1 do time_mask[0,i]=fix(total(time_mask[*,i]))
time_mask=where(reform(time_mask[0,*]) ge 1)

tmp_fn='c:\convert_tmp.txt'
if (strlowcase(!version.os_family) eq 'unix') then tmp_fn='~/convert_tmp.txt'
t=data.geolocation_fields.time._data[time_mask]
dt=fltarr(n_elements(t)-1)
for i=0L,n_elements(t)-2 do dt[i]=t[i+1]-t[i]
dum=where(dt gt 10.*60.*2.,count)
i0=0L
if count gt 0 then i0=dum[0]
tai2cal,t[i0],year, month, day, hour, minute, seconds, STATUS=status
tai2cal,t[0],y0,m0,d0,h0,min0,s0
tai2cal,t[n_elements(t)-1],y1,m1,d1,h1,min1,s1
start_time=[y0,m0,d0,h0,min0,fix(s0)]
end_time=[y1,m1,d1,h1,min1,fix(s1)]

openw,1,tmp_fn & printf,1,format='(4i02)',month,day,hour,minute & close,1
dumregel=''
openr,1,tmp_fn & readf,1,dumregel & close,1
file_delete,tmp_fn,/allow_nonexistent

datestring=strmid(st(year),2,2)+dumregel

if (file_test(fn) eq 1) then begin
  fid=hdf_open(fn,/rdwr)
endif else begin
  print,'The output file doesnt exist'
  stop
endelse

if fid le 0 then begin
  print,'Error openening file ',fn
  return
endif

if keyword_set(no2) then begin
  target=0 & swathname='NO2_'+datestring & class='FRESCO_NO2_SCD'
endif
if keyword_set(geo) then begin
  target=1 & swathname='GEO_'+datestring & class='GEOLOCATION DATA'
endif
if keyword_set(anc) then begin
  target=2 & swathname='ANC_'+datestring & class='ANCILLARY DATA'
endif

datafieldslocation=fill_datacontainer()
index_array=where(datafieldslocation.target eq target,count)
if count le 0 then begin
  print,'No valid datatypes'
  return
endif
nfields=count

swid=hdf_vd_attach(fid,-1)
hdf_vd_setinfo,swid,name=swathname,/full,class=class

n_bytes=0L
fieldnames=''
vector_length=0
for i=0,nfields-1 do begin
    dumdata=0
    hdf_name=strlowcase(datafieldslocation[index_array[i]].hdf_name)
    hdfeos_name=strlowcase(datafieldslocation[index_array[i]].hdfeos_name)
    pos=where(datafieldnames eq hdfeos_name,count)

    if (datafieldslocation[index_array[i]].origin eq -1) then count=1

    use_geo=0
    if count le 0 then begin
      pos=where(geolocationnames eq hdfeos_name,count)
      use_geo=1
    endif
    if count le 0 then continue ; this fieldname cannot be found
    datafieldslocation[index_array[i]].use_geo=use_geo
    if (use_geo eq 0) then begin
      if pos ge 0 then begin
        dumdata=data.data_fields.(pos)._data
        dum=data.data_fields.(pos).scalefactor._data &  scalefactor=dum[0]
        if scalefactor lt .1 then begin
          dumdata=float(dumdata)*scalefactor
;          print,hdfeos_name,' ',st(scalefactor)
        endif
; exception for the InstrumentConfigurationID field
        if strlowcase(hdf_name) EQ 'ssc' then begin
          tmpdata=dumdata & dumdata=pixelnr & dumdata[*]=-1
          for ii=0,n_elements(dumdata[*,0])-1 do dumdata[ii,*]=tmpdata
          tmpdata=0
        endif
      endif else begin
; exception for the sigvcds field
        if finite(datafieldslocation[index_array[i]].fixedvalue) then dumdata=replicate(datafieldslocation[index_array[i]].fixedvalue,tot_count)
; exception for the pixelnr field
        if strlowcase(hdf_name) EQ 'pixelnr' then dumdata=pixelnr
; exception for the imagenr field
        if strlowcase(hdf_name) EQ 'imagenr' then dumdata=imagenr
      endelse
    endif else begin
      if ((strlowcase(hdf_name) eq 'date') OR (strlowcase(hdf_name) eq 'time')) then begin
        dumdata=timedata_orig
      endif else begin
        if pos ge 0 then begin
          dumdata=data.geolocation_fields.(pos)._data
        endif else begin
          if finite(datafieldslocation[index_array[i]].fixedvalue) then dumdata=replicate(datafieldslocation[index_array[i]].fixedvalue,tot_count)
        endelse
      endelse
    endelse

    if ((strlowcase(hdf_name) eq 'cltpres') OR (strlowcase(hdf_name) eq 'psurf')) then dumdata=dumdata*100.
    if strlowcase(hdf_name) eq 'fltrop' then begin
      dumdata=fix(dumdata)
      dum=where(dumdata eq 255,count)
      if count gt 0 then dumdata[dum]=-1
    endif

    fieldname=hdf_name & hdf_order=1
    if ((strlowcase(fieldname) eq 'time') OR (strlowcase(fieldname) eq 'date') ) then begin
      tai2cal,dumdata,year, month, day, hour, minute, seconds, STATUS=status
      timedata=hour*10000L + minute*100L + long(seconds)
      datedata=year*10000L + month*100L + day
      dumdata=datedata
      if (strlowcase(fieldname) eq 'time') then dumdata=timedata
    endif

    if n_elements(dumdata) le 5 then continue

    datafieldslocation[index_array[i]].datatype=size(dumdata,/type)
    dum=size(dumdata) 
    if (dum[0] gt 2) then begin
      hdf_order=dum[3]
      tmp_array=intarr(dum[3],n_elements(data_index_array))
      case size(dumdata,/type) of
         1:tmp_array=byte(tmp_array)
         3:tmp_array=long(tmp_array)
         4:tmp_array=float(tmp_array)
         5:tmp_array=double(tmp_array)
        12:tmp_array=uint(tmp_array)
        13:tmp_array=ulong(tmp_array)
        14:tmp_array=long64(tmp_array)
        15:tmp_array=ulong64(tmp_array)
      endcase
      for ii=0L,n_elements(data_index_array)-1 do tmp_array[*,ii]=dumdata[pixelnr[data_index_array[ii]]-1,imagenr[data_index_array[ii]]-1,*]
      dumdata=temporary(tmp_array)
    endif else begin
      dumdata=dumdata[data_index_array]
    endelse
;    vector_length = vector_length > n_elements(dumdata)
    vector_length = vector_length > n_elements(data_index_array)
    ptr_free,datafieldslocation[index_array[i]].data[0]
    datafieldslocation[index_array[i]].data[0]=ptr_new(dumdata)

    case datafieldslocation[index_array[i]].datatype of
      1:begin  & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/byte   & n_bytes=n_bytes+1*hdf_order & end
      2:begin  & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/int    & n_bytes=n_bytes+2*hdf_order & end
      3:begin  & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/long   & n_bytes=n_bytes+4*hdf_order & end
      4:begin  & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/float  & n_bytes=n_bytes+4*hdf_order & end
      5:begin  & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/double & n_bytes=n_bytes+4*hdf_order & end
      12:begin & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/uint   & n_bytes=n_bytes+2*hdf_order & end
      13:begin & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/ulong  & n_bytes=n_bytes+4*hdf_order & end
      14:begin & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/dlong  & n_bytes=n_bytes+4*hdf_order & end
      15:begin & hdf_vd_fdefine,swid,fieldname,order=hdf_order,/dulong & n_bytes=n_bytes+4*hdf_order & end
    endcase
  endfor ; nfields

  dumdata=0
  total_data=bytarr(n_bytes,vector_length)
  n_bytes=0L

  for i=0,nfields-1 do begin
    tmp_data=0
    fieldname=datafieldslocation[index_array[i]].hdf_name
    if (ptr_valid(datafieldslocation[index_array[i]].data[0]) ne 1) then continue
    dumdata=*datafieldslocation[index_array[i]].data[0]
    ptr_free,datafieldslocation[index_array[i]].data[0]
    dum=size(dumdata) & hdf_order=1
    if (dum[0] ge 2) then hdf_order=dum[1]
    fieldnames=fieldnames+','+fieldname
    if (i eq 0) then fieldnames=fieldname
    hdf_type=type_lut[datafieldslocation[index_array[i]].datatype]
    tmp_data=hdf_packdata(dumdata,hdf_order=hdf_order,hdf_type=hdf_type)
    dum=size(tmp_data)
    total_data[n_bytes:n_bytes+dum[1]-1,*]=tmp_data
    n_bytes=n_bytes+dum[1]
  endfor ; nfields
  dumdata=0 & tmp_data=0

  fieldnames=strtrim(fieldnames,2)
  hdf_vd_write,swid,fieldnames,total_data,/full

  if keyword_set(no2) then begin
    name='orbit_identifier'     & hdf_vd_attrset,swid,-1,name,datestring,/string
    pos=strpos(data._file,'-o') & ddata=long(strmid(data._file,pos+2,5)) & name='orbit_number' & hdf_vd_attrset,swid,-1,name,ddata,/dfnt_int32
    name='start_time'           & hdf_vd_attrset,swid,-1,name,start_time,/dfnt_int32
    name='end_time'             & hdf_vd_attrset,swid,-1,name,end_time,/dfnt_int32
  endif

  hdf_vd_detach,swid
  hdf_close,fid

end ; write_swath

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

pro write_pressure_grid,fn,data

  fid=hdf_open(fn,/create)

  a=data.data_fields.tm4pressurelevela._data
  b=data.data_fields.tm4pressurelevelb._data
  total_data=fltarr(2,n_elements(a))
  total_data[0,*]=a & total_data[1,*]=b

  swathname='pressure_grid'
  swid=hdf_vd_attach(fid,-1,/write)
  hdf_vd_setinfo,swid,name=swathname,/full;,class=  ; nog doen!!

  name='Unit_of_pressure' & text_string='Pascal' & hdf_vd_attrset,swid,-1,name,text_string
  name='Pressure_level_formula' & text_string='a_lev + p_surf*b_lev' & hdf_vd_attrset,swid,-1,name,text_string
  name='Number_of_pressure_levels' & hdf_vd_attrset,swid,-1,name,n_elements(a)

  fieldname='a_lev' & hdf_vd_fdefine,swid,fieldname,order=1,/float
  fieldname='b_lev' & hdf_vd_fdefine,swid,fieldname,order=1,/float
  fieldnames='a_lev,b_lev'
  hdf_vd_write,swid,fieldnames,total_data,/full

  hdf_vd_detach,swid
  hdf_close,fid

end

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

function tai2julday, tai, STATUS=status

  status = -1 
  if ( size( tai, /TYPE ) ne 5 ) then return, -1
  if ( min(tai) le 189302405.0D ) then return, -1

  leapsec = 5.0
  jul_tai0 = JULDAY(1, 1, 1993, 0, 0, 0)
  jul = jul_tai0 + (tai- leapsec) / ( 86400D )    ; convert tai to dayfrac and add offset

  status = 0
  return, jul
end

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

pro tai2cal, tai, year, month, day, hour, minute, seconds, STATUS=status

  status = -1 
  if ( size( tai, /TYPE ) ne 5 ) then return
  if ( min(tai) le 189302405.0D ) then return   ; error if date is before 1/1/1999
  jul = tai2julday( tai, STATUS=status )

  if (status ne 0) then return
  CALDAT, jul, month , day , year , hour , minute , seconds

  status = 0
  return
end

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

pro convert_hdfeos_hdf

;-------------------------------------------------------------------------------
;
; Convert DOMINO hdf-eos5 files back to old hdf4 format.
; 
;     Ruud Dirksen, KNMI, 24 April 2008.
;-------------------------------------------------------------------------------


!except=0
DOMINOSWATHNAME='DominoNO2'

path=user_inputs()

time_struct={year:0,month:0,day:0,hour:0,minute:0,julday:0D}
attributes_struct={name:'',data:ptr_new(),count:0,type:'dfnt_char88'}
attributes=replicate(attributes_struct,150)

if (hdf_exists() eq 0 ) then begin
  print,'The HDF library is not installed on your system'
  stop
endif

version=float(!version.release)
if (version lt 6.1) then begin
  print,'You need idl version 6.1 or higher to run this program'
  stop
endif

files=dialog_pickfile(filter="*.he5",path=path,/multiple_files)
n_files=n_elements(files)
if n_files lt 1 then begin
  print,'No files selected'
  stop
endif
files=files[sort(files)]

dayindexdata=ptrarr(n_elements(files))
dumstring=strmid(files,/reverse,46,14)
measurement_dates=replicate(time_struct,n_elements(files))
measurement_dates.year=fix(strmid(dumstring,0,4))
measurement_dates.month=fix(strmid(dumstring,5,2))
measurement_dates.day=fix(strmid(dumstring,7,2))
measurement_dates.hour=fix(strmid(dumstring,10,2))
measurement_dates.minute=fix(strmid(dumstring,12,2))
m=measurement_dates
measurement_dates.julday=julday(m.month,m.day,m.year,m.hour,m.minute,0.)
m=measurement_dates

ii=0 & daynumber=0
while (ii lt n_elements(measurement_dates)) do begin
  m0=measurement_dates[ii]
  jd0=julday(m0.month,m0.day,m0.year,0,0,0)
  dum=where(((m.julday - jd0) lt 1. ) AND ((m.julday - jd0) gt 0.),count)
  if count le 0 then begin
    dayindexdata[daynumber]=ptr_new(ii)
    ii=ii+1
  endif else begin
    ii=ii+count
    dayindexdata[daynumber]=ptr_new(dum)
  endelse
  daynumber=daynumber+1
endwhile
dum=where(ptr_valid(dayindexdata) eq 1,count)
if count le 0 then begin
  print,'No valid files/days'
  stop
endif
dayindexdata=dayindexdata[dum]

for dayindex=0L,n_elements(dayindexdata)-1 do begin
  output_fn=path+'no2track_'+st(dayindex)+'.hdf'
  file_indices=*dayindexdata[dayindex]
  n_files=n_elements(file_indices)
  first_time=0
  attributes.name=''
  added_files=0
  for file_index=0L,n_files-1 do begin
    fn=files[file_indices[file_index]]
    if (h5f_is_hdf5(fn) ne 1) then begin
      print,'file '+fn+' is not a valid hdf-eos5 file, skipping to next file'
      continue
    endif

    print,'Processing '+fn
    res=h5_parse(fn)
    dum=where(strlowcase(tag_names(res.hdfeos.swaths)) eq strlowcase(DOMINOSWATHNAME),count)
    if count le 0 then begin
      print,'file '+fn+' doesnt contain tropospheric NO2 swath, skipping to next file'
      continue
    endif
    res=h5_parse(fn,/read_data)
    data=res.hdfeos.swaths.(dum[0])
    if first_time eq 0 then write_pressure_grid,output_fn,data
    write_swath,output_fn,data,datestring,/no2
    if first_time eq 0 then begin
      pad=strmid(output_fn,0,strpos(output_fn,'/',/reverse_search)+1)
      new_output_fn=pad+'no2track20'+strmid(datestring,0,6)+'.hdf'
      file_move,output_fn,new_output_fn,/overwrite
      output_fn=new_output_fn
    endif
    first_time=1
    write_swath,output_fn,data,datestring,/geo
    write_swath,output_fn,data,datestring,/anc
    read_global_attributes,data,attributes
    added_files=added_files+1
  endfor ; file_index

  dum=where(ptr_valid(attributes.data) eq 1)
  global_attributes=attributes[dum]
  sort_attributes,global_attributes
  write_global_attributes,output_fn,global_attributes
  ptr_free,attributes.data & ptr_free,global_attributes.data
  global_attributes=0
  print,st(added_files)+' orbits written to '+output_fn
  print,''
endfor ; dayindex
!except=1
print,'done!'
end ; convert_hdfeos_hdf (main program)