Changeset 3000

Timestamp:

2011-10-26T15:44:20+02:00 (13 years ago)

Author:

djlea

Message:

Updated obstools. Addition of headers to programs which explain what each utility does and how to run it. All the programs now build using the naketools utility.

Location:

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src

Files:

• bufr2fb.F90 (deleted)
• bufrdata.F90 (deleted)
• convmerge.F90 (modified) (3 diffs)
• coords.F90 (added)
• corio2fb.F90 (modified) (1 diff)
• ctl_stop.h90 (modified) (1 diff)
• date_utils.F90 (added)
• ddatetoymdhms.h90 (added)
• distance.h90 (added)
• dom_oce.F90 (added)
• enact2fb.F90 (modified) (1 diff)
• fbaccdata.F90 (modified) (1 diff)
• fbacctype.F90 (added)
• fbcomb.F90 (modified) (11 diffs)
• fbcp.F90 (deleted)
• fbmatchup.F90 (modified) (21 diffs)
• fbprint.F90 (modified) (5 diffs)
• fbsel.F90 (modified) (11 diffs)
• fbstat.F90 (modified) (19 diffs)
• fbstatncio.F90 (added)
• fbthin.F90 (modified) (3 diffs)
• greg2jul.h90 (added)
• grt_cir_dis.h90 (added)
• grt_cir_dis_saa.h90 (added)
• in_out_manager.F90 (added)
• ioncdf.F90 (added)
• lbcnfd.F90 (added)
• lib_mpp.F90 (added)
• nclatlon.F90 (deleted)
• nctools.F90 (added)
• obs_conv.F90 (added)
• obs_conv_functions.h90 (added)
• obs_fbm.F90 (added)
• obs_prof_io.F90 (added)
• obs_sla_io.F90 (added)
• obs_sla_types.F90 (added)
• obs_utils.F90 (added)
• obs_vel_io.F90 (added)
• obsprof_io.h90 (added)
• obssla_io.h90 (added)
• obssla_types.h90 (added)
• obsvel_io.h90 (added)
• omonainfo.F90 (added)
• par_kind.F90 (added)
• par_oce.F90 (added)
• proftools.F90 (added)
• sla2fb.F90 (modified) (4 diffs)
• str_c_to_for.h90 (added)
• vel2fb.F90 (modified) (1 diff)
• wrk_nemo.F90 (added)

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/convmerge.F90

@@ -3 +3 @@
    USE toolspar_kind
    USE obs_fbm
+   USE obs_utils
    IMPLICIT NONE
 
@@ -18 +19 @@
       !!  ** Action  :
       !!
+      !!   Optional :
+      !!     namelist = namobs.in    to select the observation range 
+      !!
       !!   History :
+      !!        ! 2010 (K. Mogensen) Initial version
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -122 +127 @@
 #include "greg2jul.h90"
 
-#include "ddatetoymdhms.h90"
+!#include "ddatetoymdhms.h90"
 
 END MODULE convmerge

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/corio2fb.F90

@@ -1 +1 @@
 PROGRAM corio2fb
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM corio2fb **
+   !!
+   !!  ** Purpose : Convert Coriolis format profiles to feedback format
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage:
+   !!     corio2fb.exe outputfile inputfile1 inputfile2 ...
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE obs_fbm
    USE obs_prof_io

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/ctl_stop.h90

@@ -10 +10 @@
       !!  ** Action  :
       !!
-      !!   History :   (08-12) K. Mogensen. NEMOVAR version
+      !!   History :   (2008-12) K. Mogensen. NEMOVAR version
       !!----------------------------------------------------------------------
       CHARACTER(len=*) :: cerr

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/enact2fb.F90

@@ -1 +1 @@
 PROGRAM enact2fb
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM corio2fb **
+   !!
+   !!  ** Purpose : Convert ENACT format profiles to feedback format
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage:
+   !!     enact2fb.exe outputfile inputfile1 inputfile2 ...
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE obs_fbm
    USE obs_prof_io

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/fbaccdata.F90

@@ -1 +1 @@
 MODULE fbaccdata
-   INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: inum
-   REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: zmean,zrms,zsdev
+   USE fbacctype
+   IMPLICIT NONE
+   INTEGER, DIMENSION(:,:,:,:,:), ALLOCATABLE :: inum,inumov,inumbv
+   INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: inuma
+   REAL, DIMENSION(:,:,:,:,:), ALLOCATABLE :: zbias,zrms,zsdev,zomean,zmmean,&
+      & zoemea,zovmea,zbemea,zbvmea
+   REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: zoamean
    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
-
+   REAL, DIMENSION(maxvars) :: zxymax, zxymin, zxystep
+   TYPE(xytype), DIMENSION(maxvars) :: xy
 END MODULE fbaccdata
 

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/fbcomb.F90

@@ -1 +1 @@
 PROGRAM fbcomb
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM fbcomb **
+   !!
+   !!  ** Purpose : Combine MPI decomposed feedback files into one file
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  : 
+   !!
+   !!   Usage:
+   !!     fbcomb.exe outputfile inputfile1 inputfile2 ...
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE toolspar_kind
    USE obs_fbm
@@ -23 +39 @@
    REAL(KIND=dp),ALLOCATABLE :: zsort(:,:)
    INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
+   INTEGER :: iwmo
    !
    ! Output data
@@ -30 +47 @@
    ! Loop variables
    !
-   INTEGER :: ia,iv,ii,ij, ist
+   INTEGER :: ia,iv,ii,ij
    !
    ! Get number of command line arguments
@@ -48 +65 @@
    ntotobs = 0
    ninfiles = nargs - 1
-   ist=-1
    DO ia=1, ninfiles
       CALL getarg( ia+1, cdinfile(ia) )
@@ -55 +71 @@
       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
@@ -62 +77 @@
    ! Check that the data is confirming
    !
-   DO ia=ist+1, ninfiles
-      IF ( obsdata(ia)%cdjuldref /= obsdata(ist)%cdjuldref ) THEN
+   DO ia=2, ninfiles
+      IF ( obsdata(ia)%cdjuldref /= obsdata(1)%cdjuldref ) THEN
          WRITE(*,*)'Different julian date reference. Aborting'
          CALL abort
       ENDIF
-      IF ( obsdata(ia)%nvar /= obsdata(ist)%nvar ) THEN
+      IF ( obsdata(ia)%nvar /= obsdata(1)%nvar ) THEN
          WRITE(*,*)'Different number of variables. Aborting'
          CALL abort
       ENDIF
-      IF  (obsdata(ia)%nadd /= obsdata(ist)%nadd ) THEN
+      IF  (obsdata(ia)%nadd /= obsdata(1)%nadd ) THEN
          WRITE(*,*)'Different number of additional entries. Aborting'
          CALL abort
       ENDIF
-      IF ( obsdata(ia)%next /= obsdata(ist)%next ) THEN
+      IF ( obsdata(ia)%next /= obsdata(1)%next ) THEN
          WRITE(*,*)'Different number of additional variables. Aborting'
          CALL abort
       ENDIF
