Index: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 10246)
+++ branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 10247)
@@ -7,12 +7,9 @@
 
    !!---------------------------------------------------------------------
-   !!   obs_pre_pro  : First level check and screening of T/S profiles
-   !!   obs_pre_sla  : First level check and screening of SLA observations
-   !!   obs_pre_sst  : First level check and screening of SLA observations
-   !!   obs_pre_seaice : First level check and screening of sea ice observations
-   !!   obs_pre_vel  : First level check and screening of velocity obs.
-   !!   obs_scr      : Basic screening of the observations
-   !!   obs_coo_tim  : Compute number of time steps to the observation time
-   !!   obs_sor      : Sort the observation arrays
+   !!   obs_pre_prof  : First level check and screening of profile observations
+   !!   obs_pre_surf  : First level check and screening of surface observations
+   !!   obs_scr       : Basic screening of the observations
+   !!   obs_coo_tim   : Compute number of time steps to the observation time
+   !!   obs_sor       : Sort the observation arrays
    !!---------------------------------------------------------------------
    !! * Modules used
@@ -27,4 +24,8 @@
    USE obs_inter_sup      ! Interpolation support
    USE obs_oper           ! Observation operators
+#if defined key_bdy
+   USE bdy_oce, ONLY : &        ! Boundary information
+      idx_bdy, nb_bdy
+#endif
    USE lib_mpp, ONLY : &
       & ctl_warn, ctl_stop
@@ -36,10 +37,7 @@
 
    PUBLIC &
-      & obs_pre_pro, &    ! First level check and screening of profiles
-      & obs_pre_sla, &    ! First level check and screening of SLA data
-      & obs_pre_sst, &    ! First level check and screening of SLA data
-      & obs_pre_seaice, & ! First level check and screening of sea ice data
-      & obs_pre_vel, &     ! First level check and screening of velocity profiles
-      & calc_month_len     ! Calculate the number of days in the months of a year  
+      & obs_pre_prof, &    ! First level check and screening of profile obs
+      & obs_pre_surf, &    ! First level check and screening of surface obs
+      & calc_month_len     ! Calculate the number of days in the months of a year
 
    !!----------------------------------------------------------------------
@@ -49,299 +47,17 @@
    !!----------------------------------------------------------------------
 
+!! * Substitutions 
+#  include "domzgr_substitute.h90"  
+
 CONTAINS
 
-   SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, &
-      &                    kdailyavtypes )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE obs_pre_pro  ***
-      !!
-      !! ** Purpose : First level check and screening of T and S profiles
-      !!
-      !! ** Method  : First level check and screening of T and S profiles
-      !!
-      !! ** Action  : 
-      !!
-      !! References :
-      !!   
-      !! History :
-      !!        !  2007-01  (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d 
-      !!        !  2007-03  (K. Mogensen) General handling of profiles
-      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE domstp              ! Domain: set the time-step
-      USE par_oce             ! Ocean parameters
-      USE dom_oce, ONLY : &   ! Geographical information
-         & glamt,   &
-         & gphit,   &
-         & gdept_1d,&
-         & tmask,   &
-         & nproc
-      !! * Arguments
-      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Full set of profile data
-      TYPE(obs_prof), INTENT(INOUT) :: prodatqc     ! Subset of profile data not failing screening
-      LOGICAL, INTENT(IN) :: ld_t3d         ! Switch for temperature
-      LOGICAL, INTENT(IN) :: ld_s3d         ! Switch for salinity
-      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land
-      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
-         & kdailyavtypes! Types for daily averages
-      !! * Local declarations   
-      INTEGER :: iyea0         ! Initial date
-      INTEGER :: imon0         !  - (year, month, day, hour, minute)
-      INTEGER :: iday0   
-      INTEGER :: ihou0
-      INTEGER :: imin0
-      INTEGER :: icycle        ! Current assimilation cycle
-                               ! Counters for observations that
-      INTEGER :: iotdobs       !  - outside time domain
-      INTEGER :: iosdtobs      !  - outside space domain (temperature)
-      INTEGER :: iosdsobs      !  - outside space domain (salinity)
-      INTEGER :: ilantobs      !  - within a model land cell (temperature)
-      INTEGER :: ilansobs      !  - within a model land cell (salinity)
-      INTEGER :: inlatobs      !  - close to land (temperature)
-      INTEGER :: inlasobs      !  - close to land (salinity)
-      INTEGER :: igrdobs       !  - fail the grid search
-                               ! Global counters for observations that
-      INTEGER :: iotdobsmpp    !  - outside time domain
-      INTEGER :: iosdtobsmpp   !  - outside space domain (temperature)
-      INTEGER :: iosdsobsmpp   !  - outside space domain (salinity)
-      INTEGER :: ilantobsmpp   !  - within a model land cell (temperature)
-      INTEGER :: ilansobsmpp   !  - within a model land cell (salinity)
-      INTEGER :: inlatobsmpp   !  - close to land (temperature)
-      INTEGER :: inlasobsmpp   !  - close to land (salinity)
-      INTEGER :: igrdobsmpp    !  - fail the grid search
-      TYPE(obs_prof_valid) ::  llvalid     ! Profile selection 
-      TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
-         & llvvalid            ! T,S selection 
-      INTEGER :: jvar          ! Variable loop variable
-      INTEGER :: jobs          ! Obs. loop variable
-      INTEGER :: jstp          ! Time loop variable
-      INTEGER :: inrc          ! Time index variable
-      
-      IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...'
-
-      ! Initial date initialization (year, month, day, hour, minute)
-      iyea0 =   ndate0 / 10000
-      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
-      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
-      ihou0 = 0
-      imin0 = 0
-
-      icycle = no     ! Assimilation cycle
-
-      ! Diagnotics counters for various failures.
-
-      iotdobs  = 0
-      igrdobs  = 0
-      iosdtobs = 0
-      iosdsobs = 0
-      ilantobs = 0
-      ilansobs = 0
-      inlatobs = 0
-      inlasobs = 0
-
-      ! -----------------------------------------------------------------------
-      ! Find time coordinate for profiles
-      ! -----------------------------------------------------------------------
-
-      IF ( PRESENT(kdailyavtypes) ) THEN
-         CALL obs_coo_tim_prof( icycle, &
-            &                iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-            &                profdata%nprof,   profdata%nyea, profdata%nmon, &
-            &                profdata%nday,    profdata%nhou, profdata%nmin, &
-            &                profdata%ntyp,    profdata%nqc,  profdata%mstp, &
-            &                iotdobs, kdailyavtypes = kdailyavtypes        )
-      ELSE
-         CALL obs_coo_tim_prof( icycle, &
-            &                iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-            &                profdata%nprof,   profdata%nyea, profdata%nmon, &
-            &                profdata%nday,    profdata%nhou, profdata%nmin, &
-            &                profdata%ntyp,    profdata%nqc,  profdata%mstp, &
-            &                iotdobs )
-      ENDIF
-      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
-      
-      ! -----------------------------------------------------------------------
-      ! Check for profiles failing the grid search
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi, profdata%mj, &
-         &              profdata%nqc,     igrdobs                         )
-
-      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
-
-      ! -----------------------------------------------------------------------
-      ! Reject all observations for profiles with nqc > 10
-      ! -----------------------------------------------------------------------
-
-      CALL obs_pro_rej( profdata )
-
-      ! -----------------------------------------------------------------------
-      ! Check for land points. This includes points below the model
-      ! bathymetry so this is done for every point in the profile
-      ! -----------------------------------------------------------------------
-
-      ! Temperature
-
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   &
-         &                 profdata%npvsta(:,1),  profdata%npvend(:,1), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          & 
-         &                 profdata%var(1)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(1)%vdep,                        &
-         &                 glamt,                 gphit,                &
-         &                 gdept_1d,              tmask,                &
-         &                 profdata%nqc,          profdata%var(1)%nvqc, &
-         &                 iosdtobs,              ilantobs,             &
-         &                 inlatobs,              ld_nea                )
-
-      CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp )
-      CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp )
-      CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp )
-
-      ! Salinity
-
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   &
-         &                 profdata%npvsta(:,2),  profdata%npvend(:,2), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          & 
-         &                 profdata%var(2)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(2)%vdep,                        &
-         &                 glamt,                 gphit,                &
-         &                 gdept_1d,              tmask,                &
-         &                 profdata%nqc,          profdata%var(2)%nvqc, &
-         &                 iosdsobs,              ilansobs,             &
-         &                 inlasobs,              ld_nea                )
-
-      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
-      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
-      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
-
-      ! -----------------------------------------------------------------------
-      ! Copy useful data from the profdata data structure to
-      ! the prodatqc data structure 
-      ! -----------------------------------------------------------------------
-
-      ! Allocate the selection arrays
-
-      ALLOCATE( llvalid%luse(profdata%nprof) )
-      DO jvar = 1,profdata%nvar
-         ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
-      END DO
-
-      ! We want all data which has qc flags <= 10
-
-      llvalid%luse(:) = ( profdata%nqc(:)  <= 10 )
-      DO jvar = 1,profdata%nvar
-         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
-      END DO
-
-      ! The actual copying
-
-      CALL obs_prof_compress( profdata,     prodatqc,       .TRUE.,  numout, &
-         &                    lvalid=llvalid, lvvalid=llvvalid )
-
-      ! Dellocate the selection arrays
-      DEALLOCATE( llvalid%luse )
-      DO jvar = 1,profdata%nvar
-         DEALLOCATE( llvvalid(jvar)%luse )
-      END DO
-
-      ! -----------------------------------------------------------------------
-      ! Print information about what observations are left after qc
-      ! -----------------------------------------------------------------------
-
-      ! Update the total observation counter array
-      
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'obs_pre_pro :'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*)
-         WRITE(numout,*) ' Profiles outside time domain                = ', &
-            &            iotdobsmpp
-         WRITE(numout,*) ' Remaining profiles that failed grid search  = ', &
-            &            igrdobsmpp
-         WRITE(numout,*) ' Remaining T data outside space domain       = ', &
-            &            iosdtobsmpp
-         WRITE(numout,*) ' Remaining T data at land points             = ', &
-            &            ilantobsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining T data near land points (removed) = ',&
-               &            inlatobsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining T data near land points (kept)    = ',&
-               &            inlatobsmpp
-         ENDIF
-         WRITE(numout,*) ' T data accepted                             = ', &
-            &            prodatqc%nvprotmpp(1)
-         WRITE(numout,*) ' Remaining S data outside space domain       = ', &
-            &            iosdsobsmpp
-         WRITE(numout,*) ' Remaining S data at land points             = ', &
-            &            ilansobsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining S data near land points (removed) = ',&
-               &            inlasobsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining S data near land points (kept)    = ',&
-               &            inlasobsmpp
-         ENDIF
-         WRITE(numout,*) ' S data accepted                             = ', &
-            &            prodatqc%nvprotmpp(2)
-
-         WRITE(numout,*)
-         WRITE(numout,*) ' Number of observations per time step :'
-         WRITE(numout,*)
-         WRITE(numout,997)
-         WRITE(numout,998)
-      ENDIF
-      
-      DO jobs = 1, prodatqc%nprof
-         inrc = prodatqc%mstp(jobs) + 2 - nit000
-         prodatqc%npstp(inrc)  = prodatqc%npstp(inrc) + 1
-         DO jvar = 1, prodatqc%nvar
-            IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
-               prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
-                  &                      ( prodatqc%npvend(jobs,jvar) - &
-                  &                        prodatqc%npvsta(jobs,jvar) + 1 )
-            ENDIF
-         END DO
-      END DO
-      
-      
-      CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
-         &                       nitend - nit000 + 2 )
-      DO jvar = 1, prodatqc%nvar
-         CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
-            &                       prodatqc%nvstpmpp(:,jvar), &
-            &                       nitend - nit000 + 2 )
-      END DO
-
-      IF ( lwp ) THEN
-         DO jstp = nit000 - 1, nitend
-            inrc = jstp - nit000 + 2
-            WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
-               &                    prodatqc%nvstpmpp(inrc,1), &
-               &                    prodatqc%nvstpmpp(inrc,2)
-         END DO
-      ENDIF
-
-997   FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity')
-998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------')
-999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
-      
-   END SUBROUTINE obs_pre_pro
-
-   SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea )
+   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, &
+                            kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_pre_sla  ***
       !!
