program main
  use constants_module
  use datastruct_module
  use read_hdfeos_module
  use metadata_module
  use write_hdf4_module
  use index_module
  implicit none


  integer                                              :: status,nxtrack,nimages,npressurelevels
  integer                                              :: ntotimages,i,ii,file_index,nfiles,ndays
  integer                                              :: Nentries,day_index
  integer,dimension(:),allocatable                     :: index_array
  character(len=STD_STRING),dimension(:),allocatable   :: filelist
  character(len=STD_STRING)                            :: filename,outputfilename,outputpath,datestring
  real,dimension(MAX_presslevels)                      :: dumA,dumB
  real,dimension(:),allocatable                        :: a_lev,b_lev
  type(Tdatastruct),dimension(max_nxtrack,max_ntimes)  :: temp_no2data
  type(Tdatastruct),dimension(:),allocatable           :: selected_data
  type(Tdatastruct),dimension(:,:),allocatable         :: no2data
  type(Tmetadata),dimension(:),allocatable             :: file_metadata,metadata
  type(Tindex),dimension(:),allocatable                :: day_index_array,tmp_index_array
  logical                                              :: first,favail
  
  !write(outputpath,'(a)') '/usr/people/dirksen/programs/fortran/DOMINO_he5_hdf4_converter/'
  !write(filename,'(a)') '/usr/people/dirksen/programs/fortran/DOMINO_he5_hdf4_converter/OMIFilelist'

  write(outputpath,'(a)') '/usr/people/boersma/omi/test_v1.01/hdf5to4/f90/'
  write(filename,'(a)') '/usr/people/boersma/omi/test_v1.01/hdf5to4/f90/OMIFilelist'

  inquire(file=filename,exist=favail)
  if (.not. favail) then
    write(*,'(3a)') 'ERROR: the file ',trim(filename),' does not exist'
    stop
  end if
      
  nfiles=Nentries(filename)
  allocate(filelist(nfiles))
  call readfilelist(filename,filelist,nfiles)
  allocate(tmp_index_array(nfiles))
  call CreateDateindex(filelist,nfiles,tmp_index_array,ndays)
  allocate(day_index_array(ndays))
  do i=1,ndays 
    day_index_array(i)%index=tmp_index_array(i)%index
    write(day_index_array(i)%datestring,'(a)') trim(tmp_index_array(i)%datestring)
  end do

  if (allocated(tmp_index_array)) deallocate(tmp_index_array)

  allocate(metadata(NMetadataFields))
  do i=1,NMetadataFields
    write(metadata(i)%fieldname,'(a)') '                                              '
    write(metadata(i)%value,'(a)') '                                                                 '
    metadata(i)%used = .false.
  end do
  
  first=.true. ; day_index=1
  do file_index=1,nfiles
    if (file_index == day_index_array(day_index)%index) first=.true.
    if (first) ii=0
    write(filename,'(a)') trim(filelist(file_index))
    call read_omi_tropno2(filename,temp_no2data,dumA,dumB,nxtrack,nimages,npressurelevels,.false.,status)
    if ( status /= 0 ) then ; write(*,'(2a)') 'Error reading data from the file ',trim(filename) ; cycle ; end if

    allocate(file_metadata(NMetadataFields)) 
    call read_omi_tropno2_metadata(filename,file_metadata,status)
    if ( status /= 0 ) then
      write(*,'(2a)') 'Error reading metadata from the file ',trim(filename)
      if (allocated(metadata)) deallocate(file_metadata)
      cycle 
    end if

    call copy_metadata(metadata,file_metadata)
    if (first) then
      allocate(a_lev(npressurelevels)) ; allocate(b_lev(npressurelevels))
      a_lev=dumA(1:npressurelevels) ; b_lev=dumB(1:npressurelevels)
    end if

    allocate(no2data(nxtrack,nimages))
    call copy_data2D(temp_no2data,no2data,nxtrack,nimages,npressurelevels)
  
    allocate(index_array(nxtrack*nimages))
    call create_image_index(no2data,index_array,nxtrack,nimages,ntotimages)
    write(*,'(2a)') 'Processing ',trim(filename)
    write(*,'(i6,a,i6,a)') ntotimages,' of ',nxtrack*nimages,' pixels selected'
    allocate(selected_data(ntotimages))
  
    call copy_data1D(no2data,selected_data,index_array,nxtrack,nimages,npressurelevels,ntotimages)

    write(outputfilename,'(4a)') trim(outputpath),'no2track',trim(day_index_array(day_index)%datestring),'.hdf'
    status=write2HDF4(outputfilename,selected_data,file_metadata,a_lev,b_lev,ntotimages,npressurelevels,first)
    ii=ii+1
    ! write metadata/swathattributes only if we are dealing with the last file of the day
    if (file_index == nfiles) then
      call writeNO2HdfMetadata(outputfilename,metadata)
      write(*,'(2a,": ",i3)') '# files written to ',trim(outputfilename),ii
    end if
    if (day_index < ndays) then
      if ((file_index+1) == day_index_array(day_index+1)%index) then
        call writeNO2HdfMetadata(outputfilename,metadata)      
        day_index=day_index+1
        write(*,'(2a,": ",i3)') '# files written to ',trim(outputfilename),ii
      end if
    end if

    if (allocated(selected_data))     deallocate(selected_data)
    if (allocated(no2data))           deallocate(no2data)
    if (allocated(file_metadata))     deallocate(file_metadata)
    if (allocated(index_array))       deallocate(index_array)
    if (first .and. allocated(a_lev)) deallocate(a_lev)
    if (first .and. allocated(b_lev)) deallocate(b_lev)
    first=.false.
  end do

  if (allocated(filelist))        deallocate(filelist)
  if (allocated(file_metadata))   deallocate(file_metadata)
  if (allocated(metadata))        deallocate(metadata)
  if (allocated(day_index_array)) deallocate(day_index_array)
