Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/bufr2fb.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/bufr2fb.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/bufr2fb.F90	(revision 2893)
@@ -0,0 +1,183 @@
+PROGRAM bufr2fb
+#ifdef HAS_BUFR
+   USE toolspar_kind
+   USE obs_fbm
+   USE bufrdata
+   USE convmerge
+   IMPLICIT NONE
+   !
+   ! Command line arguments for output file and input files
+   !
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdtsfilem, cdtsfile1, cduvfilem, cduvfile1
+   CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
+   !
+   ! Input data
+   !
+   TYPE(obfbdata), POINTER :: bufrtsm(:), bufrts1(:), bufruvm(:), bufruv1(:)
+   INTEGER :: ninfiles
+   INTEGER :: ntottsm, nmaxlevts, ntotts1
+   INTEGER :: ntotuvm, nmaxlevuv, ntotuv1
+   !
+   ! Loop variables
+   !
+   INTEGER :: ia
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF (nargs < 4) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'bufr2fb TS_outputfile_prof TS_outputfile_surf '//&
+         &          'UV_outputfile_prof UV_outputfile_surf '//&
+         &          'inputfile1 inputfile2 ...'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdtsfilem)
+   CALL getarg(2,cdtsfile1)
+   CALL getarg(3,cduvfilem)
+   CALL getarg(4,cduvfile1)
+   !
+   ! Get input data
+   !
+   ALLOCATE( bufrtsm(MAX(nargs-4,1)) )
+   ALLOCATE( bufrts1(MAX(nargs-4,1)) )
+   ALLOCATE( bufruvm(MAX(nargs-4,1)) )
+   ALLOCATE( bufruv1(MAX(nargs-4,1)) )
+   ALLOCATE( cdinfile(MAX(nargs-4,1)) )
+   ntottsm = 0
+   ntotts1 = 0
+   ntotuvm = 0
+   ntotuv1 = 0
+   ninfiles  = nargs - 4
+   DO ia=1,ninfiles
+      CALL getarg( ia + 4, cdinfile(ia) )
+      CALL read_bufrfile( TRIM(cdinfile(ia)), &
+         &                bufrtsm(ia), bufrts1(ia), &
+         &                bufruvm(ia), bufruv1(ia) )
+      WRITE(*,'(2A)')'File = ',TRIM(cdinfile(ia))
+      WRITE(*,'(A,I9,A)')'has',bufrtsm(ia)%nobs,' TS profiles'
+      WRITE(*,'(A,I9,A)')'and',bufruvm(ia)%nobs,' UV profiles'
+      WRITE(*,'(A,I9,A)')'and',bufrts1(ia)%nobs,' TS single level'
+      WRITE(*,'(A,I9,A)')'and',bufruv1(ia)%nobs,' UV single level'
+      ntottsm = ntottsm + bufrtsm(ia)%nobs
+      ntotuvm = ntotuvm + bufruvm(ia)%nobs
+      ntotts1 = ntotts1 + bufrts1(ia)%nobs
+      ntotuv1 = ntotuv1 + bufruv1(ia)%nobs
+   ENDDO
+   IF (ninfiles==0) THEN
+      CALL init_obfbdata( bufrtsm(1) )
+      CALL alloc_obfbdata( bufrtsm(1), 2, 0, 1, 0, 1, .FALSE. )
+      bufrtsm(1)%cname(1) = 'POTM'
+      bufrtsm(1)%cname(2) = 'PSAL'
+      bufrtsm(1)%coblong(1) = 'Potential temperature'
+      bufrtsm(1)%coblong(2) = 'Practical salinity'
+      bufrtsm(1)%cobunit(1) = 'Degrees Celsius'
+      bufrtsm(1)%cobunit(2) = 'PSU'
+      bufrtsm(1)%cextname(1) = 'TEMP'
+      bufrtsm(1)%cextlong(1) = 'Insitu temperature'
+      bufrtsm(1)%cextunit(1) = 'Degrees Celsius'
+      bufrtsm(1)%cdjuldref = '19500101000000'
+      CALL init_obfbdata( bufrts1(1) )
+      CALL alloc_obfbdata( bufrts1(1), 2, 0, 1, 0, 1, .FALSE. )
+      bufrts1(1)%cname(1) = 'POTM'
+      bufrts1(1)%cname(2) = 'PSAL'
+      bufrts1(1)%coblong(1) = 'Potential temperature'
+      bufrts1(1)%coblong(2) = 'Practical salinity'
+      bufrts1(1)%cobunit(1) = 'Degrees Celsius'
+      bufrts1(1)%cobunit(2) = 'PSU'
+      bufrts1(1)%cextname(1) = 'TEMP'
+      bufrts1(1)%cextlong(1) = 'Insitu temperature'
+      bufrts1(1)%cextunit(1) = 'Degrees Celsius'
+      bufrts1(1)%cdjuldref = '19500101000000'
+      CALL init_obfbdata( bufruvm(1) )
+      CALL alloc_obfbdata( bufruvm(1), 2, 0, 1, 0, 0, .FALSE. )
+      bufruvm(1)%cname(1) = 'UVEL'
+      bufruvm(1)%cname(2) = 'VVEL'
+      bufruvm(1)%coblong(1) = 'Zonal current'
+      bufruvm(1)%coblong(2) = 'Meridional current'
+      bufruvm(1)%cobunit(1) = 'Meters per second'
+      bufruvm(1)%cobunit(2) = 'Meters per second'
+      bufruvm(1)%cdjuldref = '19500101000000'
+      CALL init_obfbdata( bufruv1(1) )
+      CALL alloc_obfbdata( bufruv1(1), 2, 0, 1, 0, 0, .FALSE. )
+      bufruv1(1)%cname(1) = 'UVEL'
+      bufruv1(1)%cname(2) = 'VVEL'
+      bufruv1(1)%coblong(1) = 'Zonal current'
+      bufruv1(1)%coblong(2) = 'Meridional current'
+      bufruv1(1)%cobunit(1) = 'Meters per second'
+      bufruv1(1)%cobunit(2) = 'Meters per second'
+      bufruv1(1)%cdjuldref = '19500101000000'
+   ENDIF
+   WRITE(*,'(A,I8)') 'Total TS profiles : ',ntottsm
+   WRITE(*,'(A,I8)') 'Total TS single   : ',ntotts1
+   WRITE(*,'(A,I8)') 'Total UV profiles : ',ntotuvm
+   WRITE(*,'(A,I8)') 'Total UV single   : ',ntotuv1
+   !
+   ! Merge and output the data.
+   !
+   IF (TRIM(cdtsfilem)/='none') THEN
+      ! Convert insitu temperature to potential temperature using the model
+      ! salinity if no potential temperature
+      DO ia = 1, ninfiles
+         call fbpotem( bufrtsm(ia) )
+      ENDDO
+      CALL conv_fbmerge( TRIM(cdtsfilem), ninfiles, bufrtsm )
+   ENDIF
+   IF (TRIM(cdtsfile1)/='none') THEN
+      DO ia = 1, ninfiles
+         call fbpotem( bufrts1(ia) )
+      ENDDO
+      CALL conv_fbmerge( TRIM(cdtsfile1), ninfiles, bufrts1 )
+   ENDIF
+   IF (TRIM(cduvfilem)/='none') THEN
+      CALL conv_fbmerge( TRIM(cduvfilem), ninfiles, bufruvm )
+   ENDIF
+   IF (TRIM(cduvfile1)/='none') THEN
+      CALL conv_fbmerge( TRIM(cduvfile1), ninfiles, bufruv1 )
+   ENDIF
+
+CONTAINS
+
+#include "obs_conv.h90"
+   
+   SUBROUTINE fbpotem( fbdata )
+      
+      ! Convert insitu temperature to potential temperature 
+
+      TYPE(obfbdata) :: fbdata
+      REAL :: zpres
+      INTEGER :: jo,jk
+
+      DO jo = 1, fbdata%nobs
+         IF ( fbdata%pphi(jo) < 9999.0 ) THEN
+            DO jk = 1, fbdata%nlev
+               IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
+                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
+                  & ( fbdata%pob(jk,jo,2) < 9999.0 ) .AND. &
+                  & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
+                  zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
+                     &              REAL(fbdata%pphi(jo),wp) )
+                  fbdata%pob(jk,jo,1) = potemp( &
+                     &                     REAL(fbdata%pob(jk,jo,2), wp),  &
+                     &                     REAL(fbdata%pext(jk,jo,1), wp), &
+                     &                     REAL(zpres,wp), 0.0_wp )
+               ENDIF
+            ENDDO
+         ENDIF
+      ENDDO
+
+   END SUBROUTINE fbpotem
+
+#else
+   
+   WRITE(*,'(A)')'bufr2fb compiled without -DHAS_BUFR so '//&
+      &          'running it is pointless'
+   CALL abort
+
+#endif   
+
+END PROGRAM bufr2fb
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/bufrdata.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/bufrdata.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/bufrdata.F90	(revision 2893)
@@ -0,0 +1,879 @@
+#ifdef IBM
+@PROCESS NOOPT
+#endif
+MODULE bufrdata
+
+#ifdef HAS_BUFR
+
+   USE toolspar_kind
+   USE obs_utils
+   USE obs_fbm
+
+   IMPLICIT NONE
+
+CONTAINS
+
+   SUBROUTINE read_bufrfile( cdfilename, inptsm, inpts1, inpuvm, inpuv1 )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_bufrfile **
+      !!
+      !!  ** Purpose : Read from file the profile BUFR observations.
+      !!
+      !!  ** Method  : The data file is a BUFR file. 
+      !!
+      !!  ** Action  :
+      !!
+      !!   History :
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(LEN=*) :: cdfilename  ! Input file.
+      TYPE(obfbdata) :: inptsm        ! Output TS data structure (prof).
+      TYPE(obfbdata) :: inpts1        ! Output TS data structure (single).
+      TYPE(obfbdata) :: inpuvm        ! Output UV data structure (prof).
+      TYPE(obfbdata) :: inpuv1        ! Output UV data structure (single).
+      !! * Local variables
+      INTEGER :: iunit, iret
+      INTEGER :: ibytesize, isubtype, itype
+      INTEGER :: iyea, imon, iday, ihou, imin
+      INTEGER :: jk
+      REAL :: zlat,zlon
+      CHARACTER(len=8) :: csign
+      INTEGER, PARAMETER :: imaxlev = 300
+      INTEGER :: ilevts, iobstsm, iobsts1, ilevtsmax
+      INTEGER :: ilevuv, iobsuvm, iobsuv1, ilevuvmax
+      INTEGER :: iobs
+      REAL, DIMENSION(imaxlev) :: zdts, ztem, zsal, zduv, zuv, zvv
+      real :: zmdi
+
+      ! Initialize data
+      
+      ibytesize = BIT_SIZE(1)/8
+      zmdi = REAL(fbrmdi)
+      
+      ! Open file
+      CALL pbopen( iunit, cdfilename, 'r', iret )
+      IF ( iret /= 0 ) THEN
+         WRITE(*,*)'Error opening file : ', TRIM(cdfilename)
+         WRITE(*,*)'Error code : ', iret
+         CALL abort
+      ENDIF
+
+      ! Count number of observations
+
+      ilevtsmax = -1
+      ilevuvmax = -1
+      iobstsm = 0
+      iobsts1 = 0
+      iobsuvm = 0
+      iobsuv1 = 0
+
+      DO
+
+         CALL decode_bufr( &
+            & iunit, ibytesize,             &
+            & iyea, imon, iday, ihou, imin, &
+            & zlat, zlon,                   &
+            & csign, isubtype,              &
+            & imaxlev, ilevts, ilevuv,      & 
+            & zdts, ztem, zsal,             &
+            & zduv, zuv, zvv,               &
+            & zmdi, iret )
+         
+         ! No more observations if iret==-1
+
+         IF ( iret == -1 ) EXIT
+
+         ! Ignore observations without levels
+
+         IF ( ilevts > 0 ) THEN 
+
+            IF ( ilevts == 1 ) THEN
+               iobsts1 = iobsts1 + 1
+            ELSE
+               iobstsm = iobstsm + 1
+               IF ( ilevts > ilevtsmax ) ilevtsmax = ilevts
+            ENDIF
+
+         ENDIF
+
+         IF ( ilevuv > 0 ) THEN 
+
+            IF ( ilevuv == 1 ) THEN
+               iobsuv1 = iobsuv1 + 1
+            ELSE
+               iobsuvm = iobsuvm + 1
+               IF ( ilevuv > ilevuvmax ) ilevuvmax = ilevuv
+            ENDIF
+
+         ENDIF
+
+      ENDDO
+
+      ! Setup data structures
+
+      CALL init_obfbdata( inptsm )
+      CALL alloc_obfbdata( inptsm, 2, iobstsm, ilevtsmax, 0, 1, .FALSE. )
+      inptsm%cname(1) = 'POTM'
+      inptsm%cname(2) = 'PSAL'
+      inptsm%coblong(1) = 'Potential temperature'
+      inptsm%coblong(2) = 'Practical salinity'
+      inptsm%cobunit(1) = 'Degrees Celsius'
+      inptsm%cobunit(2) = 'PSU'
+      inptsm%cextname(1) = 'TEMP'
+      inptsm%cextlong(1) = 'Insitu temperature'
+      inptsm%cextunit(1) = 'Degrees Celsius'
+      CALL init_obfbdata( inpts1 )
+      CALL alloc_obfbdata( inpts1, 2, iobsts1, 1, 0, 1, .FALSE. )
+      inpts1%cname(1) = 'POTM'
+      inpts1%cname(2) = 'PSAL'
+      inpts1%coblong(1) = 'Potential temperature'
+      inpts1%coblong(2) = 'Practical salinity'
+      inpts1%cobunit(1) = 'Degrees Celsius'
+      inpts1%cobunit(2) = 'PSU'
+      inpts1%cextname(1) = 'TEMP'
+      inpts1%cextlong(1) = 'Insitu temperature'
+      inpts1%cextunit(1) = 'Degrees Celsius'
+      CALL init_obfbdata( inpuvm )
+      CALL alloc_obfbdata( inpuvm, 2, iobsuvm, ilevuvmax, 0, 0, .FALSE. )
+      inpuvm%cname(1) = 'UVEL'
+      inpuvm%cname(2) = 'VVEL'
+      inpuvm%coblong(1) = 'Zonal current'
+      inpuvm%coblong(2) = 'Meridional current'
+      inpuvm%cobunit(1) = 'Meters per second'
+      inpuvm%cobunit(2) = 'Meters per second'
+      CALL init_obfbdata( inpuv1 )
+      CALL alloc_obfbdata( inpuv1, 2, iobsuv1, 1, 0, 0, .FALSE. )
+      inpuv1%cname(1) = 'UVEL'
+      inpuv1%cname(2) = 'VVEL'
+      inpuv1%coblong(1) = 'Zonal current'
+      inpuv1%coblong(2) = 'Meridional current'
+      inpuv1%cobunit(1) = 'Meters per second'
+      inpuv1%cobunit(2) = 'Meters per second'
+
+      ! Rewind the file
+
+      CALL pbseek( iunit, 0, 0, iret )
+      IF ( iret /= 0 ) THEN
+         WRITE(*,*)'Error rewinding file : ',TRIM(cdfilename)
+         WRITE(*,*)'Error code : ',iret
+         CALL abort
+      ENDIF
+      
+      ! Retrieve the observations and put them into the data structure
+
+      iobstsm = 0
+      iobsts1 = 0
+      iobsuvm = 0
+      iobsuv1 = 0
+      iret = 0
+      inptsm%cdjuldref = '19500101000000'
+      inpts1%cdjuldref = '19500101000000'
+      inpuvm%cdjuldref = '19500101000000'
+      inpuv1%cdjuldref = '19500101000000'
+      iobs = 0
+      DO
+
+         CALL decode_bufr( &
+            & iunit, ibytesize,             &
+            & iyea, imon, iday, ihou, imin, &
+            & zlat, zlon,                   &
+            & csign, isubtype,              &
+            & imaxlev, ilevts, ilevuv,      & 
+            & zdts, ztem, zsal,             &
+            & zduv, zuv, zvv,               &
+            & zmdi, iret )
+
+         ! No more observations if iret==-1
+
+         IF ( iret == -1 ) EXIT
+         iobs = iobs + 1
+         
+         ! Convert isubtype to something EN3 like.
+
+         IF ( isubtype == 131 ) THEN
+            itype = 820
+         ELSEIF ( isubtype == 132 ) THEN
+            itype = 401
+         ELSEIF ( isubtype == 133 ) THEN
+            IF (LEN_TRIM(CSIGN)>=7) THEN
+               itype = 831
+            ELSE
+               itype = 741
+            ENDIF
+         ELSE
+            WRITE(*,*)'Unknown subtype = ',isubtype
+            itype = isubtype
+         ENDIF
+
+         ! TS data
+
+         ! Ignore observations without levels
+
+         IF ( ilevts > 0 ) THEN
+
+            IF ( ilevts==1 ) THEN
+
+               iobsts1 = iobsts1 + 1
+
+               ! Position and time
+               
+               inpts1%kindex(iobsts1) = iobs
+               inpts1%pphi(iobsts1) = zlat
+               inpts1%plam(iobsts1) = zlon
+               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
+                  &           inpts1%ptim(iobsts1) )
+
+               ! Call sign and type
+               
+               inpts1%cdwmo(iobsts1) = csign
+               WRITE(inpts1%cdtyp(iobsts1),'(I4.4)') itype
+
+               ! Depth, Salinity and Insitu Temperature
+               
+               DO jk = 1, ilevts
+                  
+                  inpts1%pdep(jk,iobsts1) = zdts(jk)
+                  inpts1%pext(jk,iobsts1,1) = ztem(jk)
+                  inpts1%pob(jk,iobsts1,2) = zsal(jk)
+                  
+                  ! If no insitu temperature set the QC flag to 4
+                  
+                  IF ( ztem(jk) == zmdi ) THEN
+                     inpts1%ivlqc(jk,iobsts1,1) = 4
+                  ELSE
+                     inpts1%ivlqc(jk,iobsts1,1) = 1
+                  ENDIF
+
+                  ! If no salinity set the QC flag to 4
+                  
+                  IF ( zsal(jk) == zmdi ) THEN
+                     inpts1%ivlqc(jk,iobsts1,2) = 4
+                  ELSE
+                     inpts1%ivlqc(jk,iobsts1,2) = 1
+                  ENDIF
+                  
+               ENDDO
+
+            ELSE
+
+               iobstsm = iobstsm + 1
+               
+               ! Position and time
+               
+               inptsm%kindex(iobstsm) = iobs
+               inptsm%pphi(iobstsm) = zlat
+               inptsm%plam(iobstsm) = zlon
+               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
+                  &           inptsm%ptim(iobstsm) )
+
+               ! Call sign and type
+               
+               inptsm%cdwmo(iobstsm) = csign
+               WRITE(inptsm%cdtyp(iobstsm),'(I4.4)') itype
+
+               ! Depth, Salinity and Insitu Temperature
+               
+               DO jk = 1, ilevts
+               
+                  inptsm%pdep(jk,iobstsm) = zdts(jk)
+                  inptsm%pext(jk,iobstsm,1) = ztem(jk)
+                  inptsm%pob(jk,iobstsm,2) = zsal(jk)
+
+                  ! If no insitu temperature set the QC flag to 4
+                  
+                  IF ( ztem(jk) == zmdi ) THEN
+                     inptsm%ivlqc(jk,iobstsm,1) = 4
+                  ELSE
+                     inptsm%ivlqc(jk,iobstsm,1) = 1
+                  ENDIF
+
+                  ! If no salinity set the QC flag to 4
+                  
+                  IF ( zsal(jk) == zmdi ) THEN
+                     inptsm%ivlqc(jk,iobstsm,2) = 4
+                  ELSE
+                     inptsm%ivlqc(jk,iobstsm,2) = 1
+                  ENDIF
+                  
+               ENDDO
+               
+            ENDIF
+
+         ENDIF
+
+         ! UV data
+
+         ! Ignore observations without levels
+
+         IF ( ilevuv > 0 ) THEN
+
+            ! Skip surface only observations (TBC).
+
+            IF ( ilevuv == 1 ) THEN
+
+               iobsuv1 = iobsuv1 + 1
+               
+               ! Position and time
+               
+               inpuv1%kindex(iobsuv1) = iobs
+               inpuv1%pphi(iobsuv1) = zlat
+               inpuv1%plam(iobsuv1) = zlon
+               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
+                  &           inpuv1%ptim(iobsuv1) )
+               
+               ! Call sign and type
+               
+               inpuv1%cdwmo(iobsuv1) = csign
+               WRITE(inpuv1%cdtyp(iobsuv1),'(I4.4)') itype
+
+               ! Depth, Salinity and Insitu Temperature
+               
+               DO jk = 1, ilevuv
+                  
+                  inpuv1%pdep(jk,iobsuv1)  =  zduv(jk)
+                  inpuv1%pob(jk,iobsuv1,1)  =  zuv(jk)
+                  inpuv1%pob(jk,iobsuv1,2)  =  zvv(jk)
+                  
+                  ! If no insitu temperature set the QC flag to 4
+                  
+                  IF ( zuv(jk) == zmdi ) THEN
+                     inpuv1%ivlqc(jk,iobsuv1,1) = 4
+                  ELSE
+                     inpuv1%ivlqc(jk,iobsuv1,1) = 1
+                  ENDIF
+                  
+                  ! If no v velocity set the QC flag to 4
+                  
+                  IF ( zvv(jk) == zmdi ) THEN
+                     inpuv1%ivlqc(jk,iobsuv1,2) = 4
+                  ELSE
+                     inpuv1%ivlqc(jk,iobsuv1,2) = 1
+                  ENDIF
+                  
+               ENDDO
+               
+            ELSE
+
+               iobsuvm = iobsuvm + 1
+               
+               ! Position and time
+               
+               inpuvm%kindex(iobsuvm) = iobs
+               inpuvm%pphi(iobsuvm) = zlat
+               inpuvm%plam(iobsuvm) = zlon
+               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
+                  &           inpuvm%ptim(iobsuvm) )
+               
+               ! Call sign and type
+               
+               inpuvm%cdwmo(iobsuvm) = csign
+               WRITE(inpuvm%cdtyp(iobsuvm),'(I4.4)') itype
+
+               ! Depth, Salinity and Insitu Temperature
+               
+               DO jk = 1, ilevuv
+                  
+                  inpuvm%pdep(jk,iobsuvm)  =  zduv(jk)
+                  inpuvm%pob(jk,iobsuvm,1)  =  zuv(jk)
+                  inpuvm%pob(jk,iobsuvm,2)  =  zvv(jk)
+                  
+                  ! If no insitu temperature set the QC flag to 4
+                  
+                  IF ( zuv(jk) == zmdi ) THEN
+                     inpuvm%ivlqc(jk,iobsuvm,1) = 4
+                  ELSE
+                     inpuvm%ivlqc(jk,iobsuvm,1) = 1
+                  ENDIF
+                  
+                  ! If no v velocity set the QC flag to 4
+                  
+                  IF ( zvv(jk) == zmdi ) THEN
+                     inpuvm%ivlqc(jk,iobsuvm,2) = 4
+                  ELSE
+                     inpuvm%ivlqc(jk,iobsuvm,2) = 1
+                  ENDIF
+                  
+               ENDDO
+
+            ENDIF
+
+
+         ENDIF
+
+
+      ENDDO
+      
+      ! Close the file
+
+      CALL pbclose( iunit, iret )
+      IF ( iret /= 0) THEN
+         WRITE(*,*)'Error closing file : ',TRIM(cdfilename)
+         WRITE(*,*)'Error code : ',iret
+         CALL abort
+      ENDIF
+      
+      ! Set the QC flags which can not be read from BUFR to 0
+
+      inptsm%ipqc(:)         = 0
+      inptsm%ipqcf(:,:)      = 0
+      inptsm%ivqc(:,:)       = 0
+      inptsm%ivqcf(:,:,:)    = 0
+      inptsm%ivqc(:,:)       = 0
+      inptsm%ivqcf(:,:,: )   = 0
+      inptsm%ioqc(:)         = 0
+      inptsm%ioqcf(:,:)      = 0
+      inptsm%idqc(:,:)       = 0
+      inptsm%idqcf(:,:,:)    = 0
+      inptsm%ivlqcf(:,:,:,:) = 0
+      inpts1%ipqc(:)         = 0
+      inpts1%ipqcf(:,:)      = 0
+      inpts1%ivqc(:,:)       = 0
+      inpts1%ivqcf(:,:,:)    = 0
+      inpts1%ivqc(:,:)       = 0
+      inpts1%ivqcf(:,:,: )   = 0
+      inpts1%ioqc(:)         = 0
+      inpts1%ioqcf(:,:)      = 0
+      inpts1%idqc(:,:)       = 0
+      inpts1%idqcf(:,:,:)    = 0
+      inpts1%ivlqcf(:,:,:,:) = 0
+      inpuvm%ipqc(:)         = 0
+      inpuvm%ipqcf(:,:)      = 0
+      inpuvm%ivqc(:,:)       = 0
+      inpuvm%ivqcf(:,:,:)    = 0
+      inpuvm%ivqc(:,:)       = 0
+      inpuvm%ivqcf(:,:,: )   = 0
+      inpuvm%ioqc(:)         = 0
+      inpuvm%ioqcf(:,:)      = 0
+      inpuvm%idqc(:,:)       = 0
+      inpuvm%idqcf(:,:,:)    = 0
+      inpuvm%ivlqcf(:,:,:,:) = 0
+      inpuv1%ipqc(:)         = 0
+      inpuv1%ipqcf(:,:)      = 0
+      inpuv1%ivqc(:,:)       = 0
+      inpuv1%ivqcf(:,:,:)    = 0
+      inpuv1%ivqc(:,:)       = 0
+      inpuv1%ivqcf(:,:,: )   = 0
+      inpuv1%ioqc(:)         = 0
+      inpuv1%ioqcf(:,:)      = 0
+      inpuv1%idqc(:,:)       = 0
+      inpuv1%idqcf(:,:,:)    = 0
+      inpuv1%ivlqcf(:,:,:,:) = 0
+
+   END SUBROUTINE read_bufrfile
+
+   SUBROUTINE decode_bufr( kunit, kwsize, &
+      & kyea, kmon, kday, khou, kmin,     &
+      & plat,plon, csign, ksty,           &
+      & kmlevs, klevts, klevve,           &
+      & pdts, ptem, psal,                 &
+      & pdve, puve, pvve,                 &
+      & pvalnull, ierr)
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE decode_bufr **
+      !!
+      !!  ** Purpose : Decode a singe BUFR record for ocean data
+      !!
+      !!  ** Method  : Call to BUFREX/BUSEL in EMOS lib.
+      !!
+      !!  ** Action  :
+      !!
+      !!   History :   (??-??) A. Vidard  . ODASYS version. 
+      !!               (08-12) K. Mogensen. NEMOVAR version
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT(IN) :: kunit      ! Bufr file unit number
+      INTEGER, INTENT(IN) :: kwsize     ! No of bytes in a word
+      INTEGER, INTENT(OUT) :: &         ! Date
+         & kyea, kmon, kday, khou, kmin
+      REAL, INTENT(OUT) :: &            ! Position
+         & plat, plon
+      CHARACTER(LEN=*), INTENT(OUT) ::  &
+         & csign                        ! Call sign
+      INTEGER, INTENT(OUT) :: ksty      ! Subtype
+      INTEGER, INTENT(IN) :: kmlevs     ! Maximum number of levels
+      INTEGER, INTENT(OUT) :: klevts    ! Actual number of levels
+      INTEGER, INTENT(OUT) :: klevve    ! Actual number of levels
+      REAL, DIMENSION(kmlevs), INTENT(OUT) :: &
+         & pdts, &                      ! Depths for t,s
+         & ptem, &                      ! Insitu temperature
+         & psal, &                      ! Salinity
+         & pdve, &                      ! Depths for currents
+         & puve, &                      ! U velocity
+         & pvve                         ! V velocity
+      REAL, INTENT(IN) :: pvalnull      ! Missing value
+      INTEGER, INTENT(INOUT) :: ierr    ! Error code
+      !! * Define BUFR parameters
+      INTEGER, PARAMETER :: imax_bufr_bytes=80000 ! Max bufr length in bytes
+      INTEGER, PARAMETER :: imax_elements=20000   ! Kax no of elements
+      INTEGER, PARAMETER :: imax_values=20000     ! Max no of values 
+      INTEGER, PARAMETER :: imax_sec0=3           ! Length of section 0
+      INTEGER, PARAMETER :: imax_sec1=40          ! Length of section 1
+      INTEGER, PARAMETER :: imax_sec2=64          ! Length of section 2
+      INTEGER, PARAMETER :: imax_sec3=4           ! Length of section 3
+      INTEGER, PARAMETER :: imax_sec4=2           ! Length of section 4
+      INTEGER, PARAMETER :: idim_ksup=9           ! Length of array ksup
+      !! * Define local variables
+      INTEGER :: ksup(idim_ksup)        ! Expanded bufr additional information
+      INTEGER :: ksec0(imax_sec0)       ! Expanded bufr section 0
+      INTEGER :: ksec1(imax_sec1)       ! Expanded bufr section 1
+      INTEGER :: ksec2(imax_sec2)       ! Expanded bufr section 2
+      INTEGER :: ksec3(imax_sec3)       ! Expanded bufr section 3
+      INTEGER :: ksec4(imax_sec4)       ! Expanded bufr section 4
+      INTEGER :: kdtlen,ktdexl          ! Bufr lengths
+      INTEGER :: &                      ! BUSEL variables
+         & ktdlst(imax_elements), ktdexp(imax_elements)
+      INTEGER :: kbuff(imax_bufr_bytes) ! Bufr message array 
+      REAL :: values(imax_values)       ! Bufr message values
+      CHARACTER(LEN=64) :: &
+         & cnames(imax_elements), &     ! Element names
+         & cunits(imax_elements), &     ! Element units
+         & cvals(imax_values)           ! Character values
+      REAL, DIMENSION((kmlevs)) :: &
+         & zspd, &                      ! Current speed
+         & zdir                         ! Current direction
+      INTEGER :: ilenbuf                ! Word len. of compressed bufr message
+      INTEGER :: kerr                   ! Error return code
+      INTEGER :: ibytes                 ! Actual no. of bytes in bufr message
+      INTEGER :: ji                     ! Loop counter
+      INTEGER :: icur                   ! Pointer to current data
+      ! Define physical constants
+      REAL, PARAMETER :: zero_celcius = 273.15
+
+      ! Initialize variables
+
+      cvals(1) = '         '
+
+      DO ji = 1, kmlevs
+         pdts(ji) = pvalnull
+         ptem(ji) = pvalnull
+         psal(ji) = pvalnull
+         puve(ji) = pvalnull
+         pvve(ji) = pvalnull
+      END DO
+      
+      DO ji = 1, imax_elements
+         ktdlst(ji) = 0
+         ktdexp(ji) = 0
+      END DO
+      values(:) = pvalnull
+
+      ! Read bufr message
+
+      kerr = 0
+      CALL pbbufr( kunit, kbuff, imax_bufr_bytes, ibytes, kerr) 
+      IF ( kerr == -1 ) THEN
+         ierr=-1
+         RETURN
+      ELSE IF ( kerr < -1 ) THEN
+         WRITE(*,*)'ERROR in PBBUFR with code ',kerr
+         CALL abort
+      ENDIF
+      
+      ! Decode and expand bufr message
+
+      ! Calculate length of bufr message in words
+      ilenbuf = ibytes/kwsize+1
+      
+      ! Expand bufr message
+
+      kerr=0
+      CALL bufrex( ilenbuf, kbuff, ksup, ksec0, ksec1, ksec2, ksec3, ksec4,   &
+         &         imax_elements, cnames, cunits, imax_values, values, cvals, &
+         &         kerr)
+      IF ( kerr /= 0 ) THEN
+         WRITE(*,*)'ERROR in BUFREX with code ',kerr
+         klevts = 0
+         klevve = 0
+         RETURN
+      ENDIF
+      CALL busel( kdtlen, ktdlst, ktdexl, ktdexp, kerr )
+      IF ( kerr /= 0 ) THEN
+         WRITE(*,*)'ERROR in BUSEL with code ',kerr
+         klevts = 0
+         klevve = 0
+         RETURN
+      ENDIF
+
+      ! Set oceanographic variables
+
+      ksty=ksec1(7)
+      kyea=ksec1(9)
+
+      ! Patch for 2 digit only date (only matter if > 2000)
+
+      IF ( kyea < 40 ) kyea = 2000 + kyea
+
+      SELECT CASE(ksty) 
+         
+      CASE(131)
+         
+         ! DRIBU
+
+         kyea = values(5)
+         kmon = values(6)
+         kday = values(7)
+         khou = values(8)
+         kmin = values(9)
+         plat = values(10)
+         plon = values(11)
+
+         ! Check that this dribu contains salinity
+
+         IF ( ktdexp(16) == 2033 ) THEN
+
+            ! Get number of levels
+            IF ( values(17) < 100000. ) THEN
+               klevts = INT( values(17) )
+            ELSE
+               ! Patch for nasty data.. without data (AV 200307)
+               klevts = 0
+            ENDIF
+
+         ELSE
+            
+            klevts = 0
+
+         ENDIF
+
+         ! Patch for missing WMO number (AV 200307)
+         IF ( values(1) >=  9999999. ) THEN
+            csign = 'XXX'
+         ELSE
+            IF ( INT( values(1) ) .LT. 100000 ) THEN
+               WRITE(csign(1:5),'(i5)') INT( values(1) )
+            ELSE
+               WRITE(csign,     '(i8)') INT( values(1) )
+            ENDIF
+         ENDIF
+
+         CALL bufr_get_data( 7062, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 18, 3, pdts )
+         
+         CALL bufr_get_data( 22043, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 19, 3, ptem )
+         
+         CALL bufr_get_data( 22062, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 20, 3, psal )
+
+         IF ( ktdexp(15) == 2031 ) THEN
+
+            klevve = INT( values(16) )
+
+            CALL bufr_get_data( 7062, imax_elements, ktdexp, &
+               &                imax_values, values, klevve, &
+               &                17, 3, pdve )
+            
+            CALL bufr_get_data( 22004, imax_elements, ktdexp, &
+               &                imax_values, values, klevve, &
+               &                18, 3, zdir )
+            
+            CALL bufr_get_data( 22031, imax_elements, ktdexp, &
+               &                imax_values, values, klevve, &
+               &                19, 3, zspd )
+
+         ELSE
+            
+            klevve = 0
+
+         ENDIF
+
+
+      CASE(132)
+         
+         ! BATHY
+
+         kyea = values(2)
+         kmon = values(3)
+         kday = values(4)
+         khou = values(5)
+         kmin = values(6)
+         plat = values(7)
+         plon = values(8)
+         
+         csign = cvals(1)
+         
+         klevts = INT( values(10) )
+
+         CALL bufr_get_data( 7062, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 11, 2, pdts )
+         
+         CALL bufr_get_data( 22042, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 12, 2, ptem )
+         
+         psal=pvalnull
+
+         klevve = 0
+
+      CASE(133)
+         
+         ! TESAC
+
+         kyea = values(2)
+         kmon = values(3)
+         kday = values(4)
+         khou = values(5)
+         kmin = values(6)
+         plat = values(7)
+         plon = values(8)
+         
+         csign = cvals(1)
+         
+         klevts = INT( values(11) )
+         
+         CALL bufr_get_data( 7062, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 12, 3, pdts )
+         
+         CALL bufr_get_data( 22043, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 13, 3, ptem )
+         
+         CALL bufr_get_data( 22062, imax_elements, ktdexp, &
+            &                imax_values, values, klevts, 14, 3, psal )
+
+         icur = 14 + 3 * (klevts - 1 ) + 1
+         
+         IF ( ktdexp(icur) == 2031 ) THEN
+
+            klevve = INT( values(icur+1) )
+
+            CALL bufr_get_data( 7062, imax_elements, ktdexp, &
+               &                imax_values, values, klevve, &
+               &                icur + 2, 3, pdve )
+            
+            CALL bufr_get_data( 22004, imax_elements, ktdexp, &
+               &                imax_values, values, klevve, &
+               &                icur + 3, 3, zdir )
+            
+            CALL bufr_get_data( 22031, imax_elements, ktdexp, &
+               &                imax_values, values, klevve, &
+               &                icur + 4, 3, zspd )
+
+         ELSE
+
+            klevve = 0
+
+         ENDIF
+         
+      CASE DEFAULT
+         
+         WRITE(*,*)'ERROR: SUBTYPE ',ksty, &
+            &           ' is not recognised as an Ocean type'
+         CALL abort
+         
+      END SELECT
+      
+      ! Convert data units
+
+      DO ji = 1, klevts
+         ! Convert temperature from Kelvin to Celcius
+         IF (( ptem(ji) > 250. ) .AND. ( ptem(ji) < 323. )) THEN
+            ptem(ji) = ptem(ji) - zero_celcius
+         ELSE
+            ptem(ji) = pvalnull
+         ENDIF
+         ! Check for physical salinities
+         IF (( psal(ji) < 0.0 ) .OR. ( psal(ji) > 50.0 )) THEN
+            psal(ji) = pvalnull
+         ENDIF
+      END DO
+
+      DO ji = 1, klevve
+         ! Convert speed and direction to u,v
+         IF (( zspd(ji) < 100 ).AND.( zdir(ji) < 361. )) THEN
+            CALL spd_to_uv( zdir(ji), zspd(ji), puve(ji), pvve(ji) )
+         ELSE
+            puve(ji) = pvalnull
+            pvve(ji) = pvalnull
+         ENDIF
+      END DO
+
+   END SUBROUTINE decode_bufr
+
+
+   SUBROUTINE bufr_get_data( ktdval, ktdlen, ktdexp, &
+      &                      kmax_values, zvalues, klevs, &
+      &                      istart, istep, zovalues )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE bufr_get_data **
+      !!
+      !!  ** Purpose : Extract data from decode buffer message
+      !!
+      !!  ** Method  : 
+      !!
+      !!  ** Action  :
+      !!
+      !!   History :   (??-??) A. Vidard  . ODASYS version. 
+      !!               (08-12) K. Mogensen. NEMOVAR version
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT(IN) :: ktdval                ! Expected descriptor
+      INTEGER, INTENT(IN) :: ktdlen                ! Size of descriptors
+      INTEGER, INTENT(IN), DIMENSION(ktdlen) :: &
+         & ktdexp                                  ! Descriptors
+      INTEGER, INTENT(IN) :: kmax_values           ! Size of input array
+      REAL, INTENT(IN) :: zvalues(kmax_values)     ! Decoded values
+      INTEGER, INTENT(IN) :: klevs                 ! Number of levels 
+      INTEGER, INTENT(IN) :: istart                ! Start in input array
+      INTEGER, INTENT(IN) :: istep                 ! Step between values
+      REAL, INTENT(OUT) :: zovalues(klevs)         ! Output values
+      !! * Local variables
+      INTEGER :: ji, ipos
+
+      DO ji = 1, klevs
+         ipos = istart + (ji-1) * istep
+         IF ( ktdexp(ipos) /= ktdval ) THEN
+            WRITE(*,*)'Problem decoding bufr data in bufr_get_data'
+            WRITE(*,*)'Expected descriptor : ', ktdval
+            WRITE(*,*)'Found descriptor : ', ktdexp(ipos)
+            CALL abort
+         ENDIF
+         zovalues(ji) = zvalues( ipos )
+      END DO
+
+   END SUBROUTINE bufr_get_data
+
+   SUBROUTINE spd_to_uv( pdir, pspd, pu, pv )
+      !!----------------------------------------------------------------------
+      !!   
+      !!                    *** ROUTINE spd_to_uv ***
+      !! 
+      !! ** Purpose : Convert ocean current direction and speed to uv.
+      !! 
+      !! ** Method  : According to WMO current directors are coded according
+      !!              to normal oceanographic convention.
+      !!
+      !! ** Action  :
+      !!
+      !! ** References:
+      !!     http://www.wmo.int/pages/prog/www/WMOCodes/OperationalCodes.html
+      !!
+      !! History
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      !! * Arguments
+      REAL, INTENT(IN) :: &
+         & pdir, &       ! Current direction
+         & pspd          ! Current speed
+      REAL, INTENT(OUT) :: &
+         & pu, &         ! Zonal corrent
+         & pv            ! Meridonal current
+      !! * Local variables
+      REAL(KIND=dp), PARAMETER ::  rad = 3.141592653589793_dp / 180.0_dp
+      REAL(KIND=dp) :: zdir
+      
+      zdir = pdir * rad
+
+      pu = pspd * SIN( zdir )
+      pv = pspd * COS( zdir )
+
+   END SUBROUTINE spd_to_uv
+
+#include "ctl_stop.h90"
+
+#include "greg2jul.h90"
+
+#endif
+
+END MODULE bufrdata
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/convmerge.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/convmerge.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/convmerge.F90	(revision 2893)
@@ -0,0 +1,126 @@
+MODULE convmerge
+
+   USE toolspar_kind
+   USE obs_fbm
+   IMPLICIT NONE
+
+CONTAINS
+
+   SUBROUTINE conv_fbmerge( cdoutfile, nfiles, fbdata )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE conv_fbmerg **
+      !!
+      !!  ** Purpose : Merge all fbfiles into a single fbfile
+      !!
+      !!  ** Method  : Use of utilities from obs_fbm.
+      !!
+      !!  ** Action  :
+      !!
+      !!   History :
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(LEN=*) :: cdoutfile                ! Input file.
+      INTEGER :: nfiles                            ! Number of files
+      TYPE(obfbdata), dimension(nfiles) :: fbdata  ! Structure to merge
+      !! * Local variables
+      type(obfbdata) :: fbmerge
+      INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
+      INTEGER :: nmaxlev
+      INTEGER :: ia,ij,ii
+      REAL(fbdp), DIMENSION(nfiles) :: djulini, djulend
+      CHARACTER(len=8) :: cl_refdate
+      INTEGER :: irefdate,iyea,imon,iday,ihou,imin,isec
+      ! Namelist variables
+      CHARACTER(len=9) :: cdnamefile='namobs.in'
+      LOGICAL :: lexists 
+      REAL(fbdp) :: dobsini,dobsend
+      NAMELIST/namobs/dobsini,dobsend
+      
+      dobsini = 0.0
+      dobsend = 99991231.235959
+
+      INQUIRE(file=cdnamefile, exist=lexists)
+      IF (lexists) THEN
+         OPEN(10,file=cdnamefile)
+         READ(10,namobs)
+         WRITE(*,namobs)
+         CLOSE(10)
+      ENDIF
+      !
+      ! Count number of data points
+      !
+      nmaxlev = 1
+      ii = 0
+      DO ia = 1,nfiles
+         IF (lexists) THEN
+            cl_refdate=fbdata(ia)%cdjuldref(1:8)
+            READ(cl_refdate,'(I8)') irefdate
+            CALL ddatetoymdhms( dobsini, iyea, imon, iday, ihou, imin, isec )
+            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(ia), &
+               &           krefdate = irefdate )
+            CALL ddatetoymdhms( dobsend, iyea, imon, iday, ihou, imin, isec )
+            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(ia), &
+            &           krefdate = irefdate )
+            DO ij = 1, fbdata(ia)%nobs
+               IF ( ( fbdata(ia)%ptim(ij) >  djulini(ia) ) .AND. &
+                  & ( fbdata(ia)%ptim(ij) <= djulend(ia) ) ) THEN
+                  ii = ii + 1
+                  nmaxlev = MAX(nmaxlev,fbdata(ia)%nlev)
+               ENDIF
+            ENDDO
+         ELSE
+            nmaxlev = MAX(nmaxlev,fbdata(ia)%nlev)
+            ii = ii + fbdata(ia)%nobs
+         ENDIF
+      ENDDO
+      !
+      ! Merge the input structures into the output structure
+      !
+      ALLOCATE( iset(ii), inum(ii), iindex(ii))
+      ii = 0
+      DO ia = 1,nfiles
+         DO ij = 1, fbdata(ia)%nobs
+            IF (lexists) THEN
+               IF ( ( fbdata(ia)%ptim(ij) >  djulini(ia) ) .AND. &
+                  & ( fbdata(ia)%ptim(ij) <= djulend(ia) ) ) THEN
+                  ii = ii + 1
+                  iset(ii)   = ia
+                  inum(ii)   = ij
+                  iindex(ii) = ii
+               ENDIF
+            ELSE
+               ii = ii + 1
+               iset(ii)   = ia
+               inum(ii)   = ij
+               iindex(ii) = ii
+            ENDIF
+         ENDDO
+      ENDDO
+      WRITE(*,*)'Output number of observations = ',ii
+      WRITE(*,*)'Output number of levels       = ',nmaxlev
+      !
+      ! Prepare fbmerge structure
+      !
+      CALL init_obfbdata( fbmerge )
+      CALL alloc_obfbdata( fbmerge, fbdata(1)%nvar, ii, nmaxlev, &
+         &                 fbdata(1)%nadd, fbdata(1)%next, fbdata(1)%lgrid )
+      CALL merge_obfbdata( nfiles, fbdata, fbmerge, iset, inum, iindex )
+      !
+      ! Write the file
+      !
+      CALL write_obfbdata( TRIM(cdoutfile), fbmerge )
+      !
+      ! Dellocate the data
+      !
+      CALL dealloc_obfbdata( fbmerge )
+
+   END SUBROUTINE conv_fbmerge
+
+#include "ctl_stop.h90"
+
+#include "greg2jul.h90"
+
+#include "ddatetoymdhms.h90"
+
+END MODULE convmerge
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/corio2fb.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/corio2fb.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/corio2fb.F90	(revision 2893)
@@ -0,0 +1,74 @@
+PROGRAM corio2fb
+   USE obs_fbm
+   USE obs_prof_io
+   USE convmerge
+   IMPLICIT NONE
+   !
+   ! Command line arguments for output file and input files
+   !
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdoutfile
+   CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
+   !
+   ! Input data
+   !
+   TYPE(obfbdata), POINTER :: coriof(:)
+   INTEGER :: ninfiles,ntotcorio,nmaxlev
+   INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
+   !
+   ! Output data
+   !
+   TYPE(obfbdata) :: fbdata
+   !
+   ! Loop variables
+   !
+   INTEGER :: ia,ii,ij
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF (nargs < 1) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'corio2fb outputfile inputfile1 inputfile2 ...'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdoutfile)
+   !
+   ! Get input data
+   !
+   ALLOCATE( coriof(MAX(nargs-1,1)) )
+   ALLOCATE( cdinfile(nargs-1) )
+   ntotcorio = 0
+   ninfiles  = nargs - 1
+   DO ia=1,ninfiles
+      CALL getarg( ia + 1, cdinfile(ia) )
+      CALL read_coriofile( TRIM(cdinfile(ia)), coriof(ia), 6, .TRUE., .FALSE. )
+      WRITE(*,'(2A)')'File = ',TRIM(cdinfile(ia))
+      WRITE(*,'(A,I9,A)')'has',coriof(ia)%nobs,' profiles'
+      ntotcorio = ntotcorio + coriof(ia)%nobs
+      nmaxlev   = MAX( nmaxlev, coriof(ia)%nlev )
+   ENDDO
+   IF (ninfiles==0) THEN
+      CALL init_obfbdata( coriof(1) )
+      CALL alloc_obfbdata( coriof(1), 2, 0, 1, 0, 1, .FALSE. )
+      coriof(1)%cname(1) = 'POTM'
+      coriof(1)%cname(2) = 'PSAL'
+      coriof(1)%coblong(1) = 'Potential temperature'
+      coriof(1)%coblong(2) = 'Practical salinity'
+      coriof(1)%cobunit(1) = 'Degrees Celsius'
+      coriof(1)%cobunit(2) = 'PSU'
+      coriof(1)%cextname(1) = 'TEMP'
+      coriof(1)%cextlong(1) = 'Insitu temperature'
+      coriof(1)%cextunit(1) = 'Degrees Celsius'
+      coriof(1)%cdjuldref = '19500101000000'
+   ENDIF
+   WRITE(*,'(A,I8)') 'Total profiles : ',ntotcorio
+   !
+   ! Merge and output the data.
+   !
+   CALL conv_fbmerge( TRIM(cdoutfile), ninfiles, coriof )
+   
+END PROGRAM corio2fb
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/ctl_stop.h90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/ctl_stop.h90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/ctl_stop.h90	(revision 2893)
@@ -0,0 +1,21 @@
+   SUBROUTINE ctl_stop( cerr )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE ctl_stop **
+      !!
+      !!  ** Purpose : Stop execution
+      !!
+      !!  ** Method  : 
+      !!
+      !!  ** Action  :
+      !!
+      !!   History :   (08-12) K. Mogensen. NEMOVAR version
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) :: cerr
+
+      WRITE(*,*)'Error:'
+      WRITE(*,*)TRIM(cerr)
+      CALL abort
+
+   END SUBROUTINE ctl_stop
+
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/enact2fb.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/enact2fb.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/enact2fb.F90	(revision 2893)
@@ -0,0 +1,74 @@
+PROGRAM enact2fb
+   USE obs_fbm
+   USE obs_prof_io
+   USE convmerge
+   IMPLICIT NONE
+   !
+   ! Command line arguments for output file and input files
+   !
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdoutfile
+   CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
+   !
+   ! Input data
+   !
+   TYPE(obfbdata), POINTER :: enactf(:)
+   INTEGER :: ninfiles,ntotenact,nmaxlev
+   INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
+   !
+   ! Output data
+   !
+   TYPE(obfbdata) :: fbdata
+   !
+   ! Loop variables
+   !
+   INTEGER :: ia,ii,ij
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF (nargs < 1) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'enact2fb outputfile inputfile1 inputfile2 ...'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdoutfile)
+   !
+   ! Get input data
+   !
+   ALLOCATE( enactf(MAX(nargs-1,1)) )
+   ALLOCATE( cdinfile(nargs-1) )
+   ntotenact = 0
+   ninfiles  = nargs - 1
+   DO ia = 1,ninfiles
+      CALL getarg( ia + 1, cdinfile(ia) )
+      CALL read_enactfile( TRIM(cdinfile(ia)), enactf(ia), 6, .TRUE., .FALSE. )
+      WRITE(*,'(2A)')'File = ',TRIM(cdinfile(ia))
+      WRITE(*,'(A,I9,A)')'has',enactf(ia)%nobs,' profiles'
+      ntotenact = ntotenact + enactf(ia)%nobs
+      nmaxlev   = MAX( nmaxlev, enactf(ia)%nlev )
+   ENDDO
+   IF (ninfiles==0) THEN
+      CALL init_obfbdata( enactf(1) )
+      CALL alloc_obfbdata( enactf(1), 2, 0, 1, 0, 1, .FALSE. )
+      enactf(1)%cname(1) = 'POTM'
+      enactf(1)%cname(2) = 'PSAL'
+      enactf(1)%coblong(1) = 'Potential temperature'
+      enactf(1)%coblong(2) = 'Practical salinity'
+      enactf(1)%cobunit(1) = 'Degrees Celsius'
+      enactf(1)%cobunit(2) = 'PSU'
+      enactf(1)%cextname(1) = 'TEMP'
+      enactf(1)%cextlong(1) = 'Insitu temperature'
+      enactf(1)%cextunit(1) = 'Degrees Celsius'
+      enactf(1)%cdjuldref = '19500101000000'
+   ENDIF
+   WRITE(*,'(A,I8)') 'Total profiles : ',ntotenact
+   !
+   ! Merge and output the data.
+   !
+   CALL conv_fbmerge( TRIM(cdoutfile), ninfiles, enactf )
+
+END PROGRAM enact2fb
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbaccdata.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbaccdata.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbaccdata.F90	(revision 2893)
@@ -0,0 +1,14 @@
+MODULE fbaccdata
+   INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: inum
+   REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: zmean,zrms,zsdev
+   INTEGER, PARAMETER :: maxvars = 10
+   REAL, DIMENSION(maxvars) :: zcheck
+   REAL, DIMENSION(maxvars) :: zhistmax, zhistmin, zhiststep
+   TYPE histtype
+      INTEGER :: npoints
+      INTEGER, POINTER, DIMENSION(:,:,:,:) :: nhist
+   END TYPE histtype
+   TYPE(histtype), DIMENSION(maxvars) :: hist
+
+END MODULE fbaccdata
+
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbcomb.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbcomb.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbcomb.F90	(revision 2893)
@@ -0,0 +1,157 @@
+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
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbcp.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbcp.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbcp.F90	(revision 2893)
@@ -0,0 +1,95 @@
+PROGRAM fbcp
+   USE obs_fbm
+   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) :: cdinfile
+   CHARACTER(len=256) :: cdtmp
+   INTEGER :: nout,ninn,nadd,next,i,j,k
+   LOGICAL :: lgrid
+   !
+   ! Feedback data
+   !
+   TYPE(obfbdata) :: fbdatain,fbdataout
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF ((nargs < 2)) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'fbcp inputfile outputfile {#outer] [include inner (0/1)] <names of extra fields>'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdinfile)
+   CALL getarg(2,cdoutfile)
+   nout=-1
+   ninn=0
+   IF ( nargs > 2 ) THEN
+      CALL getarg(3,cdtmp)
+      READ(cdtmp,'(I8)')nout
+      IF ( nargs > 3 ) THEN
+         CALL getarg(4,cdtmp)
+         READ(cdtmp,'(I8)')ninn
+         IF ((ninn<0).OR.(ninn>1)) THEN
+            WRITE(*,*)'Inner loop switch should be either 0 or 1'
+            CALL abort
+         ENDIF
+      ENDIF
+   ENDIF
+   !
+   ! Initialize feedback data
+   !
+   CALL init_obfbdata( fbdatain )
+   CALL init_obfbdata( fbdataout )
+   !
+   ! Read the file
+   !
+   CALL read_obfbdata(  TRIM(cdinfile), fbdatain )
+   !
+   ! Copy the file
+   !
+   IF ( nout >= 0 ) THEN
+      nadd = nout * ( 1 + ninn) + 1 + fbdatain%nadd
+   ELSE
+      nadd = fbdatain%nadd
+   ENDIF
+   IF ( nargs > 4 ) THEN
+      next = fbdatain%next + ( nargs - 4)
+   ELSE
+      next = fbdatain%next
+   ENDIF
+   CALL copy_obfbdata( fbdatain, fbdataout, kadd = nadd, kext = next, &
+      &                lgrid = (nout >= 0) .OR. fbdatain%lgrid )
+   !
+   ! Set additional names
+   !
+   IF ( nout>=0 ) THEN
+      k=fbdatain%nadd
+      DO j=0,nout
+         k=k+1
+         WRITE(fbdataout%caddname(k),'(A,I2.2)')'Hx',j
+         IF (j>0) THEN
+            DO i=1,ninn
+               k=k+1
+               WRITE(fbdataout%caddname(k),'(A,I2.2)')'Hxa',j
+            ENDDO
+         ENDIF
+      ENDDO
+   ENDIF
+   IF ( nargs > 4 ) THEN
+      DO i = 1, nargs - 4
+         CALL getarg(i+4,fbdataout%cextname(i+fbdatain%next))
+      ENDDO
+   ENDIF
+   !
+   ! Write the file
+   !
+   CALL write_obfbdata(  TRIM(cdoutfile), fbdataout )
+
+END PROGRAM fbcp
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbmatchup.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbmatchup.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbmatchup.F90	(revision 2893)
@@ -0,0 +1,394 @@
+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
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbprint.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbprint.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbprint.F90	(revision 2893)
@@ -0,0 +1,330 @@
+PROGRAM fbprint
+   USE toolspar_kind 
+   USE obs_fbm
+   USE index_sort
+   IMPLICIT NONE
+   !
+   ! Command line arguments input file
+   !
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdinfile,cdbrief
+   LOGICAL :: lbrief, lqcflags, lstat
+   CHARACTER(len=8) :: cdstat
+   INTEGER :: nqc
+   !
+   ! Input data
+   !
+   TYPE(obfbdata) :: obsdata
+   !
+   ! Loop variables
+   !
+   INTEGER :: ii
+   !
+   ! Sorting
+   !
+   REAL(KIND=dp),ALLOCATABLE :: zsort(:,:)
+   INTEGER,ALLOCATABLE  :: iindex(:)
+   !
+   ! Get number of command line arguments
+   !
+   nargs  = IARGC()
+   lbrief = .FALSE.
+   lstat = .FALSE.
+   cdstat = 'XXXXXXX'
+   nqc = 0
+   IF ( (nargs < 1) .OR. (nargs > 3) ) THEN
+      CALL usage()
+   ENDIF
+   IF (( nargs == 2 )) THEN
+      CALL getarg( 1, cdbrief )
+      IF ( cdbrief == '-b' ) THEN
+         lbrief = .TRUE.
+      ELSEIF( cdbrief == '-q' ) THEN
+         lqcflags = .TRUE.
+         nqc=1
+      ELSEIF( cdbrief == '-Q' ) THEN
+         lqcflags = .TRUE.
+         nqc=2
+      ELSEIF( cdbrief == '-B' ) THEN
+         lqcflags = .TRUE.
+         nqc=3
+      ELSE
+         CALL usage
+      ENDIF
+   ENDIF
+   IF (( nargs == 3 )) THEN
+      CALL getarg( 1, cdbrief )
+      IF ( cdbrief == '-s' ) THEN
+         lstat = .TRUE.
+         CALL getarg( 2, cdstat )
+      ELSE
+         CALL usage
+      ENDIF
+   ENDIF
+   CALL getarg( nargs, cdinfile )
+   !
+   ! Get input data
+   !
+   CALL read_obfbdata( TRIM(cdinfile), obsdata )
+   WRITE(*,'(2A,I9,A,I9,A)')TRIM(cdinfile), ' has ', obsdata%nobs ,&
+      & ' observations and a maximum of ', obsdata%nlev, ' levels'
+   !
+   ! Sort the data
+   !   
+   ALLOCATE(zsort(3,obsdata%nobs),iindex(obsdata%nobs))
+   DO ii=1,obsdata%nobs
+      zsort(1,ii)=obsdata%ptim(ii)
+      zsort(2,ii)=obsdata%pphi(ii)
+      zsort(3,ii)=obsdata%plam(ii)
+   ENDDO
+   CALL index_sort_dp_n(zsort,3,iindex,obsdata%nobs)
+   !
+   ! Print the sorted list
+   !   
+   DO ii=1,obsdata%nobs
+      IF (lstat) THEN
+         IF (cdstat /= obsdata%cdwmo(ii)) CYCLE
+      ENDIF
+      IF (lqcflags) THEN
+         CALL print_obs_qc(obsdata,iindex(ii),nqc)
+      ELSE
+         CALL print_obs(obsdata,iindex(ii),lbrief,lqcflags,nqc)
+      ENDIF
+   ENDDO
+
+END PROGRAM fbprint
+
+SUBROUTINE usage
+   WRITE(*,'(A)')'Usage:'
+   WRITE(*,'(A)')'fbprint [-b] [-q] inputfile'
+   WRITE(*,'(A)')'where -b selects brief output'
+   WRITE(*,'(A)')'  and -q selects qc flags rather than extra fields'
+   CALL abort()
+END SUBROUTINE usage
+
+SUBROUTINE print_obs(obsdata,iindex,lshort)
+   USE obs_fbm
+   USE date_utils
+   IMPLICIT NONE
+   TYPE(obfbdata) :: obsdata
+   INTEGER :: iindex
+   LOGICAL :: lshort
+   INTEGER :: jv,ja,je,jk
+   INTEGER :: kj,iyr,imon,iday,ihou,imin,isec
+   LOGICAL :: lskip
+   CHARACTER(len=1024) :: cdfmt1,cdfmt2
+   CHARACTER(len=16) :: cdtmp
+
+   WRITE(*,*)'Fileindex           = ',obsdata%kindex(iindex)
+   WRITE(*,*)'Station identifier  = ',obsdata%cdwmo(iindex)
+   WRITE(*,*)'Station type        = ',obsdata%cdtyp(iindex)
+   WRITE(*,*)'Latitude            = ',obsdata%pphi(iindex)
+   WRITE(*,*)'Longtude            = ',obsdata%plam(iindex)
+   WRITE(*,*)'Position QC         = ',obsdata%ipqc(iindex)
+   WRITE(*,*)'Observation QC      = ',obsdata%ioqc(iindex)
+   WRITE(*,*)'Julian date         = ',obsdata%ptim(iindex)
+   CALL jul2greg(isec,imin,ihou,iday,imon,iyr,obsdata%ptim(iindex))
+   WRITE(*,'(1X,A,I4,2I2.2)') &
+      &      'Gregorian date      = ',iyr,imon,iday
+   WRITE(*,'(1X,A,I2.2,A1,I2.2,A1,I2.2)') &
+      &      'Time                = ',ihou,':',imin,':',isec
+   IF (.NOT.lshort) THEN
+      DO jv = 1,obsdata%nvar
+         WRITE(*,*)'Variable name       = ',obsdata%cname(jv)
+         WRITE(*,*)'Variable QC         = ',obsdata%ivqc(iindex,jv)
+         IF (obsdata%lgrid) THEN
+            WRITE(*,*)'Grid I              = ',obsdata%iobsi(iindex,jv)
+            WRITE(*,*)'Grid J              = ',obsdata%iobsj(iindex,jv)
+         ENDIF
+      ENDDO
+      cdfmt1='(1X,A8,1X,A8'
+      cdfmt2='(1X,F8.2,1X,I8'
+      DO jv=1, obsdata%nvar
+         cdfmt1 = TRIM(cdfmt1)//',1X,A15,1X,A8'
+         cdfmt2 = TRIM(cdfmt2)//',1X,E15.9,1X,I8'
+         IF (obsdata%nadd>0) THEN
+            WRITE(cdtmp,'(I10)')obsdata%nadd
+            cdfmt1 = TRIM(cdfmt1)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(1X,A15)'
+            cdfmt2 = TRIM(cdfmt2)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(1X,E15.9)'
+         ENDIF
+         IF (obsdata%lgrid) THEN
+            cdfmt1 = TRIM(cdfmt1)//',1X,A10'
+            cdfmt2 = TRIM(cdfmt2)//',1X,I10'
+         ENDIF
+      ENDDO
+      IF (obsdata%next>0) THEN
+         WRITE(cdtmp,'(I10)')obsdata%next
+         cdfmt1 = TRIM(cdfmt1)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(1X,A15)'
+         cdfmt2 = TRIM(cdfmt2)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(1X,E15.9)'
+      ENDIF
+      cdfmt1=TRIM(cdfmt1)//')'
+      cdfmt2=TRIM(cdfmt2)//')'
+      IF (obsdata%lgrid) THEN
+         WRITE(*,FMT=cdfmt1)&
+            & 'DEPTH', 'DEP_QC', &
+            & (TRIM(obsdata%cname(jv))//'_OBS', &
+            & TRIM(obsdata%cname(jv))//'_QC' , &
+            & (TRIM(obsdata%cname(jv))//'_'//TRIM(obsdata%caddname(ja)),&
+            & ja = 1, obsdata%nadd ), &
+            & TRIM(obsdata%cname(jv))//'_K' , &
+            & jv = 1, obsdata%nvar ), &
+            & ( TRIM(obsdata%cextname(ja)),&
+            & ja = 1, obsdata%next )
+         DO kj=1,obsdata%nlev
+            IF (obsdata%pdep(kj,iindex)<99999.0) THEN
+               WRITE (*,FMT=cdfmt2) &
+                  & obsdata%pdep(kj,iindex),   &
+                  & obsdata%idqc(kj,iindex),   &
+                  & ( obsdata%pob(kj,iindex,jv), obsdata%ivlqc(kj,iindex,jv), &
+                  & ( obsdata%padd(kj,iindex,ja,jv) , ja=1, obsdata%nadd ), &
+                  & obsdata%iobsk(kj,iindex,jv), &
+                  & jv = 1, obsdata%nvar ), &
+                  & ( obsdata%pext(kj,iindex,ja), ja=1, obsdata%next )
+            ENDIF
+         ENDDO
+      ELSE
+         cdfmt1=TRIM(cdfmt1)//')'
+         cdfmt2=TRIM(cdfmt2)//')'
+         WRITE(*,FMT=cdfmt1)&
+            & 'DEPTH', 'DEP_QC', &
+            & (TRIM(obsdata%cname(jv))//'_OBS', &
+            & TRIM(obsdata%cname(jv))//'_QC' , &
+            & (TRIM(obsdata%cname(jv))//TRIM(obsdata%caddname(ja)),&
+            & ja = 1, obsdata%nadd ), &
+            & jv = 1, obsdata%nvar ), &
+            & ( TRIM(obsdata%cextname(ja)),&
+            & ja = 1, obsdata%next )
+         DO kj=1,obsdata%nlev
+            IF (obsdata%pdep(kj,iindex)<99999.0) THEN
+               WRITE (*,FMT=cdfmt2) &
+                  & obsdata%pdep(kj,iindex),   &
+                  & obsdata%idqc(kj,iindex),   &
+                  & ( obsdata%pob(kj,iindex,jv), obsdata%ivlqc(kj,iindex,jv), &
+                  & ( obsdata%padd(kj,iindex,ja,jv) , ja=1, obsdata%nadd ), &
+                  & jv = 1, obsdata%nvar ), &
+                     & ( obsdata%pext(kj,iindex,ja), ja=1, obsdata%next )
+            ENDIF
+         ENDDO
+      ENDIF
+   ENDIF
+   WRITE(*,*)
+END SUBROUTINE print_obs
+
+SUBROUTINE print_obs_qc(obsdata,iindex,kqc)
+   USE obs_fbm
+   USE date_utils
+   IMPLICIT NONE
+   TYPE(obfbdata) :: obsdata
+   INTEGER :: iindex
+   LOGICAL :: lqc
+   INTEGER :: kqc
+   INTEGER :: jv,ja,je,jk
+   INTEGER :: kj,iyr,imon,iday,ihou,imin,isec
+   LOGICAL :: lskip
+   CHARACTER(len=1024) :: cdfmt1,cdfmt2
+   CHARACTER(len=16) :: cdtmp
+   integer :: iqcf
+
+   IF (kqc==2) THEN
+      lskip=.TRUE.
+      IF (obsdata%ipqc(iindex)>1) lskip=.FALSE.
+      IF (obsdata%ioqc(iindex)>1) lskip=.FALSE.
+      DO jv = 1,obsdata%nvar
+         IF (obsdata%ivqc(iindex,jv)>1) lskip=.FALSE.
+      ENDDO
+      DO kj=1,obsdata%nlev
+         IF (obsdata%pdep(kj,iindex)<99999.0) THEN
+            IF (obsdata%idqc(kj,iindex)>1) lskip=.FALSE.
+            DO jv = 1, obsdata%nvar
+               IF (obsdata%ivlqc(kj,iindex,jv)>1) lskip=.FALSE.
+            ENDDO
+         ENDIF
+      ENDDO
+      IF (lskip) RETURN
+   ELSEIF (kqc==3) THEN
+      lskip=.TRUE.
+      DO kj=1,obsdata%nlev
+         IF (obsdata%pdep(kj,iindex)<99999.0) THEN
+            iqcf=0
+            DO jv = 1, obsdata%nvar
+               IF (obsdata%ivlqc(kj,iindex,jv)>1) iqcf=iqcf+1
+               IF (iqcf==obsdata%nvar) lskip=.FALSE.
+            ENDDO
+         ENDIF
+      ENDDO
+      IF (lskip) RETURN
+   ENDIF
+   WRITE(*,*)'Fileindex           = ',obsdata%kindex(iindex)
+   WRITE(*,*)'Station identifier  = ',obsdata%cdwmo(iindex)
+   WRITE(*,*)'Station type        = ',obsdata%cdtyp(iindex)
+   WRITE(*,*)'Latitude            = ',obsdata%pphi(iindex)
+   WRITE(*,*)'Longtude            = ',obsdata%plam(iindex)
+   WRITE(*,*)'Position QC         = ',obsdata%ipqc(iindex)
+   WRITE(*,*)'Position QC flags   = ',obsdata%ipqcf(:,iindex)
+   WRITE(*,*)'Observation QC      = ',obsdata%ioqc(iindex)
+   WRITE(*,*)'Observation QC flags= ',obsdata%ioqcf(:,iindex)
+   WRITE(*,*)'Julian date         = ',obsdata%ptim(iindex)
+   CALL jul2greg(isec,imin,ihou,iday,imon,iyr,obsdata%ptim(iindex))
+   WRITE(*,'(1X,A,I4,2I2.2)') &
+      &      'Gregorian date      = ',iyr,imon,iday
+   WRITE(*,'(1X,A,I2.2,A1,I2.2,A1,I2.2)') &
+      &      'Time                = ',ihou,':',imin,':',isec
+   DO jv = 1,obsdata%nvar
+      WRITE(*,*)'Variable name       = ',obsdata%cname(jv)
+      WRITE(*,*)'Variable QC         = ',obsdata%ivqc(iindex,jv)
+      WRITE(*,*)'Variable QC flags   = ',obsdata%ivqcf(:,iindex,jv)
+   ENDDO
+   cdfmt1='(1X,A8,1X,A8'
+   cdfmt2='(1X,F8.2,1X,I8'
+   WRITE(cdtmp,'(I10)')obsdata%nqcf
+   cdfmt1 = TRIM(cdfmt1)//',1X,A18'
+   cdfmt2 = TRIM(cdfmt2)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(I9)'
+   DO jv=1, obsdata%nvar
+      cdfmt1 = TRIM(cdfmt1)//',1X,A15,1X,A8'
+      cdfmt2 = TRIM(cdfmt2)//',1X,E15.9,1X,I8'
+      WRITE(cdtmp,'(I10)')obsdata%nqcf
+      cdfmt1 = TRIM(cdfmt1)//',1X,A18'
+      cdfmt2 = TRIM(cdfmt2)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(I9)'
+   ENDDO
+   IF (obsdata%next>0) THEN
+      WRITE(cdtmp,'(I10)')obsdata%next
+      cdfmt1 = TRIM(cdfmt1)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(1X,A15)'
+      cdfmt2 = TRIM(cdfmt2)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(1X,E15.9)'
+   ENDIF
+   cdfmt1=TRIM(cdfmt1)//')'
+   cdfmt2=TRIM(cdfmt2)//')'
+   WRITE(*,FMT=cdfmt1)&
+      & 'DEPTH', 'DEP_QC', 'DEP_QC_FLAGS', &
+      & (TRIM(obsdata%cname(jv))//'_OBS', &
+      & TRIM(obsdata%cname(jv))//'_QC' , &
+      & TRIM(obsdata%cname(jv))//'_QC_FLAGS',&
+      & jv = 1, obsdata%nvar ), &
+      & ( TRIM(obsdata%cextname(ja)),&
+      & ja = 1, obsdata%next )
+   DO kj=1,obsdata%nlev
+      IF (kqc>=2)  THEN
+         lskip=.TRUE.
+         IF (obsdata%idqc(kj,iindex)>1) lskip=.FALSE.
+         DO jv = 1, obsdata%nvar
+            IF (obsdata%ivlqc(kj,iindex,jv)>1) lskip=.FALSE.
+         ENDDO
+         IF (lskip) CYCLE
+      ENDIF
+      IF (obsdata%pdep(kj,iindex)<99999.0) THEN
+         WRITE (*,FMT=cdfmt2) &
+            & obsdata%pdep(kj,iindex),   &
+            & obsdata%idqc(kj,iindex),   &
+            & ( obsdata%idqcf(ja,kj,iindex), ja = 1, obsdata%nqcf ), &
+            & ( obsdata%pob(kj,iindex,jv), obsdata%ivlqc(kj,iindex,jv), &
+            & ( obsdata%ivlqcf(ja,kj,iindex,jv) , ja=1, obsdata%nqcf ), &
+            & jv = 1, obsdata%nvar ), &
+            & ( obsdata%pext(kj,iindex,ja), ja=1, obsdata%next )
+      ENDIF
+   ENDDO
+   WRITE(*,*)
+
+END SUBROUTINE print_obs_qc
+   
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbsel.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbsel.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbsel.F90	(revision 2893)
@@ -0,0 +1,172 @@
+PROGRAM fbsel
+   USE obs_fbm
+   USE date_utils
+   IMPLICIT NONE
+   TYPE(obfbdata) :: fbdatain,fbdataout
+   CHARACTER(len=256) :: filenamein,filenameout,filenametmp,cnameout
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER,PARAMETER :: maxtyp=1023
+   INTEGER,PARAMETER :: maxdates=20
+   INTEGER :: nqc,ntyp,ndates,ninidate(maxdates),nenddate(maxdates)
+   LOGICAL :: lsplitqc,lsplittyp
+   INTEGER :: iqc,ityp,idate
+   REAL :: maxlat,minlat,maxlon,minlon
+   NAMELIST/namsel/nqc,ntyp,ndates,ninidate,nenddate,lsplitqc,lsplittyp, &
+      &            maxlat,minlat,maxlon,minlon
+
+   IF (iargc()/=2) THEN
+      WRITE(*,*)'Usage:'
+      WRITE(*,*)'fbsel <input filename> <output filename>'
+      CALL abort
+   ENDIF
+
+   CALL getarg(1,filenamein)
+   CALL getarg(2,filenameout)
+   
+   nqc=-1
+   ntyp=-1
+   ndates=1
+   ninidate=19500101
+   nenddate=21000101
+
+   lsplitqc=.FALSE.
+   lsplittyp=.FALSE.
+   maxlat=1e+38
+   minlat=-1e+38
+   maxlon=1e+38
+   minlon=-1e+38
+   OPEN(10,file='namsel.in')
+   READ(10,namsel)
+   CLOSE(10)
+   WRITE(*,namsel)
+
+   CALL init_obfbdata(fbdatain)
+   CALL init_obfbdata(fbdataout)
+
+   WRITE(*,*)'Reading file : ',TRIM(filenamein)
+   CALL read_obfbdata(TRIM(filenamein),fbdatain)
+   WRITE(*,*)'Number of observations before selection = ',fbdatain%nobs
+   DO idate=1,ndates
+      IF (ndates==1) THEN
+         cnameout=filenameout
+      ELSE
+         WRITE(cnameout,'(I2.2,2A)')idate,'_',TRIM(filenameout)
+      ENDIF
+      IF (lsplitqc) THEN
+         DO iqc=1,3
+            CALL fb_sel(fbdatain,fbdataout,iqc,ntyp, &
+               &        ninidate(idate),nenddate(idate), &
+               &        maxlat,minlat,maxlon,minlon)
+            WRITE(filenametmp,'(A,I2.2,A,A)')'qc_',iqc,'_',TRIM(cnameout)
+            IF (fbdataout%nobs>0) THEN
+               WRITE(*,*)'QC selected = ',iqc
+               WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
+               WRITE(*,*)'Writing file : ',TRIM(filenametmp)
+               CALL write_obfbdata(TRIM(filenametmp),fbdataout)
+            ENDIF
+            CALL dealloc_obfbdata(fbdataout)
+         ENDDO
+      ELSEIF (lsplittyp) THEN
+         DO ityp=0,maxtyp
+            CALL fb_sel(fbdatain,fbdataout,nqc,ityp, &
+               &        ninidate(idate),nenddate(idate), &
+               &        maxlat,minlat,maxlon,minlon)
+            WRITE(filenametmp,'(A,I4.4,A,A)')'typ_',ityp,'_',TRIM(cnameout)
+            IF (fbdataout%nobs>0) THEN
+               WRITE(*,*)'Type = ',ityp
+               WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
+               WRITE(*,*)'Writing file : ',TRIM(filenametmp)
+               CALL write_obfbdata(TRIM(filenametmp),fbdataout)
+            ENDIF
+            CALL dealloc_obfbdata(fbdataout)
+         ENDDO
+      ELSE
+         CALL fb_sel(fbdatain,fbdataout,nqc,ntyp, & 
+            &        ninidate(idate),nenddate(idate), &
+            &        maxlat,minlat,maxlon,minlon)
+         WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
+         WRITE(*,*)'Writing file : ',TRIM(cnameout)
+         CALL write_obfbdata(TRIM(cnameout),fbdataout)
+         CALL dealloc_obfbdata(fbdataout)
+      ENDIF
+   ENDDO
+
+CONTAINS
+
+   SUBROUTINE fb_sel(fbdatain,fbdataout,nqc,ntyp,ninidate,nenddate,&
+      &              maxlat,minlat,maxlon,minlon)
+      TYPE(obfbdata) :: fbdatain,fbdataout
+      INTEGER :: nqc,ntyp,ninidate,nenddate
+      REAL :: maxlat,minlat,maxlon,minlon
+      INTEGER :: jobs
+      INTEGER :: iqc,ityp
+      LOGICAL :: llvalid(fbdatain%nobs)
+      INTEGER :: iyea,imon,iday
+      REAL(KIND=dp) :: zjini,zjend
+
+      iyea=ninidate/10000
+      imon=ninidate/100-iyea*100
+      iday=ninidate-iyea*10000-imon*100
+      CALL greg2jul(0,0,0,iday,imon,iyea,zjini)
+      iyea=nenddate/10000
+      imon=nenddate/100-iyea*100
+      iday=nenddate-iyea*10000-imon*100
+      CALL greg2jul(0,0,0,iday,imon,iyea,zjend)
+      DO jobs = 1, fbdatain%nobs
+         llvalid(jobs)=.TRUE.
+         IF (nqc/=-1) THEN
+            CALL check_prof(fbdatain,jobs,iqc)
+            llvalid(jobs)=(iqc==nqc).AND.llvalid(jobs)
+         ENDIF
+         IF (ntyp/=-1) THEN
+            READ(fbdatain%cdtyp(jobs),'(I4)')ityp
+            llvalid(jobs)=(ityp==ntyp).AND.llvalid(jobs)
+         ENDIF
+         IF (ninidate/=-1) THEN
+            llvalid(jobs)=(fbdatain%ptim(jobs)>zjini).AND.llvalid(jobs)
+         ENDIF
+         IF (nenddate/=-1) THEN
+            llvalid(jobs)=(fbdatain%ptim(jobs)<=zjend).AND.llvalid(jobs)
+         ENDIF
+         llvalid(jobs)=(fbdatain%pphi(jobs)<=maxlat).AND. &
+            &          (fbdatain%pphi(jobs)>=minlat).AND. &
+            &          (fbdatain%plam(jobs)<=maxlon).AND. &
+            &          (fbdatain%plam(jobs)>=minlon).AND.llvalid(jobs)
+         ! Add more checks here...
+      ENDDO
+
+      CALL subsamp_obfbdata(fbdatain,fbdataout,llvalid)
+
+   END SUBROUTINE fb_sel
+
+   SUBROUTINE check_prof(fbdata,iobs,iqc)
+      
+      TYPE(obfbdata) :: fbdata
+      INTEGER :: iobs,iqc
+      INTEGER :: i,ivar
+      
+      LOGICAL :: lpart,lfull
+
+      lpart=.false.
+      lfull=.true.
+      DO ivar=1,fbdata%nvar
+         DO i=1,fbdata%nlev
+            IF ((fbdata%ivlqc(i,iobs,ivar)>2).AND.&
+               &(fbdata%ivlqc(i,iobs,ivar)<9)) lpart = .TRUE.
+            IF (fbdata%ivlqc(i,iobs,ivar)<=2) lfull = .FALSE.
+         ENDDO
+      ENDDO
+
+      IF(lfull) THEN
+         iqc=3
+      ELSEIF (lpart) then
+         iqc=2
+      ELSE
+         iqc=1
+      ENDIF
+
+   END SUBROUTINE check_prof
+
+END PROGRAM fbsel
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbstat.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbstat.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbstat.F90	(revision 2893)
@@ -0,0 +1,606 @@
+PROGRAM fbstat
+   USE obs_fbm
+   USE fbaccdata
+   USE coords
+   USE omonainfo
+   IMPLICIT NONE
+   TYPE(obfbdata) :: fbdata
+   CHARACTER(len=256) :: filename
+   INTEGER :: jfile,jbox,jlev,jfirst,jvar,jadd,ji
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   REAL,DIMENSION(:),ALLOCATABLE :: zlev
+   INTEGER :: nmlev, nfiles
+   LOGICAL :: lexists,lomona,ltext
+   REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zdat3d
+   REAL, ALLOCATABLE, DIMENSION(:,:) :: zdat2d
+   INTEGER,DIMENSION(1) :: itime
+   INTEGER :: inidate,icurdate,loopno,ityp,iloopno
+   INTEGER :: nvar,nadd
+   CHARACTER(len=4) :: expver
+   CHARACTER(len=7) :: cltyp
+   CHARACTER(len=128) :: cdfmthead,cdfmtbody
+   LOGICAL :: lnear,linner,linnerp,linnerini,lpassive,lhistogram,lfound
+   INTEGER :: nqc
+   CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: cname,caddname
+   CHARACTER(len=20) :: carea
+   INTEGER :: jstabox, jendbox
+   NAMELIST/namfbstat/ltext,lomona,nmlev,inidate,icurdate,loopno,&
+      &               expver,lnear,linner,lpassive,lhistogram, &
+      &               zhistmax,zhistmin,zhiststep,zcheck,carea
+
+   ltext=.TRUE.
+   lomona=.FALSE.
+   nmlev=31
+   inidate=19010101
+   icurdate=19010116
+   loopno=-1
+   expver='test'
+   lnear=.TRUE.
+   linner=.FALSE.
+   lpassive=.FALSE.
+   lhistogram=.FALSE.
+   zhistmin(:)=-10.0
+   zhistmax(:)=10.0
+   zhiststep(:)=0.1
+   zcheck(:)=10.0
+   nqc=2
+   carea='all'
+   INQUIRE(file='namfbstat.in',exist=lexists)
+   IF (lexists) THEN
+      OPEN(10,file='namfbstat.in')
+      READ(10,namfbstat)
+      CLOSE(10)
+      WRITE(*,namfbstat)
+   ENDIF
+   IF (iargc()==0) THEN
+      WRITE(*,*)'Usage:'
+      WRITE(*,*)'fbstat [-nmlev] <filenames>'
+      CALL abort
+   ENDIF
+   jfirst=1
+   DO ji=1,2
+      CALL getarg(jfirst,filename)
+      IF (filename=='-42') THEN
+         nmlev=42
+         jfirst=jfirst+1
+      ELSEIF(filename=='-31') THEN
+         nmlev=31
+         jfirst=jfirst+1
+      ELSEIF(filename=='-1') THEN
+         nmlev=1
+         lnear=.TRUE.
+         jfirst=jfirst+1
+      ELSEIF(filename=='-q') THEN
+         jfirst=jfirst+1
+         CALL getarg(jfirst,filename)
+         READ(filename,'(I4)')nqc
+         IF ((nqc<0).OR.(nqc>4)) THEN
+            WRITE(*,*)'Quality control option (-q) should be 1 to 4'
+            CALL abort
+         ENDIF
+         jfirst=jfirst+1
+      ENDIF
+   END DO
+   nfiles=iargc()
+   
+   IF (carea/='all') THEN
+      IF (lomona) THEN
+         WRITE(*,*)'For omona files carea has to be all'
+         CALL abort
+      ENDIF
+      lfound=.FALSE.
+      DO jbox=1,nbox
+         IF (TRIM(carea)==TRIM(cl_boxes(jbox))) THEN
+            jstabox=jbox
+            jendbox=jbox
+            lfound=.TRUE.
+         ENDIF
+      ENDDO
+      IF (.NOT.lfound) THEN
+         WRITE(*,*)'Area not found'
+         CALL abort
+      ENDIF
+   ELSE
+      jstabox=1
+      jendbox=nbox
+   ENDIF
+   PRINT *,'jstabox,jendbox=',jstabox,jendbox
+
+   IF (.NOT.lnear) nmlev=nmlev-1
+
+   ALLOCATE(&
+      & zlev(nmlev) &
+      & )
+   IF(lnear) THEN
+      CALL getlevs(nmlev,zlev)
+   ELSE
+      CALL getlevsmean(nmlev,zlev)
+   ENDIF
+
+   DO jfile=jfirst, nfiles
+      CALL getarg(jfile,filename)
+      WRITE(*,*)'Handling file : ',TRIM(filename)
+      CALL read_obfbdata(TRIM(filename),fbdata)
+      IF (jfile==jfirst) THEN
+         nvar=fbdata%nvar
+         nadd=fbdata%nadd
+         IF (lhistogram) THEN
+            IF (nvar>maxvars) THEN
+               WRITE(*,*)'fbstat.F90: Increase maxvars to ',nvar
+               WRITE(*,*)'if you want histograms'
+               CALL abort
+            ENDIF
+            DO jvar = 1, nvar
+               hist(jvar)%npoints=(zhistmax(jvar)-zhistmin(jvar))&
+                  &              /zhiststep(jvar)+1
+               WRITE(*,*)'Number of points in histogram      = ',&
+                  &  hist(jvar)%npoints
+               WRITE(*,*)'Size of histogram array (elements) = ',&
+                  &  hist(jvar)%npoints*nmlev*nbox*nadd
+               ALLOCATE(&
+                  & hist(jvar)%nhist(hist(jvar)%npoints,nmlev,nbox,nadd) &
+                  & )
+               hist(jvar)%nhist(:,:,:,:)=0
+            ENDDO
+         ENDIF
+         ALLOCATE(&
+            & inum(nmlev,nbox,nadd,nvar),  &
+            & zmean(nmlev,nbox,nadd,nvar), &
+            & zrms(nmlev,nbox,nadd,nvar),  &
+            & zsdev(nmlev,nbox,nadd,nvar), &
+            & cname(nvar),                 &
+            & caddname(nadd)               &
+            & )
+         DO jvar = 1, nvar
+            cname(jvar) = fbdata%cname(jvar)
+         END DO
+         DO jadd = 1, nadd
+            caddname(jadd) = fbdata%caddname(jadd)
+         END DO
+         inum(:,:,:,:)=0
+         zmean(:,:,:,:)=0.0
+         zrms(:,:,:,:)=0.0
+         zsdev(:,:,:,:)=0.0
+      ELSE
+         IF (fbdata%nvar/=nvar) THEN
+            WRITE(*,*)'Different number of nvar ',fbdata%nvar,' in ',&
+               & TRIM(filename)
+            CALL abort
+         ENDIF
+         IF (fbdata%nadd/=nadd) THEN
+            WRITE(*,*)'Different number of nadd ',fbdata%nadd,' in ',&
+               & TRIM(filename)
+            CALL abort
+         ENDIF
+      ENDIF
+      CALL fb_stat(fbdata,jstabox,jendbox,nmlev,zlev,lnear,nqc,lhistogram)
+      CALL dealloc_obfbdata(fbdata)
+   ENDDO
+
+   DO jvar=1, nvar
+      DO jadd=1, nadd
+         DO jbox=jstabox, jendbox
+            DO jlev = 1, nmlev
+               IF ( inum(jlev,jbox,jadd,jvar) > 0 ) THEN
+                  zmean(jlev,jbox,jadd,jvar) = &
+                     & zmean(jlev,jbox,jadd,jvar)/inum(jlev,jbox,jadd,jvar)
+                  zrms(jlev,jbox,jadd,jvar) = &
+                     & SQRT(zrms(jlev,jbox,jadd,jvar)/inum(jlev,jbox,jadd,jvar))
+                  zsdev(jlev,jbox,jadd,jvar) = &
+                     & SQRT(MAX(zrms(jlev,jbox,jadd,jvar)**2-zmean(jlev,jbox,jadd,jvar)**2,0.0))
+               ELSE
+                  zmean(jlev,jbox,jadd,jvar) = fbrmdi
+                  zrms(jlev,jbox,jadd,jvar) = fbrmdi
+                  zsdev(jlev,jbox,jadd,jvar) = fbrmdi
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+   ENDDO
+
+   IF (ltext) THEN
+      DO jvar=1, nvar
+         DO jadd=1, nadd
+            DO jbox=jstabox, jendbox
+               WRITE(filename,'(5A)')TRIM(cname(jvar)),&
+                  &                  TRIM(caddname(jadd)),'_',&
+                  &                  TRIM(cl_boxes(jbox)),'.dat'
+               OPEN(10,file=TRIM(filename))
+               DO jlev = 1, nmlev
+                  WRITE(10,'(F16.7,2I12,3F17.10)') zlev(jlev), &
+                     & jlev, inum(jlev,jbox,jadd,jvar), &
+                     & zmean(jlev,jbox,jadd,jvar), &
+                     & zrms(jlev,jbox,jadd,jvar), &
+                     & zsdev(jlev,jbox,jadd,jvar)
+               ENDDO
+               CLOSE(10)
+            ENDDO
+         ENDDO
+      ENDDO
+   ENDIF
+
+   IF (lhistogram) THEN
+      DO jvar=1, nvar
+         DO jadd=1, nadd
+            DO jbox=jstabox, jendbox
+               WRITE(filename,'(5A)')TRIM(cname(jvar)),&
+                  &                  TRIM(caddname(jadd)),'_',&
+                  &                  TRIM(cl_boxes(jbox)),'_histogram.dat'
+               OPEN(10,file=TRIM(filename))
+               WRITE(10,'(A10,100F10.2)')'#',(zlev(jlev),jlev=1,nmlev)
+               DO ji=1,hist(jvar)%npoints
+                  WRITE(10,'(F10.2,100I10)') &
+                     & zhistmin(jvar)+(ji-1)*zhiststep(jvar), &
+                     & (hist(jvar)%nhist(ji,jlev,jbox,jadd),jlev=1,nmlev)
+               ENDDO
+               CLOSE(10)
+            ENDDO
+         ENDDO
+      ENDDO
+   ENDIF
+
+   IF (lomona) THEN
+
+      IF (nmlev>1) THEN
+         ALLOCATE(zdat3d(nmlev,nbox,1))
+      ELSE
+         ALLOCATE(zdat2d(nbox,1))
+      ENDIF
+
+      cl_expnam=expver
+      WRITE(cl_date,'(I8.8)')inidate
+      i_dp = nmlev
+      itime=icurdate
+      linnerp=.TRUE.
+      iloopno = loopno
+      linnerini = linner
+      DO jvar = 1, nvar
+         linner = linnerini
+         loopno = iloopno
+         SELECT CASE (TRIM(cname(jvar)))
+         CASE('POTM')
+            cl_var = 'votemper'
+         CASE('PSAL')
+            cl_var='vosaline'
+         CASE('SLA')
+            cl_var='sossheig'
+         CASE('SST')
+            cl_var='sosstsst'
+         END SELECT
+         DO jadd = 1, nadd
+            linner = (caddname(jadd)(1:3)=='Hxa').OR.linner
+            IF (lpassive) THEN
+               ityp=145
+            ELSE
+               IF (linner) THEN
+                  linnerp=.TRUE.
+                  ityp=144
+                  IF (jadd>1) loopno=loopno+1
+               ELSE
+                  ityp=142
+                  IF (.NOT.linnerp) THEN
+                     IF (jadd>1) loopno=loopno+1
+                  ENDIF
+               ENDIF
+            ENDIF
+            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
+            CALL obs_variable_att(cltyp)
+            IF (nmlev>1) THEN
+               zdat3d(:,:,1) = zmean(:,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
+            ELSE
+               zdat2d(:,1) = zmean(1,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+            ENDIF
+            IF (lpassive) THEN
+               ityp=245
+            ELSE
+               IF (linner) THEN
+                  ityp=244
+               ELSE
+                  ityp=242
+               ENDIF
+            ENDIF
+            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
+            CALL obs_variable_att(cltyp)
+            IF (nmlev>1) THEN
+               zdat3d(:,:,1) = zrms(:,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
+            ELSE
+               zdat2d(:,1) = zrms(1,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+            ENDIF
+            IF (lpassive) THEN
+               ityp=345
+            ELSE
+               IF (linner) THEN
+                  ityp=344
+               ELSE
+                  ityp=342
+               ENDIF
+            ENDIF
+            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
+            CALL obs_variable_att(cltyp)
+            IF (nmlev>1) THEN
+               zdat3d(:,:,1) = zsdev(:,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
+            ELSE
+               zdat2d(:,1) = zsdev(1,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+            ENDIF
+            IF (lpassive) THEN
+               ityp=445
+            ELSE
+               IF (linner) THEN
+                  ityp=444
+               ELSE
+                  ityp=442
+               ENDIF
+            ENDIF
+            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
+            CALL obs_variable_att(cltyp)
+            IF (nmlev>1) THEN
+               zdat3d(:,:,1) = inum(:,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
+            ELSE
+               zdat2d(:,1) = inum(1,:,jadd,jvar)
+               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                  &                    cl_boxes,REAL(fbrmdi))
+            ENDIF
+            linner=.FALSE.
+         ENDDO
+      ENDDO
+
+      IF (nmlev>1) THEN
+         DEALLOCATE(zdat3d)
+      ELSE
+         DEALLOCATE(zdat2d)
+      ENDIF
+      
+   ENDIF
+
+CONTAINS
+
+   SUBROUTINE fb_stat(fbdata,nstabox, nendbox, nmlev,zlev,lnear,kqc,lhist)
+      USE fbaccdata
+      USE coords
+      TYPE(obfbdata) :: fbdata
+      INTEGER :: nstabox,nendbox
+      INTEGER :: nmlev
+      REAL :: zlev(nmlev)
+      LOGICAL :: lnear
+      INTEGER :: kqc
+      LOGICAL :: lhist
+      INTEGER :: jlev, jobs, jvar, klev,jlev2,ih
+      REAL :: zarea(4),zlat,zlon,zdiff,zdiff2
+      
+      DO jbox = nstabox, nendbox
+         CALL coord_area(cl_boxes(jbox),zarea)
+         DO jobs = 1, fbdata%nobs
+            zlat = fbdata%pphi(jobs)
+            zlon = fbdata%plam(jobs)
+            IF (zlon<0) zlon=zlon+360
+            IF (zlon>360) zlon=zlon-360
+            IF ( ( zlat .GE. zarea(3) ) .AND. &
+               & ( zlat .LE. zarea(4) ) .AND. &
+               & ( ( ( zlon .GE. zarea(1) ) .AND. &
+               &     ( zlon .LE. zarea(2) ) ) .OR. &
+               &   ( ( zarea(2) .LE. zarea(1) ) .AND. &
+               &     ( zlon .GE. zarea(1) ) .AND. &
+               &     ( zlon .LE. 360        ) ) .OR. &
+               &   ( ( zarea(2)   .LE. zarea(1) ) .AND. &
+               &     ( zlon .GE. 0          ) .AND. &
+               &     ( zlon .LE. zarea(2) ) ) ) ) THEN
+            
+               DO jlev = 1, fbdata%nlev
+                  DO jvar = 1, fbdata%nvar
+                     IF (nmlev==1) THEN
+                        klev=1
+                     ELSE
+                        IF (lnear) THEN
+                           zdiff=ABS(fbdata%pdep(jlev,jobs)-zlev(1))
+                           klev=1
+                           DO jlev2=2,nmlev
+                              zdiff2=ABS(fbdata%pdep(jlev,jobs)-zlev(jlev2))
+                              IF (zdiff2<zdiff) THEN
+                                 klev=jlev2
+                                 zdiff=zdiff2
+                              ENDIF
+                           ENDDO
+                        ELSE
+                           klev = fbdata%iobsk(jlev,jobs,jvar)-1
+                        ENDIF
+                        IF ( klev > nmlev ) THEN
+                           DO jadd = 1, fbdata%nvar
+                              IF ( ABS(fbdata%padd(jlev,jobs,jadd,jvar))<9000 ) THEN
+                                 WRITE(*,*)'Error in fb_stat'
+                                 WRITE(*,*)'Increase nmlev to at least ',klev
+                                 klev=nmlev
+                                 CALL abort
+                              ENDIF
+                           ENDDO
+                        ENDIF
+                     ENDIF
+                     IF ( fbdata%ivlqc(jlev,jobs,jvar) > kqc ) CYCLE 
+                     IF (( klev > 0 ).AND. &
+                        &(ABS(fbdata%pob(jlev,jobs,jvar)) < 9000 )) THEN
+                        DO jadd = 1, fbdata%nadd
+                           IF ( ABS(fbdata%padd(jlev,jobs,jadd,jvar)) < 9000 ) THEN
+                              zdiff = ( fbdata%padd(jlev,jobs,jadd,jvar) - &
+                                 &      fbdata%pob(jlev,jobs,jvar) )
+                              inum(klev,jbox,jadd,jvar) = inum(klev,jbox,jadd,jvar) + 1
+                              zmean(klev,jbox,jadd,jvar) = zmean(klev,jbox,jadd,jvar) + &
+                                 & zdiff
+                              zrms(klev,jbox,jadd,jvar) = zrms(klev,jbox,jadd,jvar) + &
+                                 & zdiff * zdiff
+                              IF (ABS(zdiff)>zcheck(jvar)) THEN
+                                 WRITE(*,*)'Departure outside check range ',&
+                                    & TRIM(fbdata%cname(jvar)),' entry ',&
+                                    & fbdata%caddname(jadd)
+                                 WRITE(*,*)'Depar = ',zdiff
+                                 WRITE(*,*)'Check = ',zcheck(jvar)
+                                 WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
+                                 WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
+                                 WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
+                                 WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
+                                 WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
+                                 WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
+                                 WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
+                                 WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
+                                 WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
+                              ENDIF
+                              IF (lhist) THEN
+                                 ih=NINT((zdiff-zhistmin(jvar))/zhiststep(jvar))+1
+                                 IF ((ih>=1).AND.(ih<=hist(jvar)%npoints)) THEN
+                                    hist(jvar)%nhist(ih,klev,jbox,jadd) = &
+                                       hist(jvar)%nhist(ih,klev,jbox,jadd) +1
+                                 ELSE
+                                    WRITE(*,*)'Histogram value outside range for ',&
+                                       & TRIM(fbdata%cname(jvar)),' entry ',&
+                                       & fbdata%caddname(jadd)
+                                    WRITE(*,*)'Value = ',zdiff
+                                    WRITE(*,*)'Range = ',zhistmin(jvar),zhistmax(jvar)
+                                    WRITE(*,*)'Step  = ',zhiststep(jvar)
+                                    WRITE(*,*)'Index = ',ih
+                                    WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
+                                    WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
+                                    WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
+                                    WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
+                                    WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
+                                    WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
+                                    WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
+                                    WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
+                                    WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
+                                 ENDIF
+                              ENDIF
+                           ENDIF
+                        ENDDO
+                     ENDIF
+                  ENDDO
+               ENDDO
+            ENDIF
+         ENDDO
+      ENDDO
+      
+   END SUBROUTINE fb_stat
+
+   SUBROUTINE getlevsmean(nmlev,zlev)
+      IMPLICIT NONE
+      INTEGER :: nmlev
+      REAL,DIMENSION(nmlev) :: zlev
+      REAL,DIMENSION(nmlev+1) :: ztmp
+      INTEGER :: i
+
+      zlev(:)=9999.9
+      CALL getlevs(nmlev+1,ztmp)
+      DO i=1,nmlev
+         zlev(i)=0.5*(ztmp(i)+ztmp(i+1))
+      ENDDO
+
+   END SUBROUTINE getlevsmean
+
+   SUBROUTINE getlevs(nmlev,zlev)
+      IMPLICIT NONE
+      INTEGER :: nmlev
+      REAL,DIMENSION(nmlev) :: zlev
+      
+      zlev(:)=9999.9
+
+      IF (nmlev==42) THEN
+         zlev(1)=5.02159
+         zlev(2)=15.07854
+         zlev(3)=25.16046
+         zlev(4)=35.27829
+         zlev(5)=45.44776
+         zlev(6)=55.69149
+         zlev(7)=66.04198
+         zlev(8)=76.54591
+         zlev(9)=87.27029
+         zlev(10)=98.31118
+         zlev(11)=109.8062
+         zlev(12)=121.9519
+         zlev(13)=135.0285
+         zlev(14)=149.4337
+         zlev(15)=165.7285
+         zlev(16)=184.6975
+         zlev(17)=207.4254
+         zlev(18)=235.3862
+         zlev(19)=270.5341
+         zlev(20)=315.3741
+         zlev(21)=372.9655
+         zlev(22)=446.8009
+         zlev(23)=540.5022
+         zlev(24)=657.3229
+         zlev(25)=799.5496
+         zlev(26)=967.9958
+         zlev(27)=1161.806
+         zlev(28)=1378.661
+         zlev(29)=1615.291
+         zlev(30)=1868.071
+         zlev(31)=2133.517
+         zlev(32)=2408.583
+         zlev(33)=2690.780
+         zlev(34)=2978.166
+         zlev(35)=3269.278
+         zlev(36)=3563.041
+         zlev(37)=3858.676
+         zlev(38)=4155.628
+         zlev(39)=4453.502
+         zlev(40)=4752.021
+         zlev(41)=5050.990
+         zlev(42)=5350.272
+      ELSEIF (nmlev==31) THEN
+         zlev(1)=4.999938
+         zlev(2)=15.00029
+         zlev(3)=25.00176
+         zlev(4)=35.00541
+         zlev(5)=45.01332
+         zlev(6)=55.0295
+         zlev(7)=65.06181
+         zlev(8)=75.12551
+         zlev(9)=85.25037
+         zlev(10)=95.49429
+         zlev(11)=105.9699
+         zlev(12)=116.8962
+         zlev(13)=128.6979
+         zlev(14)=142.1953
+         zlev(15)=158.9606
+         zlev(16)=181.9628
+         zlev(17)=216.6479
+         zlev(18)=272.4767
+         zlev(19)=364.303
+         zlev(20)=511.5348
+         zlev(21)=732.2009
+         zlev(22)=1033.217
+         zlev(23)=1405.698
+         zlev(24)=1830.885
+         zlev(25)=2289.768
+         zlev(26)=2768.242
+         zlev(27)=3257.479
+         zlev(28)=3752.442
+         zlev(29)=4250.401
+         zlev(30)=4749.913
+         zlev(31)=5250.227
+      ELSEIF (nmlev==1) THEN
+         zlev(1)=0.0
+      ELSE
+         WRITE(*,*) 'Unknown number of levels'
+         CALL abort
+      ENDIF
+
+   END SUBROUTINE getlevs
+
+END PROGRAM fbstat
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbthin.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbthin.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbthin.F90	(revision 2893)
@@ -0,0 +1,197 @@
+PROGRAM fbthin
+   USE obs_fbm
+   IMPLICIT NONE
+   !
+   ! Command line arguments for output file and input file
+   !
+#ifndef NOIARGTHINROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdoutfile
+   CHARACTER(len=256) :: cdinfile
+   CHARACTER(len=256) :: cdtmp
+   INTEGER :: nout,ninn,nadd,next,i,j,k
+   LOGICAL :: lgrid
+   !
+   ! Feedback data
+   !
+   TYPE(obfbdata) :: fbdatain
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF ((nargs /= 2)) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'fbthin inputfile outputfile'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdinfile)
+   CALL getarg(2,cdoutfile)
+   !
+   ! Initialize feedback data
+   !
+   CALL init_obfbdata( fbdatain )
+   !
+   ! Read the file
+   !
+   CALL read_obfbdata( TRIM(cdinfile), fbdatain )
+   !
+   ! Do the thining
+   !
+   CALL fb_thin( fbdatain )
+   !
+   ! Write the file
+   !
+   CALL write_obfbdata( TRIM(cdoutfile), fbdatain )
+
+CONTAINS
+   
+   SUBROUTINE fb_thin( fbdata )
+      !
+      ! Observation thinning stub
+      !
+      IMPLICIT NONE
+      TYPE(obfbdata) :: fbdata
+      REAL :: dlat = 1.0
+      REAL :: dlon = 1.0
+      LOGICAL :: lred = .TRUE.
+      INTEGER, ALLOCATABLE :: thinobs(:,:), thinqc(:,:)
+      REAL, ALLOCATABLE :: thindist(:,:), thingridlat(:), thingridlon(:,:)
+      INTEGER :: nlon,nlat
+      INTEGER :: ii,ij,iv,ilon,ilat,iqc,iobs,irej
+      REAL :: zlon,zdlon,zrpi,zdist
+      LOGICAL :: lexists
+      NAMELIST/namthin/dlat,dlon,lred
+      
+      INQUIRE(file='namthin.in',exist=lexists)
+      IF (lexists) THEN
+         OPEN(10,file='namthin.in')
+         READ(10,namthin)
+         CLOSE(10)
+      ENDIF
+      WRITE(*,namthin)
+
+      nlon=360/dlon
+      nlat=180/dlat+1
+      ALLOCATE(thinobs(nlon,nlat))
+      ALLOCATE(thinqc(nlon,nlat))
+      ALLOCATE(thindist(nlon,nlat))
+      ALLOCATE(thingridlat(nlat))
+      ALLOCATE(thingridlon(nlon,nlat))
+      thinobs(:,:)=0
+      thinqc(:,:)=0
+      thindist(:,:)=0
+      zrpi=ATAN(1.0)/45.
+      DO ij=1,nlat
+         thingridlat(ij)=(ij-1)*dlat-90.0
+         IF (lred) THEN
+            zdlon=MIN(dlon/MAX(COS(ABS(zrpi*thingridlat(ij))),0.0001),360.0)
+         ELSE
+            zdlon=dlon
+         ENDIF
+         DO ii=1,nlon
+            thingridlon(ii,ij)=(ii-1)*zdlon
+         ENDDO
+      ENDDO
+      irej=0
+      iobs=0
+
+      DO ii=1,fbdata%nobs
+
+         ! Skip data with missing lon and lat.
+         IF (fbdata%plam(ii)==fbrmdi) CYCLE
+         IF (fbdata%pphi(ii)==fbrmdi) CYCLE
+
+         ! Count number of observations with qcflag<=2.
+         iqc=0
+         DO iv=1,fbdata%nvar
+            DO ij=1,fbdata%nlev
+               IF(fbdata%ivlqc(ij,ii,iv)<=2) iqc=iqc+1
+            ENDDO
+         ENDDO
+
+         ! Find ilat,ilon positions
+         zlon=fbdata%plam(ii)
+         IF (zlon>=360.0) zlon=zlon-360
+         IF (zlon<0.0) zlon=zlon+360
+         ! If reduced grid change zdlon accordingly
+         IF (lred) THEN
+            zdlon=MIN(dlon/MAX(COS(ABS(zrpi*fbdata%pphi(ii))),0.0001),360.0)
+         ELSE
+            zdlon=dlon
+         ENDIF
+         ilon=INT(zlon/zdlon)+1
+         ilat=INT((fbdata%pphi(ii)+90.0)/dlat)+1
+         iobs=iobs+1
+         ! Check which observation to pick based on number of valid obs.
+         IF (thinobs(ilon,ilat)>0) THEN
+            irej=irej+1
+            IF (iqc>=thinqc(ilon,ilat)) THEN
+               zdist=distance(zrpi,zlon,fbdata%pphi(ii),&
+                  & thingridlon(ilon,ilat),thingridlat(ilat))
+               IF(iqc==thinqc(ilon,ilat)) THEN
+                  IF (zdist<thindist(ilon,ilat)) THEN
+                     fbdata%ivlqc(:,thinobs(ilon,ilat),:)=4
+                     thinobs(ilon,ilat)=ii
+                     thinqc(ilon,ilat)=iqc
+                     thindist(ilon,ilat)=zdist
+                  ELSE
+                     fbdata%ivlqc(:,ii,:)=4
+                  ENDIF
+               ELSE
+                  fbdata%ivlqc(:,thinobs(ilon,ilat),:)=4
+                  thinobs(ilon,ilat)=ii
+                  thinqc(ilon,ilat)=iqc
+                  thindist(ilon,ilat)=zdist
+               ENDIF
+            ELSE
+               fbdata%ivlqc(:,ii,:)=4
+            ENDIF
+         ELSE
+            thinobs(ilon,ilat)=ii
+            thinqc(ilon,ilat)=iqc
+            thindist(ilon,ilat)=distance(zrpi,zlon,fbdata%pphi(ii),&
+               & thingridlon(ilon,ilat),thingridlat(ilat))
+         ENDIF
+      ENDDO
+      WRITE(*,*)'Observations considered = ',iobs
+      WRITE(*,*)'Observations rejected   = ',irej
+
+      DEALLOCATE(thinobs)
+      DEALLOCATE(thinqc)
+      DEALLOCATE(thindist)
+      DEALLOCATE(thingridlat)
+      DEALLOCATE(thingridlon)
+      
+   END SUBROUTINE fb_thin
+
+   REAL FUNCTION distance( prpi, plon, plat, gridlon, gridlat )
+      REAL :: prpi,plon,plat,gridlat,gridlon
+      REAL :: zplat,zplon,zgridlat,zgridlon
+      REAL :: za1,za2,zb1,zb2,zc1,zc2,zcos1,zcos2
+
+      zplon=plon
+      zgridlon=gridlon
+      IF (zplon<-180) zplon=zplon+360.0
+      IF (zplon>=180) zplon=zplon-360.0
+      IF (zgridlon<-180) zgridlon=zgridlon+360.0
+      IF (zgridlon>=180) zgridlon=zgridlon-360.0
+      zplon=zplon*prpi
+      zgridlon=zgridlon*prpi
+      zplat=plat*prpi
+      zgridlat=gridlat*prpi
+      zcos1=COS(zplat)
+      zcos2=COS(zgridlat)
+      za1=SIN(zplat)
+      za2=SIN(zgridlat)
+      zb1=zcos1*COS(zplon)
+      zb2=zcos2*COS(zgridlon)
+      zc1=zcos1*SIN(zplon)
+      zc2=zcos2*SIN(zgridlon)
+
+      distance=ASIN(SQRT(1.0-(za1*za2+zb1*zb2+zc1*zc2)**2))
+      
+    END FUNCTION distance
+
+END PROGRAM fbthin
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/index_sort.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/index_sort.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/index_sort.F90	(revision 2893)
@@ -0,0 +1,281 @@
+MODULE index_sort
+
+CONTAINS
+
+   LOGICAL FUNCTION lessn(a,b,n)
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE lessn  ***
+      !!          
+      !! ** Purpose : Compare two array and return true if the first
+      !!              element of array "a" different from the corresponding 
+      !!              array "b" element is less than the this element
+      !!
+      !! ** Method  : 
+      !!
+      !! ** Action  : 
+      !!
+      !! References : 
+      !!
+      !! History :
+      !!        !  08-02  (K. Mogensen)  Original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      USE toolspar_kind
+      IMPLICIT NONE
+      INTEGER :: n
+      REAL(KIND=dp), DIMENSION(n) :: a,b
+      INTEGER :: i,j
+
+      DO i=1,n
+         IF (a(i)/=b(i)) THEN
+            IF (a(i)<b(i)) THEN
+               lessn=.TRUE.
+            ELSE
+               lessn=.FALSE.
+            ENDIF
+            EXIT
+         ENDIF
+      ENDDO
+
+   END FUNCTION lessn
+
+   SUBROUTINE  index_sort_dp_n(pval, n, kindx, kvals)
+      USE toolspar_kind
+      IMPLICIT NONE
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE index_sort  ***
+      !!          
+      !! ** Purpose : Get indicies for ascending order for a
+      !!              double precision array 2D
+      !!
+      !! ** Method  : Heapsort with call to lessn for comparision
+      !!
+      !! ** Action  : 
+      !!
+      !! References : http://en.wikipedia.org/wiki/Heapsort
+      !!
+      !! History :
+      !!        !  08-02  (K. Mogensen)  Original code based on index_sort_dp
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & n                               ! Number of keys
+      INTEGER, INTENT(IN) :: &
+         & kvals                           ! Number of values
+      REAL(KIND=dp),DIMENSION(n,kvals),INTENT(IN) :: &
+         & pval                            ! Array to be sorted
+      INTEGER,DIMENSION(kvals),INTENT(INOUT) :: &
+         & kindx                           ! Indicies for ordering
+      !! * Local variables
+      INTEGER :: ji, jj, jt, jn, jparent, jchild
+
+      DO ji = 1, kvals
+         kindx(ji) = ji
+      END DO
+
+      IF (kvals > 1) THEN
+
+         ji = kvals/2 + 1
+         jn = kvals
+
+         main_loop : DO
+
+            IF ( ji > 1 ) THEN
+               ji = ji-1
+               jt = kindx(ji)
+            ELSE
+               jt = kindx(jn)
+               kindx(jn) = kindx(1)
+               jn = jn-1
+               IF ( jn == 1 ) THEN
+                  kindx(1) = jt
+                  EXIT main_loop
+               ENDIF
+            ENDIF
+
+            jparent = ji
+            jchild =  2*ji
+
+            inner_loop : DO
+               IF ( jchild > jn ) EXIT inner_loop
+               IF ( jchild < jn ) THEN
+                  IF ( lessn(pval(:,kindx(jchild)),pval(:,kindx(jchild+1)),n) ) THEN
+                     jchild = jchild+1
+                  ENDIF
+               ENDIF
+               IF  ( lessn(pval(:,jt),pval(:,kindx(jchild)),n) ) THEN
+                  kindx(jparent) = kindx(jchild)
+                  jparent = jchild
+                  jchild = jchild*2
+               ELSE 
+                  jchild = jn + 1 
+               ENDIF
+            ENDDO inner_loop
+
+            kindx(jparent) = jt
+
+         END DO  main_loop
+      ENDIF
+
+   END SUBROUTINE index_sort_dp_n
+
+   SUBROUTINE  index_sort_dp(pval, kindx, kvals)
+      USE toolspar_kind
+      IMPLICIT NONE
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE index_sort  ***
+      !!          
+      !! ** Purpose : Get indicies for ascending order for a
+      !!              double precision array
+      !!
+      !! ** Method  : Heapsort
+      !!
+      !! ** Action  : 
+      !!
+      !! References : http://en.wikipedia.org/wiki/Heapsort
+      !!
+      !! History :
+      !!        !  06-05  (K. Mogensen)  Original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(KIND=dp),DIMENSION(*),INTENT(IN) :: &
+         & pval                            ! Array to be sorted
+      INTEGER,DIMENSION(*),INTENT(INOUT) :: &
+         & kindx                           ! Indicies for ordering
+      INTEGER, INTENT(IN) :: &
+         & kvals                           ! Number of values
+
+      !! * Local variables
+      INTEGER :: ji, jj, jt, jn, jparent, jchild
+
+      DO ji = 1, kvals
+         kindx(ji) = ji
+      END DO
+
+      IF (kvals > 1) THEN
+
+         ji = kvals/2 + 1
+         jn = kvals
+
+         main_loop : DO
+
+            IF ( ji > 1 ) THEN
+               ji = ji-1
+               jt = kindx(ji)
+            ELSE
+               jt = kindx(jn)
+               kindx(jn) = kindx(1)
+               jn = jn-1
+               IF ( jn == 1 ) THEN
+                  kindx(1) = jt
+                  EXIT main_loop
+               ENDIF
+            ENDIF
+
+            jparent = ji
+            jchild =  2*ji
+
+            inner_loop : DO
+               IF ( jchild > jn ) EXIT inner_loop
+               IF ( jchild < jn ) THEN
+                  IF ( pval(kindx(jchild)) <  pval(kindx(jchild+1)) ) THEN
+                     jchild = jchild+1
+                  ENDIF
+               ENDIF
+               IF  ( pval(jt) < pval(kindx(jchild))) THEN
+                  kindx(jparent) = kindx(jchild)
+                  jparent = jchild
+                  jchild = jchild*2
+               ELSE 
+                  jchild = jn + 1 
+               ENDIF
+            ENDDO inner_loop
+
+            kindx(jparent) = jt
+
+         END DO  main_loop
+
+      ENDIF
+   END SUBROUTINE index_sort_dp
+
+   SUBROUTINE  index_sort_int(kval, kindx, kvals)
+      USE toolspar_kind
+      IMPLICIT NONE
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE index_sort  ***
+      !!          
+      !! ** Purpose : Get indicies for ascending order for an
+      !!              integer array
+      !!
+      !! ** Method  : Heapsort
+      !!
+      !! ** Action  : 
+      !!
+      !! References : http://en.wikipedia.org/wiki/Heapsort
+      !!
+      !! History :
+      !!        !  06-05  (K. Mogensen)  Original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER,DIMENSION(*),INTENT(IN) :: &
+         & kval                            ! Array to be sorted
+      INTEGER,DIMENSION(*),INTENT(INOUT) :: &
+         & kindx                           ! Indicies for ordering
+      INTEGER, INTENT(IN) :: &
+         & kvals                           ! Number of values
+
+      !! * Local variables
+      INTEGER :: ji, jj, jt, jn, jparent, jchild
+
+      DO ji = 1, kvals
+         kindx(ji) = ji
+      END DO
+
+      IF (kvals > 1 ) THEN
+
+         ji = kvals/2 + 1
+         jn = kvals
+
+         main_loop : DO
+
+            IF ( ji > 1 ) THEN
+               ji = ji-1
+               jt = kindx(ji)
+            ELSE
+               jt = kindx(jn)
+               kindx(jn) = kindx(1)
+               jn = jn-1
+               IF ( jn == 1 ) THEN
+                  kindx(1) = jt
+                  EXIT main_loop
+               ENDIF
+            ENDIF
+
+            jparent = ji
+            jchild =  2*ji
+
+            inner_loop : DO
+               IF ( jchild > jn ) EXIT inner_loop
+               IF ( jchild < jn ) THEN
+                  IF ( kval(kindx(jchild)) <  kval(kindx(jchild+1)) ) THEN
+                     jchild = jchild+1
+                  ENDIF
+               ENDIF
+               IF  ( kval(jt) < kval(kindx(jchild))) THEN
+                  kindx(jparent) = kindx(jchild)
+                  jparent = jchild
+                  jchild = jchild*2
+               ELSE 
+                  jchild = jn + 1 
+               ENDIF
+            ENDDO inner_loop
+
+            kindx(jparent) = jt
+
+         END DO  main_loop
+
+      ENDIF
+
+   END SUBROUTINE index_sort_int
+
+END MODULE index_sort
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/nclatlon.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/nclatlon.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/nclatlon.F90	(revision 2893)
@@ -0,0 +1,198 @@
+PROGRAM nclatlon
+   USE netcdf
+   IMPLICIT NONE
+   !
+   ! Command line arguments for output file and input file
+   !
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdinfile
+   CHARACTER(len=64) :: cdlatvar
+   CHARACTER(len=64) :: cdlonvar
+   CHARACTER(len=256) :: cdoutfile
+   CHARACTER(len=32) :: cdscale
+   !
+   ! Input data
+   !
+   REAL, ALLOCATABLE, DIMENSION(:) :: zlat, zlon
+   REAL :: zscalat,zscalon
+   INTEGER :: nopoint
+   !
+   ! NetCDF ids
+   !
+   INTEGER :: idfile, idlat, idlon
+   INTEGER, ALLOCATABLE, DIMENSION(:) :: idimlats, idimlons, idimlens
+   INTEGER :: ndimlat, ndimlon
+   !
+   ! Loop variables
+   !
+   INTEGER :: i, isize
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF ((nargs<4).OR.(nargs>5)) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'nclatlon inputfile latvar lonvar outputfile [scale]'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdinfile)
+   CALL getarg(2,cdlatvar)
+   CALL getarg(3,cdlonvar)
+   CALL getarg(4,cdoutfile)
+   WRITE(*,*)'cdinfile  = ',TRIM(cdinfile)
+   WRITE(*,*)'cdlatvar  = ',TRIM(cdlatvar)
+   WRITE(*,*)'cdlonvar  = ',TRIM(cdlonvar)
+   WRITE(*,*)'cdoutfile = ',TRIM(cdoutfile)
+   !
+   IF (nargs==5) THEN
+      CALL getarg(5,cdscale)
+      WRITE(*,*)'cdscale = ',TRIM(cdscale)
+      READ(cdscale,*)zscalat
+      zscalon=zscalat
+   ELSE
+      zscalat=1.0
+      zscalon=1.0
+   ENDIF
+   !
+   ! Open netCDF file
+   ! 
+   IF (nf90_open(TRIM(cdinfile),nf90_nowrite,idfile)/=nf90_noerr) THEN
+      WRITE(*,*)'Error opening input file'
+      CALL abort
+   ENDIF
+   !
+   ! Get latitude id
+   !
+   IF (nf90_inq_varid(idfile,TRIM(cdlatvar),idlat)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting lat id'
+      CALL abort
+   ENDIF
+   !
+   ! Get number of latitude dimensions
+   !
+   IF (nf90_inquire_variable(idfile,idlat,ndims=ndimlat)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting number of latitude dimensions'
+      CALL abort
+   ENDIF
+   WRITE(*,*)'ndimlat=',ndimlat
+   !
+   ! Get latitude dimensions
+   !
+   ALLOCATE(idimlats(ndimlat))
+   IF (nf90_inquire_variable(idfile,idlat,dimids=idimlats)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting latitude dimensions'
+      CALL abort
+   ENDIF
+   WRITE(*,*)'idimlats=',idimlats
+   !
+   ! Get longitude id
+   !
+   IF (nf90_inq_varid(idfile,TRIM(cdlonvar),idlon)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting lon id'
+      CALL abort
+   ENDIF
+   !
+   ! Get number of longitude dimensions
+   !
+   IF (nf90_inquire_variable(idfile,idlon,ndims=ndimlon)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting number of longitude dimensions'
+      CALL abort
+   ENDIF
+   WRITE(*,*)'ndimlon=',ndimlon
+   !
+   ! Get longitude dimensions
+   !
+   ALLOCATE(idimlons(ndimlon))
+   IF (nf90_inquire_variable(idfile,idlon,dimids=idimlons)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting longitude dimensions'
+      CALL abort
+   ENDIF
+   WRITE(*,*)'idimlons=',idimlons
+   !
+   ! Check that latitude and longitude has the same shape
+   !
+   IF (ndimlon/=ndimlat) THEN
+      WRITE(*,*)'Different number of dimensions for latitude and longitude'
+      CALL abort
+   ENDIF
+   DO i=1,ndimlat
+      IF (idimlons(i)/=idimlats(i)) THEN
+         WRITE(*,*)'Different dimension for latitude and longitude'
+         CALL abort
+      ENDIF
+   ENDDO
+   !
+   ! Get dimensions
+   !
+   ALLOCATE(idimlens(ndimlat))
+   DO i=1,ndimlat
+      IF (nf90_inquire_dimension(idfile,idimlats(i),len=idimlens(i))/=nf90_noerr) THEN
+         WRITE(*,*)'Error getting dimension length'
+         CALL abort
+      ENDIF
+   ENDDO
+   WRITE(*,*)'idimlens=',idimlens
+   !
+   ! Get the data
+   !
+   isize=1
+   DO i=1,ndimlat
+      isize=isize*idimlens(i)
+   ENDDO
+   ALLOCATE(zlat(isize),zlon(isize))
+   IF (nf90_get_var(idfile,idlat,zlat)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting latitude'
+      CALL abort
+   ENDIF
+   IF (nf90_get_var(idfile,idlon,zlon)/=nf90_noerr) THEN
+      WRITE(*,*)'Error getting longitude'
+      CALL abort
+   ENDIF
+   ! 
+   ! Get the scale_factor for latitude (if it exists)
+   !
+   IF (zscalat==1.0) THEN
+      IF (nf90_inquire_attribute(idfile,idlat,&
+         &                     "scale_factor")==nf90_noerr) THEN
+         IF (nf90_get_att(idfile,idlat,"scale_factor",zscalat) &
+            & /=nf90_noerr) THEN
+            WRITE(*,*)'Error getting latitude scale_factor'
+            CALL abort
+         ENDIF
+         WRITE(*,*)'Latitude scale factor',zscalat
+      ENDIF
+   ENDIF
+   ! 
+   ! Get the scale_factor for longitude (if it exists)
+   !
+   IF (zscalon==1.0) THEN
+      IF (nf90_inquire_attribute(idfile,idlon,&
+         &                     "scale_factor")==nf90_noerr) THEN
+         IF (nf90_get_att(idfile,idlon,"scale_factor",zscalon) &
+            & /=nf90_noerr) THEN
+            WRITE(*,*)'Error getting longitude scale_factor'
+            CALL abort
+         ENDIF
+         WRITE(*,*)'Longitude scale factor',zscalon
+      ENDIF
+   ENDIF
+   !
+   ! Close netCDF file
+   ! 
+   IF (nf90_close(idfile)/=nf90_noerr) THEN
+      WRITE(*,*)'Error closing input file'
+      CALL abort
+   ENDIF
+   !
+   ! Write the data
+   !
+   OPEN(10,file=cdoutfile)
+   DO i=1,isize
+      WRITE(10,'(2F12.6)')zscalat*zlat(i),zscalon*zlon(i)
+   ENDDO
+   CLOSE(10)
+   !
+END PROGRAM nclatlon
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/sla2fb.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/sla2fb.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/sla2fb.F90	(revision 2893)
@@ -0,0 +1,64 @@
+PROGRAM sla2fb
+   USE obs_fbm
+   USE obs_sla_io
+   USE convmerge
+   IMPLICIT NONE
+   !
+   ! Command line arguments for output file and input files
+   !
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdoutfile
+   CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
+   !
+   ! Input data
+   !
+   TYPE(obfbdata), POINTER :: slaf(:)
+   INTEGER :: ninfiles,ntotobs
+   !
+   ! Output data
+   !
+   TYPE(obfbdata) :: fbdata
+   !
+   ! Loop variables
+   !
+   INTEGER :: ip,ia,ji,jk
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF (nargs < 1) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'sla2fb outputfile inputfile1 inputfile2 ...'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdoutfile)
+   !
+   ! Get input data
+   !
+   ALLOCATE( slaf(MAX(nargs-1,1)) )
+   ALLOCATE( cdinfile(nargs-1) )
+   ntotobs = 0
+   ninfiles  = nargs - 1
+   DO ia=1,ninfiles
+      CALL getarg( ia + 1, cdinfile(ia) )
+      WRITE(*,'(2A)')'File = ',TRIM(cdinfile(ia))
+      CALL read_avisofile( TRIM(cdinfile(ia)), slaf(ia), 6, .TRUE., .FALSE. )
+      WRITE(*,'(A,I9,A)')'has',slaf(ia)%nobs,' observations'
+      ntotobs = ntotobs + slaf(ia)%nobs
+   ENDDO
+   IF (ninfiles==0) THEN
+      CALL init_obfbdata( slaf(1) )
+      CALL alloc_obfbdata( slaf(1), 1, 0, 1, 0, 0, .FALSE. )
+      slaf(1)%cname(1) = 'SLA'
+      slaf(1)%cdjuldref = '19500101000000'
+   ENDIF
+   WRITE(*,'(A,I8)') 'Total observations : ',ntotobs
+   !
+   ! Merge and output the data.
+   !
+   CALL conv_fbmerge( TRIM(cdoutfile), ninfiles, slaf )
+
+END PROGRAM sla2fb
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/toolspar_kind.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/toolspar_kind.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/toolspar_kind.F90	(revision 2893)
@@ -0,0 +1,8 @@
+MODULE toolspar_kind
+  IMPLICIT NONE
+  INTEGER, PUBLIC, PARAMETER ::          &  !: Floating point section
+       sp = SELECTED_REAL_KIND( 6, 37),  &  !: single precision (real 4)
+       dp = SELECTED_REAL_KIND(12,307),  &  !: double precision (real 8)
+       wp = SELECTED_REAL_KIND(12,307),  &  !: double precision (real 8)
+       ik = SELECTED_INT_KIND(6)            !: integer precision 
+END MODULE toolspar_kind
Index: /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/vel2fb.F90
===================================================================
--- /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/vel2fb.F90	(revision 2893)
+++ /branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/vel2fb.F90	(revision 2893)
@@ -0,0 +1,70 @@
+PROGRAM vel2fb
+   USE obs_fbm
+   USE obs_vel_io
+   USE convmerge
+   IMPLICIT NONE
+   !
+   ! Command line arguments for output file and input files
+   !
+#ifndef NOIARGCPROTO
+   INTEGER,EXTERNAL :: iargc
+#endif
+   INTEGER :: nargs
+   CHARACTER(len=256) :: cdoutfile
+   CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
+   !
+   ! Input data
+   !
+   TYPE(obfbdata), POINTER :: velf(:)
+   INTEGER :: ninfiles,ntotvel,nmaxlev
+   INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
+   !
+   ! Output data
+   !
+   TYPE(obfbdata) :: fbdata
+   !
+   ! Loop variables
+   !
+   INTEGER :: ia,ii,ij
+   !
+   ! Get number of command line arguments
+   !
+   nargs=IARGC()
+   IF (nargs < 1) THEN
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'vel2fb outputfile inputfile1 inputfile2 ...'
+      CALL abort()
+   ENDIF
+   CALL getarg(1,cdoutfile)
+   !
+   ! Get input data
+   !
+   ALLOCATE( velf(MAX(nargs-1,1)) )
+   ALLOCATE( cdinfile(nargs-1) )
+   ntotvel = 0
+   ninfiles  = nargs - 1
+   DO ia = 1,ninfiles
+      CALL getarg( ia + 1, cdinfile(ia) )
+      CALL read_taondbc( TRIM(cdinfile(ia)), velf(ia), 6, .TRUE., .FALSE. )
+      WRITE(*,'(2A)')'File = ',TRIM(cdinfile(ia))
+      WRITE(*,'(A,I9,A)')'has',velf(ia)%nobs,' profiles'
+      ntotvel = ntotvel + velf(ia)%nobs
+      nmaxlev   = MAX( nmaxlev, velf(ia)%nlev )
+   ENDDO
+   IF (ninfiles==0) THEN
+      CALL init_obfbdata( velf(1) )
+      CALL alloc_obfbdata( velf(1), 2, 0, 1, 0, 1, .FALSE. )
+      velf(1)%cname(1) = 'UVEL'
+      velf(1)%cname(2) = 'VVEL'
+      velf(1)%coblong(1) = 'Zonal current'
+      velf(1)%coblong(2) = 'Meridional current'
+      velf(1)%cobunit(1) = 'Meters per second'
+      velf(1)%cobunit(2) = 'Meters per second'
+   ENDIF
+   WRITE(*,'(A,I8)') 'Total profiles : ',ntotvel
+   !
+   ! Merge and output the data.
+   !
+   CALL conv_fbmerge( TRIM(cdoutfile), ninfiles, velf )
+
+END PROGRAM vel2fb