-      !! ** Purpose : First level check and screening of SLA observations
-      !!
-      !! ** Method  : First level check and screening of SLA observations
+      !! ** Purpose : First level check and screening of surface observations
+      !!
+      !! ** Method  : First level check and screening of surface observations
       !!
       !! ** Action  : 
@@ -352,4 +68,5 @@
       !!        !  2007-03  (A. Weaver, K. Mogensen) Original
       !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
+      !!        !  2015-02  (M. Martin) Combined routine for surface types.
       !!----------------------------------------------------------------------
       !! * Modules used
@@ -362,9 +79,11 @@
          & nproc
       !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: sladata    ! Full set of SLA data
-      TYPE(obs_surf), INTENT(INOUT) :: sladatqc   ! Subset of SLA data not failing screening
-      LOGICAL, INTENT(IN) :: ld_sla         ! Switch for SLA data
-      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land
+      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data
+      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc  ! Subset of surface data not failing screening
+      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land
+      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary
+      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
       !! * Local declarations
+      INTEGER :: iqc_cutoff = 255   ! cut off for QC value
       INTEGER :: iyea0        ! Initial date
       INTEGER :: imon0        !  - (year, month, day, hour, minute)
@@ -379,4 +98,5 @@
       INTEGER :: inlasobs     !  - close to land
       INTEGER :: igrdobs      !  - fail the grid search
+      INTEGER :: ibdysobs     !  - close to open boundary
                               ! Global counters for observations that
       INTEGER :: iotdobsmpp     !  - outside time domain
@@ -385,4 +105,5 @@
       INTEGER :: inlasobsmpp    !  - close to land
       INTEGER :: igrdobsmpp     !  - fail the grid search
+      INTEGER :: ibdysobsmpp  !  - close to open boundary
       LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
          & llvalid            ! SLA data selection
@@ -391,6 +112,7 @@
       INTEGER :: inrc         ! Time index variable
 
-      IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...'
-
+      IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      
       ! Initial date initialization (year, month, day, hour, minute)
       iyea0 =   ndate0 / 10000
@@ -409,23 +131,27 @@
       ilansobs = 0
       inlasobs = 0
-
-      ! -----------------------------------------------------------------------
-      ! Find time coordinate for SLA data
+      ibdysobs = 0 
+
+      ! Set QC cutoff to optional value if provided
+      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
+
+      ! -----------------------------------------------------------------------
+      ! Find time coordinate for surface data
       ! -----------------------------------------------------------------------
 
       CALL obs_coo_tim( icycle, &
          &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-         &              sladata%nsurf,   sladata%nyea, sladata%nmon, &
-         &              sladata%nday,    sladata%nhou, sladata%nmin, &
-         &              sladata%nqc,     sladata%mstp, iotdobs        )
+         &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, &
+         &              surfdata%nday,    surfdata%nhou, surfdata%nmin, &
+         &              surfdata%nqc,     surfdata%mstp, iotdobs        )
 
       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
       
       ! -----------------------------------------------------------------------
-      ! Check for SLA data failing the grid search
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_grd( sladata%nsurf,   sladata%mi, sladata%mj, &
-         &              sladata%nqc,     igrdobs                         )
+      ! Check for surface data failing the grid search
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_grd( surfdata%nsurf,   surfdata%mi, surfdata%mj, &
+         &              surfdata%nqc,     igrdobs                         )
 
       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
