module write_hdf4_module
  !
  ! Write output data file of the NO2 retrieval-assimilation 
  ! in HDF4 format 
  !
  use datetime_module, only: date2tau,tau2date,Hdatetime2idate
  use datastruct_module
  use constants_Module  
  use metadata_module
  use hdf4_module

  implicit none
!  private
!  public :: writeNO2HdfHeader, writeNO2HdfTrack, writeNO2HdfMetadata


contains

  integer function write2HDF4(filename,no2data,metadata,a_lev,b_lev,N_images,npressurelevels,first)  
    implicit none
    
    type(Tdatastruct),dimension(N_images),intent(in)      :: no2data
    type(Tmetadata),dimension(NMetadataFields),intent(in) :: metadata
    integer,intent(in)                                    :: N_images,npressurelevels
    real,dimension(npressurelevels),intent(in)            :: a_lev, b_lev
    integer                                               :: retstatus
    character(len=STD_STRING)                             :: filename
    logical,intent(in)                                    :: first
    
    retstatus = -1
    if (first) call writeNO2HdfHeader(filename,a_lev,b_lev,npressurelevels)
    call writeNO2HdfTrack(filename,no2data,metadata,N_images,npressurelevels)

    retstatus = 0
    write2HDF4=retstatus

  end function write2HDF4
  