-      IF ( obsdata(ia)%lgrid .NEQV. obsdata(ist)%lgrid ) THEN
+      IF ( obsdata(ia)%lgrid .NEQV. obsdata(1)%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
+         IF ( obsdata(ia)%cname(iv) /= obsdata(1)%cname(iv) ) THEN
             WRITE(*,*)'Variable name ', TRIM(obsdata(ia)%cname(iv)), &
-               &      ' is different from ', TRIM(obsdata(ist)%cname(iv)), &
+               &      ' is different from ', TRIM(obsdata(1)%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
+         IF ( obsdata(1)%lgrid ) THEN
+            IF ( obsdata(ia)%cgrid(iv) /= obsdata(1)%cgrid(iv) ) THEN
+               IF (obsdata(1)%nobs==0) THEN
+                  obsdata(1)%cgrid(iv) = obsdata(ia)%cgrid(iv)
+               ELSE
+                  IF (obsdata(ia)%nobs>0) THEN
+                     WRITE(*,*)'Grid name ', TRIM(obsdata(ia)%cgrid(iv)), &
+                        &      ' is different from ', &
+                        &      TRIM(obsdata(1)%cgrid(iv)), '. Aborting'
+                     CALL abort
+                  ENDIF
+               ENDIF
             ENDIF
          ENDIF
       ENDDO
       DO iv=1,obsdata(ia)%nadd
-         IF ( obsdata(ia)%caddname(iv) /= obsdata(ist)%caddname(iv) ) THEN
+         IF ( obsdata(ia)%caddname(iv) /= obsdata(1)%caddname(iv) ) THEN
             WRITE(*,*)'Additional name ', TRIM(obsdata(ia)%caddname(iv)), &
-               &      ' is different from ', TRIM(obsdata(ist)%caddname(iv)), &
+               &      ' is different from ', TRIM(obsdata(1)%caddname(iv)), &
                &      '. Aborting'
             CALL abort
@@ -108 +129 @@
       ENDDO
       DO iv=1,obsdata(ia)%next
-         IF ( obsdata(ia)%cextname(iv) /= obsdata(ist)%cextname(iv) ) THEN
+         IF ( obsdata(ia)%cextname(iv) /= obsdata(1)%cextname(iv) ) THEN
             WRITE(*,*)'Extra name ', TRIM(obsdata(ia)%cextname(iv)), &
-               &      ' is different from ', TRIM(obsdata(ist)%cextname(iv)), &
+               &      ' is different from ', TRIM(obsdata(1)%cextname(iv)), &
                &      '. Aborting'
             CALL abort
@@ -119 +140 @@
    ! Construct sorting arrays
    !
-   ALLOCATE( zsort(3,ntotobs), iset(ntotobs), &
+   ALLOCATE( zsort(5,ntotobs), iset(ntotobs), &
       & inum(ntotobs), iindex(ntotobs))
    ii = 0
@@ -128 +149 @@
          zsort(2,ii) = obsdata(ia)%pphi(ij)
          zsort(3,ii) = obsdata(ia)%plam(ij)
+         iwmo = TRANSFER( obsdata(ia)%cdwmo(ij)(1:4), iwmo )
+         zsort(4,ii) = iwmo
+         iwmo = TRANSFER( obsdata(ia)%cdwmo(ij)(5:8), iwmo )
+         zsort(5,ii) = iwmo
          iset(ii) = ia
          inum(ii) = ij
@@ -135 +160 @@
    ! Get indexes for time sorting.
    !
-   CALL index_sort_dp_n(zsort,3,iindex,ntotobs)
+   CALL index_sort_dp_n(zsort,5,iindex,ntotobs)
    !
    ! Allocate output data
@@ -144 +169 @@
    ENDDO
    CALL init_obfbdata( obsoutdata )
-   CALL alloc_obfbdata( obsoutdata, obsdata(ist)%nvar, ntotobs, nlev, &
-      &                 obsdata(ist)%nadd, obsdata(ist)%next, obsdata(ist)%lgrid )
+   CALL alloc_obfbdata( obsoutdata, obsdata(1)%nvar, ntotobs, nlev, &
+      &                 obsdata(1)%nadd, obsdata(1)%next, obsdata(1)%lgrid )
    !
    ! Copy input data into output data

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/fbmatchup.F90

@@ -1 +1 @@
 PROGRAM fbmatchup
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM fbmatchup **
+   !!
+   !!  ** Purpose : Find matching obs in feedback files
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !! 
+   !!   Usage:
+   !!     fbmatchup outputfile inputfile1 varname1 inputfile2 varname2 ...
+   !!
+   !!   Optional: 
+   !!     namelist = namfbmatchup.in       to set ldaily820
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE toolspar_kind
    USE obs_fbm
@@ -14 +33 @@
    CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
    CHARACTER(len=ilenname),ALLOCATABLE :: cdnames(:)
-   INTEGER :: nqc
+   CHARACTER(len=2*ilenname) :: cdtmp
    LOGICAL :: ldaily820
-   NAMELIST/namfbmatchup/nqc,ldaily820
+   NAMELIST/namfbmatchup/ldaily820
    !
    ! Input data
@@ -22 +41 @@
    TYPE(obfbdata)         :: obsdatatmp(1)
    TYPE(obfbdata),POINTER :: obsdata(:)
-   INTEGER :: ninfiles,ntotobs,nlev
+   INTEGER :: ninfiles,ntotobs,nlev,nadd,next
    !
    ! Time sorting arrays
@@ -35 +54 @@
    REAL(KIND=fbsp) :: ztim,zphi,zlam
    INTEGER(KIND=SELECTED_INT_KIND(12)) :: iwmo
+   LOGICAL, ALLOCATABLE :: ltaken(:)
    !
    ! Output data
@@ -40 +60 @@
    TYPE(obfbdata) :: obsoutdata
    !
+   ! Storage for extra to search for unique.
+   !
+   CHARACTER(len=ilenname), ALLOCATABLE :: cexttmp(:)
+   TYPE extout
+      LOGICAL, POINTER, DIMENSION(:) :: luse
+      INTEGER, POINTER, DIMENSION(:) :: ipos
+   END TYPE extout
+   TYPE(extout), POINTER, DIMENSION(:) :: pextinf
+   !
    ! Loop variables
    !
-   INTEGER :: ia,ip,i1,ii,ij,ik1,ik2,iv,ist
+   INTEGER :: ifi,ia,ip,i1,ii,ij,ik1,ik2,iv,ist,iadd,ie,iext
    LOGICAL :: llfound
-   LOGICAL :: lexists
+   LOGICAL :: lexists,lnotobs
    INTEGER :: ityp
    !
@@ -59 +88 @@
    ! Read namelist if present
    !
-   nqc=1
    ldaily820=.FALSE.
    INQUIRE(file='namfbmatchup.in',exist=lexists)
@@ -71 +99 @@
    ! Get input data
    !
-   ninfiles = ( nargs -1 )/ 2
+   ninfiles = ( nargs -1 ) / 2
    ALLOCATE( obsdata( ninfiles ) )
    ALLOCATE( cdinfile( ninfiles ) )
    ALLOCATE( cdnames( ninfiles ) )
    ip = 1
-   DO ia=1, ninfiles
+   DO ifi = 1, ninfiles
       !
       ! Read the unsorted file
       !
       ip = ip + 1
-      CALL getarg( ip, cdinfile(ia) )
+      CALL getarg( ip, cdinfile(ifi) )
       ip = ip + 1
-      CALL getarg( ip, cdnames(ia) )
+      CALL getarg( ip, cdnames(ifi) )
       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
+      CALL read_obfbdata( TRIM(cdinfile(ifi)), obsdatatmp(1) )
+      !
+      ! Check if we have fewer levels than in the first file
+      !
+      IF ( ifi > 1 )  THEN
          IF ( obsdatatmp(1)%nlev > obsdata(1)%nlev ) THEN
             WRITE(*,*)'Warning. More levels in file than the first file'
@@ -111 +127 @@
       ENDIF
       !
-      ! Check that we have few observations than in the first file
-      !
-      IF ( ia > 1 )  THEN
+      ! Check if we have fewer observations than in the first file
+      !
+      IF ( ifi > 1 )  THEN
          IF ( obsdatatmp(1)%nobs > obsdata(1)%nobs ) THEN
             WRITE(*,*)'Warning. More obs in file than the first file'
@@ -125 +141 @@
       ! Check that we have the same number of variables
       !
-      IF ( ia > 1 )  THEN
+      IF ( ifi > 1 )  THEN
          IF ( obsdatatmp(1)%nvar /= obsdata(1)%nvar ) THEN
             WRITE(*,*)'Error. Different number of variables.'
@@ -136 +152 @@
       ! Check reference datas
       !
-      IF ( ia > 1 )  THEN
+      IF ( ifi > 1 )  THEN
          IF ( obsdatatmp(1)%cdjuldref /= obsdata(1)%cdjuldref ) THEN
             WRITE(*,*)'Different reference dates'
@@ -145 +161 @@
       ! Special fix for daily average MRB data (820) for the first file
       !
-      IF (ldaily820.AND.(ia==1)) THEN
+      IF (ldaily820.AND.(ifi==1)) THEN
          DO ij = 1,obsdatatmp(1)%nobs
             READ(obsdatatmp(1)%cdtyp(ij),'(I5)')ityp
@@ -171 +187 @@
       !
       CALL index_sort_dp_n(zsort,3,iindex,obsdatatmp(1)%nobs)
-      CALL init_obfbdata( obsdata(ia) )
-      CALL alloc_obfbdata( obsdata(ia), &
+      CALL init_obfbdata( obsdata(ifi) )
+      CALL alloc_obfbdata( obsdata(ifi), &
          &                 obsdatatmp(1)%nvar, obsdatatmp(1)%nobs, &
          &                 obsdatatmp(1)%nlev, obsdatatmp(1)%nadd, &
@@ -179 +195 @@
       ! Copy input data into output data
       !
-      CALL merge_obfbdata( 1, obsdatatmp, obsdata(ia), iset, inum, iindex )
+      CALL merge_obfbdata( 1, obsdatatmp, obsdata(ifi), iset, inum, iindex )
       CALL dealloc_obfbdata( obsdatatmp(1) )
       
-      WRITE(*,'(2A)')'File = ', TRIM(cdinfile(ia))
-      WRITE(*,'(A,I9,A)')'has', obsdata(ia)%nobs, ' observations'
+      WRITE(*,'(2A)')'File = ', TRIM(cdinfile(ifi))
+      WRITE(*,'(A,I9,A)')'has', obsdata(ifi)%nobs, ' observations'
 
       DEALLOCATE( zsort, iset, inum, iindex )
@@ -190 +206 @@
    !
    ! Prepare output data
-   !
+   ! 
    CALL init_obfbdata( obsoutdata )
    !
+   ! Count number of additional fields
+   !
+   nadd = 0
+   DO ifi = 1, ninfiles
+      nadd = nadd + obsdata(ifi)%nadd
+   ENDDO
+   !
+   ! Count number of unique extra fields
+   !
+   ! First the maximum to construct list
+   next = 0
+   DO ifi = 1, ninfiles
+      next = next + obsdata(ifi)%next
+   ENDDO
+   ALLOCATE( &
+      & cexttmp(next) &
+      & )
+   ! Setup pextinf structure and search for unique extra fields
+   ALLOCATE( &
+      & pextinf(ninfiles) &
+      & )
+   next = 0
+   DO ifi = 1, ninfiles
+      ALLOCATE( &
+         & pextinf(ifi)%luse(obsdata(ifi)%next), &
+         & pextinf(ifi)%ipos(obsdata(ifi)%next)  &
+         & )
+      DO ie = 1, obsdata(ifi)%next
+         llfound = .FALSE.
+         DO ii = 1, next
+            IF ( cexttmp(ii) == obsdata(ifi)%cextname(ie) ) THEN
+               llfound = .TRUE.
+            ENDIF
+         ENDDO
+         IF (llfound) THEN
+            pextinf(ifi)%luse(ie) = .FALSE.
+            pextinf(ifi)%ipos(ie) = -1
+         ELSE
+            next = next + 1
+            pextinf(ifi)%luse(ie) = .TRUE.
+            pextinf(ifi)%ipos(ie) = next
+            cexttmp(next) = obsdata(ifi)%cextname(ie)
+         ENDIF
+      ENDDO
+   ENDDO
+   !
    ! 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,:)
+      &                kadd = nadd, kext = next )
+   ALLOCATE( ltaken(obsoutdata%nlev) )
+   iadd = 0
+   DO ifi = 1, ninfiles
+      DO ia = 1, obsdata(ifi)%nadd
+         cdtmp = TRIM(obsdata(ifi)%caddname(ia))//TRIM(cdnames(ifi))
+         obsoutdata%caddname(iadd+ia) = cdtmp(1:ilenname)
+         DO iv = 1, obsdata(ifi)%nvar
+            obsoutdata%caddlong(iadd+ia,iv) = obsdata(ifi)%caddlong(ia,iv)
+            obsoutdata%caddunit(iadd+ia,iv) = obsdata(ifi)%caddunit(ia,iv)
+         ENDDO
+      ENDDO
+      DO ie = 1, obsdata(ifi)%next
+         IF ( pextinf(ifi)%luse(ie) ) THEN
+            obsoutdata%cextname(pextinf(ifi)%ipos(ie)) = &
+               & obsdata(ifi)%cextname(ie)
+            obsoutdata%cextlong(pextinf(ifi)%ipos(ie)) = &
+               & obsdata(ifi)%cextlong(ie)
+            obsoutdata%cextunit(pextinf(ifi)%ipos(ie)) = &
+               & obsdata(ifi)%cextunit(ie)
+         ENDIF
+      ENDDO
+      iadd = iadd + obsdata(ifi)%nadd
    ENDDO
    !
@@ -220 +300 @@
    ! Merge extra data into output data
    !
-   DO ia = 2, ninfiles
+   iadd = obsdata(1)%nadd
+   DO ifi = 2, ninfiles
       ist = 1
-      DO ii = 1, obsdata(ia)%nobs
+      DO ii = 1, obsdata(ifi)%nobs
          IF (MOD(ii,10000)==1) THEN
-            WRITE(*,*)'Handling observation no ',ii,' for file no ',ia
+            WRITE(*,*)'Handling observation no ',ii,' for file no ',ifi
          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 ) 
+         iwmo = TRANSFER( obsdata(ifi)%cdwmo(ii), iwmo )
+         ztim = REAL( obsdata(ifi)%ptim(ii), fbsp )
+         zphi = REAL( obsdata(ifi)%pphi(ii), fbsp )
+         zlam = REAL( obsdata(ifi)%plam(ii), fbsp ) 
          ! Check if the the same index is the right one.
          IF ( iwmo == irwmo(ii) ) THEN
@@ -237 +318 @@
                   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
+                     i1 = ii
                   ENDIF
                ENDIF
@@ -261 +332 @@
                         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
+                           EXIT
                         ENDIF
                      ENDIF
@@ -289 +347 @@
                         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
+                           EXIT
                         ENDIF
                      ENDIF
@@ -311 +356 @@
          ! 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%ioqc(i1) = & 
+               & MAX( obsoutdata%ioqc(i1), obsdata(ifi)%ioqc(ii) )
+            obsoutdata%ipqc(i1) = &
+               & MAX( obsoutdata%ipqc(i1), obsdata(ifi)%ipqc(ii) )
+            obsoutdata%itqc(i1) = &
+               & MAX( obsoutdata%itqc(i1), obsdata(ifi)%itqc(ii) )
+            obsoutdata%ivqc(i1,:) = &
+               & MAX( obsoutdata%ivqc(i1,:), obsdata(ifi)%ivqc(ii,:) )
+            obsoutdata%ioqcf(:,i1)   = IOR( obsdata(ifi)%ioqcf(:,ii), &
                &                            obsoutdata%ioqcf(:,i1) )
-            obsoutdata%ipqcf(:,i1)   = IOR( obsdata(ia)%ipqcf(:,ii), &
+            obsoutdata%ipqcf(:,i1)   = IOR( obsdata(ifi)%ipqcf(:,ii), &
                &                            obsoutdata%ipqcf(:,i1) )
-            obsoutdata%itqcf(:,i1)   = IOR( obsdata(ia)%itqcf(:,ii), &
+            obsoutdata%itqcf(:,i1)   = IOR( obsdata(ifi)%itqcf(:,ii), &
                &                            obsoutdata%itqcf(:,i1) )
-            obsoutdata%ivqcf(:,i1,:) = IOR( obsdata(ia)%ivqcf(:,ii,:), &
+            obsoutdata%ivqcf(:,i1,:) = IOR( obsdata(ifi)%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 ) ==  &
+            ! Search for levels 
+            ltaken(:) = .FALSE.
+            DO ik1 = 1, obsdata(ifi)%nlev
+               levloop : DO ik2 = 1, obsoutdata%nlev
+                  IF ( REAL( obsdata(ifi)%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
+                     lnotobs=.TRUE.
+                     IF (ltaken(ik2)) CYCLE
+                     DO iv = 1, obsdata(ifi)%nvar
+                        IF ( REAL( obsdata(ifi)%pob(ik1,ii,iv), fbsp ) == &
+                           & REAL( obsoutdata%pob(ik2,i1,iv), fbsp ) ) THEN
+                           lnotobs=.FALSE.
+                        ENDIF
+                     ENDDO
+                     IF (lnotobs) CYCLE levloop
+                     ltaken(ik2)=.TRUE.
+                     DO ia = 1, obsdata(ifi)%nadd
+                        obsoutdata%padd(ik2,i1,iadd+ia,:) = &
+                           & obsdata(ifi)%padd(ik1,ii,ia,:)
+                     ENDDO
+                     DO ie = 1, obsdata(ifi)%next
+                        IF ( pextinf(ifi)%luse(ie) ) THEN
+                           obsoutdata%pext(ik2,i1,pextinf(ifi)%ipos(ie)) = &
+                              & obsdata(ifi)%pext(ik1,ii,ie)
+                        ENDIF
+                     ENDDO
+                     obsoutdata%idqc(ik2,i1) = &
+                        & MAX( obsoutdata%idqc(ik2,i1), obsdata(ifi)%idqc(ik1,ii) )
+                     obsoutdata%ivlqc(ik2,i1,:) = &
+                        & MAX( obsoutdata%ivlqc(ik2,i1,:), obsdata(ifi)%ivlqc(ik1,ii,:) )
                      obsoutdata%idqcf(:,ik2,i1)    = &
-                        &                IOR( obsdata(ia)%idqcf(:,ik1,ii), &
+                        &                IOR( obsdata(ifi)%idqcf(:,ik1,ii), &
                         &                     obsoutdata%idqcf(:,ik2,i1) )
                      obsoutdata%ivlqcf(:,ik2,i1,:) = &
-                        &                IOR( obsdata(ia)%ivlqcf(:,ik1,ii,:), &
+                        &                IOR( obsdata(ifi)%ivlqcf(:,ik1,ii,:), &
                         &                     obsoutdata%ivlqcf(:,ik2,i1,:) )
                      llfound = .TRUE.
                      EXIT
                   ENDIF
-               ENDDO
+               ENDDO levloop
                ! Write warning if level not found
-               IF (.NOT.llfound.AND.(obsdata(ia)%pdep(ik1,ii)/=fbrmdi)) THEN
+               IF (.NOT.llfound.AND.(obsdata(ifi)%pdep(ik1,ii)/=fbrmdi)) THEN
                   WRITE(*,*)'Level not found in first file : ',&
                      &      TRIM( cdinfile(1)  )
                   WRITE(*,*)'Data file                     : ',&
-                     &      TRIM( cdinfile(ia) )
+                     &      TRIM( cdinfile(ifi) )
                   WRITE(*,*)'Identifier                    : ',&
-                     &      obsdata(ia)%cdwmo(ii)
-                  WRITE(*,*)'Julian date                   : ',&
-                     &      obsdata(ia)%ptim(ii)
+                     &      obsdata(ifi)%cdwmo(ii)
+                  WRITE(*,*)'Julifin date                   : ',&
+                     &      obsdata(ifi)%ptim(ii)
                   WRITE(*,*)'Latitude                      : ',&
-                     &      obsdata(ia)%pphi(ii)
+                     &      obsdata(ifi)%pphi(ii)
                   WRITE(*,*)'Longitude                     : ',&
-                     &      obsdata(ia)%plam(ii)
+                     &      obsdata(ifi)%plam(ii)
                   WRITE(*,*)'Depth                         : ',&
-                     &      obsdata(ia)%pdep(ik1,ii)
+                     &      obsdata(ifi)%pdep(ik1,ii)
                ENDIF
             ENDDO
@@ -371 +437 @@
                &      TRIM( cdinfile(1)  )
             WRITE(*,*)'Data file                                : ',&
-               &      TRIM( cdinfile(ia) )
+               &      TRIM( cdinfile(ifi) )
             WRITE(*,*)'Identifier                               : ',&
-               &      obsdata(ia)%cdwmo(ii)
-            WRITE(*,*)'Julian date                              : ',&
-               &      obsdata(ia)%ptim(ii)
+               &      obsdata(ifi)%cdwmo(ii)
+            WRITE(*,*)'Julifin date                              : ',&
+               &      obsdata(ifi)%ptim(ii)
             WRITE(*,*)'Latitude                                 : ',&
-               &      obsdata(ia)%pphi(ii)
+               &      obsdata(ifi)%pphi(ii)
             WRITE(*,*)'Longitude                                : ',&
-               &      obsdata(ia)%plam(ii)
+               &      obsdata(ifi)%plam(ii)
             ist = 1
          ENDIF
       ENDDO
-      IF (obsdata(ia)%nobs>0) THEN
-         WRITE(*,*)'Handled last obs. no    ',ii,' for file no ',ia
-      ENDIF
+      IF (obsdata(ifi)%nobs>0) THEN
+         WRITE(*,*)'Handled last obs. no    ',ii,' for file no ',ifi
+      ENDIF
+      iadd = iadd + obsdata(ifi)%nadd
    ENDDO
    !
@@ -391 +458 @@
    !
    CALL write_obfbdata( TRIM(cdoutfile), obsoutdata )
+   !
+   ! Deallocate temporary data
+   ! 
+   DEALLOCATE(zrtim,zrphi,zrlam,irwmo )
+   DEALLOCATE( &
+      & cexttmp &
+      & )
+   DO ifi = 1, ninfiles
+      DEALLOCATE( &
+         & pextinf(ifi)%luse, &
+         & pextinf(ifi)%ipos  &
+         & )
+   ENDDO
+   DEALLOCATE( &
+      & pextinf &
+      & )
 
 END PROGRAM fbmatchup

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/fbprint.F90

@@ -1 +1 @@
 PROGRAM fbprint
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM fbprint **
+   !!
+   !!  ** Purpose : Print feedback file contents as text
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage :
+   !!     fbprint [options] inputfile
+   !!   Options :
+   !!     -b            shorter output
+   !!     -q            QC flags (nqc=1) select observations based on QC flags
+   !!     -Q            QC flags (nqc=2) select observations based on QC flags
+   !!     -B            QC flags (nqc=3) select observations based on QC flags
+   !!     -u            unsorted
+   !!     -s ID         select station ID  
+   !!     -t TYPE       select observation type
+   !!     -v NUM1-NUM2  select variable range to print by number (default all)
+   !!     -a NUM1-NUM2  select additional variable range to print by number (default all)
+   !!     -e NUM1-NUM2  select extra variable range to print by number (default all)
+   !!     -d            output date range
+   !!     -D            print depths
+   !!     -z            use zipped files
+   !! 
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE toolspar_kind 
    USE obs_fbm
    USE index_sort
+   USE date_utils
+   USE proftools
+
    IMPLICIT NONE
    !
@@ -11 +44 @@
 #endif
    INTEGER :: nargs
-   CHARACTER(len=256) :: cdinfile,cdbrief
-   LOGICAL :: lbrief, lqcflags, lstat
-   CHARACTER(len=8) :: cdstat
+   CHARACTER(len=256) :: cdinfile, cdbrief
+   LOGICAL :: lbrief, lqcflags, lstat, ltyp, lsort, ldaterange, lzinp, ldepths
+   CHARACTER(len=ilenwmo) :: cdstat
+   CHARACTER(len=ilentyp) :: cdtyp
    INTEGER :: nqc
+   INTEGER :: nvar1, nvar2, nadd1, nadd2, next1, next2
    !
    ! Input data
@@ -22 +57 @@
    ! Loop variables
    !
-   INTEGER :: ii
+   INTEGER :: ii, iarg, ip
    !
    ! Sorting
    !
+   INTEGER :: iwmo
    REAL(KIND=dp),ALLOCATABLE :: zsort(:,:)
    INTEGER,ALLOCATABLE  :: iindex(:)
@@ -34 +70 @@
    lbrief = .FALSE.
    lstat = .FALSE.
+   ltyp = .FALSE.
+   ldaterange = .FALSE.
+   ldepths = .FALSE.
    cdstat = 'XXXXXXX'
+   cdtyp = 'XXXX'
+   lsort = .TRUE.
    nqc = 0
-   IF ( (nargs < 1) .OR. (nargs > 3) ) THEN
+   nvar1 = -1
+   nvar2 = -1
+   nadd1 = -1
+   nadd2 = -1
+   next1 = -1
+   next2 = -1
+   lzinp = .FALSE.
+   IF ( nargs < 1 ) THEN
       CALL usage()
    ENDIF
-   IF (( nargs == 2 )) THEN
-      CALL getarg( 1, cdbrief )
+   iarg = 1
+   DO
+      IF ( iarg == nargs ) EXIT
+      CALL getarg( iarg, cdbrief )
       IF ( cdbrief == '-b' ) THEN
          lbrief = .TRUE.
+         iarg = iarg + 1
       ELSEIF( cdbrief == '-q' ) THEN
          lqcflags = .TRUE.
          nqc=1
+         iarg = iarg + 1
       ELSEIF( cdbrief == '-Q' ) THEN
          lqcflags = .TRUE.
          nqc=2
+         iarg = iarg + 1
       ELSEIF( cdbrief == '-B' ) THEN
          lqcflags = .TRUE.
          nqc=3
+         iarg = iarg + 1
+      ELSEIF( cdbrief == '-u' ) THEN
+         lsort = .FALSE.
+         iarg = iarg + 1
+      ELSEIF ( cdbrief == '-s' ) THEN
+         lstat = .TRUE.
+         CALL getarg( iarg + 1, cdstat )
+         iarg = iarg + 2 
+      ELSEIF ( cdbrief == '-t' ) THEN
+         ltyp = .TRUE.
+         CALL getarg( iarg + 1, cdtyp )
+         iarg = iarg + 2 
+      ELSEIF ( cdbrief == '-v' ) THEN
+         CALL getarg( iarg + 1, cdbrief )
+         ip=INDEX(cdbrief,'-')
+         IF (ip==0) THEN
+            READ(cdbrief,'(I10)') nvar1
+            IF (nvar1==0) THEN
+               nvar2=-1
+            ELSE
+               nvar2 = nvar1
+            ENDIF
+         ELSEIF(ip==1) THEN
+            nvar1=1
+            READ(cdbrief(ip+1:),'(I10)') nvar2
+         ELSE
+            READ(cdbrief(1:ip-1),'(I10)') nvar1
+            READ(cdbrief(ip+1:),'(I10)') nvar2
+         ENDIF
+         iarg = iarg + 2
+      ELSEIF ( cdbrief == '-a' ) THEN
+         CALL getarg( iarg + 1, cdbrief )
+         ip=INDEX(cdbrief,'-')
+         IF (ip==0) THEN
+            READ(cdbrief,'(I10)') nadd1
+            IF (nadd1==0) THEN
+               nadd2=-1
+            ELSE
+               nadd2 = nadd1
+            ENDIF
+         ELSEIF(ip==1) THEN
+            nadd1=1
+            READ(cdbrief(ip+1:),'(I10)') nadd2
+         ELSE
+            READ(cdbrief(1:ip-1),'(I10)') nadd1
+            READ(cdbrief(ip+1:),'(I10)') nadd2
+         ENDIF
+         iarg = iarg + 2
+      ELSEIF ( cdbrief == '-e' ) THEN
+         CALL getarg( iarg + 1, cdbrief )
+         ip=INDEX(cdbrief,'-')
+         IF (ip==0) THEN
+            READ(cdbrief,'(I10)') next1
+            IF (next1==0) THEN
+               next2=-1
+            ELSE
+               next2 = next1
+            ENDIF
+         ELSEIF(ip==1) THEN
+            next1=1
+            READ(cdbrief(ip+1:),'(I10)') next2
+         ELSE
+            READ(cdbrief(1:ip-1),'(I10)') next1
+            READ(cdbrief(ip+1:),'(I10)') next2
+         ENDIF
+         iarg = iarg + 2
+      ELSEIF ( cdbrief == '-d' ) THEN
+         ldaterange=.TRUE.
+         iarg = iarg + 1
+      ELSEIF ( cdbrief == '-D' ) THEN
+         ldepths=.TRUE.
+         iarg = iarg + 1
+      ELSEIF ( cdbrief == '-z' ) THEN
+         lzinp=.TRUE.
+         iarg = iarg + 1
       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
+   ENDDO
    CALL getarg( nargs, cdinfile )
    !
    ! Get input data
    !
-   CALL read_obfbdata( TRIM(cdinfile), obsdata )
+   IF (lzinp) THEN
+#if defined NOSYSTEM
+      WRITE(*,*)'Compressed files need the system subroutine call'
+      CALL abort
+#else
+      CALL system( 'cp '//TRIM(cdinfile)//' fbprint_tmp.nc.gz' )
+      CALL system( 'gzip -df fbprint_tmp.nc.gz' )
+      CALL read_obfbdata( 'fbprint_tmp.nc', obsdata )
+      CALL system( 'rm -f fbprint_tmp.nc' )
+#endif
+   ELSE
+      CALL read_obfbdata( TRIM(cdinfile), obsdata )
+   ENDIF
+   CALL sealsfromargo( obsdata )
    WRITE(*,'(2A,I9,A,I9,A)')TRIM(cdinfile), ' has ', obsdata%nobs ,&
       & ' observations and a maximum of ', obsdata%nlev, ' levels'
+   IF (nvar1<0) THEN
+      nvar1 = 1
+      nvar2 = obsdata%nvar
+   ENDIF
+   IF (nadd1<0) THEN
+      nadd1 = 1
+      nadd2 = obsdata%nadd
+   ENDIF
+   IF (next1<0) THEN
+      next1 = 1
+      next2 = obsdata%next
+   ENDIF
    !
    ! 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
+   ALLOCATE(zsort(5,obsdata%nobs),iindex(obsdata%nobs))
+   IF (lsort) THEN
+      DO ii=1,obsdata%nobs
+         zsort(1,ii)=obsdata%ptim(ii)
+         zsort(2,ii)=obsdata%pphi(ii)
+         zsort(3,ii)=obsdata%plam(ii)
+         iwmo = TRANSFER( obsdata%cdwmo(ii)(1:4), iwmo )
+         zsort(4,ii) = iwmo
+         iwmo = TRANSFER( obsdata%cdwmo(ii)(5:8), iwmo )
+         zsort(5,ii) = iwmo
+      ENDDO
+      CALL index_sort_dp_n(zsort,5,iindex,obsdata%nobs)
+   ELSE
+      DO ii=1,obsdata%nobs
+         iindex(ii)=ii
+      ENDDO
+   ENDIF
+   IF (ldaterange) THEN
+      IF (obsdata%nobs>0) THEN
+         WRITE(*,'(A)')'First observation'
+         CALL print_time(obsdata%ptim(1))
+         WRITE(*,'(A)')'Last observation'
+         CALL print_time(obsdata%ptim(obsdata%nobs))
       ENDIF
-      IF (lqcflags) THEN
-         CALL print_obs_qc(obsdata,iindex(ii),nqc)
-      ELSE
-         CALL print_obs(obsdata,iindex(ii),lbrief,lqcflags,nqc)
+   ELSE
+      !
+      ! Print the sorted list
+      !   
+      DO ii=1,obsdata%nobs
+         IF (lstat) THEN
+            IF (TRIM(ADJUSTL(cdstat)) /= &
+               &TRIM(ADJUSTL(obsdata%cdwmo(iindex(ii))))) CYCLE
+         ENDIF
+         IF (ltyp) THEN
+            IF (TRIM(ADJUSTL(cdtyp)) /= &
+               &TRIM(ADJUSTL(obsdata%cdtyp(iindex(ii))))) CYCLE
+         ENDIF
+         IF (ldepths) THEN
+            CALL print_depths(obsdata,iindex(ii))
+         ELSE
+            IF (lqcflags) THEN
+               CALL print_obs_qc(obsdata,iindex(ii),nqc,nvar1,nvar2)
+            ELSE
+               CALL print_obs(obsdata,iindex(ii),lbrief,&
+                  &           nvar1,nvar2,nadd1,nadd2,next1,next2)
+            ENDIF
+         ENDIF
+      ENDDO
+
+   ENDIF
+
+CONTAINS
+
+   SUBROUTINE usage
+      WRITE(*,'(A)')'Usage:'
+      WRITE(*,'(A)')'fbprint [options] inputfile'
+      CALL abort()
+   END SUBROUTINE usage
+   
+   SUBROUTINE print_depths(obsdata,iindex)
+      IMPLICIT NONE
+      TYPE(obfbdata) :: obsdata
+      INTEGER :: iindex
+      INTEGER :: kj
+      REAL :: mindep,maxdep
+
+      mindep=10000
+      maxdep=0
+      DO kj=1,obsdata%nlev
+         IF (obsdata%pdep(kj,iindex)<99999.0) THEN
+            IF (obsdata%pdep(kj,iindex)>maxdep) maxdep=obsdata%pdep(kj,iindex)
+            IF (obsdata%pdep(kj,iindex)<mindep) mindep=obsdata%pdep(kj,iindex)
+         ENDIF
+      ENDDO
+      
+      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)
+      CALL print_time( obsdata%ptim(iindex) )
+      WRITE(*,*)'Position QC         = ',obsdata%ipqc(iindex)
+      WRITE(*,*)'Observation QC      = ',obsdata%ioqc(iindex)
+      WRITE(*,*)'Minimum obs. depth  = ',mindep
+      WRITE(*,*)'Maximum obs. depth  = ',maxdep
+      WRITE(*,*)
+
+   END SUBROUTINE print_depths
+
+   SUBROUTINE print_obs(obsdata,iindex,lshort,&
+      &                 kvar1,kvar2,kadd1,kadd2,kext1,kext2)
+      IMPLICIT NONE
+      TYPE(obfbdata) :: obsdata
+      INTEGER :: iindex
+      LOGICAL :: lshort
+      INTEGER :: kvar1,kvar2,kadd1,kadd2,kext1,kext2
+      INTEGER :: jv,ja,je,jk
+      INTEGER :: kj
+      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)
+      CALL print_time( obsdata%ptim(iindex) )
+      WRITE(*,*)'Position QC         = ',obsdata%ipqc(iindex)
+      WRITE(*,*)'Observation QC      = ',obsdata%ioqc(iindex)
+      IF (.NOT.lshort) THEN
+         DO jv = kvar1, kvar2
+            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 = kvar1, kvar2
+            cdfmt1 = TRIM(cdfmt1)//',1X,A15,1X,A8'
+            cdfmt2 = TRIM(cdfmt2)//',1X,E15.9,1X,I8'
+            IF (kadd2-kadd1+1>0) THEN
+               WRITE(cdtmp,'(I10)')kadd2-kadd1+1
+               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 (kext2-kext1+1>0) THEN
+            WRITE(cdtmp,'(I10)')kext2-kext1+1
+            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 = kadd1, kadd2 ), &
+               & TRIM(obsdata%cname(jv))//'_K' , &
+               & jv = kvar1, kvar2 ), &
+               & ( TRIM(obsdata%cextname(ja)),&
+               & ja = kext1, kext2 )
+            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 = kadd1, kadd2 ), &
+                     & obsdata%iobsk(kj,iindex,jv), &
+                     & jv = kvar1, kvar2 ), &
+                     & ( obsdata%pext(kj,iindex,ja), ja = kext1, kext2 )
+               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 = kadd1, kadd2 ), &
+               & jv = kvar1, kvar2 ), &
+               & ( TRIM(obsdata%cextname(ja)),&
+               & ja = kext1, kext2 )
+            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 = kadd1, kadd2 ), &
+                     & jv = kvar1, kvar2 ), &
+                     & ( obsdata%pext(kj,iindex,ja), ja = kext1, kext2 )
+               ENDIF
+            ENDDO
+         ENDIF
       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(*,*)