@@ -435,33 +161,36 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_spc_2d( sladata%nsurf,              &
+      CALL obs_coo_spc_2d( surfdata%nsurf,              &
          &                 jpi,          jpj,          &
-         &                 sladata%mi,   sladata%mj,   & 
-         &                 sladata%rlam, sladata%rphi, &
+         &                 surfdata%mi,   surfdata%mj,   & 
+         &                 surfdata%rlam, surfdata%rphi, &
          &                 glamt,        gphit,        &
-         &                 tmask(:,:,1), sladata%nqc,  &
+         &                 tmask(:,:,1), surfdata%nqc,  &
          &                 iosdsobs,     ilansobs,     &
-         &                 inlasobs,     ld_nea        )
+         &                 inlasobs,     ld_nea,       &
+         &                 ibdysobs,     ld_bound_reject, &
+         &                 iqc_cutoff                     )
 
       CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
       CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
       CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
-
-      ! -----------------------------------------------------------------------
-      ! Copy useful data from the sladata data structure to
-      ! the sladatqc data structure 
+      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Copy useful data from the surfdata data structure to
+      ! the surfdataqc data structure 
       ! -----------------------------------------------------------------------
 
       ! Allocate the selection arrays
 
-      ALLOCATE( llvalid(sladata%nsurf) )
-      
-      ! We want all data which has qc flags <= 10
-
-      llvalid(:)  = ( sladata%nqc(:)  <= 10 )
+      ALLOCATE( llvalid(surfdata%nsurf) )
+      
+      ! We want all data which has qc flags <= iqc_cutoff
+
+      llvalid(:)  = ( surfdata%nqc(:)  <= iqc_cutoff )
 
       ! The actual copying
 
-      CALL obs_surf_compress( sladata,     sladatqc,       .TRUE.,  numout, &
+      CALL obs_surf_compress( surfdata,     surfdataqc,       .TRUE.,  numout, &
          &                    lvalid=llvalid )
 
@@ -477,38 +206,38 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'obs_pre_sla :'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*)
-         WRITE(numout,*) ' SLA data outside time domain                  = ', &
+         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain                  = ', &
             &            iotdobsmpp
-         WRITE(numout,*) ' Remaining SLA data that failed grid search    = ', &
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search    = ', &
             &            igrdobsmpp
-         WRITE(numout,*) ' Remaining SLA data outside space domain       = ', &
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain       = ', &
             &            iosdsobsmpp
-         WRITE(numout,*) ' Remaining SLA data at land points             = ', &
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points             = ', &
             &            ilansobsmpp
          IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining SLA data near land points (removed) = ', &
+            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', &
                &            inlasobsmpp
          ELSE
-            WRITE(numout,*) ' Remaining SLA data near land points (kept)    = ', &
+            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept)    = ', &
                &            inlasobsmpp
          ENDIF
-         WRITE(numout,*) ' SLA data accepted                             = ', &
-            &            sladatqc%nsurfmpp
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', &
+            &            ibdysobsmpp  
+         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', &
+            &            surfdataqc%nsurfmpp
 
          WRITE(numout,*)
          WRITE(numout,*) ' Number of observations per time step :'
          WRITE(numout,*)
-         WRITE(numout,1997)
-         WRITE(numout,1998)
+         WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1)
+         WRITE(numout,'(10X,A,5X,A)')'---------','-----------------'
+         CALL FLUSH(numout)
       ENDIF
       
-      DO jobs = 1, sladatqc%nsurf
-         inrc = sladatqc%mstp(jobs) + 2 - nit000
-         sladatqc%nsstp(inrc)  = sladatqc%nsstp(inrc) + 1
-      END DO
-      
-      CALL obs_mpp_sum_integers( sladatqc%nsstp, sladatqc%nsstpmpp, &
+      DO jobs = 1, surfdataqc%nsurf
+         inrc = surfdataqc%mstp(jobs) + 2 - nit000
+         surfdataqc%nsstp(inrc)  = surfdataqc%nsstp(inrc) + 1
+      END DO
+      
+      CALL obs_mpp_sum_integers( surfdataqc%nsstp, surfdataqc%nsstpmpp, &
          &                       nitend - nit000 + 2 )
 
@@ -516,28 +245,32 @@
          DO jstp = nit000 - 1, nitend
             inrc = jstp - nit000 + 2
-            WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc)
+            WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc)
+            CALL FLUSH(numout)
          END DO
       ENDIF
 
-1997  FORMAT(10X,'Time step',5X,'Sea level anomaly')
-1998  FORMAT(10X,'---------',5X,'-----------------')
 1999  FORMAT(10X,I9,5X,I17)
 
-   END SUBROUTINE obs_pre_sla
-
-   SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE obs_pre_sst  ***
-      !!
-      !! ** Purpose : First level check and screening of SST observations
-      !!
-      !! ** Method  : First level check and screening of SST observations
-      !!
-      !! ** Action  : 
-      !!
-      !! References :
-      !!   
+   END SUBROUTINE obs_pre_surf
+
+
+   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, &
+      &                     kpi, kpj, kpk, &
+      &                     zmask, pglam, pgphi,  &
+      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff )
+
+!!----------------------------------------------------------------------
+      !!                    ***  ROUTINE obs_pre_prof  ***
+      !!
+      !! ** Purpose : First level check and screening of profiles
+      !!
+      !! ** Method  : First level check and screening of profiles
+      !!
       !! History :
-      !!        !  2007-03  (S. Ricci) SST data preparation 
+      !!        !  2007-06  (K. Mogensen) original : T and S profile data
+      !!        !  2008-09  (M. Valdivieso) : TAO velocity data
+      !!        !  2009-01  (K. Mogensen) : New feedback stricture
+      !!        !  2015-02  (M. Martin) : Combined profile routine.
+      !!
       !!----------------------------------------------------------------------
       !! * Modules used
@@ -545,39 +278,62 @@
       USE par_oce             ! Ocean parameters
       USE dom_oce, ONLY : &   ! Geographical information
-         & glamt,   &
-         & gphit,   &
-         & tmask,   &
+         & gdept_1d,             &
          & nproc
+
       !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: sstdata     ! Full set of SST data
-      TYPE(obs_surf), INTENT(INOUT) :: sstdatqc    ! Subset of SST data not failing screening
-      LOGICAL :: ld_sst             ! Switch for SST data
-      LOGICAL :: ld_nea             ! Switch for rejecting observation near land
+      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data
+      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening
+      LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: &
+         & ld_var                                 ! Observed variables switches
+      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land
+      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary
+      INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes
+      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
+         & kdailyavtypes                          ! Types for daily averages
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: &
+         & zmask
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: &
+         & pglam, &
+         & pgphi
+      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
+
       !! * Local declarations
+      INTEGER :: iqc_cutoff = 255   ! cut off for QC value
       INTEGER :: iyea0        ! Initial date
       INTEGER :: imon0        !  - (year, month, day, hour, minute)
-      INTEGER :: iday0   
+      INTEGER :: iday0    
       INTEGER :: ihou0    
       INTEGER :: imin0
       INTEGER :: icycle       ! Current assimilation cycle
-                              ! Counters for observations that
-      INTEGER :: iotdobs      !  - outside time domain
-      INTEGER :: iosdsobs     !  - outside space domain
-      INTEGER :: ilansobs     !  - within a model land cell
-      INTEGER :: inlasobs     !  - close to land
-      INTEGER :: igrdobs      !  - fail the grid search
-                              ! Global counters for observations that
-      INTEGER :: iotdobsmpp   !  - outside time domain
-      INTEGER :: iosdsobsmpp  !  - outside space domain
-      INTEGER :: ilansobsmpp  !  - within a model land cell
-      INTEGER :: inlasobsmpp  !  - close to land
-      INTEGER :: igrdobsmpp   !  - fail the grid search
-      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
-         & llvalid            ! SST data selection
+                                                       ! Counters for observations that are
+      INTEGER                           :: iotdobs     !  - outside time domain
+      INTEGER, DIMENSION(profdata%nvar) :: iosdvobs    !  - outside space domain
+      INTEGER, DIMENSION(profdata%nvar) :: ilanvobs    !  - within a model land cell
+      INTEGER, DIMENSION(profdata%nvar) :: inlavobs    !  - close to land
+      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs    !  - boundary   
+      INTEGER                           :: igrdobs     !  - fail the grid search
+      INTEGER                           :: iuvchku     !  - reject UVEL if VVEL rejected
+      INTEGER                           :: iuvchkv     !  - reject VVEL if UVEL rejected
+                                                       ! Global counters for observations that are
+      INTEGER                           :: iotdobsmpp  !  - outside time domain
+      INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp !  - outside space domain
+      INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp !  - within a model land cell
+      INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp !  - close to land
+      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp !  - boundary
+      INTEGER :: igrdobsmpp                            !  - fail the grid search
+      INTEGER :: iuvchkumpp                            !  - reject UVEL if VVEL rejected
+      INTEGER :: iuvchkvmpp                            !  - reject VVEL if UVEL rejected
+      TYPE(obs_prof_valid) ::  llvalid      ! Profile selection 
+      TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
+         & llvvalid           ! vars selection 
+      INTEGER :: jvar         ! Variable loop variable
       INTEGER :: jobs         ! Obs. loop variable
       INTEGER :: jstp         ! Time loop variable
       INTEGER :: inrc         ! Time index variable
