PROGRAM fbmatchup
   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(:)
   CHARACTER(len=ilenname),ALLOCATABLE :: cdnames(:)
   INTEGER :: nqc
   LOGICAL :: ldaily820
   NAMELIST/namfbmatchup/nqc,ldaily820
   !
   ! Input data
   !
   TYPE(obfbdata)         :: obsdatatmp(1)
   TYPE(obfbdata),POINTER :: obsdata(:)
   INTEGER :: ninfiles,ntotobs,nlev
   !
   ! Time sorting arrays
   !
   REAL(KIND=dp),ALLOCATABLE :: zsort(:,:)
   INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
   !
   ! Comparison arrays and scalars
   !
   REAL(KIND=fbsp), ALLOCATABLE :: zrtim(:),zrphi(:),zrlam(:)
   INTEGER(KIND=SELECTED_INT_KIND(12)), ALLOCATABLE :: irwmo(:)
   REAL(KIND=fbsp) :: ztim,zphi,zlam
   INTEGER(KIND=SELECTED_INT_KIND(12)) :: iwmo
   !
   ! Output data
   !
   TYPE(obfbdata) :: obsoutdata
   !
   ! Loop variables
   !
   INTEGER :: ia,ip,i1,ii,ij,ik1,ik2,iv,ist
   LOGICAL :: llfound
   LOGICAL :: lexists
   INTEGER :: ityp
   !
   ! Get number of command line arguments
   !
   nargs = IARGC()
   IF ( ( MOD(nargs,2) /= 1 ) .OR. ( nargs == 0 )  ) THEN
      WRITE(*,'(A)')'Usage:'
      WRITE(*,'(A)')'fbmatchup outputfile inputfile1 varname1 inputfile2 varname2 ...'
      CALL abort()
   ENDIF
   CALL getarg( 1, cdoutfile )
   !
   ! Read namelist if present
   !
   nqc=1
   ldaily820=.FALSE.
   INQUIRE(file='namfbmatchup.in',exist=lexists)
   IF (lexists) THEN
      OPEN(10,file='namfbmatchup.in')
      READ(10,namfbmatchup)
      CLOSE(10)
      WRITE(*,namfbmatchup)
   ENDIF
   !
   ! Get input data
   !
   ninfiles = ( nargs -1 )/ 2
   ALLOCATE( obsdata( ninfiles ) )
   ALLOCATE( cdinfile( ninfiles ) )
   ALLOCATE( cdnames( ninfiles ) )
   ip = 1
   DO ia=1, ninfiles
      !
      ! Read the unsorted file
      !
      ip = ip + 1
      CALL getarg( ip, cdinfile(ia) )
      ip = ip + 1
      CALL getarg( ip, cdnames(ia) )
      CALL init_obfbdata( obsdatatmp(1) )
      CALL read_obfbdata( TRIM(cdinfile(ia)), obsdatatmp(1) )
      !
      ! Check that only one variable present in input file
      !
      IF ( obsdatatmp(1)%nadd > 1 ) THEN
         WRITE(*,*)'Warning. More than one variable in input file'
         WRITE(*,*)'Number of variables = ', obsdatatmp(1)%nadd
         WRITE(*,*)'Only the first one is going to be stored'
      ENDIF
      IF ( obsdatatmp(1)%nadd < 1 ) THEN
         WRITE(*,*)'Error. less than one variable in input file'
         CALL abort
      ENDIF
      !
      ! Check that we have few levels than in the first file
      !
      IF ( ia > 1 )  THEN
         IF ( obsdatatmp(1)%nlev > obsdata(1)%nlev ) THEN
            WRITE(*,*)'Warning. More levels in file than the first file'
            WRITE(*,*)'Number of levels in current file = ', obsdatatmp(1)%nlev
            WRITE(*,*)'Number of levels in first file   = ', obsdata(1)%nlev
            WRITE(*,*)'Only the number of levels in the first'//&
               &' file will be used'
         ENDIF
      ENDIF
      !
      ! Check that we have few observations than in the first file
      !
      IF ( ia > 1 )  THEN
         IF ( obsdatatmp(1)%nobs > obsdata(1)%nobs ) THEN
            WRITE(*,*)'Warning. More obs in file than the first file'
            WRITE(*,*)'Number of obs in current file = ', obsdatatmp(1)%nobs
            WRITE(*,*)'Number of obs in first file   = ', obsdata(1)%nobs
            WRITE(*,*)'Only the observations in the first'//&
               &' file will be stored'
         ENDIF
      ENDIF
      !
      ! Check that we have the same number of variables
      !
      IF ( ia > 1 )  THEN
         IF ( obsdatatmp(1)%nvar /= obsdata(1)%nvar ) THEN
            WRITE(*,*)'Error. Different number of variables.'
            WRITE(*,*)'Number of var in current file = ', obsdatatmp(1)%nvar
            WRITE(*,*)'Number of var in first file   = ', obsdata(1)%nvar
            CALL abort
         ENDIF
      ENDIF
      !
      ! Check reference datas
      !
      IF ( ia > 1 )  THEN
         IF ( obsdatatmp(1)%cdjuldref /= obsdata(1)%cdjuldref ) THEN
            WRITE(*,*)'Different reference dates'
            CALL abort
         ENDIF
      ENDIF
      !
      ! Special fix for daily average MRB data (820) for the first file
      !
      IF (ldaily820.AND.(ia==1)) THEN
         DO ij = 1,obsdatatmp(1)%nobs
            READ(obsdatatmp(1)%cdtyp(ij),'(I5)')ityp
            IF (ityp==820) THEN
               obsdatatmp(1)%ptim(ij)=INT(obsdatatmp(1)%ptim(ij))+1.0
            ENDIF
         ENDDO
      ENDIF
      !
      ! Construct sorting arrays
      !
      ALLOCATE( zsort(3,obsdatatmp(1)%nobs), iset(obsdatatmp(1)%nobs), &
         & inum(obsdatatmp(1)%nobs), iindex(obsdatatmp(1)%nobs))
      ii = 0
      DO ij = 1,obsdatatmp(1)%nobs
         ii = ii+1
         zsort(1,ii) = obsdatatmp(1)%ptim(ij)
         zsort(2,ii) = obsdatatmp(1)%pphi(ij)
         zsort(3,ii) = obsdatatmp(1)%plam(ij)
         iset(ii) = 1
         inum(ii) = ij
      ENDDO
      !
      ! Get indexes for time sorting.
      !
      CALL index_sort_dp_n(zsort,3,iindex,obsdatatmp(1)%nobs)
      CALL init_obfbdata( obsdata(ia) )
      CALL alloc_obfbdata( obsdata(ia), &
         &                 obsdatatmp(1)%nvar, obsdatatmp(1)%nobs, &
         &                 obsdatatmp(1)%nlev, obsdatatmp(1)%nadd, &
         &                 obsdatatmp(1)%next, obsdatatmp(1)%lgrid )
      !
      ! Copy input data into output data
      !
      CALL merge_obfbdata( 1, obsdatatmp, obsdata(ia), iset, inum, iindex )
      CALL dealloc_obfbdata( obsdatatmp(1) )
      
      WRITE(*,'(2A)')'File = ', TRIM(cdinfile(ia))
      WRITE(*,'(A,I9,A)')'has', obsdata(ia)%nobs, ' observations'

      DEALLOCATE( zsort, iset, inum, iindex )

   ENDDO
   !
   ! Prepare output data
   !
   CALL init_obfbdata( obsoutdata )
   !
   ! Copy the first input data to output data
   !
   CALL copy_obfbdata( obsdata(1), obsoutdata, &
      &                kadd = ninfiles )
   DO ia = 1, ninfiles
      obsoutdata%caddname(ia) = cdnames(ia)
      obsoutdata%caddlong(ia,:) = obsdata(ia)%caddlong(1,:)
      obsoutdata%caddunit(ia,:) = obsdata(ia)%caddunit(1,:)
   ENDDO
   !
   ! Allocate comparison arrays and file them
   ! 
   IF (ilenwmo>8) THEN
      WRITE(*,*)'Fix fbmatchup to allow string length > 8'
      CALL abort
   ENDIF
   ALLOCATE(zrtim(obsoutdata%nobs),zrphi(obsoutdata%nobs), &
      &     zrlam(obsoutdata%nobs),irwmo(obsoutdata%nobs))
   DO i1 = 1, obsoutdata%nobs
      irwmo(i1) = TRANSFER( obsoutdata%cdwmo(i1), irwmo(i1) )
      zrtim(i1) = REAL( obsoutdata%ptim(i1), fbsp )
      zrphi(i1) = REAL( obsoutdata%pphi(i1), fbsp )
      zrlam(i1) = REAL( obsoutdata%plam(i1), fbsp ) 
   ENDDO
   !
   ! Merge extra data into output data
   !
   DO ia = 2, ninfiles
      ist = 1
      DO ii = 1, obsdata(ia)%nobs
         IF (MOD(ii,10000)==1) THEN
            WRITE(*,*)'Handling observation no ',ii,' for file no ',ia
         ENDIF
         llfound = .FALSE.
         iwmo = TRANSFER( obsdata(ia)%cdwmo(ii), iwmo )
         ztim = REAL( obsdata(ia)%ptim(ii), fbsp )
         zphi = REAL( obsdata(ia)%pphi(ii), fbsp )
         zlam = REAL( obsdata(ia)%plam(ii), fbsp ) 
         ! Check if the the same index is the right one.
         IF ( iwmo == irwmo(ii) ) THEN
            IF ( ztim == zrtim(ii) ) THEN
               IF (  zphi == zrphi(ii) ) THEN
                  IF ( zlam == zrlam(ii) ) THEN
                     llfound = .TRUE.
                     DO iv = 1, obsdata(ia)%nvar
                        ! Since the inner loop don't change the
                        ! qc decisions use this to ensure match
                        ! for duplicate observations
                        IF ( obsdata(ia)%ivqc(ii,iv) /= &
                           & obsoutdata%ivqc(ii,iv)) THEN
                           llfound = .FALSE.
                           CYCLE
                        ENDIF
                     ENDDO
                     IF (llfound) i1 = ii
                  ENDIF
               ENDIF
            ENDIF
         ENDIF
         ! Search for position in from previous found position
         ! if not the same index
         IF (.NOT.llfound) THEN
            DO i1 = ist, obsoutdata%nobs
               IF ( iwmo == irwmo(i1) ) THEN
                  IF ( ztim == zrtim(i1) ) THEN
                     IF (  zphi == zrphi(i1) ) THEN
                        IF ( zlam == zrlam(i1) ) THEN
                           llfound = .TRUE.
                           DO iv = 1, obsdata(ia)%nvar
                              ! Since the inner loop don't change the
                              ! qc decisions use this to ensure match
                              ! for duplicate observations
                              IF ( obsdata(ia)%ivqc(ii,iv) /= &
                                 & obsoutdata%ivqc(i1,iv)) THEN
                                 llfound = .FALSE.
                                 WRITE (*,*)'QC flags different for ',&
                                    &  TRIM(obsdata(ia)%cdwmo(ii)),' at ', &
                                    & obsdata(ia)%ptim(ii)
                                 CYCLE
                              ENDIF
                           ENDDO
                           IF (llfound) EXIT
                        ENDIF
                     ENDIF
                  ENDIF
               ENDIF
            ENDDO
         ENDIF
         ! If not fount try agan from the beginnning
         IF ( .NOT.llfound ) THEN
            DO i1 = 1, obsoutdata%nobs
               IF ( iwmo == irwmo(i1) ) THEN
                  IF ( ztim == zrtim(i1) ) THEN
                     IF (  zphi == zrphi(i1) ) THEN
                        IF ( zlam == zrlam(i1) ) THEN
                           llfound = .TRUE.
                           DO iv = 1, obsdata(ia)%nvar
                              ! Since the inner loop don't change the
                              ! qc decisions use this to ensure match
                              ! for duplicate observations
                              IF ( obsdata(ia)%ivqc(ii,iv) /= &
                                 & obsoutdata%ivqc(i1,iv)) THEN
                                 llfound = .FALSE.
                                 WRITE (*,*)'QC flags different for ',&
                                    &  TRIM(obsdata(ia)%cdwmo(ii)),' at ', &
                                    & obsdata(ia)%ptim(ii)
                                 CYCLE
                              ENDIF
                           ENDDO
                           IF (llfound) EXIT
                        ENDIF
                     ENDIF
                  ENDIF
               ENDIF
            ENDDO
         ENDIF
         ! If found put the data into the common structure
         IF (llfound) THEN
            IF ( nqc == ia ) THEN
               obsoutdata%ioqc(i1)   = obsdata(ia)%ioqc(ii)
               obsoutdata%ipqc(i1)   = obsdata(ia)%ipqc(ii)
               obsoutdata%itqc(i1)   = obsdata(ia)%itqc(ii)
               obsoutdata%ivqc(i1,:) = obsdata(ia)%ivqc(ii,:)
            ENDIF
            obsoutdata%ioqcf(:,i1)   = IOR( obsdata(ia)%ioqcf(:,ii), &
               &                            obsoutdata%ioqcf(:,i1) )
            obsoutdata%ipqcf(:,i1)   = IOR( obsdata(ia)%ipqcf(:,ii), &
               &                            obsoutdata%ipqcf(:,i1) )
            obsoutdata%itqcf(:,i1)   = IOR( obsdata(ia)%itqcf(:,ii), &
               &                            obsoutdata%itqcf(:,i1) )
            obsoutdata%ivqcf(:,i1,:) = IOR( obsdata(ia)%ivqcf(:,ii,:), &
               &                            obsoutdata%ivqcf(:,i1,:) )
            llfound = .FALSE.
            ! Search for levels
            DO ik1 = 1, obsdata(ia)%nlev
               DO ik2 = 1, obsoutdata%nlev
                  IF ( REAL( obsdata(ia)%pdep(ik1,ii), fbsp ) ==  &
                     & REAL( obsoutdata%pdep(ik2,i1),  fbsp ) ) THEN
                     obsoutdata%padd(ik2,i1,ia,:) = &
                        & obsdata(ia)%padd(ik1,ii,1,:)
                     IF ( nqc == ia ) THEN
                        obsoutdata%idqc(ik2,i1)    = obsdata(ia)%idqc(ik1,ii)
                        obsoutdata%ivlqc(ik2,i1,:) = obsdata(ia)%ivlqc(ik1,ii,:)
                     ENDIF
                     obsoutdata%idqcf(:,ik2,i1)    = &
                        &                IOR( obsdata(ia)%idqcf(:,ik1,ii), &
                        &                     obsoutdata%idqcf(:,ik2,i1) )
                     obsoutdata%ivlqcf(:,ik2,i1,:) = &
                        &                IOR( obsdata(ia)%ivlqcf(:,ik1,ii,:), &
                        &                     obsoutdata%ivlqcf(:,ik2,i1,:) )
                     llfound = .TRUE.
                     EXIT
                  ENDIF
               ENDDO
               ! Write warning if level not found
               IF (.NOT.llfound.AND.(obsdata(ia)%pdep(ik1,ii)/=fbrmdi)) THEN
                  WRITE(*,*)'Level not found in first file : ',&
                     &      TRIM( cdinfile(1)  )
                  WRITE(*,*)'Data file                     : ',&
                     &      TRIM( cdinfile(ia) )
                  WRITE(*,*)'Identifier                    : ',&
                     &      obsdata(ia)%cdwmo(ii)
                  WRITE(*,*)'Julian date                   : ',&
                     &      obsdata(ia)%ptim(ii)
                  WRITE(*,*)'Latitude                      : ',&
                     &      obsdata(ia)%pphi(ii)
                  WRITE(*,*)'Longitude                     : ',&
                     &      obsdata(ia)%plam(ii)
                  WRITE(*,*)'Depth                         : ',&
                     &      obsdata(ia)%pdep(ik1,ii)
               ENDIF
            ENDDO
            ist = i1
         ELSE
            ! Write warning if observation not found
            WRITE(*,*)'Observation not found in first data file : ',&
               &      TRIM( cdinfile(1)  )
            WRITE(*,*)'Data file                                : ',&
               &      TRIM( cdinfile(ia) )
            WRITE(*,*)'Identifier                               : ',&
               &      obsdata(ia)%cdwmo(ii)
            WRITE(*,*)'Julian date                              : ',&
               &      obsdata(ia)%ptim(ii)
            WRITE(*,*)'Latitude                                 : ',&
               &      obsdata(ia)%pphi(ii)
            WRITE(*,*)'Longitude                                : ',&
               &      obsdata(ia)%plam(ii)
            ist = 1
         ENDIF
      ENDDO
      IF (obsdata(ia)%nobs>0) THEN
         WRITE(*,*)'Handled last obs. no    ',ii,' for file no ',ia
      ENDIF
   ENDDO
   !
   ! Write output file
   !
   CALL write_obfbdata( TRIM(cdoutfile), obsoutdata )

END PROGRAM fbmatchup