+   END SUBROUTINE print_obs
+   
+   SUBROUTINE print_obs_qc(obsdata,iindex,kqc,kvar1,kvar2)
+      IMPLICIT NONE
+      TYPE(obfbdata) :: obsdata
+      INTEGER :: iindex
+      LOGICAL :: lqc
+      INTEGER :: kqc
+      INTEGER :: kvar1,kvar2
+      INTEGER :: jv,ja,je,jk
+      INTEGER :: kj
+      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 = kvar1, kvar2
+            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 = kvar1, kvar2
+                  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 = kvar1, kvar2
+                  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)
+      CALL print_time( obsdata%ptim(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)
+      DO jv = kvar1, kvar2
          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
+         WRITE(*,*)'Variable QC flags   = ',obsdata%ivqcf(:,iindex,jv)
       ENDDO
       cdfmt1='(1X,A8,1X,A8'
       cdfmt2='(1X,F8.2,1X,I8'
-      DO jv=1, obsdata%nvar
+      WRITE(cdtmp,'(I10)')obsdata%nqcf
+      cdfmt1 = TRIM(cdfmt1)//',1X,A18'
+      cdfmt2 = TRIM(cdfmt2)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(I9)'
+      DO jv = kvar1, kvar2
          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
+         WRITE(cdtmp,'(I10)')obsdata%nqcf
+         cdfmt1 = TRIM(cdfmt1)//',1X,A18'
+         cdfmt2 = TRIM(cdfmt2)//',1X,'//TRIM(ADJUSTL(cdtmp))//'(I9)'
       ENDDO
       IF (obsdata%next>0) THEN
@@ -163 +488 @@
       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
+      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 = kvar1, kvar2 ), &
+         & ( TRIM(obsdata%cextname(ja)),&
+         & ja = 1, obsdata%next )
       DO kj=1,obsdata%nlev