-
-      IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...'
+      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line
+      CHARACTER(LEN=256) :: cout2  ! Diagnostic output line
+
+      IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
 
       ! Initial date initialization (year, month, day, hour, minute)
@@ -590,429 +346,58 @@
       icycle = no     ! Assimilation cycle
 
-      ! Diagnotics counters for various failures.
-
-      iotdobs  = 0
-      igrdobs  = 0
-      iosdsobs = 0
-      ilansobs = 0
-      inlasobs = 0
-
-      ! -----------------------------------------------------------------------
-      ! Find time coordinate for SST data
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_tim( icycle, &
-         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-         &              sstdata%nsurf,   sstdata%nyea, sstdata%nmon, &
-         &              sstdata%nday,    sstdata%nhou, sstdata%nmin, &
-         &              sstdata%nqc,     sstdata%mstp, iotdobs        )
+      ! Diagnostics counters for various failures.
+
+      iotdobs     = 0
+      igrdobs     = 0
+      iosdvobs(:) = 0
+      ilanvobs(:) = 0
+      inlavobs(:) = 0
+      ibdyvobs(:) = 0
+      iuvchku     = 0
+      iuvchkv     = 0
+
+
+      ! Set QC cutoff to optional value if provided
+      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
+
+      ! -----------------------------------------------------------------------
+      ! Find time coordinate for profiles
+      ! -----------------------------------------------------------------------
+
+      IF ( PRESENT(kdailyavtypes) ) THEN
+         CALL obs_coo_tim_prof( icycle, &
+            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+            &              profdata%nprof,   profdata%nyea, profdata%nmon, &
+            &              profdata%nday,    profdata%nhou, profdata%nmin, &
+            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
+            &              iotdobs, kdailyavtypes = kdailyavtypes,         &
+            &              kqc_cutoff = iqc_cutoff )
+      ELSE
+         CALL obs_coo_tim_prof( icycle, &
+            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+            &              profdata%nprof,   profdata%nyea, profdata%nmon, &
+            &              profdata%nday,    profdata%nhou, profdata%nmin, &
+            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
+            &              iotdobs,          kqc_cutoff = iqc_cutoff )
+      ENDIF
+
       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
-      ! -----------------------------------------------------------------------
-      ! Check for SST data failing the grid search
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_grd( sstdata%nsurf,   sstdata%mi, sstdata%mj, &
-         &              sstdata%nqc,     igrdobs                         )
+      
+      ! -----------------------------------------------------------------------
+      ! Check for profiles failing the grid search
+      ! -----------------------------------------------------------------------
+
+      DO jvar = 1, profdata%nvar
+         CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,jvar), profdata%mj(:,jvar), &
+            &              profdata%nqc,     igrdobs )
+      END DO
+
       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
 
       ! -----------------------------------------------------------------------