!**********************************************************************  
  
  subroutine writeNO2HdfMetadata(filename,metadata)  
    implicit none
    character(len=STD_STRING),intent(in)       :: filename
    type(Tmetadata),dimension(NMetadataFields) :: metadata
    integer                                    :: i,j,io,status
    character(len=STD_STRING)                  :: VersionStr,Authors,Email,Affil,proc_mode
    character(len=SWATHNAME_LENGTH)            :: attrfieldname

    character(8)                               :: cdate
    character(10)                              :: ctime
    character(5)                               :: czone
    integer,dimension(8)                       :: dvalues
    !
    character(11)                              :: creation_date
    character(3),dimension(12)                 :: monthnames = &
         (/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec' /)


!    write(*,'(2a)') 'writeNO2HdfMetadata: hdf file = ', trim(filename)
    do i=1,NMetadataFields
      if (index(metadata(i)%fieldname,'Authors') > 0) then ; write(Authors,'(a)') trim(metadata(i)%value) ; metadata(i)%used=.false. ; end if
      if (index(metadata(i)%fieldname,'Email') > 0)   then ; write(Email,'(a)') trim(metadata(i)%value) ; metadata(i)%used=.false. ; end if
      if (index(metadata(i)%fieldname,'Affiliation') > 0) then ; write(Affil,'(a)') trim(metadata(i)%value) ; metadata(i)%used=.false. ; end if
      if (index(metadata(i)%fieldname,'Processing_mode') > 0) then ; write(proc_mode,'(a)') trim(metadata(i)%value) ; metadata(i)%used=.false. ; end if
      if (index(metadata(i)%fieldname,'PGE_name') > 0) then
        do j=1,NMetadataFields ; if (index(metadata(j)%fieldname,'PGE_version')>0) exit ; end do
        write(VersionStr,'(3a)') trim(metadata(i)%value),' Version ',trim(metadata(j)%value)
        metadata(i)%used=.false. ; metadata(j)%used=.false.
      end if
    end do
            
    ! ----------------------------------------------
    ! write header
    ! ----------------------------------------------

    call date_and_time(cdate,ctime,czone,dvalues)
    write(creation_date,'(i2,x,a3,i5)') dvalues(3),monthnames(dvalues(2)),dvalues(1)

    io = sfstart(filename,DFACC_WRITE)
    status = sfscatt(io, 'Authors',dfnt_char,         len(trim(Authors)),   trim(Authors) )
    status = sfscatt(io, 'Affiliation',dfnt_char,     len(trim(Affil)),     trim(Affil) )
    status = sfscatt(io, 'Email',dfnt_char,           len(trim(Email)),     trim(Email) )
    status = sfscatt(io, 'Data_created_by',dfnt_char, len(trim(VersionStr)),trim(VersionStr) )
    status = sfscatt(io, 'Creation_date',dfnt_char, 11, creation_date )
    status = sfscatt(io, 'Processing_mode',dfnt_char,len(trim(proc_mode)),trim(proc_mode))
    status = sfscatt(io, 'Unit_of_NO2_column',dfnt_char, 18, '1e15 molecules/cm2')

    do i=1,NMetadataFields
      if (index(metadata(i)%fieldname,'Input_pointer') < 0) cycle
      if (.not. metadata(i)%used) cycle
      if (len(trim(metadata(i)%value)) < 1 ) cycle
      write(attrfieldname,'(2a)') trim(metadata(i)%fieldname)
      status = sfscatt(io, trim(attrfieldname),dfnt_char,len(trim(metadata(i)%value)),trim(metadata(i)%value))
    end do
    status = sfend(io)	  
  
  end subroutine writeNO2HdfMetadata
    
!**********************************************************************    
  
  subroutine writeNO2HdfHeader(filename,a_lev,b_lev,n_lev)
    !------------------------------------------------------
    ! Write the header of the NO2 data file
    !------------------------------------------------------

    implicit none
    !
    character(len=STD_STRING), intent(in)      :: filename
    integer,intent(in)                         :: n_lev
    real,dimension(n_lev),intent(in)           :: a_lev, b_lev

    !
    integer                                    :: io, fid, vd_ref, vd_id, status
    integer                                    :: num_of_records, buf_size
    integer,dimension(2*n_lev+1)               :: databuf

    io = sfstart(filename,DFACC_CREATE)    
    status = sfend(io)
    
    ! ----------------------------------------------
    ! write pressure vdata table
    ! ----------------------------------------------

    ! Re-open HDF file
    fid = hopen(filename, dfacc_write, 0)

    ! Start vdata interface
    status = vfstart(fid)

    ! Create new vdata and attach to it
    vd_ref = -1
    vd_id = vsfatch(fid, vd_ref, 'w')    
    status = vsfsnam(vd_id, 'pressure_grid')
    status = vsfscls(vd_id, 'PRESSURE')
   
    ! Add attributes
    status = vsfsnat(vd_id, -1, 'Unit_of_pressure',dfnt_char, 6, 'Pascal')
    status = vsfsnat(vd_id, -1, 'Pressure_level_formula',dfnt_char, 24, 'p = a_lev + p_surf*b_lev')
    status = vsfsnat(vd_id, -1, 'Number_of_pressure_levels',dfnt_int32, 1, n_lev )

    ! Initialize fields
    status = vsffdef(vd_id, 'a_lev', dfnt_float32, 1)
    status = vsffdef(vd_id, 'b_lev', dfnt_float32, 1)

    ! set field list
    status = vsfsfld(vd_id, 'a_lev,b_lev')

    ! Pack the data
    buf_size = 4*2*n_lev
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size,n_lev, 'a_lev', a_lev )
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size,n_lev, 'b_lev', b_lev )

    ! Write the data
    num_of_records = vsfwrit(vd_id, databuf, n_lev, full_interlace)

    ! Detach from vdata
    status = vsfdtch(vd_id)

    ! End vdata interace 
    status = vfend(fid)

    ! close the HDF file
    status = hclose(fid)

  end subroutine writeNO2HdfHeader

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

  subroutine writeNO2HdfTrack(filename, no2data,metadata,N_images, npressurelevels)
    
    implicit none
    !
    character(len=STD_STRING), intent(in)            :: filename
    integer,intent(in)                               :: npressurelevels,N_images
    type(Tdatastruct),dimension(N_images),intent(in) :: no2data
    type(Tmetadata),dimension(NMetadataFields)       :: metadata
    integer, parameter                               :: nfields_1 = 16
    integer                                          :: total_order_1 
    
    character(len=9), dimension(nfields_1) :: fld_names_1 = &
         (/ 'date     ','time     ','lon      ','lat      ','vcd      ', &
            'sigvcd   ','vcdtrop  ','sigvcdt  ','vcdstrat ','sigvcds  ', &
            'fltrop   ','psurf    ','sigvcdak ','sigvcdtak','kernel   ','ghostcol ' /)
	    
    integer, parameter :: nfields_2 = 12 , total_order_2 = 1+1+1+1+1+4+4+1+1+1+1+1
    character(len=7), dimension(nfields_2) :: fld_names_2 = &
         (/ 'sza    ','vza    ','saa    ','raa    ','ssc    ','loncorn','latcorn','pixelnr','imagenr',&
            'gpxflag','theight','mheight' /)

    integer, parameter :: nfields_3 = 12 , total_order_3 = 1+1+1+1+1+1+1+1+1+1+1+1
    character(len=7), dimension(nfields_3) :: fld_names_3 = &
         (/ 'scd    ','sigscd ','amf    ','amftrop','amfgeo ','scdstr ', &
            'clfrac ','cltpres','albclr ','crfrac ','ltropo ','fcvcdtr'  /)

    integer :: fid, vd_ref, vd_id, status, num_of_records
    integer :: buf_size,ifield,tracknumber

    integer,dimension(:),allocatable :: databuf
    
    real,dimension(:,:),allocatable  :: latitudecornerpoints, longitudecornerpoints,averagingkernel 
    integer,dimension(6)             :: idate6
    integer                          :: i,pos

    character(len=8)                 :: wridate
    character(len=6)                 :: writime
    character(len=18)                :: vd_name
    character(len=nfields_1*9)       :: fn_list
    character(len=10)                :: trackid
    character(len=5)                 :: c_orbitnumber


    call Hdatetime2idate(no2data(1)%hdate,no2data(1)%hTime,idate6)
    write(trackid,'(5i2.2)') mod(idate6(1),10),idate6(2:5)
        
    do i=1,NMetadataFields
      if (index(metadata(i)%value,'OMI-Aura') > 0)  then
        pos=index(metadata(i)%value,'-o')
        write(c_orbitnumber,'(a)') metadata(i)%value(pos+2:pos+6)
        read(c_orbitnumber,'(i6)') tracknumber
        exit
      end if
    end do

    if (tracknumber < 0) then
      write(*,'(a)') 'writeNO2HdfTrack: Error could not read the orbitnumber'
      return
    end if

    total_order_1=nfields_1 + npressurelevels - 1
    status=0