-         IF (obsdata%pdep(kj,iindex)<99999.0) THEN
+         IF (kqc>=2)  THEN
+            lskip=.TRUE.
             IF (obsdata%idqc(kj,iindex)>1) lskip=.FALSE.
-            DO jv = 1, obsdata%nvar
+            DO jv = kvar1, kvar2
                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 (lskip) CYCLE
+         ENDIF
          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
+            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 = kvar1, kvar2 ), &
+               & ( obsdata%pext(kj,iindex,ja), ja=1, obsdata%next )
+         ENDIF
+      ENDDO
+      WRITE(*,*)
+      
+   END SUBROUTINE print_obs_qc
+
+   SUBROUTINE print_time(ptim)
+      IMPLICIT NONE
+      REAL(fbdp) :: ptim
+      INTEGER:: iyr,imon,iday,ihou,imin,isec      
+      WRITE(*,*)'Julian date         = ',ptim
+      CALL jul2greg(isec,imin,ihou,iday,imon,iyr,ptim)
+      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
+   END  SUBROUTINE print_time
+
+END PROGRAM fbprint
+
    

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/fbsel.F90

@@ -1 +1 @@
 PROGRAM fbsel
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM fbsel **
+   !!
+   !!  ** Purpose : Select or subsample observations
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage:
+   !!     fbsel <input filename> <output filename>
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE obs_fbm
    USE date_utils
@@ -11 +27 @@
    INTEGER,PARAMETER :: maxdates=20
    INTEGER :: nqc,ntyp,ndates,ninidate(maxdates),nenddate(maxdates)
-   LOGICAL :: lsplitqc,lsplittyp
-   INTEGER :: iqc,ityp,idate
+   LOGICAL :: lsplitqc,lsplittyp,lsplitstat
+   INTEGER :: iqc,ityp,idate,istat
    REAL :: maxlat,minlat,maxlon,minlon
+   CHARACTER(len=ilenwmo) :: cdwmo,cdwmobeg,cdwmoend
+   CHARACTER(len=ilenwmo), DIMENSION(:), POINTER :: clstatids
+   INTEGER :: nstat
    NAMELIST/namsel/nqc,ntyp,ndates,ninidate,nenddate,lsplitqc,lsplittyp, &
-      &            maxlat,minlat,maxlon,minlon
+      &            lsplitstat,maxlat,minlat,maxlon,minlon,cdwmo,&
+      &            cdwmobeg,cdwmoend
 
    IF (iargc()/=2) THEN
@@ -34 +54 @@
    lsplitqc=.FALSE.
    lsplittyp=.FALSE.
+   lsplitstat=.FALSE.
+   cdwmo=REPEAT('X',ilenwmo)
+   cdwmobeg=cdwmo
+   cdwmoend=cdwmo
    maxlat=1e+38
    minlat=-1e+38
@@ -41 +65 @@
    READ(10,namsel)
    CLOSE(10)
+   IF (cdwmobeg==REPEAT('X',ilenwmo)) cdwmobeg=cdwmo
+   IF (cdwmoend==REPEAT('X',ilenwmo)) cdwmoend=cdwmo
    WRITE(*,namsel)
 
@@ -59 +85 @@
             CALL fb_sel(fbdatain,fbdataout,iqc,ntyp, &
                &        ninidate(idate),nenddate(idate), &
-               &        maxlat,minlat,maxlon,minlon)
+               &        maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
             WRITE(filenametmp,'(A,I2.2,A,A)')'qc_',iqc,'_',TRIM(cnameout)
             IF (fbdataout%nobs>0) THEN
@@ -73 +99 @@
             CALL fb_sel(fbdatain,fbdataout,nqc,ityp, &
                &        ninidate(idate),nenddate(idate), &
-               &        maxlat,minlat,maxlon,minlon)
+               &        maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
             WRITE(filenametmp,'(A,I4.4,A,A)')'typ_',ityp,'_',TRIM(cnameout)
             IF (fbdataout%nobs>0) THEN
@@ -83 +109 @@
             CALL dealloc_obfbdata(fbdataout)
          ENDDO
+      ELSEIF (lsplitstat) THEN
+         CALL fb_sel_uniqueids(fbdatain,clstatids,nstat)
+         DO istat=1,nstat
+            CALL fb_sel(fbdatain,fbdataout,nqc,ntyp, &
+               &        ninidate(idate),nenddate(idate), &
+               &        maxlat,minlat,maxlon,minlon,clstatids(istat),clstatids(istat))
+            WRITE(filenametmp,'(4A)')'statid_', &
+               & TRIM(clstatids(istat)),'_',TRIM(cnameout)
+            IF (fbdataout%nobs>0) THEN
+               WRITE(*,*)'Station = ',clstatids(istat)
+               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)
+            &        maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
          WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
          WRITE(*,*)'Writing file : ',TRIM(cnameout)
@@ -97 +139 @@
 
    SUBROUTINE fb_sel(fbdatain,fbdataout,nqc,ntyp,ninidate,nenddate,&
-      &              maxlat,minlat,maxlon,minlon)
+      &              maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
       TYPE(obfbdata) :: fbdatain,fbdataout
       INTEGER :: nqc,ntyp,ninidate,nenddate
       REAL :: maxlat,minlat,maxlon,minlon
+      CHARACTER(len=ilenwmo) :: cdwmobeg,cdwmoend
       INTEGER :: jobs
       INTEGER :: iqc,ityp
@@ -106 +149 @@
       INTEGER :: iyea,imon,iday
       REAL(KIND=dp) :: zjini,zjend
-
+      LOGICAL :: lcheckwmo
+
+      lcheckwmo=(cdwmobeg/=REPEAT('X',ilenwmo)).OR.&
+         &      (cdwmoend/=REPEAT('X',ilenwmo))
       iyea=ninidate/10000
       imon=ninidate/100-iyea*100
@@ -131 +177 @@
             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)
+         llvalid(jobs)=(fbdatain%pphi(jobs)<=maxlat).AND.       &
+            &          (fbdatain%pphi(jobs)>=minlat).AND.       &
+            &          (((fbdatain%plam(jobs)<=maxlon).AND.     &
+            &            (fbdatain%plam(jobs)>=minlon)).OR.     &
+            &           ((fbdatain%plam(jobs)+360<=maxlon).AND. &
+            &            (fbdatain%plam(jobs)+360>=minlon)).OR. &
+            &           ((fbdatain%plam(jobs)-360<=maxlon).AND. &
+            &            (fbdatain%plam(jobs)-360>=minlon))).AND.llvalid(jobs) 
+         IF (lcheckwmo) THEN
+            llvalid(jobs)=LGE(TRIM(fbdatain%cdwmo(jobs)),TRIM(cdwmobeg)) &
+               &    .AND. LLE(TRIM(fbdatain%cdwmo(jobs)),TRIM(cdwmoend)) &
+               &    .AND. llvalid(jobs)
+         ENDIF
          ! Add more checks here...
       ENDDO
@@ -142 +197 @@
    END SUBROUTINE fb_sel
 