-      ! Check for land points. 
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_spc_2d( sstdata%nsurf,              &
-         &                 jpi,          jpj,          &
-         &                 sstdata%mi,   sstdata%mj,   & 
-         &                 sstdata%rlam, sstdata%rphi, &
-         &                 glamt,        gphit,        &
-         &                 tmask(:,:,1), sstdata%nqc,  &
-         &                 iosdsobs,     ilansobs,     &
-         &                 inlasobs,     ld_nea        )
-
-      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
-      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
-      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
-
-      ! -----------------------------------------------------------------------
-      ! Copy useful data from the sstdata data structure to
-      ! the sstdatqc data structure 
-      ! -----------------------------------------------------------------------
-
-      ! Allocate the selection arrays
-
-      ALLOCATE( llvalid(sstdata%nsurf) )
-      
-      ! We want all data which has qc flags <= 0
-
-      llvalid(:)  = ( sstdata%nqc(:)  <= 10 )
-
-      ! The actual copying
-
-      CALL obs_surf_compress( sstdata,     sstdatqc,       .TRUE.,  numout, &
-         &                    lvalid=llvalid )
-
-      ! Dellocate the selection arrays
-      DEALLOCATE( llvalid )
-
-      ! -----------------------------------------------------------------------
-      ! Print information about what observations are left after qc
-      ! -----------------------------------------------------------------------
-
-      ! Update the total observation counter array
-      
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'obs_pre_sst :'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*)
-         WRITE(numout,*) ' SST data outside time domain                  = ', &
-            &            iotdobsmpp
-         WRITE(numout,*) ' Remaining SST data that failed grid search    = ', &
-            &            igrdobsmpp
-         WRITE(numout,*) ' Remaining SST data outside space domain       = ', &
-            &            iosdsobsmpp
-         WRITE(numout,*) ' Remaining SST data at land points             = ', &
-            &            ilansobsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining SST data near land points (removed) = ', &
-               &            inlasobsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining SST data near land points (kept)    = ', &
-               &            inlasobsmpp
-         ENDIF
-         WRITE(numout,*) ' SST data accepted                             = ', &
-            &            sstdatqc%nsurfmpp
-
-         WRITE(numout,*)
-         WRITE(numout,*) ' Number of observations per time step :'
-         WRITE(numout,*)
-         WRITE(numout,1997)
-         WRITE(numout,1998)
-      ENDIF
-      
-      DO jobs = 1, sstdatqc%nsurf
-         inrc = sstdatqc%mstp(jobs) + 2 - nit000
-         sstdatqc%nsstp(inrc)  = sstdatqc%nsstp(inrc) + 1
-      END DO
-      
-      CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, &
-         &                       nitend - nit000 + 2 )
-
-      IF ( lwp ) THEN
-         DO jstp = nit000 - 1, nitend
-            inrc = jstp - nit000 + 2
-            WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc)
-         END DO
-      ENDIF
-
-1997  FORMAT(10X,'Time step',5X,'Sea surface temperature')
-1998  FORMAT(10X,'---------',5X,'-----------------')
-1999  FORMAT(10X,I9,5X,I17)
-      
-   END SUBROUTINE obs_pre_sst
-
-   SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE obs_pre_seaice  ***
-      !!
-      !! ** Purpose : First level check and screening of Sea Ice observations
-      !!
-      !! ** Method  : First level check and screening of Sea Ice observations
-      !!
-      !! ** Action  : 
-      !!
-      !! References :
-      !!   
-      !! History :
-      !!        !  2007-11 (D. Lea) based on obs_pre_sst
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE domstp              ! Domain: set the time-step
-      USE par_oce             ! Ocean parameters
-      USE dom_oce, ONLY : &   ! Geographical information
-         & glamt,   &
-         & gphit,   &
-         & tmask,   &
-         & nproc
-      !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: seaicedata     ! Full set of Sea Ice data
-      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc    ! Subset of sea ice data not failing screening
-      LOGICAL :: ld_seaice     ! Switch for sea ice data
-      LOGICAL :: ld_nea        ! Switch for rejecting observation near land
-      !! * Local declarations
-      INTEGER :: iyea0         ! Initial date
-      INTEGER :: imon0         !  - (year, month, day, hour, minute)
-      INTEGER :: iday0    
-      INTEGER :: ihou0    
-      INTEGER :: imin0
-      INTEGER :: icycle       ! Current assimilation cycle
-                              ! Counters for observations that
-      INTEGER :: iotdobs      !  - outside time domain
-      INTEGER :: iosdsobs     !  - outside space domain
-      INTEGER :: ilansobs     !  - within a model land cell
-      INTEGER :: inlasobs     !  - close to land
-      INTEGER :: igrdobs      !  - fail the grid search
-                              ! Global counters for observations that
-      INTEGER :: iotdobsmpp   !  - outside time domain
-      INTEGER :: iosdsobsmpp  !  - outside space domain
-      INTEGER :: ilansobsmpp  !  - within a model land cell
-      INTEGER :: inlasobsmpp  !  - close to land
-      INTEGER :: igrdobsmpp   !  - fail the grid search
-      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
-         & llvalid            ! data selection
-      INTEGER :: jobs         ! Obs. loop variable
-      INTEGER :: jstp         ! Time loop variable
-      INTEGER :: inrc         ! Time index variable
-
-      IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...'
-
-      ! Initial date initialization (year, month, day, hour, minute)
-      iyea0 =   ndate0 / 10000
-      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
-      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
-      ihou0 = 0
-      imin0 = 0
-
-      icycle = no     ! Assimilation cycle
-
-      ! Diagnotics counters for various failures.
-
-      iotdobs  = 0
-      igrdobs  = 0
-      iosdsobs = 0
-      ilansobs = 0
-      inlasobs = 0
-
-      ! -----------------------------------------------------------------------
-      ! Find time coordinate for sea ice data
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_tim( icycle, &
-         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-         &              seaicedata%nsurf,   seaicedata%nyea, seaicedata%nmon, &
-         &              seaicedata%nday,    seaicedata%nhou, seaicedata%nmin, &
-         &              seaicedata%nqc,     seaicedata%mstp, iotdobs        )
-      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
-      ! -----------------------------------------------------------------------
-      ! Check for sea ice data failing the grid search
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_grd( seaicedata%nsurf,   seaicedata%mi, seaicedata%mj, &
-         &              seaicedata%nqc,     igrdobs                         )
-      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
-
-      ! -----------------------------------------------------------------------
-      ! Check for land points. 
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_spc_2d( seaicedata%nsurf,                 &
-         &                 jpi,             jpj,             &
-         &                 seaicedata%mi,   seaicedata%mj,   & 
-         &                 seaicedata%rlam, seaicedata%rphi, &
-         &                 glamt,           gphit,           &
-         &                 tmask(:,:,1),    seaicedata%nqc,  &
-         &                 iosdsobs,        ilansobs,        &
-         &                 inlasobs,        ld_nea           )
-
-      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
-      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
-      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
-
-      ! -----------------------------------------------------------------------
-      ! Copy useful data from the seaicedata data structure to
-      ! the seaicedatqc data structure 
-      ! -----------------------------------------------------------------------
-
-      ! Allocate the selection arrays
-
-      ALLOCATE( llvalid(seaicedata%nsurf) )
-      
-      ! We want all data which has qc flags <= 0
-
-      llvalid(:)  = ( seaicedata%nqc(:)  <= 10 )
-
-      ! The actual copying
-
-      CALL obs_surf_compress( seaicedata,     seaicedatqc,       .TRUE.,  numout, &
-         &                    lvalid=llvalid )
-
-      ! Dellocate the selection arrays
-      DEALLOCATE( llvalid )
-
-      ! -----------------------------------------------------------------------
-      ! Print information about what observations are left after qc
-      ! -----------------------------------------------------------------------
-
-      ! Update the total observation counter array
-      
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'obs_pre_seaice :'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*)
-         WRITE(numout,*) ' Sea ice data outside time domain                  = ', &
-            &            iotdobsmpp
-         WRITE(numout,*) ' Remaining sea ice data that failed grid search    = ', &
-            &            igrdobsmpp
-         WRITE(numout,*) ' Remaining sea ice data outside space domain       = ', &
-            &            iosdsobsmpp
-         WRITE(numout,*) ' Remaining sea ice data at land points             = ', &
-            &            ilansobsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', &
-               &            inlasobsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining sea ice data near land points (kept)    = ', &
-               &            inlasobsmpp
-         ENDIF
-         WRITE(numout,*) ' Sea ice data accepted                             = ', &
-            &            seaicedatqc%nsurfmpp
-
-         WRITE(numout,*)
-         WRITE(numout,*) ' Number of observations per time step :'
-         WRITE(numout,*)
-         WRITE(numout,1997)
-         WRITE(numout,1998)
-      ENDIF
-      
-      DO jobs = 1, seaicedatqc%nsurf
-         inrc = seaicedatqc%mstp(jobs) + 2 - nit000
-         seaicedatqc%nsstp(inrc)  = seaicedatqc%nsstp(inrc) + 1
-      END DO
-      
-      CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, &
-         &                       nitend - nit000 + 2 )
-
-      IF ( lwp ) THEN
-         DO jstp = nit000 - 1, nitend
-            inrc = jstp - nit000 + 2
-            WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc)
-         END DO
-      ENDIF
-
-1997  FORMAT(10X,'Time step',5X,'Sea ice data           ')
-1998  FORMAT(10X,'---------',5X,'-----------------')
-1999  FORMAT(10X,I9,5X,I17)
-      
-   END SUBROUTINE obs_pre_seaice
-
-   SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE obs_pre_taovel  ***
-      !!
-      !! ** Purpose : First level check and screening of U and V profiles
-      !!
-      !! ** Method  : First level check and screening of U and V profiles
-      !!
-      !! History :
-      !!        !  2007-06  (K. Mogensen) original : T and S profile data
-      !!        !  2008-09  (M. Valdivieso) : TAO velocity data
-      !!        !  2009-01  (K. Mogensen) : New feedback strictuer
-      !!
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE domstp              ! Domain: set the time-step
-      USE par_oce             ! Ocean parameters
-      USE dom_oce, ONLY : &   ! Geographical information
-         & glamt, glamu, glamv,    &
-         & gphit, gphiu, gphiv,    &
-         & gdept_1d,             &
-         & tmask, umask, vmask,  &
-         & nproc
-      !! * Arguments
-      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data
-      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening
-      LOGICAL, INTENT(IN) :: ld_vel3d      ! Switch for zonal and meridional velocity components
-      LOGICAL, INTENT(IN) :: ld_nea        ! Switch for rejecting observation near land
-      LOGICAL, INTENT(IN) :: ld_dailyav    ! Switch for daily average data
-      !! * Local declarations
-      INTEGER :: iyea0        ! Initial date
-      INTEGER :: imon0        !  - (year, month, day, hour, minute)
-      INTEGER :: iday0    
-      INTEGER :: ihou0    
-      INTEGER :: imin0
-      INTEGER :: icycle       ! Current assimilation cycle
-                              ! Counters for observations that
-      INTEGER :: iotdobs      !  - outside time domain
-      INTEGER :: iosduobs     !  - outside space domain (zonal velocity component)
-      INTEGER :: iosdvobs     !  - outside space domain (meridional velocity component)
-      INTEGER :: ilanuobs     !  - within a model land cell (zonal velocity component)
-      INTEGER :: ilanvobs     !  - within a model land cell (meridional velocity component)
-      INTEGER :: inlauobs     !  - close to land (zonal velocity component)
-      INTEGER :: inlavobs     !  - close to land (meridional velocity component)
-      INTEGER :: igrdobs      !  - fail the grid search
-      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa
-      INTEGER :: iuvchkv      !
-                              ! Global counters for observations that
-      INTEGER :: iotdobsmpp   !  - outside time domain
-      INTEGER :: iosduobsmpp  !  - outside space domain (zonal velocity component)
-      INTEGER :: iosdvobsmpp  !  - outside space domain (meridional velocity component)
-      INTEGER :: ilanuobsmpp  !  - within a model land cell (zonal velocity component)
-      INTEGER :: ilanvobsmpp  !  - within a model land cell (meridional velocity component)
-      INTEGER :: inlauobsmpp  !  - close to land (zonal velocity component)
-      INTEGER :: inlavobsmpp  !  - close to land (meridional velocity component)
-      INTEGER :: igrdobsmpp   !  - fail the grid search
-      INTEGER :: iuvchkumpp   !  - reject u if v rejected and vice versa
-      INTEGER :: iuvchkvmpp   !
-      TYPE(obs_prof_valid) ::  llvalid      ! Profile selection 
-      TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
-         & llvvalid           ! U,V selection 
-      INTEGER :: jvar         ! Variable loop variable
-      INTEGER :: jobs         ! Obs. loop variable
-      INTEGER :: jstp         ! Time loop variable
-      INTEGER :: inrc         ! Time index variable
-
-      IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data'
-
-      ! Initial date initialization (year, month, day, hour, minute)
-      iyea0 =   ndate0 / 10000
-      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
-      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
-      ihou0 = 0
-      imin0 = 0
-
-      icycle = no     ! Assimilation cycle
-
-      ! Diagnotics counters for various failures.
-
-      iotdobs  = 0
-      igrdobs  = 0
-      iosduobs = 0
-      iosdvobs = 0
-      ilanuobs = 0
-      ilanvobs = 0
-      inlauobs = 0
-      inlavobs = 0
-      iuvchku  = 0
-      iuvchkv = 0
-
-      ! -----------------------------------------------------------------------
-      ! Find time coordinate for profiles
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_tim_prof( icycle, &
-         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-         &              profdata%nprof,   profdata%nyea, profdata%nmon, &
-         &              profdata%nday,    profdata%nhou, profdata%nmin, &
-         &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
-         &              iotdobs, ld_dailyav = ld_dailyav        )
-	 
-      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
-      
-      ! -----------------------------------------------------------------------
-      ! Check for profiles failing the grid search
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), &
-         &              profdata%nqc,     igrdobs                         )
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), &
-         &              profdata%nqc,     igrdobs                         )
-
-      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
-
-      ! -----------------------------------------------------------------------
-      ! Reject all observations for profiles with nqc > 10
-      ! -----------------------------------------------------------------------
-
-      CALL obs_pro_rej( profdata )
+      ! Reject all observations for profiles with nqc > iqc_cutoff
+      ! -----------------------------------------------------------------------
+
+      CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff )
 
       ! -----------------------------------------------------------------------