end program main

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

subroutine create_image_index(no2data,index_array,nxtrack,nimages,ntotimages)
  use datastruct_module
  implicit none

  integer,intent(in)                           :: nxtrack,nimages
  integer,intent(out)                          :: ntotimages
  integer,dimension(nxtrack*nimages)           :: index_array
  integer                                      :: ii,i,j
  type(Tdatastruct),dimension(nxtrack,nimages),intent(in) :: no2data

  index_array = 0 ; ntotimages = 0
  
  do j=1,nimages
    do i=1,nxtrack
      ii = j*nxtrack + i-1
      if (no2data(i,j)%SlantColumnAmountNO2 < -1.e10) cycle
      if (no2data(i,j)%TroposphericVerticalColumn < -1.e10) cycle
      ntotimages = ntotimages + 1
      index_array(ntotimages) = ii
    end do
  end do

end subroutine create_image_index

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

subroutine printdata2D(no2data,nxtrack,ntimes,npressurelevels)
  use datastruct_module
  implicit none
  integer,intent(in)                                      :: nxtrack,ntimes,npressurelevels
  type(Tdatastruct),dimension(nxtrack,ntimes),intent(in)  :: no2data
  integer                                                 :: i

  i=npressurelevels
  write(*,'(a,t30,e,t45,f)') 'SurfaceAlbedo',             minval(no2data%SurfaceAlbedo),maxval(no2data%SurfaceAlbedo)
  write(*,'(a,t30,e,t45,f)') 'CloudRadianceFraction',     minval(no2data%CloudRadianceFraction),maxval(no2data%CloudRadianceFraction)
  write(*,'(a,t30,e,t45,f)') 'TroposphericVerticalColumn', minval(no2data%TroposphericVerticalColumn),maxval(no2data%TroposphericVerticalColumn)
  write(*,'(a,t30,e,t45,f)') 'SlantColumnAmountNO2'     , minval(no2data%SlantColumnAmountNO2),maxval(no2data%SlantColumnAmountNO2)
  write(*,'(a,t30,i6,t40,i6)')'InstrumentConfigurationId', minval(no2data%InstrumentConfigurationId),maxval(no2data%InstrumentConfigurationId)
  write(*,'(a,t30,i6,t40,i6)')'GroundPixelQualityFlags',   minval(no2data%GroundPixelQualityFlags),maxval(no2data%GroundPixelQualityFlags)
  write(*,'(a,t30,i6,t40,i6)')'MeasurementQualityFlags',   minval(no2data%MeasurementQualityFlags),maxval(no2data%MeasurementQualityFlags)
  write(*,'(a,t30,i6,t40,i6)')'TM4TropoPauseLevel',        minval(no2data%TM4TropoPauseLevel),     maxval(no2data%TM4TropoPauseLevel)
  write(*,'(a,t30,i6,t40,i6)')'TM4TerrainHeight',          minval(no2data%TM4TerrainHeight),       maxval(no2data%TM4TerrainHeight)
  write(*,'(a,t30,i6,t40,i6)')'TroposphericColumnFlag',    minval(no2data%TroposphericColumnFlag), maxval(no2data%TroposphericColumnFlag)
  write(*,'(a,t30,i6,t40,i6)') 'TerrainHeight',            minval(no2data%TerrainHeight), maxval(no2data%TerrainHeight)
  write(*,'(a,t30,f,t45,f)') 'AveragingKernel', minval(no2data(22,555)%averagingkernel),maxval(no2data(22,555)%averagingkernel)  
end subroutine printdata2D

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