+   SUBROUTINE fb_sel_uniqueids(fbdatain,clstatids,nstat)
+      TYPE(obfbdata) :: fbdatain
+      CHARACTER(len=ilenwmo), DIMENSION(:), POINTER :: clstatids
+      INTEGER :: nstat
+      INTEGER :: jobs,kobs
+      LOGICAL, DIMENSION(fbdatain%nobs) :: lunique
+
+      lunique(:)=.TRUE.
+      DO jobs=1,fbdatain%nobs
+         IF (.NOT.lunique(jobs)) CYCLE
+         DO kobs=jobs+1,fbdatain%nobs
+            IF (.NOT.lunique(kobs)) CYCLE
+            IF (fbdatain%cdwmo(jobs)==fbdatain%cdwmo(kobs)) THEN
+               lunique(kobs)=.FALSE.
+            ENDIF
+         ENDDO
+      ENDDO
+      nstat=COUNT(lunique)
+      ALLOCATE(clstatids(nstat))
+      kobs=0
+      DO jobs=1,fbdatain%nobs
+         IF (lunique(jobs)) THEN
+            kobs=kobs+1
+            clstatids(kobs)=fbdatain%cdwmo(jobs)
+         ENDIF
+      ENDDO
+      WRITE(*,*)'Unique station ids'
+      DO jobs=1,nstat
+         WRITE(*,'(I5,1X,A)')jobs,clstatids(jobs)
+      ENDDO
+
+   END SUBROUTINE fb_sel_uniqueids
+   
    SUBROUTINE check_prof(fbdata,iobs,iqc)
       

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/fbstat.F90

@@ -1 +1 @@
 PROGRAM fbstat
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM fbstat **
+   !!
+   !!  ** Purpose : Output feedback file summary info/statistics
+   !!         into a number of .dat files for different areas
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage:
+   !!     fbstat [-nmlev] <filenames>
+   !!   Optional:
+   !!     namelist = namfbstat.in
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE obs_fbm
    USE fbaccdata
    USE coords
    USE omonainfo
+   USE fbstatncio
+   USE proftools
    IMPLICIT NONE
    TYPE(obfbdata) :: fbdata
-   CHARACTER(len=256) :: filename
-   INTEGER :: jfile,jbox,jlev,jfirst,jvar,jadd,ji
+   CHARACTER(len=256) :: filename,outfilename
+   INTEGER :: jfile,jbox,jlev,jfirst,jvar,jadd,ji,ja,jt,jboxl
 #ifndef NOIARGCPROTO
    INTEGER,EXTERNAL :: iargc
@@ -13 +34 @@
    REAL,DIMENSION(:),ALLOCATABLE :: zlev
    INTEGER :: nmlev, nfiles
-   LOGICAL :: lexists,lomona,ltext
+   LOGICAL :: lexists,lomona,ltext,lnetcdf,lzinp
    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zdat3d
    REAL, ALLOCATABLE, DIMENSION(:,:) :: zdat2d
    INTEGER,DIMENSION(1) :: itime
    INTEGER :: inidate,icurdate,loopno,ityp,iloopno
-   INTEGER :: nvar,nadd
+   INTEGER :: nvar,nadd,noberr,nbgerr
    CHARACTER(len=4) :: expver
-   CHARACTER(len=7) :: cltyp
+   CHARACTER(len=20) :: 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
+   LOGICAL :: lxyplot,lrmmean
+   INTEGER :: nqc,nqco
+   REAL :: rlspc,rlmax
+   CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: cname,caddname,&
+      & cobename,cbgename
+   INTEGER, PARAMETER :: nmaxareas = 20
+   CHARACTER(len=20), DIMENSION(nmaxareas) :: carea
+   LOGICAL, DIMENSION(:), ALLOCATABLE :: lskipbox
+   INTEGER, parameter :: maxtyp = 10
+   CHARACTER(len=ilentyp), DIMENSION(maxtyp) :: ctyp
+   INTEGER :: ntyp,nboxuserl,ipdcst
+   REAL :: mindcst
+   NAMELIST/namfbstat/ltext,lomona,lnetcdf,nmlev,inidate,icurdate,loopno,&
+      &               expver,lnear,linner,lpassive,lhistogram,&
+      &               zhistmax,zhistmin,zhiststep,zcheck,carea,nmlev,&
+      &               nqc,nqco, &
+      &               rlspc,rlmax,ntyp,ctyp,&
+      &               lxyplot,zxymin,zxymax,zxystep,lzinp,lrmmean,mindcst
 
    ltext=.TRUE.
+   lnetcdf=.TRUE.
    lomona=.FALSE.
    nmlev=31
    inidate=19010101
    icurdate=19010116
-   loopno=-1
+   loopno=0
    expver='test'
    lnear=.TRUE.
@@ -45 +78 @@
    zhistmax(:)=10.0
    zhiststep(:)=0.1
-   zcheck(:)=10.0
+   zcheck(:)=1000.0
    nqc=2
-   carea='all'
+   nqco=2
+   carea(:)='all'
+   rlmax=5000.0
+   rlspc=-0.1
+   ntyp=1
+   ctyp(:)='all'
+   lxyplot=.FALSE.
+   zxymin(:)=-5.0
+   zxymax(:)=45.0
+   zxystep(:)=0.5
+   lzinp=.FALSE.
+   lrmmean=.FALSE.
+   mindcst=-1.0
+
    INQUIRE(file='namfbstat.in',exist=lexists)
    IF (lexists) THEN
@@ -55 +101 @@
       WRITE(*,namfbstat)
    ENDIF
+   mindcst=mindcst*1000.0 !From km to m.
    IF (iargc()==0) THEN
       WRITE(*,*)'Usage:'
@@ -85 +132 @@
    END DO
    nfiles=iargc()
-   
-   IF (carea/='all') THEN
+
+   CALL coord_user_init('o')
+
+   ALLOCATE(lskipbox(nboxuser))
+   lskipbox(:)=.FALSE.
+
+   IF (carea(1)/='all') THEN
       IF (lomona) THEN
-         WRITE(*,*)'For omona files carea has to be all'
+         WRITE(*,*)'For omona files carea(1) 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.
+      lskipbox(:)=.TRUE.
+      DO ji=1,nmaxareas
+         IF (carea(ji)/='all') THEN
+            lfound=.FALSE.
+            DO jbox=1,nboxuser
+               IF (TRIM(carea(ji))==TRIM(cl_boxes_user(jbox))) THEN
+                  lskipbox(jbox)=.FALSE.
+                  lfound=.TRUE.
+               ENDIF
+            ENDDO
+            IF (.NOT.lfound) THEN
+               WRITE(*,*)'Area ',TRIM(carea(ji)),' not found'
+               CALL abort
+            ENDIF
          ENDIF
       ENDDO
-      IF (.NOT.lfound) THEN
-         WRITE(*,*)'Area not found'
+      nboxuserl=0
+      DO ji=1,nboxuser
+         WRITE(*,*)'Area ',TRIM(cl_boxes_user(ji)),' is set to ',lskipbox(ji)
+         IF (.NOT.lskipbox(ji)) nboxuserl=nboxuserl+1
+      ENDDO
+      WRITE(*,*)'Total areas for statistics = ',nboxuserl
+      IF (lomona.AND.(nboxuserl/=nboxuser)) THEN
+         WRITE(*,*)'Omona files only possible if all areas'
          CALL abort
       ENDIF
    ELSE
-      jstabox=1
-      jendbox=nbox
+      nboxuserl=nboxuser
    ENDIF
-   PRINT *,'jstabox,jendbox=',jstabox,jendbox
-
-   IF (.NOT.lnear) nmlev=nmlev-1
-
-   ALLOCATE(&
-      & zlev(nmlev) &
-      & )
-   IF(lnear) THEN
-      CALL getlevs(nmlev,zlev)
+
+   IF (rlspc>0.0) THEN
+      lnear=.TRUE.
+      nmlev=rlmax/rlspc+1
+      ALLOCATE(zlev(nmlev))
+      DO ji=1,nmlev
+         zlev(ji)=(ji-1)*rlspc
+      ENDDO
    ELSE
-      CALL getlevsmean(nmlev,zlev)
+      IF (.NOT.lnear) nmlev=nmlev-1
+      ALLOCATE(&
+         & zlev(nmlev) &
+         & )
+      IF(lnear) THEN
+         CALL getlevs(nmlev,zlev)
+      ELSE
+         CALL getlevsmean(nmlev,zlev)
+      ENDIF
    ENDIF
 
@@ -123 +195 @@
       CALL getarg(jfile,filename)
       WRITE(*,*)'Handling file : ',TRIM(filename)
-      CALL read_obfbdata(TRIM(filename),fbdata)
+      CALL flush(6)
+      IF (lzinp) THEN
+#if defined NOSYSTEM
+         WRITE(*,*)'Compressed files need the system subroutine call'
+         CALL abort
+#else
+         CALL system('cp '//TRIM(filename)//' fbstat_tmp.nc.gz')
+         CALL system('gzip -df fbstat_tmp.nc.gz')
+         CALL read_obfbdata('fbstat_tmp.nc',fbdata)
+         CALL system('rm -f fbstat_tmp.nc')
+#endif
+      ELSE
+         CALL read_obfbdata(TRIM(filename),fbdata)
+      ENDIF
+      CALL sealsfromargo( fbdata )
       IF (jfile==jfirst) THEN
          nvar=fbdata%nvar
-         nadd=fbdata%nadd
+         nadd=0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:2)=='Hx') nadd=nadd+1
+         ENDDO
+         noberr=0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:3)=='OBE') noberr=noberr+1
+         ENDDO
+         nbgerr=0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:3)=='BGE') nbgerr=nbgerr+1
+         ENDDO
          IF (lhistogram) THEN
             IF (nvar>maxvars) THEN
@@ -139 +236 @@
                   &  hist(jvar)%npoints
                WRITE(*,*)'Size of histogram array (elements) = ',&
-                  &  hist(jvar)%npoints*nmlev*nbox*nadd
+                  &  hist(jvar)%npoints*nmlev*nboxuserl*nadd
                ALLOCATE(&
-                  & hist(jvar)%nhist(hist(jvar)%npoints,nmlev,nbox,nadd) &
+                  & hist(jvar)%nhist(hist(jvar)%npoints,nmlev,nboxuserl,nadd,ntyp) &
                   & )
-               hist(jvar)%nhist(:,:,:,:)=0
+               hist(jvar)%nhist(:,:,:,:,:)=0
+            ENDDO
+         ENDIF
+         ipdcst=0
+         IF (mindcst>0) THEN
+            DO ja= 1, fbdata%next
+               IF (TRIM(fbdata%cextname(ja))=='DCST') THEN
+                  ipdcst=ja
+                  EXIT 
+               ENDIF
+            ENDDO
+            IF (ipdcst==0) THEN
+               WRITE(*,*)'Distance to coast not found in file, but mindcst>0'
+               WRITE(*,*)'Extra variables:'
+               DO ja= 1, fbdata%next
+                  WRITE(*,*)ja,TRIM(fbdata%cextname(ja))
+               ENDDO
+               CALL abort
+            ENDIF
+         ENDIF
+         IF (lxyplot) THEN
+            IF (nvar>maxvars) THEN
+               WRITE(*,*)'fbstat.F90: Increase maxvars to ',nvar
+               WRITE(*,*)'if you want xyplots'
+               CALL abort
+            ENDIF
+            DO jvar = 1, nvar
+               xy(jvar)%npoints=(zxymax(jvar)-zxymin(jvar))&
+                  &              /zxystep(jvar)+1
+               WRITE(*,*)'Number of points in x and y for xyplots = ',&
+                  &  xy(jvar)%npoints
+               WRITE(*,*)'Size of xyplot array (elements)         = ',&
+                  &  xy(jvar)%npoints*xy(jvar)%npoints*nmlev*nboxuserl*nadd
+               ALLOCATE(&
+                  & xy(jvar)%nxy(xy(jvar)%npoints,xy(jvar)%npoints,&
+                  &              nmlev,nboxuserl,nadd,ntyp) &
+                  & )
+               xy(jvar)%nxy(:,:,:,:,:,:)=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)               &
+            & inum(nmlev,nboxuserl,nadd,nvar,ntyp),  &
+            & inumov(nmlev,nboxuserl,noberr,nvar,ntyp),  &
+            & inumbv(nmlev,nboxuserl,nbgerr,nvar,ntyp),  &
+            & inuma(nmlev,nboxuserl,nvar,ntyp),      &
+            & zbias(nmlev,nboxuserl,nadd,nvar,ntyp), &
+            & zrms(nmlev,nboxuserl,nadd,nvar,ntyp),  &
+            & zsdev(nmlev,nboxuserl,nadd,nvar,ntyp), &
+            & zomean(nmlev,nboxuserl,nadd,nvar,ntyp),&
+            & zmmean(nmlev,nboxuserl,nadd,nvar,ntyp),&
+            & zoemea(nmlev,nboxuserl,noberr,nvar,ntyp),&
+            & zovmea(nmlev,nboxuserl,noberr,nvar,ntyp),&
+            & zbemea(nmlev,nboxuserl,nbgerr,nvar,ntyp),&
+            & zbvmea(nmlev,nboxuserl,nbgerr,nvar,ntyp),&
+            & zoamean(nmlev,nboxuserl,nvar,ntyp),    &
+            & cname(nvar),                           &
+            & caddname(nadd),                        &
+            & cobename(noberr),                      &
+            & cbgename(nbgerr)                       &
             & )
          DO jvar = 1, nvar
             cname(jvar) = fbdata%cname(jvar)
          END DO