@@ -1021,43 +406,26 @@
       ! -----------------------------------------------------------------------
 
-      ! Zonal Velocity Component
-
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   &
-         &                 profdata%npvsta(:,1),  profdata%npvend(:,1), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          & 
-         &                 profdata%var(1)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(1)%vdep,                        &
-         &                 glamu,                 gphiu,                &
-         &                 gdept_1d,              umask,                &
-         &                 profdata%nqc,          profdata%var(1)%nvqc, &
-         &                 iosduobs,              ilanuobs,             &
-         &                 inlauobs,              ld_nea                )
-
-      CALL obs_mpp_sum_integer( iosduobs, iosduobsmpp )
-      CALL obs_mpp_sum_integer( ilanuobs, ilanuobsmpp )
-      CALL obs_mpp_sum_integer( inlauobs, inlauobsmpp )
-
-      ! Meridional Velocity Component
-
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   &
-         &                 profdata%npvsta(:,2),  profdata%npvend(:,2), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          & 
-         &                 profdata%var(2)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(2)%vdep,                        &
-         &                 glamv,                 gphiv,                &
-         &                 gdept_1d,              vmask,                &
-         &                 profdata%nqc,          profdata%var(2)%nvqc, &
-         &                 iosdvobs,              ilanvobs,             &
-         &                 inlavobs,              ld_nea                )
-
-      CALL obs_mpp_sum_integer( iosdvobs, iosdvobsmpp )
-      CALL obs_mpp_sum_integer( ilanvobs, ilanvobsmpp )
-      CALL obs_mpp_sum_integer( inlavobs, inlavobsmpp )
+      DO jvar = 1, profdata%nvar
+         CALL obs_coo_spc_3d( profdata%nprof,          profdata%nvprot(jvar),   &
+            &                 profdata%npvsta(:,jvar), profdata%npvend(:,jvar), &
+            &                 jpi,                     jpj,                     &
+            &                 jpk,                                              &
+            &                 profdata%mi,             profdata%mj,             &
+            &                 profdata%var(jvar)%mvk,                           &
+            &                 profdata%rlam,           profdata%rphi,           &
+            &                 profdata%var(jvar)%vdep,                          &
+            &                 pglam(:,:,jvar),         pgphi(:,:,jvar),         &
+            &                 gdept_1d,                zmask(:,:,:,jvar),       &
+            &                 profdata%nqc,            profdata%var(jvar)%nvqc, &
+            &                 iosdvobs(jvar),          ilanvobs(jvar),          &
+            &                 inlavobs(jvar),          ld_nea,                  &
+            &                 ibdyvobs(jvar),          ld_bound_reject,         &
+            &                 iqc_cutoff       )
+
+         CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) )
+      END DO
 
       ! -----------------------------------------------------------------------
@@ -1065,7 +433,9 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_uv_rej( profdata, iuvchku, iuvchkv )
-      CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
-      CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
+      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
+         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff )
+         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
+         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
+      ENDIF
 
       ! -----------------------------------------------------------------------
@@ -1081,9 +451,9 @@
       END DO
 
-      ! We want all data which has qc flags = 0
-
-      llvalid%luse(:) = ( profdata%nqc(:)  <= 10 )
+      ! We want all data which has qc flags <= iqc_cutoff
+
+      llvalid%luse(:) = ( profdata%nqc(:)  <= iqc_cutoff )
       DO jvar = 1,profdata%nvar
-         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
+         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff )
       END DO
 
@@ -1106,48 +476,46 @@
       
       IF(lwp) THEN
+      
          WRITE(numout,*)
-         WRITE(numout,*) 'obs_pre_vel :'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*)
-         WRITE(numout,*) ' Profiles outside time domain                = ', &
+         WRITE(numout,*) ' Profiles outside time domain                       = ', &
             &            iotdobsmpp
-         WRITE(numout,*) ' Remaining profiles that failed grid search  = ', &
+         WRITE(numout,*) ' Remaining profiles that failed grid search         = ', &
             &            igrdobsmpp
-         WRITE(numout,*) ' Remaining U data outside space domain       = ', &
-            &            iosduobsmpp
-         WRITE(numout,*) ' Remaining U data at land points             = ', &
-            &            ilanuobsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining U data near land points (removed) = ',&
-               &            inlauobsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining U data near land points (kept)    = ',&
-               &            inlauobsmpp
-         ENDIF
-         WRITE(numout,*) ' U observation rejected since V rejected     = ', &
-            &            iuvchku     
-         WRITE(numout,*) ' U data accepted                             = ', &
-            &            prodatqc%nvprotmpp(1)
-         WRITE(numout,*) ' Remaining V data outside space domain       = ', &
-            &            iosdvobsmpp
-         WRITE(numout,*) ' Remaining V data at land points             = ', &
-            &            ilanvobsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining V data near land points (removed) = ',&
-               &            inlavobsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining V data near land points (kept)    = ',&
-               &            inlavobsmpp
-         ENDIF
-         WRITE(numout,*) ' V observation rejected since U rejected     = ', &
-            &            iuvchkv     
-         WRITE(numout,*) ' V data accepted                             = ', &
-            &            prodatqc%nvprotmpp(2)
+         DO jvar = 1, profdata%nvar
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain       = ', &
+               &            iosdvobsmpp(jvar)
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points             = ', &
+               &            ilanvobsmpp(jvar)
+            IF (ld_nea) THEN
+               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ',&
+                  &            inlavobsmpp(jvar)
+            ELSE
+               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept)    = ',&
+                  &            inlavobsmpp(jvar)
+            ENDIF
+            IF ( TRIM(profdata%cvars(jvar)) == 'UVEL' ) THEN
+               WRITE(numout,*) ' U observation rejected since V rejected     = ', &
+                  &            iuvchku
+            ELSE IF ( TRIM(profdata%cvars(jvar)) == 'VVEL' ) THEN
+               WRITE(numout,*) ' V observation rejected since U rejected     = ', &
+                  &            iuvchkv
+            ENDIF
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ',&
+                  &            ibdyvobsmpp(jvar)
+            WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted                             = ', &
+               &            prodatqc%nvprotmpp(jvar)
+         END DO
 
          WRITE(numout,*)
          WRITE(numout,*) ' Number of observations per time step :'
          WRITE(numout,*)