!    allocate(databuf(total_order_1*obsMaxpix+1),stat=status)
    allocate(databuf(4*2*total_order_1*N_images+1),stat=status)
    if (status /= 0) then
      write(*,*) 'writeNO2HdfTrack: ERROR cannot allocate array databuf'
      stop
    end if

    allocate(averagingkernel(npressurelevels,N_images))
    allocate(latitudecornerpoints(4,N_images))
    allocate(longitudecornerpoints(4,N_images))
    
    do i=1,N_images
      longitudecornerpoints(:,i) = no2data(i)%longitudecornerpoints(:)
      latitudecornerpoints(:,i)  = no2data(i)%latitudecornerpoints(:)
      averagingkernel(:,i)       = no2data(i)%averagingkernel(:)
    end do
    
    ! Open HDF file
    fid = hopen(filename, dfacc_write, 0)

    ! ----------------------------------------------
    ! write NO2 vdata table
    ! ----------------------------------------------

    ! Start vdata interface
    status = vfstart(fid)

    ! Create new vdata and attach to it
    vd_ref = -1
    vd_id = vsfatch(fid, vd_ref, 'w')    

    ! Generate vdata name using starting date-time
    write(wridate,'(i8)')   no2data(1)%hdate
    write(writime,'(i6.6)') no2data(1)%htime
    vd_name = 'NO2_'//trackid(1:10)

!    write(*,*) 'Making vdata :    ', vd_name
    status = vsfsnam(vd_id, vd_name)
    status = vsfscls(vd_id, 'FRESCO_NO2_SCD')

    ! Initialize fields
    status = vsffdef(vd_id, fld_names_1(1), dfnt_int32, 1)
    status = vsffdef(vd_id, fld_names_1(2), dfnt_int32, 1)
    status = vsffdef(vd_id, fld_names_1(3), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(4), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(5), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(6), dfnt_float32, 1)      
    status = vsffdef(vd_id, fld_names_1(7), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(8), dfnt_float32, 1)      
    status = vsffdef(vd_id, fld_names_1(9), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(10), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(11), dfnt_int32, 1)
    status = vsffdef(vd_id, fld_names_1(12), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(13), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(14), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_1(15), dfnt_float32, npressurelevels)
    status = vsffdef(vd_id, fld_names_1(16), dfnt_float32, 1)
    
    ! Set field list
    fn_list=trim(fld_names_1(1))
    do ifield = 2, nfields_1
       fn_list=trim(fn_list)//','//trim(fld_names_1(ifield))
    enddo