-         DO jadd = 1, nadd
-            caddname(jadd) = fbdata%caddname(jadd)
+         jadd = 0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:2)=='Hx') THEN
+               jadd=jadd+1
+               caddname(jadd) = fbdata%caddname(ja)
+            ENDIF
          END DO
-         inum(:,:,:,:)=0
-         zmean(:,:,:,:)=0.0
-         zrms(:,:,:,:)=0.0
-         zsdev(:,:,:,:)=0.0
+         jadd = 0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:3)=='OBE') THEN
+               jadd=jadd+1
+               cobename(jadd) = fbdata%caddname(ja)
+            ENDIF
+         END DO
+         jadd = 0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:3)=='BGE') THEN
+               jadd=jadd+1
+               cbgename(jadd) = fbdata%caddname(ja)
+            ENDIF
+         END DO
+         IF (nadd>0) THEN
+            inum(:,:,:,:,:)=0
+            zbias(:,:,:,:,:)=0.0
+            zrms(:,:,:,:,:)=0.0
+            zsdev(:,:,:,:,:)=0.0
+            zomean(:,:,:,:,:)=0.0 
+            zmmean(:,:,:,:,:)=0.0
+         ENDIF
+         IF (noberr>0) THEN
+            inumov(:,:,:,:,:)=0
+            zoemea(:,:,:,:,:)=0
+            zovmea(:,:,:,:,:)=0
+         ENDIF
+         IF (nbgerr>0) THEN
+            inumbv(:,:,:,:,:)=0
+            zbemea(:,:,:,:,:)=0
+            zbvmea(:,:,:,:,:)=0
+         ENDIF
+         inuma(:,:,:,:)=0
+         zoamean(:,:,:,:)=0.0
       ELSE
          IF (fbdata%nvar/=nvar) THEN
@@ -170 +350 @@
             CALL abort
          ENDIF
-         IF (fbdata%nadd/=nadd) THEN
-            WRITE(*,*)'Different number of nadd ',fbdata%nadd,' in ',&
+         jadd = 0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:2)=='Hx') THEN
+               jadd=jadd+1
+            ENDIF
+         END DO
+         IF (jadd/=nadd) THEN
+            WRITE(*,*)'Different number of nadd ',jadd,' in ',&
                & TRIM(filename)
             CALL abort
          ENDIF
-      ENDIF
-      CALL fb_stat(fbdata,jstabox,jendbox,nmlev,zlev,lnear,nqc,lhistogram)
+         jadd = 0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:3)=='OBE') THEN
+               jadd=jadd+1
+            ENDIF
+         END DO
+         IF (jadd/=noberr) THEN
+            WRITE(*,*)'Different number of noberr ',jadd,' in ',&
+               & TRIM(filename)
+            CALL abort
+         ENDIF
+         jadd = 0
+         DO ja= 1, fbdata%nadd
+            IF (fbdata%caddname(ja)(1:3)=='BGE') THEN
+               jadd=jadd+1
+            ENDIF
+         END DO
+         IF (jadd/=nbgerr) THEN
+            WRITE(*,*)'Different number of nbgerr ',jadd,' in ',&
+               & TRIM(filename)
+            CALL abort
+         ENDIF
+         IF (ipdcst>0) THEN
+            IF (ipdcst>fbdata%next) THEN
+               WRITE(*,*)'Distrance to coast in file not compatible with first file'
+               CALL abort
+            ENDIF
+            IF (TRIM(fbdata%cextname(ipdcst))/='DCST') THEN
+               WRITE(*,*)'Distrance to coast in file not compatible with first file'
+               CALL abort
+            ENDIF
+         ENDIF
+      ENDIF
+      IF (lrmmean) THEN
+         CALL fb_rmmean(fbdata)
+      ENDIF
+      CALL fb_stat(fbdata,lskipbox,nmlev,zlev,lnear,nqc,nqco,&
+         &         lhistogram,lxyplot,ntyp,ctyp,ipdcst,mindcst)
       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