-         WRITE(numout,997)
-         WRITE(numout,998)
+         WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles'
+         WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------'
+         DO jvar = 1, prodatqc%nvar
+            WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar))
+            WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------'
+         END DO
+         WRITE(numout,*) cout1
+         WRITE(numout,*) cout2
       ENDIF
       
@@ -1176,15 +544,13 @@
          DO jstp = nit000 - 1, nitend
             inrc = jstp - nit000 + 2
-            WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
-               &                    prodatqc%nvstpmpp(inrc,1), &
-               &                    prodatqc%nvstpmpp(inrc,2)
+            WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc)
+            DO jvar = 1, prodatqc%nvar
+               WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar)
+            END DO
+            WRITE(numout,*) cout1
          END DO
       ENDIF
 
-997   FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')
-998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
-999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
-
-   END SUBROUTINE obs_pre_vel
+   END SUBROUTINE obs_pre_prof
 
    SUBROUTINE obs_coo_tim( kcycle, &
@@ -1293,5 +659,5 @@
             &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
             kobsstp(jobs) = -1
-            kobsqc(jobs)  = kobsqc(jobs) + 11
+            kobsqc(jobs)  = IBSET(kobsqc(jobs),13)
             kotdobs       = kotdobs + 1
             CYCLE
@@ -1344,5 +710,5 @@
          IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
             & .OR.( kobsstp(jobs) > nitend ) ) THEN
-            kobsqc(jobs) = kobsqc(jobs) + 12
+            kobsqc(jobs) = IBSET(kobsqc(jobs),13)
             kotdobs = kotdobs + 1
             CYCLE
@@ -1389,5 +755,5 @@
       &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
       &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, &
-      &                    ld_dailyav )
+      &                    kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_tim ***
@@ -1433,7 +799,9 @@
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes    ! Types for daily averages
-      LOGICAL, OPTIONAL :: ld_dailyav    ! All types are daily averages
+      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
+
       !! * Local declarations
       INTEGER :: jobs
+      INTEGER :: iqc_cutoff=255
 
       !-----------------------------------------------------------------------
@@ -1454,9 +822,9 @@
          DO jobs = 1, kobsno
             
-            IF ( kobsqc(jobs) <= 10 ) THEN
+            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
                
                IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
                   & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
-                  kobsqc(jobs) = kobsqc(jobs) + 14
+                  kobsqc(jobs) = IBSET(kobsqc(jobs),13)
                   kotdobs      = kotdobs + 1
                   CYCLE
@@ -1467,24 +835,4 @@
       ENDIF
 
-      !------------------------------------------------------------------------
-      ! If ld_dailyav is set then all data assumed to be daily averaged
-      !------------------------------------------------------------------------
-      
-      IF ( PRESENT( ld_dailyav) ) THEN
-         IF (ld_dailyav) THEN
-            DO jobs = 1, kobsno
-               
-               IF ( kobsqc(jobs) <= 10 ) THEN
-                  
-                  IF ( kobsstp(jobs) == (nit000 - 1) ) THEN
-                     kobsqc(jobs) = kobsqc(jobs) + 14
-                     kotdobs      = kotdobs + 1
-                     CYCLE
-                  ENDIF
-                  
-               ENDIF
-            END DO
-         ENDIF
-      ENDIF
 
    END SUBROUTINE obs_coo_tim_prof
@@ -1521,5 +869,5 @@
       DO jobs = 1, kobsno
          IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
-            kobsqc(jobs) = kobsqc(jobs) + 18
+            kobsqc(jobs) = IBSET(kobsqc(jobs),12)
             kgrdobs = kgrdobs + 1
          ENDIF
@@ -1532,5 +880,7 @@
       &                       plam,   pphi,    pmask,            &
       &                       kobsqc, kosdobs, klanobs,          &
-      &                       knlaobs,ld_nea                     )
+      &                       knlaobs,ld_nea,                    &
+      &                       kbdyobs,ld_bound_reject,           &
+      &                       kqc_cutoff                         )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_spc_2d  ***
@@ -1565,11 +915,20 @@
       INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
          & kobsqc             ! Observation quality control
-      INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain
-      INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell
-      INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land
-      LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land
+      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain
+      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell
+      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land
+      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary
+      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land
+      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary 
+      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value
+
       !! * Local declarations
       REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
          & zgmsk              ! Grid mask
+#if defined key_bdy 
+      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
+         & zbmsk              ! Boundary mask
+      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
+#endif 
       REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
          & zglam, &           ! Model longitude at grid points
@@ -1588,5 +947,5 @@
          ! For invalid points use 2,2
 
-         IF ( kobsqc(jobs) >= 10 ) THEN
+         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
 
             igrdi(1,1,jobs) = 1
@@ -1613,13 +972,27 @@
 
       END DO
-      
-      CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pmask, zgmsk )
-      CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, plam, zglam )
-      CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pphi, zgphi )
+
+#if defined key_bdy             
+      ! Create a mask grid points in boundary rim
+      IF (ld_bound_reject) THEN
+         zbdymask(:,:) = 1.0_wp
+         DO ji = 1, nb_bdy
+            DO jj = 1, idx_bdy(ji)%nblen(1)
+               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
+            ENDDO
+         ENDDO
+ 
+         CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
+      ENDIF
+#endif       
+      
+      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk )
+      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam )
+      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
 
       DO jobs = 1, kobsno
 
          ! Skip bad observations
-         IF ( kobsqc(jobs) >= 10 ) CYCLE
+         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
 
          ! Flag if the observation falls outside the model spatial domain
@@ -1628,5 +1001,5 @@
             &  .OR. ( pobsphi(jobs) <  -90. ) &
             &  .OR. ( pobsphi(jobs) >   90. ) ) THEN
-            kobsqc(jobs) = kobsqc(jobs) + 11
+            kobsqc(jobs) = IBSET(kobsqc(jobs),11)
             kosdobs = kosdobs + 1
             CYCLE
@@ -1635,5 +1008,5 @@
          ! Flag if the observation falls with a model land cell
          IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
-            kobsqc(jobs) = kobsqc(jobs)  + 12
+            kobsqc(jobs) = IBSET(kobsqc(jobs),10)
             klanobs = klanobs + 1
             CYCLE
@@ -1649,6 +1022,6 @@
                IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
                   & .AND. &
-                  & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
-                  & ) THEN
+                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) )  &
+                  & < 1.0e-6_wp ) ) THEN
                   lgridobs = .TRUE.
                   iig = ji
@@ -1657,22 +1030,41 @@
             END DO
          END DO
-  
-         ! For observations on the grid reject them if their are at
-         ! a masked point
-         
+ 
          IF (lgridobs) THEN
             IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
-               kobsqc(jobs) = kobsqc(jobs) + 12
+               kobsqc(jobs) = IBSET(kobsqc(jobs),10)
                klanobs = klanobs + 1
                CYCLE
             ENDIF
          ENDIF
-                      
+
+ 
          ! Flag if the observation falls is close to land
          IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
-            IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14
             knlaobs = knlaobs + 1
-            CYCLE
+            IF (ld_nea) THEN
+               kobsqc(jobs) = IBSET(kobsqc(jobs),9)
+               CYCLE
+            ENDIF
          ENDIF
+
+#if defined key_bdy
+         ! Flag if the observation falls close to the boundary rim
+         IF (ld_bound_reject) THEN
+            IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
+               kobsqc(jobs) = IBSET(kobsqc(jobs),8)
+               kbdyobs = kbdyobs + 1
+               CYCLE
+            ENDIF
+            ! for observations on the grid...
+            IF (lgridobs) THEN
+               IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
+                  kobsqc(jobs) = IBSET(kobsqc(jobs),8)
+                  kbdyobs = kbdyobs + 1
+                  CYCLE
+               ENDIF
+            ENDIF
+         ENDIF
+#endif 
             
       END DO