!    write(*,*) 'with fields :   ', trim(fn_list)
    status = vsfsfld(vd_id, trim(fn_list))

    ! Pack the data
    buf_size = 4*total_order_1*N_images

    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(1), no2data%hdate) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(2), no2data%htime) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(3), no2data%Longitude)                   
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(4), no2data%Latitude)                     
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(5), no2data%TotalVerticalColumn )         
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(6), no2data%TotalVerticalColumnError )    
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(7), no2data%TroposphericVerticalColumn )  
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(8), no2data%TroposphericVerticalColumnError )  
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(9), no2data%AssStratVerticalColumn )         
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(10), no2data%AssStratVerticalColumnError ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(11), int(no2data%TroposphericColumnFlag,kind=4) )    
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(12), no2data%TM4SurfacePressure*100. )        
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(13), no2data%VCDErrorUsingAvKernel )      
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(14), no2data%VCDTropErrorUsingAvKernel )  
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(15), averagingkernel )           
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_1(16), no2data%ghostcolumn )               

    ! Write the data
    num_of_records = vsfwrit(vd_id, databuf, N_images, full_interlace)

    ! Add attributes: track identifier string (e.g. 0503010050) and
    !                 start + end time (these will be separate vdata's)
    status = vsfsnat(vd_id, -1, 'orbit_identifier', dfnt_char, 10, trackid )
    status = vsfsnat(vd_id, -1, 'orbit_number', dfnt_int32, 1, trackNumber )
    call Hdatetime2idate(no2data(1)%hdate,no2data(1)%hTime,idate6)
    status = vsfsnat(vd_id, -1, 'start_time', dfnt_int32, 6, idate6)
    call Hdatetime2idate(no2data(N_images)%hdate,no2data(N_images)%hTime,idate6)
    status = vsfsnat(vd_id, -1, 'end_time', dfnt_int32, 6, idate6)

    ! Detach from vdata
    status = vsfdtch(vd_id)

    ! ----------------------------------------------
    ! write GEOLOCATION vdata table
    ! ----------------------------------------------

    ! Create new vdata and attach to it
    vd_ref = -1
    vd_id = vsfatch(fid, vd_ref, 'w')    

    ! Generate vdata name using starting date-time
    vd_name = 'GEO_'//trackid(1:10)

!    write(*,*) 'Making vdata :    ', vd_name
    status = vsfsnam(vd_id, vd_name)
    status = vsfscls(vd_id, 'GEOLOCATION DATA')

    ! Initialize fields
    status = vsffdef(vd_id, fld_names_2(1), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_2(2), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_2(3), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_2(4), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_2(5), dfnt_int32  , 1)
    status = vsffdef(vd_id, fld_names_2(6), dfnt_float32, 4)       
    status = vsffdef(vd_id, fld_names_2(7), dfnt_float32, 4)
    status = vsffdef(vd_id, fld_names_2(8), dfnt_int32, 1)
    status = vsffdef(vd_id, fld_names_2(9), dfnt_int32, 1)
    status = vsffdef(vd_id, fld_names_2(10), dfnt_int32, 1)  ! groundpixel quality flags
    status = vsffdef(vd_id, fld_names_2(11), dfnt_int32, 1)  ! terrain height
    status = vsffdef(vd_id, fld_names_2(12), dfnt_int32, 1)  ! model terrain height
	
    ! Set field list
    fn_list=trim(fld_names_2(1))
    do ifield = 2, nfields_2
       fn_list=trim(fn_list)//','//trim(fld_names_2(ifield))
    enddo
!    write(*,*) 'with fields :   ', trim(fn_list)
    status = vsfsfld(vd_id, trim(fn_list))

    ! Pack the data
    buf_size = 4*total_order_2*N_images
    i=1

    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(1), no2data%solarZenithAngle )
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(2), no2data%ViewingZenithAngle )
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(3), no2data%solarAzimuthAngle )
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(4), no2data%ViewingAzimuthAngle )
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(5), int(no2data%InstrumentConfigurationId,kind=4) ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(6), LongitudeCornerpoints )                         
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(7), LatitudeCornerpoints )                          
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(8), int(no2data%pixelnr,kind=4))                    
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(9), int(no2data%imagenr,kind=4))                    
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(10), no2data%GroundPixelQualityFlags)               
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(11), int(no2data%terrainheight,kind=4))             
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_2(12), int(no2data%TM4TerrainHeight,kind=4))          

    ! Write the data
    num_of_records = vsfwrit(vd_id, databuf, N_images, full_interlace)

    ! Detach from vdata
    status = vsfdtch(vd_id)

    ! ----------------------------------------------
    ! write ANCILLARY vdata table
    ! ----------------------------------------------

    ! Create new vdata and attach to it
    vd_ref = -1
    vd_id = vsfatch(fid, vd_ref, 'w')    

    ! Generate vdata name using starting date-time
    vd_name = 'ANC_'//trackid(1:10)