+   DO jt=1,ntyp
+      DO jvar=1,nvar
+         DO jadd=1,nadd
+            jboxl=0
+            DO jbox=1,nboxuser
+               IF (lskipbox(jbox)) CYCLE
+               jboxl=jboxl+1
+               DO jlev = 1, nmlev
+                  IF ( inum(jlev,jboxl,jadd,jvar,jt) > 0 ) THEN
+                     zbias(jlev,jboxl,jadd,jvar,jt) = &
+                        & zbias(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt)
+                     zrms(jlev,jboxl,jadd,jvar,jt) = &
+                        & SQRT(zrms(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt))
+                     zsdev(jlev,jboxl,jadd,jvar,jt) = &
+                        & SQRT(MAX(zrms(jlev,jboxl,jadd,jvar,jt)**2-zbias(jlev,jboxl,jadd,jvar,jt)**2,0.0))
+                     zomean(jlev,jboxl,jadd,jvar,jt) = &
+                        & zomean(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt)
+                     zmmean(jlev,jboxl,jadd,jvar,jt) = &
+                        & zmmean(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt)
+                  ELSE
+                     zbias(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                     zrms(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                     zsdev(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                     zomean(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                     zmmean(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDDO
+         DO jadd=1,noberr
+            jboxl=0
+            DO jbox=1,nboxuser
+               IF (lskipbox(jbox)) CYCLE
+               jboxl=jboxl+1
+               DO jlev = 1, nmlev
+                  IF ( inumov(jlev,jboxl,jadd,jvar,jt) > 0 ) THEN
+                     zoemea(jlev,jboxl,jadd,jvar,jt) = &
+                        & zoemea(jlev,jboxl,jadd,jvar,jt)/inumov(jlev,jboxl,jadd,jvar,jt)
+                     zovmea(jlev,jboxl,jadd,jvar,jt) = &
+                        & zovmea(jlev,jboxl,jadd,jvar,jt)/inumov(jlev,jboxl,jadd,jvar,jt)
+                  ELSE
+                     zoemea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                     zovmea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDDO
+         DO jadd=1,nbgerr
+            jboxl=0
+            DO jbox=1,nboxuser
+               IF (lskipbox(jbox)) CYCLE
+               jboxl=jboxl+1
+               DO jlev = 1, nmlev
+                  IF ( inumbv(jlev,jboxl,jadd,jvar,jt) > 0 ) THEN
+                     zbemea(jlev,jboxl,jadd,jvar,jt) = &
+                        & zbemea(jlev,jboxl,jadd,jvar,jt)/inumbv(jlev,jboxl,jadd,jvar,jt)
+                     zbvmea(jlev,jboxl,jadd,jvar,jt) = &
+                        & zbvmea(jlev,jboxl,jadd,jvar,jt)/inumbv(jlev,jboxl,jadd,jvar,jt)
+                  ELSE
+                     zbemea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                     zbvmea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
+                  ENDIF
+               ENDDO
             ENDDO
          ENDDO
       ENDDO
    ENDDO
+   DO jt=1,ntyp
+      DO jvar=1,nvar
+         jboxl=0
+         DO jbox=1,nboxuser
+            IF (lskipbox(jbox)) CYCLE
+            jboxl=jboxl+1
+            DO jlev = 1, nmlev
+               IF ( inuma(jlev,jboxl,jvar,jt) > 0 ) THEN
+                  zoamean(jlev,jboxl,jvar,jt) = &
+                     & zoamean(jlev,jboxl,jvar,jt)/inuma(jlev,jboxl,jvar,jt)
+               ELSE
+                  zoamean(jlev,jboxl,jvar,jt) = 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'
+
+      DO jt=1,ntyp
+         DO jvar=1,nvar
+            DO jadd=1,nadd
+               jboxl=0
+               DO jbox=1,nboxuser
+                  IF (lskipbox(jbox)) CYCLE
+                  jboxl=jboxl+1
+                  WRITE(filename,'(7A)')TRIM(cname(jvar)),&
+                     &                  TRIM(caddname(jadd)),'_',&
+                     &                  TRIM(cl_boxes_user(jbox)),'_',&
+                     &                  TRIM(ADJUSTL(ctyp(jt))),'.dat'
+                  OPEN(10,file=TRIM(filename))
+                  DO jlev = 1, nmlev
+                     WRITE(10,'(F16.7,2I12,5F17.10)') zlev(jlev), &
+                        & jlev, inum(jlev,jboxl,jadd,jvar,jt), &
+                        & zbias(jlev,jboxl,jadd,jvar,jt), &
+                        & zrms(jlev,jboxl,jadd,jvar,jt), &
+                        & zsdev(jlev,jboxl,jadd,jvar,jt), &
+                        & zomean(jlev,jboxl,jadd,jvar,jt), &
+                        & zmmean(jlev,jboxl,jadd,jvar,jt)
+                  ENDDO
+                  CLOSE(10)
+               ENDDO
+            ENDDO
+            DO jadd=1,noberr
+               jboxl=0
+               DO jbox=1,nboxuser
+                  IF (lskipbox(jbox)) CYCLE
+                  jboxl=jboxl+1
+                  WRITE(filename,'(7A)')TRIM(cname(jvar)),&
+                     &                  TRIM(cobename(jadd)),'_',&
+                     &                  TRIM(cl_boxes_user(jbox)),'_',&
+                     &                  TRIM(ADJUSTL(ctyp(jt))),'.dat'
+                  OPEN(10,file=TRIM(filename))
+                  DO jlev = 1, nmlev
+                     WRITE(10,'(F16.7,2I12,5F17.10)') zlev(jlev), &
+                        & jlev, inumov(jlev,jboxl,jadd,jvar,jt), &
+                        & zoemea(jlev,jboxl,jadd,jvar,jt), &
+                        & zovmea(jlev,jboxl,jadd,jvar,jt)
+                  ENDDO
+                  CLOSE(10)
+               ENDDO
+            ENDDO
+            DO jadd=1,nbgerr
+               jboxl=0
+               DO jbox=1,nboxuser
+                  IF (lskipbox(jbox)) CYCLE
+                  jboxl=jboxl+1
+                  WRITE(filename,'(7A)')TRIM(cname(jvar)),&
+                     &                  TRIM(cbgename(jadd)),'_',&
+                     &                  TRIM(cl_boxes_user(jbox)),'_',&
+                     &                  TRIM(ADJUSTL(ctyp(jt))),'.dat'
+                  OPEN(10,file=TRIM(filename))
+                  DO jlev = 1, nmlev
+                     WRITE(10,'(F16.7,2I12,5F17.10)') zlev(jlev), &
+                        & jlev, inumbv(jlev,jboxl,jadd,jvar,jt), &
+                        & zbemea(jlev,jboxl,jadd,jvar,jt), &
+                        & zbvmea(jlev,jboxl,jadd,jvar,jt)
+                  ENDDO
+                  CLOSE(10)
+               ENDDO
+            ENDDO
+         ENDDO
+      ENDDO
+
+      DO jt=1,ntyp
+         DO jvar=1,nvar
+            jboxl=0
+            DO jbox=1,nboxuser
+               IF (lskipbox(jbox)) CYCLE
+               jboxl=jboxl+1
+               WRITE(filename,'(7A)')TRIM(cname(jvar)),'_',&
+                  &                  TRIM(cl_boxes_user(jbox)),'_',&
+                  &                  TRIM(ADJUSTL(ctyp(jt))),'.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)
+                  WRITE(10,'(F16.7,2I12,F17.10)') zlev(jlev), &
+                     & jlev, inuma(jlev,jboxl,jvar,jt), &
+                     & zoamean(jlev,jboxl,jvar,jt)
                ENDDO
                CLOSE(10)
@@ -220 +573 @@
          ENDDO
       ENDDO
+
+      IF (lhistogram) THEN
+         DO jt=1,ntyp
+            DO jvar=1,nvar
+               DO jadd=1,nadd
+                  jboxl=0
+                  DO jbox=1,nboxuser
+                     IF (lskipbox(jbox)) CYCLE
+                     jboxl=jboxl+1
+                     WRITE(filename,'(7A)')TRIM(cname(jvar)),&
+                        &                  TRIM(caddname(jadd)),'_',&
+                        &                  TRIM(cl_boxes_user(jbox)),'_',&
+                        &                  TRIM(ADJUSTL(ctyp(jt))),&
+                        &                  '_histogram.dat'
+                     OPEN(10,file=TRIM(filename))
+                     WRITE(10,'(A10,1000F10.2)')'#',(zlev(jlev),jlev=1,nmlev)
+                     DO ji=1,hist(jvar)%npoints
+                        WRITE(10,'(F10.2,1000I10)') &
+                           & zhistmin(jvar)+(ji-1)*zhiststep(jvar), &
+                           & (hist(jvar)%nhist(ji,jlev,jboxl,jadd,jt),jlev=1,nmlev)
+                     ENDDO
+                     CLOSE(10)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDDO
+      ENDIF
+      
    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
+   IF (lnetcdf) THEN
+
+      IF (nadd>0) THEN
+         DO jt=1,ntyp
+            WRITE(outfilename,'(3A)')'fbstat_',TRIM(ADJUSTL(ctyp(jt))),'.nc'
+            CALL fbstat_ncwrite(TRIM(outfilename),&
+               & nvar,cname,nadd,caddname,noberr,cobename,nbgerr,cbgename,&
+               & nboxuser,nboxuserl,20,cl_boxes_user,lskipbox,nmlev,zlev,&
+               & inum(:,:,:,:,jt),zbias(:,:,:,:,jt),zrms(:,:,:,:,jt), &
+               & zsdev(:,:,:,:,jt),zomean(:,:,:,:,jt),zmmean(:,:,:,:,jt),&
+               & inuma(:,:,:,jt),zoamean(:,:,:,jt), &
+               & inumov(:,:,:,:,jt),zoemea(:,:,:,:,jt),zovmea(:,:,:,:,jt), &
+               & inumbv(:,:,:,:,jt),zbemea(:,:,:,:,jt),zbvmea(:,:,:,:,jt) )
+            IF (lhistogram) THEN
+               WRITE(outfilename,'(3A)')'fbstat_hist_',TRIM(ADJUSTL(ctyp(jt))),'.nc'
+               CALL fbstat_ncwrite_hist(TRIM(outfilename),&
+                  & nvar,cname,nadd,caddname,&
+                  & nboxuser,20,cl_boxes_user,lskipbox,nmlev,zlev,&
+                  & hist,zhistmin,zhiststep,jt)
+            ENDIF
+            IF (lxyplot) THEN
+               WRITE(outfilename,'(3A)')'fbstat_xyplot_',TRIM(ADJUSTL(ctyp(jt))),'.nc'
+               CALL fbstat_ncwrite_xy(TRIM(outfilename),&
+                  & nvar,cname,nadd,caddname,&
+                  & nboxuser,20,cl_boxes_user,lskipbox,nmlev,zlev,&
+                  & xy,zxymin,zxystep,jt)
+            ENDIF
          ENDDO
-      ENDDO
+      ENDIF
    ENDIF
 
    IF (lomona) THEN
 
+      IF (ntyp>1) THEN
+         WRITE(*,*)'Omona file only supported for the first type which is : ',TRIM(ctyp(1))
+      ENDIF
       IF (nmlev>1) THEN
-         ALLOCATE(zdat3d(nmlev,nbox,1))
+         ALLOCATE(zdat3d(nmlev,nboxuser,1))
       ELSE
-         ALLOCATE(zdat2d(nbox,1))
+         ALLOCATE(zdat2d(nboxuser,1))
       ENDIF
 
@@ -257 +652 @@
       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
+      i_fill=0
+
+      DO jt=1,ntyp
+         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'
+            CASE('UVEL')
+               cl_var='vozocrtx'
+            CASE('VVEL')
+               cl_var='vomecrty'
+            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,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zbias(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zbias(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               IF (lpassive) THEN
+                  ityp=245
+               ELSE
+                  IF (linner) THEN
+                     ityp=244
+                  ELSE
+                     ityp=242
+                  ENDIF
+               ENDIF
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zrms(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zrms(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               IF (lpassive) THEN
+                  ityp=345
+               ELSE
+                  IF (linner) THEN
+                     ityp=344
+                  ELSE
+                     ityp=342
+                  ENDIF
+               ENDIF
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+               & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zsdev(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zsdev(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               IF (lpassive) THEN
+                  ityp=445
+               ELSE
+                  IF (linner) THEN
+                     ityp=444
+                  ELSE
+                     ityp=442
+                  ENDIF
+               ENDIF
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = inum(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = inum(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               IF (lpassive) THEN
+                  ityp=545
+               ELSE
+                  IF (linner) THEN
+                     ityp=544
+                  ELSE
+                     ityp=542
+                  ENDIF
+               ENDIF
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zomean(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zomean(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               IF (lpassive) THEN
+                  ityp=645
+               ELSE
+                  IF (linner) THEN
+                     ityp=644
+                  ELSE
+                     ityp=642
+                  ENDIF
+               ENDIF
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zmmean(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zmmean(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               linner=.FALSE.
+            ENDDO
+            loopno = iloopno
+            DO jadd = 1, noberr
+               linner = .TRUE.
+               ityp = 139
+               IF (jadd>1) loopno=loopno+1
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zoemea(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zoemea(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               ityp = 239
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zovmea(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zovmea(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+            ENDDO
+            loopno = iloopno
+            DO jadd = 1, nbgerr
+               linner = .TRUE.
+               ityp = 141
+               IF (jadd>1) loopno=loopno+1
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zbemea(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zbemea(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+               ityp = 241
+               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
+                  & TRIM(ADJUSTL(ctyp(jt)))
+               CALL obs_variable_att(cltyp)
+               IF (nmlev>1) THEN
+                  zdat3d(:,:,1) = zbvmea(:,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
+               ELSE
+                  zdat2d(:,1) = zbvmea(1,:,jadd,jvar,jt)
+                  i_fill=0
+                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
+                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               ENDIF
+            ENDDO
             IF (lpassive) THEN
-               ityp=145
+               ityp=745
             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
+               ityp=742
+            ENDIF
+            WRITE(cltyp,'(I3.3,A1,A)')ityp,'_',&
+               & TRIM(ADJUSTL(ctyp(jt)))
             CALL obs_variable_att(cltyp)
             IF (nmlev>1) THEN
-               zdat3d(:,:,1) = zmean(:,:,jadd,jvar)
+               zdat3d(:,:,1) = inuma(:,:,jvar,jt)
+               i_fill=0
                CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
-                  &                    cl_boxes,REAL(fbrmdi))
-               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
+                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
             ELSE
-               zdat2d(:,1) = zmean(1,:,jadd,jvar)
+               zdat2d(:,1) = inuma(1,:,jvar,jt)
+               i_fill=0
                CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
-                  &                    cl_boxes,REAL(fbrmdi))
+                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
             ENDIF
             IF (lpassive) THEN
-               ityp=245
+               ityp=845
             ELSE
-               IF (linner) THEN
-                  ityp=244
-               ELSE
-                  ityp=242
-               ENDIF
-            ENDIF
-            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
+               ityp=842
+            ENDIF
+            WRITE(cltyp,'(I3.3,A1,A)')ityp,'_',&
+               & TRIM(ADJUSTL(ctyp(jt)))
             CALL obs_variable_att(cltyp)
             IF (nmlev>1) THEN
-               zdat3d(:,:,1) = zrms(:,:,jadd,jvar)
+               zdat3d(:,:,1) = zoamean(:,:,jvar,jt)
+               i_fill=0
                CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
-                  &                    cl_boxes,REAL(fbrmdi))
-               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
+                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+               CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
             ELSE
-               zdat2d(:,1) = zrms(1,:,jadd,jvar)
+               zdat2d(:,1) = zoamean(1,:,jvar,jt)
+               i_fill=0
                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.
+                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
+            ENDIF
          ENDDO
       ENDDO
@@ -375 +952 @@
 CONTAINS
 
-   SUBROUTINE fb_stat(fbdata,nstabox, nendbox, nmlev,zlev,lnear,kqc,lhist)
+   SUBROUTINE fb_stat(fbdata,lskipbox,nmlev,zlev,lnear,kqc,kqco,lhist,lxyplot,&
+      &               ntyp,ctyp,ipdcst,mindcst)
       USE fbaccdata
       USE coords
       TYPE(obfbdata) :: fbdata
-      INTEGER :: nstabox,nendbox
+      LOGICAL, DIMENSION(nboxuser) :: lskipbox
       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
+      INTEGER :: kqc,kqco
+      LOGICAL :: lhist,lxyplot
+      INTEGER :: ntyp
+      CHARACTER(len=ilentyp), DIMENSION(ntyp) :: ctyp
+      INTEGER :: ipdcst
+      REAL :: mindcst
+      INTEGER, DIMENSION(nboxuser) :: jlboxnum
+      INTEGER :: jlev, jobs, jvar, klev,jlev2,ih,ja,jadd,jbox,jt,ix,iy,jboxl
+      REAL :: zarea(4),zlat,zlon,zdiff,zdiff2,zvar
       
-      DO jbox = nstabox, nendbox
-         CALL coord_area(cl_boxes(jbox),zarea)
+      jboxl=0
+      jlboxnum=-1
+      DO jbox = 1, nboxuser
+         IF (lskipbox(jbox)) CYCLE
+         jboxl=jboxl+1
+         jlboxnum(jbox)=jboxl
+      ENDDO
+
+      !$omp parallel do default(shared) private(jlev,jobs,jvar,klev,jlev2,ih,ja,jadd,jbox,jt,ix,iy,jboxl,zarea,zlat,zlon,zdiff,zdiff2)
+      DO jbox = 1, nboxuser
+         IF (lskipbox(jbox)) CYCLE
+         jboxl=jlboxnum(jbox)
+         CALL coord_area_user(cl_boxes_user(jbox),zarea)
          DO jobs = 1, fbdata%nobs
+            ! Reject observations with observation, position or time flag rejections
+            IF (fbdata%ioqc(jobs)>kqco) CYCLE
+            IF (fbdata%ipqc(jobs)>kqco) CYCLE
+            IF (fbdata%itqc(jobs)>kqco) CYCLE
             zlat = fbdata%pphi(jobs)
             zlon = fbdata%plam(jobs)
@@ -407 +1005 @@
             
                DO jlev = 1, fbdata%nlev
+                  IF (ipdcst>0) THEN
+                     IF (fbdata%pext(jlev,jobs,ipdcst)==fbrmdi) CYCLE
+                     IF (fbdata%pext(jlev,jobs,ipdcst)<mindcst) CYCLE
+                  ENDIF
                   DO jvar = 1, fbdata%nvar
                      IF (nmlev==1) THEN
@@ -425 +1027 @@
                         ENDIF
                         IF ( klev > nmlev ) THEN
-                           DO jadd = 1, fbdata%nvar
-                              IF ( ABS(fbdata%padd(jlev,jobs,jadd,jvar))<9000 ) THEN
+                           DO ja = 1, fbdata%nadd
+                              IF ( fbdata%caddname(ja)(1:2) /= 'Hx' ) CYCLE
+                              IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar))<9000 ) THEN
                                  WRITE(*,*)'Error in fb_stat'
                                  WRITE(*,*)'Increase nmlev to at least ',klev
@@ -435 +1038 @@
                         ENDIF
                      ENDIF
+                     IF (( klev > 0 ).AND. &
+                        &(ABS(fbdata%pob(jlev,jobs,jvar)) < 9000 )) THEN
+                        DO jt=1,ntyp
+                           IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
+                              IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
+                           ENDIF
+                           inuma(klev,jboxl,jvar,jt) = inuma(klev,jboxl,jvar,jt) + 1
+                           zoamean(klev,jboxl,jvar,jt) = zoamean(klev,jboxl,jvar,jt) + &
+                              & fbdata%pob(jlev,jobs,jvar)
+                        ENDDO
+                     ENDIF
+                     IF ( fbdata%ivlqc(jlev,jobs,jvar) < 0 ) CYCLE 
                      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) - &
+                        jadd = 0
+                        DO ja = 1, fbdata%nadd
+                           IF ( fbdata%caddname(ja)(1:2) /= 'Hx' ) CYCLE
+                           jadd = jadd + 1
+                           IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar)) < 9000 ) THEN
+                              zdiff = ( fbdata%padd(jlev,jobs,ja,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
+                              DO jt=1,ntyp
+                                 IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
+                                    IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
+                                 ENDIF
+                                 inum(klev,jboxl,jadd,jvar,jt) = inum(klev,jboxl,jadd,jvar,jt) + 1
+                                 zbias(klev,jboxl,jadd,jvar,jt) = zbias(klev,jboxl,jadd,jvar,jt) + &
+                                    & zdiff
+                                 zrms(klev,jboxl,jadd,jvar,jt) = zrms(klev,jboxl,jadd,jvar,jt) + &
+                                    & zdiff * zdiff
+                                 zomean(klev,jboxl,jadd,jvar,jt) = zomean(klev,jboxl,jadd,jvar,jt) + &
+                                    & fbdata%pob(jlev,jobs,jvar)
+                                 zmmean(klev,jboxl,jadd,jvar,jt) = zmmean(klev,jboxl,jadd,jvar,jt) + &
+                                    & fbdata%padd(jlev,jobs,ja,jvar)
+                              ENDDO
                               IF (ABS(zdiff)>zcheck(jvar)) THEN
                                  WRITE(*,*)'Departure outside check range ',&
@@ -459 +1086 @@
                                  WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
                                  WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
-                                 WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
+                                 WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,ja,jvar)
                                  WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
                                  WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
@@ -466 +1093 @@
                                  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
+                                    DO jt=1,ntyp
+                                       IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
+                                          IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
+                                       ENDIF
+                                       hist(jvar)%nhist(ih,klev,jboxl,jadd,jt) = &
+                                          hist(jvar)%nhist(ih,klev,jboxl,jadd,jt) +1
+                                    ENDDO
                                  ELSE
                                     WRITE(*,*)'Histogram value outside range for ',&
@@ -476 +1108 @@
                                     WRITE(*,*)'Step  = ',zhiststep(jvar)
                                     WRITE(*,*)'Index = ',ih
-                                    WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
+                                    WRITE(*,*)'Id    = ',TRIM(fbdata%cdwmo(jobs))
+                                    WRITE(*,*)'Typ   = ',TRIM(fbdata%cdtyp(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
+                              IF (lxyplot) THEN
+                                 ix=NINT((fbdata%pob(jlev,jobs,jvar)-zxymin(jvar))/&
+                                    &    zxystep(jvar))+1
+                                 iy=NINT((fbdata%padd(jlev,jobs,ja,jvar)-zxymin(jvar))/&
+                                    &    zxystep(jvar))+1
+                                 IF ((ix>=1).AND.(ix<=xy(jvar)%npoints).AND. &
+                                    &(iy>=1).AND.(iy<=xy(jvar)%npoints)) THEN
+                                    DO jt=1,ntyp
+                                       IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
+                                          IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
+                                       ENDIF
+                                       xy(jvar)%nxy(ix,iy,klev,jboxl,jadd,jt) = &
+                                          xy(jvar)%nxy(ix,iy,klev,jboxl,jadd,jt) +1
+                                    ENDDO
+                                 ELSE
+                                    WRITE(*,*)'xy plot values outside range for ',&
+                                       & TRIM(fbdata%cname(jvar)),' entry ',&
+                                       & fbdata%caddname(jadd)
+                                    WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
+                                    WRITE(*,*)'Model = ',fbdata%padd(jlev,jobs,ja,jvar)
+                                    WRITE(*,*)'Range = ',zxymin(jvar),zxymax(jvar)
+                                    WRITE(*,*)'Step  = ',zxystep(jvar)
+                                    WRITE(*,*)'Index = ',ih
+                                    WRITE(*,*)'Id    = ',TRIM(fbdata%cdwmo(jobs))
+                                    WRITE(*,*)'Typ   = ',TRIM(fbdata%cdtyp(jobs))
                                     WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
                                     WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
@@ -489 +1157 @@
                            ENDIF
                         ENDDO
+                        jadd = 0
+                        DO ja = 1, fbdata%nadd
+                           IF ( fbdata%caddname(ja)(1:3) /= 'OBE' ) CYCLE
+                           jadd = jadd + 1
+                           IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar)) < 9000 ) THEN
+                              zvar = fbdata%padd(jlev,jobs,ja,jvar)*fbdata%padd(jlev,jobs,ja,jvar)
+                              DO jt=1,ntyp
+                                 IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
+                                    IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
+                                 ENDIF
+                                 inumov(klev,jboxl,jadd,jvar,jt) = inumov(klev,jboxl,jadd,jvar,jt) + 1
+                                 zoemea(klev,jboxl,jadd,jvar,jt) = zoemea(klev,jboxl,jadd,jvar,jt) + &
+                                    & fbdata%padd(jlev,jobs,ja,jvar)                                 
+                                 zovmea(klev,jboxl,jadd,jvar,jt) = zovmea(klev,jboxl,jadd,jvar,jt) + zvar
+                              ENDDO
+                           ENDIF
+                        ENDDO
+                        jadd = 0
+                        DO ja = 1, fbdata%nadd
+                           IF ( fbdata%caddname(ja)(1:3) /= 'BGE' ) CYCLE
+                           jadd = jadd + 1
+                           IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar)) < 9000 ) THEN
+                              zvar = fbdata%padd(jlev,jobs,ja,jvar)*fbdata%padd(jlev,jobs,ja,jvar)
+                              DO jt=1,ntyp
+                                 IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
+                                    IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
+                                 ENDIF
+                                 inumbv(klev,jboxl,jadd,jvar,jt) = inumbv(klev,jboxl,jadd,jvar,jt) + 1
+                                 zbemea(klev,jboxl,jadd,jvar,jt) = zbemea(klev,jboxl,jadd,jvar,jt) + &
+                                    & fbdata%padd(jlev,jobs,ja,jvar)
+                                 zbvmea(klev,jboxl,jadd,jvar,jt) = zbvmea(klev,jboxl,jadd,jvar,jt) + zvar
+                              ENDDO
+                           ENDIF
+                        ENDDO
                      ENDIF
                   ENDDO
@@ -495 +1197 @@
          ENDDO
       ENDDO
+      !$omp end parallel do
       
    END SUBROUTINE fb_stat
+
+   SUBROUTINE fb_rmmean(fbdata)
+      TYPE(obfbdata) :: fbdata
+      INTEGER :: jadd,jmean
+
+      jmean=0
+      DO jadd=1,fbdata%nadd
+         IF (fbdata%caddname(jadd)(1:4)=='MEAN') THEN
+            jmean=jadd
+            EXIT
+         ENDIF
+      ENDDO
+      IF (jmean==0) THEN
+         WRITE(*,*)'Warning: MEAN additional variable not found'
+         RETURN
+      ENDIF
+      IF (fbdata%nobs>0) THEN
+         DO jadd=1,fbdata%nadd 
+            IF (fbdata%caddname(jadd)(1:2)=='Hx') THEN
+               fbdata%padd(:,:,jadd,:)=fbdata%padd(:,:,jadd,:)&
+                  &                   +fbdata%padd(:,:,jmean,:)
+            ENDIF
+         ENDDO
+      ENDIF
+
+   END SUBROUTINE fb_rmmean
 
    SUBROUTINE getlevsmean(nmlev,zlev)

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/fbthin.F90

@@ -1 +1 @@
 PROGRAM fbthin
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM fbthin **
+   !!
+   !!  ** Purpose : Thin the data to 1 degree resolution
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage:
+   !!     fbthin inputfile outputfile
+   !!
+   !!   Required:
+   !!     namelist = namthin.in
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE obs_fbm
    IMPLICIT NONE
@@ -5 +24 @@
    ! Command line arguments for output file and input file
    !
-#ifndef NOIARGTHINROTO
+#ifndef NOIARGCPROTO
    INTEGER,EXTERNAL :: iargc
 #endif
@@ -50 +69 @@
    SUBROUTINE fb_thin( fbdata )
       !
-      ! Observation thinning stub
+      ! Observation thinning 
       !
       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
+      ! Namelist parameters
+      INTEGER, PARAMETER :: nmaxtypes = 10
+      CHARACTER(len=ilentyp), DIMENSION(nmaxtypes) :: thintypes
+      REAL, DIMENSION(nmaxtypes) :: thindists, thindtime
+      ! Local variables
+      NAMELIST/namthin/thintypes, thindists, thindtime
+      INTEGER :: it,ii,ij,iv,iobs,irej
+      REAL :: zdist
+
+      ! Get namelist
+      thintypes(:) = 'XXXX'
+      ! Distance in km
+      thindists(:) = 100.0
+      ! Time difference in days
+      thindtime(:) = 0.99999999
+      OPEN(10,file='namthin.in')
+      READ(10,namthin)
+      CLOSE(10)
+      WRITE(*,namthin)
       
-      INQUIRE(file='namthin.in',exist=lexists)
-      IF (lexists) THEN
-         OPEN(10,file='namthin.in')
-         READ(10,namthin)
-         CLOSE(10)
-      ENDIF
-      WRITE(*,namthin)
+      ! Convert to meters
+      thindists(:) = thindists(:) * 1000.0
 
-      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
+      DO it = 1, nmaxtypes
+
+         IF ( TRIM(thintypes(it)) == 'XXXX' ) CYCLE
+
+         iobs = 0 
+         irej = 0
+
+         master_loop: DO ii= 1, fbdata%nobs 
+            
+            IF ( TRIM(ADJUSTL(thintypes(it))) /= 'all' ) THEN
+               IF ( TRIM(ADJUSTL(fbdata%cdtyp(ii))) /= &
+                  & TRIM(ADJUSTL(thintypes(it))) ) CYCLE
+            ENDIF
+            
+            iobs = iobs + 1
+
+            ! Skip data with missing lon and lat and observation flag rejected.
+
+            IF (fbdata%plam(ii)==fbrmdi) CYCLE
+            IF (fbdata%pphi(ii)==fbrmdi) CYCLE
+            IF (fbdata%ioqc(ii)>2) CYCLE
+            
+            DO ij=ii+1, fbdata%nobs
+
+               ! Skip data with missing lon and lat and observation flag rejected.
+               
+               IF (fbdata%plam(ij)==fbrmdi) CYCLE
+               IF (fbdata%pphi(ij)==fbrmdi) CYCLE
+               IF (fbdata%ioqc(ij)>2) CYCLE
+               
+               ! Skip different type unless thintypes is 'all'
+
+               IF ( TRIM(ADJUSTL(thintypes(it))) /= 'all' ) THEN
+                  IF ( TRIM(ADJUSTL(fbdata%cdtyp(ij))) /= &
+                     & TRIM(ADJUSTL(thintypes(it))) ) CYCLE
+               ENDIF
+
+               IF ( ABS( fbdata%ptim(ij) - fbdata%ptim(ii) ) &
+                  & >= thindtime(it) ) CYCLE
+
+               zdist = distance( fbdata%plam(ii), fbdata%pphi(ii), &
+                  &              fbdata%plam(ij), fbdata%pphi(ij) )
+
+               IF ( zdist < thindists(it) ) THEN
+
+                  irej = irej + 1
+                  fbdata%ioqc(ij)    = 4
+                  fbdata%ioqcf(2,ij) = fbdata%ioqcf(2,ij) + 32
+
+               ENDIF
+            ENDDO
+            
+         ENDDO master_loop
+
+         WRITE(*,*)'For type                = ',TRIM(thintypes(it))
+         WRITE(*,*)'Observations considered = ',iobs
+         WRITE(*,*)'Observations rejected   = ',irej
+
       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
+#include "distance.h90"
 
 END PROGRAM fbthin

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/sla2fb.F90

@@ -1 +1 @@
 PROGRAM sla2fb
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM sla2fb **
+   !!
+   !!  ** Purpose : Convert AVISO SLA format to feedback format
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage:
+   !!     sla2fb [-s type] outputfile inputfile1 inputfile2 ...
+   !!   Option:
+   !!     -s            Select altimeter data_source
+   !!
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE obs_fbm
    USE obs_sla_io
@@ -13 +31 @@
    CHARACTER(len=256) :: cdoutfile
    CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
+   CHARACTER(len=256) :: cdtmp
+   CHARACTER(len=5) :: cdsource
    !
    ! Input data
@@ -25 +45 @@
    ! Loop variables
    !
-   INTEGER :: ip,ia,ji,jk
+   INTEGER :: ip,ia,ji,jk,noff
    !
    ! Get number of command line arguments
@@ -32 +52 @@
    IF (nargs < 1) THEN
       WRITE(*,'(A)')'Usage:'
-      WRITE(*,'(A)')'sla2fb outputfile inputfile1 inputfile2 ...'
+      WRITE(*,'(A)')'sla2fb [-s type] outputfile inputfile1 inputfile2 ...'
       CALL abort()
    ENDIF
-   CALL getarg(1,cdoutfile)
+   cdsource=''
    !
    ! Get input data
    !
-   ALLOCATE( slaf(MAX(nargs-1,1)) )
-   ALLOCATE( cdinfile(nargs-1) )
+   noff=1
+   IF ( nargs > 1 ) THEN
+      CALL getarg(1,cdtmp)
+      IF (TRIM(cdtmp)=='-s') THEN
+         IF ( nargs < 3 ) THEN
+            WRITE(*,*)'Missing arguments to -s <datasource>'
+            CALL abort
+         ENDIF
+         CALL getarg(2,cdsource)
+         noff=3
+      ENDIF
+   ENDIF
+   CALL getarg(noff,cdoutfile)
+   ninfiles = nargs - noff
+   ALLOCATE( slaf(MAX(nargs-noff,1)) )
+   ALLOCATE( cdinfile(nargs-noff) )
    ntotobs = 0
-   ninfiles  = nargs - 1
    DO ia=1,ninfiles
-      CALL getarg( ia + 1, cdinfile(ia) )
+      CALL getarg( ia + noff, 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'
+      IF (LEN_TRIM(cdsource)>0) THEN
+         DO ji=1,slaf(ia)%nobs
+            slaf(ia)%cdwmo(ji)=TRIM(slaf(ia)%cdwmo(ji))//'_'//TRIM(cdsource)
+         ENDDO
+      ENDIF
       ntotobs = ntotobs + slaf(ia)%nobs
    ENDDO

branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/src/vel2fb.F90

@@ -1 +1 @@
 PROGRAM vel2fb
+   !!---------------------------------------------------------------------
+   !!
+   !!                     ** PROGRAM vel2fb **
+   !!
+   !!  ** Purpose : Convert TAO/PIRATA/RAMA currents to feedback format
+   !!
+   !!  ** Method  : Use of utilities from obs_fbm.
+   !!
+   !!  ** Action  :
+   !!
+   !!   Usage:
+   !!     vel2fb outputfile inputfile1 inputfile2 ...
+   !! 
+   !!   History :
+   !!        ! 2010 (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
    USE obs_fbm
    USE obs_vel_io