subroutine printdata1D(no2data,nimages)
  use datastruct_module
  implicit none
  integer,intent(in)                               :: nimages
  type(Tdatastruct),dimension(nimages),intent(in)  :: no2data

  write(*,'(a,t30,e,t45,f)') 'SurfaceAlbedo',             minval(no2data%SurfaceAlbedo),maxval(no2data%SurfaceAlbedo)
  write(*,'(a,t30,e,t45,f)') 'CloudRadianceFraction',     minval(no2data%CloudRadianceFraction),maxval(no2data%CloudRadianceFraction)
  write(*,'(a,t30,e,t45,f)') 'TroposphericVerticalColumn', minval(no2data%TroposphericVerticalColumn),maxval(no2data%TroposphericVerticalColumn)
  write(*,'(a,t30,e,t45,f)') 'SlantColumnAmountNO2'     , minval(no2data%SlantColumnAmountNO2),maxval(no2data%SlantColumnAmountNO2)
  write(*,'(a,t30,i6,t40,i6)')'InstrumentConfigurationId', minval(no2data%InstrumentConfigurationId),maxval(no2data%InstrumentConfigurationId)
  write(*,'(a,t30,i6,t40,i6)')'GroundPixelQualityFlags',   minval(no2data%GroundPixelQualityFlags),maxval(no2data%GroundPixelQualityFlags)
  write(*,'(a,t30,i6,t40,i6)')'MeasurementQualityFlags',   minval(no2data%MeasurementQualityFlags),maxval(no2data%MeasurementQualityFlags)
  write(*,'(a,t30,i6,t40,i6)')'TM4TropoPauseLevel',        minval(no2data%TM4TropoPauseLevel),     maxval(no2data%TM4TropoPauseLevel)
  write(*,'(a,t30,i6,t40,i6)')'TM4TerrainHeight',          minval(no2data%TM4TerrainHeight),       maxval(no2data%TM4TerrainHeight)
  write(*,'(a,t30,i6,t40,i6)')'TroposphericColumnFlag',    minval(no2data%TroposphericColumnFlag), maxval(no2data%TroposphericColumnFlag)
  write(*,'(a,t30,i6,t40,i6)') 'TerrainHeight',            minval(no2data%TerrainHeight), maxval(no2data%TerrainHeight)
  write(*,'(a,t30,f,t45,f)') 'AveragingKernel', minval(no2data(555)%averagingkernel),maxval(no2data(555)%averagingkernel)  
end subroutine printdata1D

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

subroutine print_swathattributes(metadata)
  use metadata_module
  implicit none
  integer                                    :: i
  type(Tmetadata),dimension(NMetadataFields) :: metadata
  
  
  write(*,'(a)') 'The swath attributes'
  write(*,*)
  do i=1,NMetadataFields
    if (.not. metadata(i)%used ) cycle
    if ( len(trim(metadata(i)%value)) == 0) then
      write(*,'(i2,2x,a,t22,i5)') i,trim(metadata(i)%fieldname),metadata(i)%numvalue
    else
      write(*,'(i2,2x,a,t22,a)') i,trim(metadata(i)%fieldname),trim(metadata(i)%value)
    end if
  end do
  write(*,*) 

end subroutine print_swathattributes

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

subroutine CreateDateindex(filelist,nfiles,index_array,ndays)
  use constants_module
  use index_module
  implicit none
  character(len=STD_STRING),dimension(nfiles) :: filelist
  integer,intent(in)                          :: nfiles
  integer                                     :: i,pos,ndays
  character(len=10)                           :: datumstring1,datumstring2
  character(len=STD_STRING)                   :: filename
  type(Tindex),dimension(nfiles)              :: index_array

  write(datumstring1,'(a)') '         '
  write(datumstring2,'(a)') datumstring1
  ndays=0
  do i=1,nfiles
    write(filename,'(a)') trim(filelist(i))
    pos=index(trim(filename),'OMDOMINO')
    write(datumstring1,'(a)') trim(filename(pos+9:pos+17))
    if (lgt(datumstring1,datumstring2)) then
      ndays=ndays+1
      index_array(ndays)%index=i
      write(index_array(ndays)%datestring,'(2a)') trim(datumstring1(1:4)),trim(datumstring1(6:9))
    end if
    write(datumstring2,'(a)') datumstring1
  end do

end subroutine CreateDateindex

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

integer function Nentries(filename)
  use constants_module
  implicit none
  integer                              :: nfiles,status
  integer,parameter                    :: kflf=37845  
  character(len=STD_STRING),intent(in) :: filename
  character(len=STD_STRING)            :: dumstr
  logical                              :: favail
  
  nfiles=0 ; open(kflf,file=filename)
  do 
    read(kflf,'(a)',iostat=status) dumstr
    if (status /= 0) exit
    inquire(file=dumstr,exist=favail)
    if ( .not. favail ) cycle
    nfiles=nfiles+1
  end do
  close(kflf)
  Nentries=nfiles
  
end function Nentries

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

subroutine readfilelist(filename,filelist,nfiles)
  use constants_module
  implicit none
  integer,parameter                           :: kflf = 37845
  integer,intent(in)                          :: nfiles
  integer                                     :: ii,i,status
  logical                                     :: favail
  character(len=STD_STRING),intent(in)        :: filename
  character(len=STD_STRING),dimension(nfiles) :: filelist
  character(len=STD_STRING)                   :: fn
    
  open(kflf,file=filename) ;  ; ii=1
  do i=1,nfiles
    read(kflf,'(a)',iostat=status) fn
    if (status /=0) exit
    inquire(file=fn,exist=favail)
    if ( .not. favail ) cycle
    write(filelist(ii),'(a)') trim(fn)
    ii=ii+1
  end do
  close(kflf)

end subroutine readfilelist
