PROGRAM fbcomb
   USE toolspar_kind
   USE obs_fbm
   USE index_sort
   IMPLICIT NONE
   !
   ! Command line arguments for output file and input file
   !
#ifndef NOIARGCPROTO
   INTEGER,EXTERNAL :: iargc
#endif
   INTEGER :: nargs
   CHARACTER(len=256) :: cdoutfile
   CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
   !
   ! Input data
   !
   TYPE(obfbdata),POINTER :: obsdata(:)
   INTEGER :: ninfiles,ntotobs,nlev
   !
   ! Time sorting arrays
   !
   REAL(KIND=dp),ALLOCATABLE :: zsort(:,:)
   INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
   !
   ! Output data
   !
   TYPE(obfbdata) :: obsoutdata
   !
   ! Loop variables
   !
   INTEGER :: ia,iv,ii,ij, ist
   !
   ! Get number of command line arguments
   !
   nargs = IARGC()
   IF ( nargs < 2 ) THEN
      WRITE(*,'(A)')'Usage:'
      WRITE(*,'(A)')'fbcomb outputfile inputfile1 inputfile2 ...'
      CALL abort()
   ENDIF
   CALL getarg( 1, cdoutfile )
   !
   ! Get input data
   !
   ALLOCATE( obsdata( nargs - 1 ) )
   ALLOCATE( cdinfile( nargs - 1 ) )
   ntotobs = 0
   ninfiles = nargs - 1
   ist=-1
   DO ia=1, ninfiles
      CALL getarg( ia+1, cdinfile(ia) )
      CALL init_obfbdata( obsdata(ia) )
      CALL read_obfbdata( TRIM(cdinfile(ia)), obsdata(ia) )
      WRITE(*,'(2A)')'File = ', TRIM(cdinfile(ia))
      WRITE(*,'(A,I9,A)')'has', obsdata(ia)%nobs, ' observations'
      IF (obsdata(ia)%nobs > 0 .AND. ist < 0) ist=ia             ! find first file with obs in
      ntotobs = ntotobs + obsdata(ia)%nobs
   ENDDO
   WRITE(*,'(A,I8)') 'Total obsfiles : ',ntotobs
   !
   ! Check that the data is confirming
   !
   DO ia=ist+1, ninfiles
      IF ( obsdata(ia)%cdjuldref /= obsdata(ist)%cdjuldref ) THEN
         WRITE(*,*)'Different julian date reference. Aborting'
         CALL abort
      ENDIF
      IF ( obsdata(ia)%nvar /= obsdata(ist)%nvar ) THEN
         WRITE(*,*)'Different number of variables. Aborting'
         CALL abort
      ENDIF
      IF  (obsdata(ia)%nadd /= obsdata(ist)%nadd ) THEN
         WRITE(*,*)'Different number of additional entries. Aborting'
         CALL abort
      ENDIF
      IF ( obsdata(ia)%next /= obsdata(ist)%next ) THEN
         WRITE(*,*)'Different number of additional variables. Aborting'
         CALL abort
      ENDIF
      IF ( obsdata(ia)%lgrid .NEQV. obsdata(ist)%lgrid ) THEN
         WRITE(*,*)'Inconsistent grid search info. Aborting'
         CALL abort
      ENDIF
      DO iv=1, obsdata(ia)%nvar
         IF ( obsdata(ia)%cname(iv) /= obsdata(ist)%cname(iv) ) THEN
            WRITE(*,*)'Variable name ', TRIM(obsdata(ia)%cname(iv)), &
               &      ' is different from ', TRIM(obsdata(ist)%cname(iv)), &
               &      '. Aborting'
            CALL abort
         ENDIF
         IF ( obsdata(ist)%lgrid .AND. obsdata(ia)%nobs > 0) THEN
            IF ( obsdata(ia)%cgrid(iv) /= obsdata(ist)%cgrid(iv) ) THEN
               WRITE(*,*)'Grid name ', TRIM(obsdata(ia)%cgrid(iv)), &
                  &      ' is different from ', TRIM(obsdata(ist)%cgrid(iv)), &
                  &      '. Aborting'
               CALL abort
            ENDIF
         ENDIF
      ENDDO
      DO iv=1,obsdata(ia)%nadd
         IF ( obsdata(ia)%caddname(iv) /= obsdata(ist)%caddname(iv) ) THEN
            WRITE(*,*)'Additional name ', TRIM(obsdata(ia)%caddname(iv)), &
               &      ' is different from ', TRIM(obsdata(ist)%caddname(iv)), &
               &      '. Aborting'
            CALL abort
         ENDIF
      ENDDO
      DO iv=1,obsdata(ia)%next
         IF ( obsdata(ia)%cextname(iv) /= obsdata(ist)%cextname(iv) ) THEN
            WRITE(*,*)'Extra name ', TRIM(obsdata(ia)%cextname(iv)), &
               &      ' is different from ', TRIM(obsdata(ist)%cextname(iv)), &
               &      '. Aborting'
            CALL abort
         ENDIF
      ENDDO
   ENDDO
   !
   ! Construct sorting arrays
   !
   ALLOCATE( zsort(3,ntotobs), iset(ntotobs), &
      & inum(ntotobs), iindex(ntotobs))
   ii = 0
   DO ia = 1,ninfiles
      DO ij = 1,obsdata(ia)%nobs
         ii = ii+1
         zsort(1,ii) = obsdata(ia)%ptim(ij)
         zsort(2,ii) = obsdata(ia)%pphi(ij)
         zsort(3,ii) = obsdata(ia)%plam(ij)
         iset(ii) = ia
         inum(ii) = ij
      ENDDO
   ENDDO
   !
   ! Get indexes for time sorting.
   !
   CALL index_sort_dp_n(zsort,3,iindex,ntotobs)
   !
   ! Allocate output data
   !   
   nlev = -1
   DO ia = 1,ninfiles
      IF ( obsdata(ia)%nlev > nlev ) nlev = obsdata(ia)%nlev
   ENDDO
   CALL init_obfbdata( obsoutdata )
   CALL alloc_obfbdata( obsoutdata, obsdata(ist)%nvar, ntotobs, nlev, &
      &                 obsdata(ist)%nadd, obsdata(ist)%next, obsdata(ist)%lgrid )
   !
   ! Copy input data into output data
   !
   CALL merge_obfbdata( ninfiles, obsdata, obsoutdata, iset, inum, iindex )
   !
   ! Save output data 
   !
   CALL write_obfbdata ( TRIM(cdoutfile), obsoutdata )
   
END PROGRAM fbcomb