!    write(*,*) 'Making vdata :    ', vd_name
    status = vsfsnam(vd_id, vd_name)
    status = vsfscls(vd_id, 'ANCILLARY DATA')

    ! Initialize fields
    status = vsffdef(vd_id, fld_names_3(1), dfnt_float32, 1)
!    status = vsffdef(vd_id, fld_names_3(2), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_3(2), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_3(3), dfnt_float32, 1)       
    status = vsffdef(vd_id, fld_names_3(4), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_3(5), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_3(6), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_3(7), dfnt_float32, 1)
    status = vsffdef(vd_id, fld_names_3(8), dfnt_float32, 1)       
    status = vsffdef(vd_id, fld_names_3(9),dfnt_float32, 1)       
    status = vsffdef(vd_id, fld_names_3(10),dfnt_float32, 1)       
    status = vsffdef(vd_id, fld_names_3(11),  dfnt_int32, 1)    
    status = vsffdef(vd_id, fld_names_3(12),dfnt_float32, 1)    
    
    ! Set field list
    fn_list=trim(fld_names_3(1))
    do ifield = 2, nfields_3
       fn_list=trim(fn_list)//','//trim(fld_names_3(ifield))
    enddo
!    write(*,*) 'with fields :   ', trim(fn_list)
    status = vsfsfld(vd_id, trim(fn_list))

    ! Pack the data
    buf_size = 4*total_order_3*N_images

    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(1), no2data%SlantColumnAmountNO2 ) 
!    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(2), no2data%OriginalSlantColumnAmountNO2 ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(2), no2data%SlantColumnAmountNO2Std ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(3), no2data%AirMassFactor )  
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(4), no2data%AirMassFactorTropospheric ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(5), no2data%AirMassFactorGeometric )    
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(6), no2data%AssStratSlantColumn ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(7), no2data%CloudFraction ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(8), no2data%CloudPressure*100. ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(9), no2data%SurfaceAlbedo ) 
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(10), no2data%CloudRadianceFraction )    
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(11), int(no2data%TM4TropoPauseLevel,kind=4 ))     
    status = vsfnpak(vd_id, hdf_vspack, ' ', databuf, buf_size, N_images, fld_names_3(12), no2data%TroposphericVerticalColumnModel ) 

    ! Write the data
    num_of_records = vsfwrit(vd_id, databuf, N_images, full_interlace)

    ! Detach from vdata
    status = vsfdtch(vd_id)

    ! -------------------------------------

    ! End vdata interace 
    status = vfend(fid)

    ! close the HDF file
    status = hclose(fid)

    if (allocated(databuf))               deallocate(databuf)
    if (allocated(latitudecornerpoints))  deallocate(latitudecornerpoints)
    if (allocated(longitudecornerpoints)) deallocate(longitudecornerpoints)
    if (allocated(averagingkernel))       deallocate(averagingkernel)

  end subroutine writeNO2HdfTrack

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