@@ -1686,5 +1078,7 @@
       &                       plam,    pphi,    pdep,    pmask, &
       &                       kpobsqc, kobsqc,  kosdobs,        &
-      &                       klanobs, knlaobs, ld_nea          )
+      &                       klanobs, knlaobs, ld_nea,         &
+      &                       kbdyobs, ld_bound_reject,         &
+      &                       kqc_cutoff                        )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_spc_3d  ***
@@ -1709,5 +1103,13 @@
       !! * Modules used
       USE dom_oce, ONLY : &       ! Geographical information
-         & gdepw_1d                        
+         & gdepw_1d,      &
+         & gdepw_0,       &                       
+#if defined key_vvl
+         & gdepw_n,       & 
+         & gdept_n,       &
+#endif
+         & ln_zco,        &
+         & ln_zps,        &
+         & lk_vvl                        
 
       !! * Arguments
@@ -1743,8 +1145,19 @@
       INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell
       INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land
+      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary
       LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land
+      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary
+      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value
+
       !! * Local declarations
       REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
          & zgmsk              ! Grid mask
+#if defined key_bdy 
+      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
+         & zbmsk              ! Boundary mask
+      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
+#endif 
+      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
+         & zgdepw         
       REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
          & zglam, &           ! Model longitude at grid points
@@ -1754,4 +1167,5 @@
          & igrdj
       LOGICAL :: lgridobs           ! Is observation on a model grid point.
+      LOGICAL :: ll_next_to_land    ! Is a profile next to land 
       INTEGER :: iig, ijg           ! i,j of observation on model grid point.
       INTEGER :: jobs, jobsp, jk, ji, jj
@@ -1763,5 +1177,5 @@
          ! For invalid points use 2,2
 
-         IF ( kpobsqc(jobs) >= 10 ) THEN
+         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
             
             igrdi(1,1,jobs) = 1
@@ -1788,13 +1202,30 @@
          
       END DO
-      
-      CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, pmask, zgmsk )
-      CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, plam, zglam )
-      CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, pphi, zgphi )
+
+#if defined key_bdy 
+      ! Create a mask grid points in boundary rim
+      IF (ld_bound_reject) THEN           
+         zbdymask(:,:) = 1.0_wp
+         DO ji = 1, nb_bdy
+            DO jj = 1, idx_bdy(ji)%nblen(1)
+               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
+            ENDDO
+         ENDDO
+      ENDIF
+ 
+      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
+#endif 
+      
+      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
+      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
+      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
+      ! Need to know the bathy depth for each observation for sco
+      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), &
+         &                  zgdepw )
 
       DO jobs = 1, kprofno
 
          ! Skip bad profiles
-         IF ( kpobsqc(jobs) >= 10 ) CYCLE
+         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
 
          ! Check if this observation is on a grid point
@@ -1807,5 +1238,5 @@
                IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
                   & .AND. &
-                  & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
+                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
                   & ) THEN
                   lgridobs = .TRUE.
@@ -1816,4 +1247,11 @@
          END DO
 
+         ! Check if next to land 
+         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN 
+            ll_next_to_land=.TRUE. 
+         ELSE 
+            ll_next_to_land=.FALSE. 
+         ENDIF 
+         
          ! Reject observations
 
@@ -1827,15 +1265,53 @@
                &  .OR. ( pobsdep(jobsp) < 0.0          )       &
                &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
-               kobsqc(jobsp) = kobsqc(jobsp) + 11
+               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
                kosdobs = kosdobs + 1
                CYCLE
             ENDIF
 
-            ! Flag if the observation falls with a model land cell
-            IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
-               &  == 0.0_wp ) THEN
-               kobsqc(jobsp) = kobsqc(jobsp) + 12
-               klanobs = klanobs + 1
-               CYCLE
+            ! To check if an observations falls within land there are two cases: 
+            ! 1: z-coordibnates, where the check uses the mask 
+            ! 2: terrain following (eg s-coordinates),  
+            !    where we use the depth of the bottom cell to mask observations 
+             
+            IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco )  ) THEN !(CASE 1) 
+                
+               ! Flag if the observation falls with a model land cell 
+               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
+                  &  == 0.0_wp ) THEN 
+                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
+                  klanobs = klanobs + 1 
+                  CYCLE 
+               ENDIF 
+             
+               ! Flag if the observation is close to land 
+               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
+                  &  0.0_wp) THEN 
+                  knlaobs = knlaobs + 1 
+                  IF (ld_nea) THEN   
+                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
+                  ENDIF  
+               ENDIF 
+             
+            ELSE ! Case 2 
+               ! Flag if the observation is deeper than the bathymetry 
+               ! Or if it is within the mask 
+               IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
+                  &     .OR. & 
+                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
+                  &  == 0.0_wp) ) THEN
+                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
+                  klanobs = klanobs + 1 
+                  CYCLE 
+               ENDIF 
+                
+               ! Flag if the observation is close to land 
+               IF ( ll_next_to_land ) THEN 
+                  knlaobs = knlaobs + 1 
+                  IF (ld_nea) THEN   
+                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
+                  ENDIF  
+               ENDIF 
+             
             ENDIF
 
@@ -1845,5 +1321,5 @@
             IF (lgridobs) THEN
                IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
-                  kobsqc(jobsp) = kobsqc(jobsp) + 12
+                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
                   klanobs = klanobs + 1
                   CYCLE
@@ -1851,11 +1327,4 @@
             ENDIF
             
-            ! Flag if the observation falls is close to land
-            IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
-               &  0.0_wp) THEN
-               IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
-               knlaobs = knlaobs + 1
-            ENDIF
-
             ! Set observation depth equal to that of the first model depth
             IF ( pobsdep(jobsp) <= pdep(1) ) THEN
@@ -1863,4 +1332,23 @@
             ENDIF
             
+#if defined key_bdy
+            ! Flag if the observation falls close to the boundary rim
+            IF (ld_bound_reject) THEN
+               IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
+                  kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
+                  kbdyobs = kbdyobs + 1
+                  CYCLE
+               ENDIF
+               ! for observations on the grid...
+               IF (lgridobs) THEN
+                  IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
+                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
+                     kbdyobs = kbdyobs + 1
+                     CYCLE
+                  ENDIF
+               ENDIF
+            ENDIF
+#endif 
+            
          END DO
       END DO
@@ -1868,5 +1356,5 @@
    END SUBROUTINE obs_coo_spc_3d
 
-   SUBROUTINE obs_pro_rej( profdata )
+   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_pro_rej ***
@@ -1886,4 +1374,6 @@
       !! * Arguments
       TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data
+      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value
+
       !! * Local declarations
       INTEGER :: jprof
@@ -1895,5 +1385,5 @@
       DO jprof = 1, profdata%nprof
 
-         IF ( profdata%nqc(jprof) > 10 ) THEN
+         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
             
             DO jvar = 1, profdata%nvar
@@ -1903,5 +1393,5 @@
                   
                   profdata%var(jvar)%nvqc(jobs) = &
-                     & profdata%var(jvar)%nvqc(jobs) + 26
+                     & IBSET(profdata%var(jvar)%nvqc(jobs),14)
 
                END DO
@@ -1915,5 +1405,5 @@
    END SUBROUTINE obs_pro_rej
 
-   SUBROUTINE obs_uv_rej( profdata, knumu, knumv )
+   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_uv_rej ***
@@ -1935,4 +1425,6 @@
       INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected
       INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected
+      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
+
       !! * Local declarations
       INTEGER :: jprof
@@ -1954,12 +1446,12 @@
          DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
             
-            IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. &
-               & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN
-               profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42
+            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. &
+               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN
+               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
                knumv = knumv + 1
             ENDIF
-            IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. &
-               & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN
-               profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42
+            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. &
+               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN
+               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
                knumu = knumu + 1
             ENDIF
