module read_hdfeos_module

use constants_Module
use hdfeos5_Module
use datastruct_module
use datetime_module

contains

subroutine read_OMI_tropNO2(filename,datastruct,nxtrack,ntimes,npressurelevels,verbose,retstatus)
  implicit none

  real,dimension(:,:),allocatable           :: SolarZenithAngle,SolarAzimuthAngle,longitude,latitude,SurfaceAlbedo
  real,dimension(:,:),allocatable           :: ViewingZenithAngle,ViewingAzimuthAngle,TM4SurfacePressure
  real,dimension(:),allocatable             :: SpacecraftAltitude,SpacecraftLatitude,SpacecraftLongitude
  real,dimension(:,:),allocatable           :: GhostColumn,AirMassFactor,AirMassFactorTropospheric,AirMassFactorGeometric
  real,dimension(:,:),allocatable           :: VCDTropErrorUsingAvKernel,VCDErrorUsingAvKernel,CloudRadianceFraction
  real,dimension(:,:),allocatable           :: OriginalSlantColumnAmountNO2,TroposphericVerticalColumnError
  real,dimension(:,:),allocatable           :: TotalVerticalColumn,TotalVerticalColumnError,TroposphericVerticalColumn
  real,dimension(:,:),allocatable           :: AssStratSlantColumn,AssStratVerticalColumn
  real,dimension(:,:),allocatable           :: AssStratVerticalColumnError,cloudfraction,cloudpressure,cloudfractionstd,cloudpressurestd
  real,dimension(:,:),allocatable           :: SlantColumnAmountNO2,SlantColumnAmountNO2Std,TroposphericVerticalColumnModel
  real,dimension(:,:,:),allocatable         :: LatitudeCornerpoints,LongitudeCornerpoints,averagingkernel
  real(flt64),   dimension(:),  allocatable :: TAITime
  integer(int32),dimension(:),allocatable   :: htime,hdate
  integer(int16),dimension(:,:),allocatable :: GroundPixelQualityFlags,int16_buffer,TerrainHeight,TM4TerrainHeight
  integer(int16),dimension(:,:,:),allocatable :: int16_buffer3d
  integer(int32),dimension(:,:),allocatable :: int32_buffer
  integer(int_8), dimension(:),allocatable  :: MeasurementQualityFlags,InstrumentConfigurationId
  integer(int_8),dimension(:,:),allocatable :: TM4TropoPauseLevel,TroposphericColumnFlag,int8_buffer
  integer(int32)                            :: itau_offset_TM_OMI
  integer(int32)                            :: numbertype

  integer,intent(out)                                   :: retstatus
  integer                                               :: fid,status,swid
  integer                                               :: rank,i
  integer,intent(in)                                    :: nxtrack,ntimes,npressurelevels
  integer(hdf_kind_type)                                :: start(3),stride(3),edge(3)
  real                                                  :: scalefactor
  logical                                               :: fexist
  logical,intent(in)                                    :: verbose
  character(len=STD_STRING)                             :: dimlist, maxdimlist
  character(len=STD_STRING),intent(in)                  :: filename
  character(len=SWATHNAME_LENGTH)                       :: swathname,fieldname
  character(len=5000)                                   :: dfieldlist,gfieldlist
  integer                                               :: ndflds,ngflds,pos,dumtime(6)
  integer,dimension(100)                                :: ntype,nrank
  type(Tdatastruct),dimension(nxtrack,ntimes)           :: datastruct

  retstatus = 1 ; fid = -1 ; swid = -1
  write(swathname,'(a)') '                                                          '
  write(swathname,'(a)') 'DominoNO2'
  inquire (file=filename,exist=fexist)
  
  if (fexist == .false.) then
    write(*,'(3a)') 'read_OMI_tropNO2: ERROR file ',trim(filename),' does not exist'
    return
  end if

  if ( (nxtrack > MAX_NXTRACK) .OR. (nxtrack <= 0) ) then
    write(*,'(a,i)') 'read_OMI_tropNO2: ERROR invalid value of nxtrack ',nxtrack
    return
  end if

  if ( (ntimes > MAX_NTIMES) .OR. (ntimes <= 0)) then
    write(*,'(a,i)') 'read_OMI_tropNO2: ERROR invalid value of ntimes ',ntimes
    return
  end if

  if ( (npressurelevels > MAX_presslevels) .OR. (npressurelevels <= 0)) then
    write(*,'(a,i)') 'read_OMI_tropNO2: ERROR invalid value of npressurelevels ',npressurelevels
    return
  end if

  ! open the file with read-only access
  fid = he5_swopen( filename, HE5F_ACC_RDONLY )         
  if (fid <= 0 ) then
    write(*,'(2a)') 'read_OMI_tropNO2: ERROR cannot open file ',trim(filename)
    return
  end if
  if (verbose) write(*,'(3a)') 'File ',trim(filename),' opened'
  
  ! attach to the swath
  swid = he5_swattach( fid, swathname )  
  if (swid <= 0 ) then
    write(*,'(2a)') 'read_OMI_tropNO2: ERROR cannot attach swath ',trim(swathname)
    status = he5_swclose(fid)
    return
  end if
  if (verbose) write(*,'(2a)') 'Attached swath ',trim(swathname)

  ndflds = he5_swinqdflds(swid,dfieldlist,nrank,ntype)
  ngflds = he5_swinqgflds(swid,gfieldlist,nrank,ntype)
  if (verbose) write(*,'("nr of data fields: ",i3," nr of geo fields: ",i3)') ndflds,ngflds
  
  pos = index(trim(gfieldlist),'Latitude')
  if (pos < 1) then
    write(*,'(a)') 'The file is invalid, it doesnt contain the Latitude field'
    return
  end if

  allocate(latitude(nxtrack,ntimes))                         ;    allocate(TotalVerticalColumn(nxtrack,ntimes))
  allocate(longitude(nxtrack,ntimes))                        ;    allocate(TotalVerticalColumnError(nxtrack,ntimes))
  allocate(SolarZenithAngle(nxtrack,ntimes))                 ;    allocate(TroposphericVerticalColumn(nxtrack,ntimes))
  allocate(SolarAzimuthAngle(nxtrack,ntimes))                ;    allocate(AssStratSlantColumn(nxtrack,ntimes))
  allocate(ViewingZenithAngle(nxtrack,ntimes))               ;    allocate(AssStratVerticalColumn(nxtrack,ntimes))
  allocate(ViewingAzimuthAngle(nxtrack,ntimes))              ;    allocate(AssStratVerticalColumnError(nxtrack,ntimes))
  allocate(cloudfraction(nxtrack,ntimes))                    ;    allocate(SlantColumnAmountNO2(nxtrack,ntimes))
  allocate(cloudpressure(nxtrack,ntimes))                    ;    allocate(SlantColumnAmountNO2Std(nxtrack,ntimes))
  allocate(terrainheight(nxtrack,ntimes))                    ;    allocate(TroposphericVerticalColumnModel(nxtrack,ntimes))
  allocate(TM4SurfacePressure(nxtrack,ntimes))               ;    allocate(LatitudeCornerpoints(nxtrack,ntimes,NCORNERPOINTS))
  allocate(TAItime(ntimes))                                  ;    allocate(LongitudeCornerpoints(nxtrack,ntimes,NCORNERPOINTS))
  allocate(htime(ntimes))                                    ;    allocate(TM4TropoPauseLevel(nxtrack,ntimes)) 
  allocate(hdate(ntimes))                                    ;    allocate(TM4TerrainHeight(nxtrack,ntimes))   
  allocate(SpacecraftAltitude(ntimes))                       ;    allocate(TroposphericColumnFlag(nxtrack,ntimes))
  allocate(SpacecraftLatitude(ntimes))                       ;    allocate(InstrumentConfigurationId(ntimes))
  allocate(SpacecraftLongitude(ntimes))                      ;    allocate(MeasurementQualityFlags(ntimes))
  allocate(GhostColumn(nxtrack,ntimes))                      ;    allocate(GroundPixelQualityFlags(nxtrack,ntimes))
  allocate(AirMassFactor(nxtrack,ntimes))                    ;    allocate(CloudPressureStd(nxtrack,ntimes))
  allocate(AirMassFactorTropospheric(nxtrack,ntimes))        ;    allocate(CloudFractionStd(nxtrack,ntimes))
  allocate(AirMassFactorGeometric(nxtrack,ntimes))           ;    allocate(int16_buffer(nxtrack,ntimes))
  allocate(SurfaceAlbedo(nxtrack,ntimes))                    ;    allocate(int32_buffer(nxtrack,ntimes))
  allocate(VCDTropErrorUsingAvKernel(nxtrack,ntimes))        ;    allocate(int8_buffer(nxtrack,ntimes))  
  allocate(VCDErrorUsingAvKernel(nxtrack,ntimes))            ;    allocate(averagingkernel(nxtrack,ntimes,npressurelevels))
  allocate(CloudRadianceFraction(nxtrack,ntimes))            ;    allocate(int16_buffer3d(nxtrack,ntimes,npressurelevels))
  allocate(OriginalSlantColumnAmountNO2(nxtrack,ntimes))     ;    allocate(TroposphericVerticalColumnError(nxtrack,ntimes))
  
  ! 1-D fields
  start(:)  = 0
  stride(:) = 1
  edge(1)   = ntimes

  taitime = -1. ; htime = -1 ; hdate = -1
  
  retstatus = 10 ; fieldname =  'Time'
  if ( index(trim(gfieldlist),trim(fieldname)) > 0 ) then 
    status = he5_swrdfld( swid, fieldname, start, stride, edge, TAItime)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid)
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  iyear0 = 1950 ; icalendo = 2 ; julday0 = julday(1,1,iyear0)
  call getTMTimeOffset ( itau_offset_TM_OMI )
  do i=1,edge(1)
    if (taitime(i) < 3.e8) cycle
    call tau2date(int(TAITime(i) - TAI_20040101UTC + itau_offset_TM_OMI),dumtime)
    htime(i) = dumtime(4)*10000 + dumtime(5)*100 + dumtime(6)
    hdate(i) = dumtime(1)*10000 + dumtime(2)*100 + dumtime(3)
  end do

  retstatus = retstatus + 1 ; fieldname =  'MeasurementQualityFlags'
  if ( index(trim(dfieldlist),trim(fieldname)) > 0 ) then 
    status = he5_swrdfld( swid, fieldname, start, stride, edge, MeasurementQualityFlags)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if
  
  retstatus = retstatus + 1 ; fieldname =  'InstrumentConfigurationId'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, InstrumentConfigurationId)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'SpacecraftAltitude'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, SpacecraftAltitude)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if
 
  retstatus = retstatus +1 ; fieldname =  'SpacecraftLatitude'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, SpacecraftLatitude)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'SpacecraftLongitude'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,SpacecraftLongitude )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  ! 2-D fields
  start(:)  = 0 ; stride(:) = 1 
  edge(1)   = nxtrack ; edge(2) = ntimes ; edge(3) = NCORNERPOINTS

  fieldname = 'Latitude' ; retstatus = retstatus + 1
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, latitude )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  fieldname = 'Longitude' ; retstatus = retstatus + 1
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, longitude )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname = 'SolarZenithAngle'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, solarZenithangle )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if
  
  retstatus = retstatus +1 ; fieldname =  'SolarAzimuthAngle'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, solarAzimuthangle )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'ViewingZenithAngle'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, viewingZenithangle )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'ViewingAzimuthAngle'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, viewingAzimuthangle )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'CloudFraction'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, int16_buffer )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
    status = he5_swrdlattr(swid,fieldname,"ScaleFactor",scalefactor)
    cloudFraction = real(int16_buffer) * scalefactor
    where(cloudfraction < (-32760.0*scalefactor)) cloudfraction = FILL_FLT32
  end if

  retstatus = retstatus +1 ; fieldname =  'CloudPressure'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, int16_buffer)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
    cloudpressure=real(int16_buffer)
    where(cloudpressure < -32760.0) cloudpressure = FILL_FLT32
  end if

  retstatus = retstatus +1 ; fieldname =  'TerrainHeight'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, terrainheight)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TM4SurfacePressure'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then  
    status = he5_swrdfld( swid, fieldname, start, stride, edge, tm4surfacepressure)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if


  retstatus = retstatus +1 ; fieldname =  'GhostColumn'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, GhostColumn)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'AirMassFactor'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,AirMassFactor )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'AirMassFactorTropospheric'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, AirMassFactorTropospheric)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if
  
  retstatus = retstatus +1 ; fieldname =  'AirMassFactorGeometric'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, AirMassFactorGeometric)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'SurfaceAlbedo'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, int16_buffer)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
    status = he5_swrdlattr(swid,fieldname,"ScaleFactor",scalefactor)
    surfacealbedo = real(int16_buffer) * scalefactor
    where(surfacealbedo < (-32760.0*scalefactor)) surfacealbedo = FILL_FLT32
  end if

  retstatus = retstatus +1 ; fieldname =  'VCDTropErrorUsingAvKernel'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, VCDTropErrorUsingAvKernel)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'VCDErrorUsingAvKernel'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, VCDErrorUsingAvKernel)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'CloudRadianceFraction'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, int16_buffer)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
    status = he5_swrdlattr(swid,fieldname,"ScaleFactor",scalefactor)
    CloudRadianceFraction = real(int16_buffer) * scalefactor    
    where(CloudRadianceFraction < (-32760.0 * scalefactor ) ) CloudRadianceFraction = FILL_FLT32
  end if

  retstatus = retstatus +1 ; fieldname =  'OriginalSlantColumnAmountNO2'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, OriginalSlantColumnAmountNO2)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TroposphericVerticalColumnError'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, TroposphericVerticalColumnError)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TotalVerticalColumn'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,TotalVerticalColumn )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TotalVerticalColumnError'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,TotalVerticalColumnError )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TroposphericVerticalColumn'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, TroposphericVerticalColumn)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'AssimilatedStratosphericSlantColumn'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, AssStratSlantColumn)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'AssimilatedStratosphericVerticalColumn'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,AssStratVerticalColumn )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'AssimilatedStratosphericVerticalColumnError'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, AssStratVerticalColumnError)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'SlantColumnAmountNO2'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, SlantColumnAmountNO2)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'SlantColumnAmountNO2Std'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, SlantColumnAmountNO2Std)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TroposphericVerticalColumnModel'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,TroposphericVerticalColumnModel )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TM4TropoPauseLevel'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, TM4TropoPauseLevel)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TM4TerrainHeight'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, TM4TerrainHeight)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'TroposphericColumnFlag'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, TroposphericColumnFlag)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'GroundPixelQualityFlags'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, GroundPixelQualityFlags)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'CloudPressureStd'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, int16_buffer)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
    CloudPressureStd=real(int16_buffer)
    where(CloudPressureStd < -32760.0)  CloudPressureStd = FILL_FLT32
  end if

  retstatus = retstatus +1 ; fieldname =  'CloudFractionStd'
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, int16_buffer)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
    status = he5_swrdlattr(swid,fieldname,"ScaleFactor",scalefactor)
    CloudFractionStd = real(int16_buffer) * scalefactor
    where(cloudfractionstd < (-32760.0*scalefactor)) cloudfractionstd = FILL_FLT32
  end if

  ! 3-D fields
  start(:)  = 0 ; stride(:) = 1 
  edge(1)   = nxtrack ; edge(2) = ntimes ; edge(3) = NCORNERPOINTS

  retstatus = retstatus +1 ; fieldname =  'LatitudeCornerpoints'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge, LatitudeCornerpoints)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'LongitudeCornerpoints'
  if (index(trim(gfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,LongitudeCornerpoints )
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
  end if

  retstatus = retstatus +1 ; fieldname =  'AveragingKernel' ; edge(3)=npressurelevels
  if (index(trim(dfieldlist),trim(fieldname)) > 0) then
    status = he5_swrdfld( swid, fieldname, start, stride, edge,int16_buffer3d)
    if (status /= 0 ) then 
      write(*,'(2a)') 'ERROR reading ',fieldname ; status = he5_swdetach(swid) 
      status = he5_swclose(fid) ; goto 100 
    end if
    status = he5_swrdlattr(swid,fieldname,"ScaleFactor",scalefactor)
    averagingkernel = real(int16_buffer3d) * scalefactor
    where(averagingkernel < (-32760.0*scalefactor)) averagingkernel = FILL_FLT32
  end if

  ! detach the swath & close the file
  status = he5_swdetach(swid) ; status = he5_swclose(fid)   

!  if (kind(TM4TropoPauseLevel) /= kind(datastruct%TM4TropoPauseLevel) )              stop 'datatype inconsistency for TM4TropoPauseLevel'
  if (kind(TM4TerrainHeight) /= kind(datastruct%TM4TerrainHeight) )                  stop 'datatype inconsistency for TM4TerrainHeight'
!  if (kind(TroposphericColumnFlag) /= kind(datastruct%TroposphericColumnFlag) )      stop 'datatype inconsistency for TroposphericColumnFlag'

  datastruct(1:nxtrack,1:ntimes)%SolarZenithAngle    = SolarZenithAngle 
  datastruct(1:nxtrack,1:ntimes)%SolarAzimuthAngle   = SolarAzimuthAngle      
  datastruct(1:nxtrack,1:ntimes)%longitude           = longitude     
  datastruct(1:nxtrack,1:ntimes)%latitude            = latitude      
  datastruct(1:nxtrack,1:ntimes)%ViewingZenithAngle  = ViewingZenithAngle         
  datastruct(1:nxtrack,1:ntimes)%ViewingAzimuthAngle = ViewingAzimuthAngle       
  datastruct(1:nxtrack,1:ntimes)%CloudFraction       = CloudFraction
  datastruct(1:nxtrack,1:ntimes)%CloudPressure       = CloudPressure
  datastruct(1:nxtrack,1:ntimes)%Terrainheight       = Terrainheight
  datastruct(1:nxtrack,1:ntimes)%tm4surfacepressure  = tm4surfacepressure  

  datastruct(1:nxtrack,1:ntimes)%GhostColumn                                  = GhostColumn                               
  datastruct(1:nxtrack,1:ntimes)%AirMassFactor                                = AirMassFactor                             
  datastruct(1:nxtrack,1:ntimes)%AirMassFactorTropospheric                    = AirMassFactorTropospheric                 
  datastruct(1:nxtrack,1:ntimes)%AirMassFactorGeometric                       = AirMassFactorGeometric                    
  datastruct(1:nxtrack,1:ntimes)%SurfaceAlbedo                                = SurfaceAlbedo                             
  datastruct(1:nxtrack,1:ntimes)%VCDTropErrorUsingAvKernel                    = VCDTropErrorUsingAvKernel                 
  datastruct(1:nxtrack,1:ntimes)%VCDErrorUsingAvKernel                        = VCDErrorUsingAvKernel                     
  datastruct(1:nxtrack,1:ntimes)%CloudRadianceFraction                        = CloudRadianceFraction                     
  datastruct(1:nxtrack,1:ntimes)%OriginalSlantColumnAmountNO2                 = OriginalSlantColumnAmountNO2              
  datastruct(1:nxtrack,1:ntimes)%TroposphericVerticalColumnError              = TroposphericVerticalColumnError           
  datastruct(1:nxtrack,1:ntimes)%TotalVerticalColumn                          = TotalVerticalColumn                       
  datastruct(1:nxtrack,1:ntimes)%TotalVerticalColumnError                     = TotalVerticalColumnError                  
  datastruct(1:nxtrack,1:ntimes)%TroposphericVerticalColumn                   = TroposphericVerticalColumn                
  datastruct(1:nxtrack,1:ntimes)%AssStratSlantColumn                          = AssStratSlantColumn       
  datastruct(1:nxtrack,1:ntimes)%AssStratVerticalColumn                       = AssStratVerticalColumn    
  datastruct(1:nxtrack,1:ntimes)%AssStratVerticalColumnError                  = AssStratVerticalColumnError
  datastruct(1:nxtrack,1:ntimes)%SlantColumnAmountNO2                         = SlantColumnAmountNO2                      
  datastruct(1:nxtrack,1:ntimes)%SlantColumnAmountNO2Std                      = SlantColumnAmountNO2Std                   
  datastruct(1:nxtrack,1:ntimes)%TroposphericVerticalColumnModel              = TroposphericVerticalColumnModel           
  datastruct(1:nxtrack,1:ntimes)%TM4TropoPauseLevel                           = TM4TropoPauseLevel
  datastruct(1:nxtrack,1:ntimes)%TM4TerrainHeight                             = TM4TerrainHeight
  datastruct(1:nxtrack,1:ntimes)%TroposphericColumnFlag                       = TroposphericColumnFlag
  datastruct(1:nxtrack,1:ntimes)%GroundPixelQualityFlags                      = GroundPixelQualityFlags
  datastruct(1:nxtrack,1:ntimes)%CloudPressureStd                             = CloudPressureStd
  datastruct(1:nxtrack,1:ntimes)%CloudFractionStd                             = CloudFractionStd
  do i=1,NCORNERPOINTS
    datastruct(1:nxtrack,1:ntimes)%LatitudeCornerpoints(i)                    = LatitudeCornerpoints(:,:,i)
    datastruct(1:nxtrack,1:ntimes)%LongitudeCornerpoints(i)                   = LongitudeCornerpoints(:,:,i)   
  end do
  do i=1,npressurelevels
    datastruct(1:nxtrack,1:ntimes)%averagingkernel(i)                         = averagingkernel(:,:,i)
  end do

  
  datastruct(1,1:ntimes)%TAItime                     = TAItime  
  datastruct(1,1:ntimes)%htime                       = htime
  datastruct(1,1:ntimes)%hdate                       = hdate
  datastruct(1,1:ntimes)%MeasurementQualityFlags     = MeasurementQualityFlags
  datastruct(1,1:ntimes)%InstrumentConfigurationId   = InstrumentConfigurationId
  datastruct(1,1:ntimes)%SpacecraftAltitude          = SpacecraftAltitude
  datastruct(1,1:ntimes)%SpacecraftLatitude          = SpacecraftLatitude
  datastruct(1,1:ntimes)%SpacecraftLongitude         = SpacecraftLongitude
      
  do i = 2,nxtrack
    datastruct(i,:)%TAItime                   = datastruct(1,:)%TAItime
    datastruct(i,:)%htime                     = datastruct(1,:)%htime
    datastruct(i,:)%hdate                     = datastruct(1,:)%hdate
    datastruct(i,:)%MeasurementQualityFlags   = datastruct(1,:)%MeasurementQualityFlags
    datastruct(i,:)%InstrumentConfigurationId = datastruct(1,:)%InstrumentConfigurationId    
    datastruct(i,:)%SpacecraftLatitude        = datastruct(1,:)%SpacecraftLatitude
    datastruct(i,:)%SpacecraftLongitude       = datastruct(1,:)%SpacecraftLongitude
    datastruct(i,:)%SpacecraftAltitude        = datastruct(1,:)%SpacecraftAltitude
  end do
  
  where(datastruct%GroundPixelQualityFlags < 0)   datastruct%GroundPixelQualityFlags   = datastruct%GroundPixelQualityFlags   + 65536
  where(datastruct%MeasurementQualityFlags < 0)   datastruct%MeasurementQualityFlags   = datastruct%MeasurementQualityFlags   + 256
  where(datastruct%InstrumentConfigurationId < 0) datastruct%InstrumentConfigurationId = datastruct%InstrumentConfigurationId + 256
  where(datastruct%TroposphericColumnFlag <0)     datastruct%TroposphericColumnFlag    = datastruct%TroposphericColumnFlag    + 256
  where(datastruct%TM4TropoPauseLevel <0)         datastruct%TM4TropoPauseLevel        = datastruct%TM4TropoPauseLevel        + 256
  
! everything succesful
  retStatus = 0

  !  deallocate temporary arrays
100 if (allocated(latitude))                                   deallocate(latitude)                                   
  if (allocated(longitude))                                    deallocate(longitude)                                  
  if (allocated(SolarZenithAngle))                             deallocate(SolarZenithAngle)                           
  if (allocated(SolarAzimuthAngle))                            deallocate(SolarAzimuthAngle)                          
  if (allocated(ViewingZenithAngle))                           deallocate(ViewingZenithAngle)                         
  if (allocated(ViewingAzimuthAngle))                          deallocate(ViewingAzimuthAngle)                        
  if (allocated(cloudfraction))                                deallocate(cloudfraction)                              
  if (allocated(cloudpressure))                                deallocate(cloudpressure)                              
  if (allocated(terrainheight))                                deallocate(terrainheight)                              
  if (allocated(TM4SurfacePressure))                           deallocate(TM4SurfacePressure)                         
  if (allocated(TAItime))                                      deallocate(TAItime)
  if (allocated(htime))                                        deallocate(htime)
  if (allocated(hdate))                                        deallocate(hdate)
  if (allocated(SpacecraftAltitude))                           deallocate(SpacecraftAltitude)                         
  if (allocated(SpacecraftLatitude))                           deallocate(SpacecraftLatitude)                         
  if (allocated(SpacecraftLongitude))                          deallocate(SpacecraftLongitude)                        
  if (allocated(GhostColumn))                                  deallocate(GhostColumn)                                
  if (allocated(AirMassFactor))                                deallocate(AirMassFactor)                              
  if (allocated(AirMassFactorTropospheric))                    deallocate(AirMassFactorTropospheric)                  
  if (allocated(AirMassFactorGeometric))                       deallocate(AirMassFactorGeometric)                     
  if (allocated(SurfaceAlbedo))                                deallocate(SurfaceAlbedo)                              
  if (allocated(VCDTropErrorUsingAvKernel))                    deallocate(VCDTropErrorUsingAvKernel)                  
  if (allocated(VCDErrorUsingAvKernel))                        deallocate(VCDErrorUsingAvKernel)                      
  if (allocated(CloudRadianceFraction))                        deallocate(CloudRadianceFraction)                      
  if (allocated(OriginalSlantColumnAmountNO2))                 deallocate(OriginalSlantColumnAmountNO2)               
  if (allocated(TroposphericVerticalColumnError))              deallocate(TroposphericVerticalColumnError)            
  if (allocated(TotalVerticalColumn))                          deallocate(TotalVerticalColumn)                        
  if (allocated(TotalVerticalColumnError))                     deallocate(TotalVerticalColumnError)                   
  if (allocated(TroposphericVerticalColumn))                   deallocate(TroposphericVerticalColumn)                 
  if (allocated(AssStratSlantColumn))                          deallocate(AssStratSlantColumn)        
  if (allocated(AssStratVerticalColumn))                       deallocate(AssStratVerticalColumn)     
  if (allocated(AssStratVerticalColumnError))                  deallocate(AssStratVerticalColumnError)
  if (allocated(SlantColumnAmountNO2))                         deallocate(SlantColumnAmountNO2)                       
  if (allocated(SlantColumnAmountNO2Std))                      deallocate(SlantColumnAmountNO2Std)                    
  if (allocated(TroposphericVerticalColumnModel))              deallocate(TroposphericVerticalColumnModel)            
  if (allocated(LatitudeCornerpoints))                         deallocate(LatitudeCornerpoints)                       
  if (allocated(LongitudeCornerpoints))                        deallocate(LongitudeCornerpoints)                      
  if (allocated(TM4TropoPauseLevel))                           deallocate(TM4TropoPauseLevel)                         
  if (allocated(TM4TerrainHeight))                             deallocate(TM4TerrainHeight)                           
  if (allocated(TroposphericColumnFlag))                       deallocate(TroposphericColumnFlag)                     
  if (allocated(InstrumentConfigurationId))                    deallocate(InstrumentConfigurationId)           
  if (allocated(MeasurementQualityFlags))                      deallocate(MeasurementQualityFlags)             
  if (allocated(GroundPixelQualityFlags))                      deallocate(GroundPixelQualityFlags)                    
  if (allocated(CloudPressureStd))                             deallocate(CloudPressureStd)
  if (allocated(CloudFractionStd))                             deallocate(CloudFractionStd)
  if (allocated(int16_buffer))                                 deallocate(int16_buffer)
  if (allocated(int32_buffer))                                 deallocate(int32_buffer)
  if (allocated(int8_buffer))                                  deallocate(int8_buffer)
  if (allocated(averagingkernel))                              deallocate(averagingkernel)
  if (allocated(int16_buffer3d))                               deallocate(int16_buffer3d)
                                                                
end subroutine read_OMI_tropNO2                               

!*********************************************

subroutine find_selected_data(no2data,lat_limits,lon_limits,index_array,nxtrack,ntimes)
  implicit none
  
  real,intent(in)                                                 :: lat_limits(2),lon_limits(2)  
  type(tdatastruct),dimension(nxtrack,ntimes),intent(in)          :: no2data
  integer,intent(in)                                              :: nxtrack,ntimes
  integer                                                         :: n_elements,l_index,img,row
  integer,dimension(max_nxtrack,max_ntimes),intent(out)           :: index_array

!  nxtrack = size(no2data,dim=1)
!  ntimes  = size(no2data,dim=2)

  index_array(:,:) = 0 ; n_elements = 0
  do row=1,nxtrack
    l_index = 2
    index_array(row,1) = 0
    do img=1,ntimes
      if ((no2data(row,img)%latitude  >= lat_limits(1)) .AND. (no2data(row,img)%latitude  <= lat_limits(2)) .AND. &
          (no2data(row,img)%longitude >= lon_limits(1)) .AND. (no2data(row,img)%longitude <= lon_limits(2))) then
        index_array(row,l_index) = img
        index_array(row,1) = index_array(row,1) + 1
        l_index = l_index + 1 ; n_elements = n_elements + 1
      end if
    end do
  end do  

  if (n_elements <= 0) then
    write(*,'(a)') 'no data selected'
  else
    write(*,'(a,i)') '# of datapoints selected ',n_elements  
  end if

end subroutine find_selected_data

!**************************************

subroutine copy_selected_data(no2data,selected_data,index_array)
  implicit none

  integer                                                         :: l_index,i,j
  integer,intent(in),dimension(max_nxtrack,max_ntimes)            :: index_array
  type(tdatastruct),dimension(max_nxtrack,max_ntimes),intent(in)  :: no2data
  type(tdatastruct),dimension(*)                                  :: selected_data
  
  l_index = 0
  do i=1,max_nxtrack
    do j=2,index_array(i,1)+1
      l_index = l_index + 1
      selected_data(l_index)%SolarZenithAngle    = no2data(i,index_array(i,j))%SolarZenithAngle   
      selected_data(l_index)%SolarAzimuthAngle   = no2data(i,index_array(i,j))%SolarAzimuthAngle  
      selected_data(l_index)%longitude           = no2data(i,index_array(i,j))%longitude          
      selected_data(l_index)%latitude            = no2data(i,index_array(i,j))%latitude           
      selected_data(l_index)%ViewingZenithAngle  = no2data(i,index_array(i,j))%ViewingZenithAngle 
      selected_data(l_index)%ViewingAzimuthAngle = no2data(i,index_array(i,j))%ViewingAzimuthAngle
      selected_data(l_index)%CloudFraction       = no2data(i,index_array(i,j))%CloudFraction      
      selected_data(l_index)%CloudPressure       = no2data(i,index_array(i,j))%CloudPressure      
      selected_data(l_index)%Terrainheight       = no2data(i,index_array(i,j))%Terrainheight  
      selected_data(l_index)%TM4surfacePressure  = no2data(i,index_array(i,j))%TM4surfacePressure    
      selected_data(l_index)%TAItime             = no2data(i,index_array(i,j))%TAItime
      selected_data(l_index)%htime               = no2data(i,index_array(i,j))%htime
      selected_data(l_index)%hdate               = no2data(i,index_array(i,j))%hdate
    end do                                                   
  end do                                                     
  
  write(*,'(i,a)') l_index,' elements copied'

end subroutine copy_selected_data

!**************************************

subroutine read_omi_tropno2_metadata(filename,metadata,N,retstatus)
  use metadata_module

  implicit none
  character(len=STD_STRING),intent(in)        :: filename
  character(len=metadata_value_length)        :: strbuf
  character(len=SWATHNAME_LENGTH)             :: swathname,attrname
  character(len=3000)                         :: attrnames
  integer,intent(in)                          :: N
  integer                                     :: i,nattributes,listlength
  integer                                     :: swid,fid,status,retstatus
  logical                                     :: fexist
  type(Tmetadata),dimension(N)                :: metadata
  
  retstatus = 1 ; fid = -1 ; swid = -1
  write(swathname,'(a)') '                                                          '
  write(swathname,'(a)') 'DominoNO2'
  inquire (file=filename,exist=fexist)
  
  if (fexist == .false.) then
    write(*,'(3a)') 'read_omi_tropno2_metadata: ERROR file ',trim(filename),' does not exist'
    return
  end if

  ! open the file with read-only access
  fid = he5_swopen( filename, HE5F_ACC_RDONLY )         
  if (fid <= 0 ) then
    write(*,'(2a)') 'read_omi_tropno2_metadata: ERROR cannot open file ',trim(filename)
    return
  end if
  
  ! attach to the swath
  swid = he5_swattach( fid, swathname )  
  if (swid <= 0 ) then
    write(*,'(2a)') 'read_omi_tropno2_metadata: ERROR cannot attach swath ',trim(swathname)
    status = he5_swclose(fid)
    return
  end if
  
  nattributes = he5_swinqattrs(swid,attrnames,listlength)

  if ( nattributes < 1) then
    write(*,'(a)') 'There are no attributes in this file'
    return
  end if

  do i=1,N
    write(strbuf,'(a)') '                                          '
    write(attrname,'(a)') trim(metadata(i)%fieldname)
    write(metadata(i)%value,'(a)') 'N/A'
    if ( index(trim(attrnames),trim(attrname)) < 1 ) cycle
    status = he5_swrdattr(swid, trim(attrname), strbuf)  
    if (status /= 0 ) then
      write(*,'(a)') 'Error reading swath attribute'
      cycle
    end if
    write(metadata(i)%value,'(a)') trim(strbuf)
  end do

  ! detach the swath & close the file
  status = he5_swdetach(swid) ; status = he5_swclose(fid)   
  
  retstatus = 0
  
end subroutine read_omi_tropno2_metadata

!*****************************************************

subroutine copy_data(data_in,data_out,nxtrack,ntimes,npressurelevels)
  implicit none
  integer,intent(in)                                             :: nxtrack,ntimes,npressurelevels
  integer                                                        :: i
  type(tdatastruct),dimension(max_nxtrack,max_ntimes),intent(in) :: data_in
  type(tdatastruct),dimension(nxtrack,ntimes)                    :: data_out


  data_out%SolarZenithAngle                            = data_in%SolarZenithAngle
  data_out%SolarAzimuthAngle                           = data_in%SolarAzimuthAngle
  data_out%longitude                                   = data_in%longitude
  data_out%latitude                                    = data_in%latitude
  data_out%ViewingZenithAngle                          = data_in%ViewingZenithAngle
  data_out%ViewingAzimuthAngle                         = data_in%ViewingAzimuthAngle
  data_out%TM4SurfacePressure                          = data_in%TM4SurfacePressure
  data_out%CloudPressure                               = data_in%CloudPressure
  data_out%CloudPressureStd                            = data_in%CloudPressureStd
  data_out%CloudFraction                               = data_in%CloudFraction
  data_out%CloudFractionStd                            = data_in%CloudFractionStd
  data_out%SpacecraftAltitude                          = data_in%SpacecraftAltitude
  data_out%SpacecraftLatitude                          = data_in%SpacecraftLatitude
  data_out%SpacecraftLongitude                         = data_in%SpacecraftLongitude
  data_out%GhostColumn                                 = data_in%GhostColumn
  data_out%AirMassFactor                               = data_in%AirMassFactor
  data_out%AirMassFactorTropospheric                   = data_in%AirMassFactorTropospheric
  data_out%AirMassFactorGeometric                      = data_in%AirMassFactorGeometric
  data_out%SurfaceAlbedo                               = data_in%SurfaceAlbedo
  data_out%VCDTropErrorUsingAvKernel                   = data_in%VCDTropErrorUsingAvKernel
  data_out%VCDErrorUsingAvKernel                       = data_in%VCDErrorUsingAvKernel
  data_out%CloudRadianceFraction                       = data_in%CloudRadianceFraction
  data_out%OriginalSlantColumnAmountNO2                = data_in%OriginalSlantColumnAmountNO2
  data_out%TotalVerticalColumn                         = data_in%TotalVerticalColumn
  data_out%TotalVerticalColumnError                    = data_in%TotalVerticalColumnError
  data_out%TroposphericVerticalColumn                  = data_in%TroposphericVerticalColumn
  data_out%TroposphericVerticalColumnError             = data_in%TroposphericVerticalColumnError          
  data_out%AssStratSlantColumn                         = data_in%AssStratSlantColumn
  data_out%AssStratVerticalColumn                      = data_in%AssStratVerticalColumn
  data_out%AssStratVerticalColumnError                 = data_in%AssStratVerticalColumnError
  data_out%SlantColumnAmountNO2                        = data_in%SlantColumnAmountNO2
  data_out%SlantColumnAmountNO2Std                     = data_in%SlantColumnAmountNO2Std
  data_out%TroposphericVerticalColumnModel             = data_in%TroposphericVerticalColumnModel
  data_out%TAItime                                     = data_in%TAItime
  data_out%htime                                       = data_in%htime
  data_out%hdate                                       = data_in%hdate
  data_out%MeasurementQualityFlags                     = data_in%MeasurementQualityFlags
  data_out%InstrumentConfigurationId                   = data_in%InstrumentConfigurationId
  data_out%GroundPixelQualityFlags                     = data_in%GroundPixelQualityFlags
  data_out%TerrainHeight                               = data_in%TerrainHeight
  data_out%TM4TerrainHeight                            = data_in%TM4TerrainHeight
  data_out%TM4TropoPauseLevel                          = data_in%TM4TropoPauseLevel
  data_out%TroposphericColumnFlag                      = data_in%TroposphericColumnFlag

  do i=1,NCORNERPOINTS
    data_out%LatitudeCornerpoints(i)   = data_in%LatitudeCornerpoints(i)
    data_out%LongitudeCornerpoints(i)  = data_in%LongitudeCornerpoints(i)
  end do

  do i=1,npressurelevels
    data_out%averagingkernel(i)  = data_in%averagingkernel(i)
  end do

end subroutine copy_data

!*****************************************************

subroutine inquire_dimensions(filename,nxtrack,ntimes,npressurelevels,verbose,retstatus)
  implicit none

  integer(int32)                            :: numbertype

  integer,intent(out)                       :: retstatus
  integer                                   :: fid,status,swid
  integer                                   :: rank,dims(3)
  integer,intent(out)                       :: nxtrack,ntimes,npressurelevels
  logical                                   :: fexist
  logical,intent(in)                        :: verbose
  character(len=STD_STRING)                 :: dimlist, maxdimlist
  character(len=STD_STRING),intent(in)      :: filename
  character(len=SWATHNAME_LENGTH)           :: swathname,fieldname
  character(len=5000)                       :: dfieldlist,gfieldlist
  integer                                   :: ndflds,ngflds,pos
  integer,dimension(100)                    :: ntype,nrank

  retstatus = 1 ; fid = -1 ; swid = -1
  write(swathname,'(a)') '                                                          '
  write(swathname,'(a)') 'DominoNO2'
  inquire (file=filename,exist=fexist)
  
  if (fexist == .false.) then
    write(*,'(3a)') 'read_OMI_tropNO2: ERROR file ',trim(filename),' does not exist'
    return
  end if

  ! open the file with read-only access
  fid = he5_swopen( filename, HE5F_ACC_RDONLY )         
  if (fid <= 0 ) then
    write(*,'(2a)') 'read_OMI_tropNO2: ERROR cannot open file ',trim(filename)
    return
  end if
  if (verbose) write(*,'(3a)') 'File ',trim(filename),' opened'
  
  ! attach to the swath
  swid = he5_swattach( fid, swathname )  
  if (swid <= 0 ) then
    write(*,'(2a)') 'read_OMI_tropNO2: ERROR cannot attach swath ',trim(swathname)
    status = he5_swclose(fid)
    return
  end if
  if (verbose) write(*,'(2a)') 'Attached swath ',trim(swathname)

  ndflds = he5_swinqdflds(swid,dfieldlist,nrank,ntype)
  ngflds = he5_swinqgflds(swid,gfieldlist,nrank,ntype)
  if (verbose) write(*,'("nr of data fields: ",i3," nr of geo fields: ",i3)') ndflds,ngflds
  
  pos = index(trim(gfieldlist),'Latitude')
  if (pos < 1) then
    write(*,'(a)') 'The file is invalid, it doesnt contain the Latitude field'
    return
  end if

  pos = index(trim(dfieldlist),'AveragingKernel')
  if (pos < 1) then
    write(*,'(a)') 'The file is invalid, it doesnt contain the AveragingKernel field'
    return
  end if

  ! get dimensions
  fieldname = 'AveragingKernel'
  status = he5_swfldinfo(swid, fieldname, rank, dims,numbertype, dimlist, maxdimlist)  
  if (status /= 0 ) then
    write(*,'(a)') 'read_OMI_tropNO2: ERROR retrieving fieldinfo'
    status = he5_swdetach(swid)
    status = he5_swclose(fid)
    return
  end if

  status = he5_swdetach(swid)
  status = he5_swclose(fid)

  nxtrack = dims(1)
  ntimes  = dims(2)
  npressurelevels = dims(3)

  if ( (nxtrack > MAX_NXTRACK) .OR. (nxtrack <= 0) ) then
    write(*,'(a,i)') 'read_OMI_tropNO2: ERROR invalid value of nxtrack ',nxtrack
    return
  end if

  if ( (ntimes > MAX_NTIMES) .OR. (ntimes <= 0)) then
    write(*,'(a,i)') 'read_OMI_tropNO2: ERROR invalid value of ntimes ',ntimes
    return
  end if

  if ( (npressurelevels > MAX_presslevels) .OR. (npressurelevels <= 0)) then
    write(*,'(a,i)') 'read_OMI_tropNO2: ERROR invalid value of npressurelevels ',npressurelevels
    return
  end if

  if (verbose) write(*,'(2a,i3,a,i5,a,i3)') trim(filename),' has dimensions nxtrack=',dims(1),' ntimes=',dims(2),' npressurelevels=',dims(3)

  retstatus = 0

end subroutine inquire_dimensions

end module read_hdfeos_module