!   subroutine writeNO2HdfMetadata
!     use omiL2_metadata, only: Tmetadata,metadata,metadataNO2HdfFileName
!     
!     implicit none
!     
!     integer              :: ii,io,status,i,j,nomno2a,pos
!     integer              :: attr_index,orbit_number,orbit_number_md    
!     integer,dimension(50):: omno2_index    
!     character(len=150)   :: buffer_string,line
!     character(len=16)    :: fieldname
!     character(len=6)     :: c_orbit_number
!     character(len=120)   :: orbit_filename,stripe_correction
!     character(len=64)    :: vdata_name2
!     
!     call sort_and_unify_metadata
!     io = sfstart(metadataNO2HdfFileName,DFACC_WRITE)    
! 
!     ii=1
!     do i=1,size(metadata)
!       if (metadata(i)%field_used) then
!         write(buffer_string,'(i2)') ii
!         write(fieldname,'(a14,a)') 'Input_pointer_',trim(adjustl(buffer_string))
! 	
!         write(buffer_string,'(a)') metadata(i)%field
! 	    write(line,'(a)') buffer_string
! 
!         write(buffer_string,'(a)') metadata(i)%value
! 	
!         if (len(trim(adjustl(buffer_string))) > 0) then
!           if (metadata(i)%field == 'Stripe_correction') then
!             write(line,'(a,x,a8,a)') trim(adjustl(line)),'method: ',trim(adjustl(buffer_string))
!           else
!   	        write(line,'(a,x,a10,a)') trim(adjustl(line)),'filename: ',trim(adjustl(buffer_string))
!           end if
!         end if
! 	
!         write(buffer_string,'(a)') metadata(i)%version
!         if (len(trim(adjustl(buffer_string))) > 0) write(line,'(a,x,a9,a)') trim(adjustl(line)),'Version: ',trim(adjustl(buffer_string))
! 	
! 	    write(line,'(a)') trim(adjustl(line))
!         status = sfscatt(io,fieldname,dfnt_char,len(line),line)
!         ii=ii+1
!       end if
!     end do
! 
!     status=sfend(io)
! 
! ! find all OMNO2A fields
!   nomno2a=0
!   do i=1,size(metadata)
!     if (metadata(i)%field_used .and. metadata(i)%field == 'OMNO2A') then
!       nomno2a=nomno2a+1
!       omno2_index(nomno2a)=i
!     end if
!   end do
! !write(*,'(a,i4)') '# OMNO2A fields:',nomno2a
! 
!   file_id = hopen(metadataNO2HdfFileName, DFACC_WRITE, 0)
!   status = vfstart(file_id)
! !write(*,*) 'filename',metadataNO2HdfFileName,'file_id',file_id,'status',status
! 
!   vd_ref = vsffnd(file_id, 'pressure_grid')
!   vd_ref_old=vd_ref
!   
!   do i=1,50
!     gotonextswath: do 
!       vd_ref     = vsfgid(file_id, vd_ref)
!       vd_id      = vsfatch(file_id, vd_ref, 'r')
!       status     = vsfgnam(vd_id, vdata_name2)
!       if (vdata_name2(1:3) == 'NO2' .or. vd_id == -1) exit gotonextswath
!       status     = vsfdtch(vd_id)
!     end do gotonextswath
!     if (vd_id == -1) exit
!     status=vsfdtch(vd_id)
!     vd_id=vsfatch(file_id,vd_ref,'w')
! 
!     attr_index=vsffdat(vd_id,-1,'orbit_number')
!     status=vsfgnat(vd_id,-1,attr_index,orbit_number)
! !write(*,*) 'Swathname',vdata_name2,'orbit number',orbit_number
! 
!     do j=1,nomno2a
!       write(orbit_filename,'(a)') metadata(omno2_index(j))%value
!       write(stripe_correction,'(a)') metadata(omno2_index(j))%stripe_correction
!       pos=index(orbit_filename,'-o')
!       c_orbit_number=orbit_filename(pos+2:pos+6)
!       read(c_orbit_number,'(i6)') orbit_number_md
!       if (orbit_number_md == orbit_number) then
!         status=vsfsnat(vd_id,-1,'input_pointer',dfnt_char,len(trim(orbit_filename)),orbit_filename)
!         status=vsfsnat(vd_id,-1,'stripe_correction',dfnt_char,len(trim(stripe_correction)),stripe_correction)
! !write(*,*) 'orbit filename',orbit_filename
!         exit
!       end if
!     end do
!     status=vsfdtch(vd_id)
!   end do
!     
!   status = sfend(file_id)
!   status = hclose(file_id)
!   
!  end subroutine writeNO2HdfMetadata


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


end module write_hdf4_module
