Index: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 10119)
+++ branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 10120)
@@ -7,37 +7,34 @@
 
    !!----------------------------------------------------------------------
-   !!   obs_pro_opt :    Compute the model counterpart of temperature and
-   !!                    salinity observations from profiles
-   !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly
-   !!                    observations
-   !!   obs_sst_opt :    Compute the model counterpart of sea surface temperature
-   !!                    observations
-   !!   obs_sss_opt :    Compute the model counterpart of sea surface salinity
-   !!                    observations
-   !!   obs_seaice_opt : Compute the model counterpart of sea ice concentration
-   !!                    observations
-   !!
-   !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional
-   !!                    components of velocity from observations.
+   !!   obs_prof_opt :    Compute the model counterpart of profile data
+   !!   obs_surf_opt :    Compute the model counterpart of surface data
    !!----------------------------------------------------------------------
 
-   !! * Modules used   
+   !! * Modules used
    USE par_kind, ONLY : &         ! Precision variables
       & wp
    USE in_out_manager             ! I/O manager
    USE obs_inter_sup              ! Interpolation support
-   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt
+   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the obs pt
       & obs_int_h2d, &
       & obs_int_h2d_init
-   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt
+   USE obs_averg_h2d, ONLY : &    ! Horizontal averaging to the obs footprint
+      & obs_avg_h2d, &
+      & obs_avg_h2d_init, &
+      & obs_max_fpsize
+   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt
       & obs_int_z1d,    &
       & obs_int_z1d_spl
-   USE obs_const,  ONLY :     &
-      & obfillflt		  ! Fillvalue   
+   USE obs_const,  ONLY :    &    ! Obs fill value
+      & obfillflt
    USE dom_oce,       ONLY : &
-      & glamt, glamu, glamv, &
-      & gphit, gphiu, gphiv
-   USE lib_mpp,       ONLY : &
+      & glamt, glamf, &
+      & gphit, gphif
+   USE lib_mpp,       ONLY : &    ! Warning and stopping routines
       & ctl_warn, ctl_stop
+   USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time
+      & sbc_dcy, nday_qsr
+   USE obs_grid,      ONLY : & 
+      & obs_level_search     
 
    IMPLICIT NONE
@@ -46,12 +43,9 @@
    PRIVATE
 
-   PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations
-      &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations
-      &   obs_sst_opt, &  ! Compute the model counterpart of SST observations
-      &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations
-      &   obs_seaice_opt, &
-      &   obs_vel_opt     ! Compute the model counterpart of velocity profile data
-
-   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
+   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs
+      &   obs_surf_opt     ! Compute the model counterpart of surface obs
+
+   INTEGER, PARAMETER, PUBLIC :: &
+      & imaxavtypes = 20   ! Max number of daily avgd obs types
 
    !!----------------------------------------------------------------------
@@ -61,9 +55,16 @@
    !!----------------------------------------------------------------------
 
+   !! * Substitutions 
+#  include "domzgr_substitute.h90" 
 CONTAINS
 
-   SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
-      &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, &
-      &                    kdailyavtypes )
+
+   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          &
+      &                     kit000, kdaystp,                      &
+      &                     pvar1, pvar2, pgdept, pgdepw,         &
+      &                     pmask1, pmask2,                       &  
+      &                     plam1, plam2, pphi1, pphi2,           &
+      &                     k1dint, k2dint, kdailyavtypes )
+
       !!-----------------------------------------------------------------------
       !!
@@ -78,5 +79,5 @@
       !!
       !!    First, a vertical profile of horizontally interpolated model
-      !!    now temperatures is computed at the obs (lon, lat) point.
+      !!    now values is computed at the obs (lon, lat) point.
       !!    Several horizontal interpolation schemes are available:
       !!        - distance-weighted (great circle) (k2dint = 0)
@@ -86,5 +87,5 @@
       !!        - polynomial (quadrilateral grid)  (k2dint = 4)
       !!
-      !!    Next, the vertical temperature profile is interpolated to the
+      !!    Next, the vertical profile is interpolated to the
       !!    data depth points. Two vertical interpolation schemes are
       !!    available:
@@ -96,14 +97,11 @@
       !!    routine.
       !!
-      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is
+      !!    If the logical is switched on, the model equivalent is
       !!    a daily mean model temperature field. So, we first compute
       !!    the mean, then interpolate only at the end of the day.
       !!
-      !!    Note: the in situ temperature observations must be converted
+      !!    Note: in situ temperature observations must be converted
       !!    to potential temperature (the model variable) prior to
       !!    assimilation. 
-      !!??????????????????????????????????????????????????????????????
-      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
-      !!??????????????????????????????????????????????????????????????
       !!
       !! ** Action  :
@@ -115,6 +113,8 @@
       !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
       !!      ! 07-03 (K. Mogensen) General handling of profiles
+      !!      ! 15-02 (M. Martin) Combined routine for all profile types
+      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes
       !!-----------------------------------------------------------------------
-  
+
       !! * Modules used
       USE obs_profiles_def ! Definition of storage space for profile obs.
@@ -123,22 +123,31 @@
 
       !! * Arguments
-      TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening
-      INTEGER, INTENT(IN) :: kt        ! Time step
-      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
+      TYPE(obs_prof), INTENT(INOUT) :: &
+         & prodatqc                  ! Subset of profile data passing QC
+      INTEGER, INTENT(IN) :: kt      ! Time step
+      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
       INTEGER, INTENT(IN) :: kpj
       INTEGER, INTENT(IN) :: kpk
-      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
-                                       !   (kit000-1 = restart time)
-      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
-      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
-      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                    
+      INTEGER, INTENT(IN) :: kit000  ! Number of the first time step
+                                     !   (kit000-1 = restart time)
+      INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header)
+      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header)
+      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
-         & ptn,    &    ! Model temperature field
-         & psn,    &    ! Model salinity field
-         & ptmask       ! Land-sea mask
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
-         & pgdept       ! Model array of depth levels
+         & pvar1,    &               ! Model field 1
+         & pvar2,    &               ! Model field 2
+         & pmask1,   &               ! Land-sea mask 1
+         & pmask2                    ! Land-sea mask 2
+      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
+         & plam1,    &               ! Model longitudes for variable 1
+         & plam2,    &               ! Model longitudes for variable 2
+         & pphi1,    &               ! Model latitudes for variable 1
+         & pphi2                     ! Model latitudes for variable 2
+      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
+         & pgdept, &                 ! Model array of depth T levels 
+         & pgdepw                    ! Model array of depth W levels 
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
-         & kdailyavtypes! Types for daily averages
+         & kdailyavtypes             ! Types for daily averages
+
       !! * Local declarations
       INTEGER ::   ji
@@ -152,130 +161,187 @@
       INTEGER ::   iend
       INTEGER ::   iobs
+      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes 
+      INTEGER ::   inum_obs
       INTEGER, DIMENSION(imaxavtypes) :: &
          & idailyavtypes
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
+         & igrdi1, &
+         & igrdi2, &
+         & igrdj1, &
+         & igrdj2
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
+
       REAL(KIND=wp) :: zlam
       REAL(KIND=wp) :: zphi
       REAL(KIND=wp) :: zdaystp
       REAL(KIND=wp), DIMENSION(kpk) :: &
-         & zobsmask, &
+         & zobsmask1, &
+         & zobsmask2, &
          & zobsk,    &
          & zobs2k
-      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
+      REAL(KIND=wp), DIMENSION(2,2,1) :: &
+         & zweig1, &
+         & zweig2, &
          & zweig
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
-         & zmask, &
-         & zintt, &
-         & zints, &
-         & zinmt, &
-         & zinms
+         & zmask1, &
+         & zmask2, &
+         & zint1,  &
+         & zint2,  &
+         & zinm1,  &
+         & zinm2,  &
+         & zgdept, & 
+         & zgdepw
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zglam, &
-         & zgphi
-      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi, &
-         & igrdj
+         & zglam1, &
+         & zglam2, &
+         & zgphi1, &
+         & zgphi2
+      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2   
+      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
+
+      LOGICAL :: ld_dailyav
 
       !------------------------------------------------------------------------
       ! Local initialization 
       !------------------------------------------------------------------------
-      ! ... Record and data counters
+      ! Record and data counters
       inrc = kt - kit000 + 2
       ipro = prodatqc%npstp(inrc)
- 
+
       ! Daily average types
+      ld_dailyav = .FALSE.
       IF ( PRESENT(kdailyavtypes) ) THEN
          idailyavtypes(:) = kdailyavtypes(:)
+         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
       ELSE
          idailyavtypes(:) = -1
       ENDIF
 
-      ! Initialize daily mean for first timestep
+      ! Daily means are calculated for values over timesteps:
+      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
       idayend = MOD( kt - kit000 + 1, kdaystp )
 
-      ! Added kt == 0 test to catch restart case 
-      IF ( idayend == 1 .OR. kt == 0) THEN
-         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
+      IF ( ld_dailyav ) THEN
+
+         ! Initialize daily mean for first timestep of the day
+         IF ( idayend == 1 .OR. kt == 0 ) THEN
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     prodatqc%vdmean(ji,jj,jk,1) = 0.0
+                     prodatqc%vdmean(ji,jj,jk,2) = 0.0
+                  END DO
+               END DO
+            END DO
+         ENDIF
+
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  prodatqc%vdmean(ji,jj,jk,1) = 0.0
-                  prodatqc%vdmean(ji,jj,jk,2) = 0.0
+                  ! Increment field 1 for computing daily mean
+                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
+                     &                        + pvar1(ji,jj,jk)
+                  ! Increment field 2 for computing daily mean
+                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
+                     &                        + pvar2(ji,jj,jk)
                END DO
             END DO
          END DO
-      ENDIF
-
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! Increment the temperature field for computing daily mean
-               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                  &                        + ptn(ji,jj,jk)
-               ! Increment the salinity field for computing daily mean
-               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                  &                        + psn(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-   
-      ! Compute the daily mean at the end of day
-      zdaystp = 1.0 / REAL( kdaystp )
-      IF ( idayend == 0 ) THEN
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                     &                        * zdaystp
-                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                  &                           * zdaystp
+
+         ! Compute the daily mean at the end of day
+         zdaystp = 1.0 / REAL( kdaystp )
+         IF ( idayend == 0 ) THEN
+            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
+            CALL FLUSH(numout)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
+                        &                        * zdaystp
+                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
+                        &                        * zdaystp
+                  END DO
                END DO
             END DO
-         END DO
+         ENDIF
+
       ENDIF
 
       ! Get the data for interpolation
       ALLOCATE( &
-         & igrdi(2,2,ipro),      &
-         & igrdj(2,2,ipro),      &
-         & zglam(2,2,ipro),      &
-         & zgphi(2,2,ipro),      &
-         & zmask(2,2,kpk,ipro),  &
-         & zintt(2,2,kpk,ipro),  &
-         & zints(2,2,kpk,ipro)   &
+         & igrdi1(2,2,ipro),      &
+         & igrdi2(2,2,ipro),      &
+         & igrdj1(2,2,ipro),      &
+         & igrdj2(2,2,ipro),      &
+         & zglam1(2,2,ipro),      &
+         & zglam2(2,2,ipro),      &
+         & zgphi1(2,2,ipro),      &
+         & zgphi2(2,2,ipro),      &
+         & zmask1(2,2,kpk,ipro),  &
+         & zmask2(2,2,kpk,ipro),  &
+         & zint1(2,2,kpk,ipro),   &
+         & zint2(2,2,kpk,ipro),   &
+         & zgdept(2,2,kpk,ipro),  & 
+         & zgdepw(2,2,kpk,ipro)   & 
          & )
 
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
          iobs = jobs - prodatqc%nprofup
-         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
-         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
-         igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
-         igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
-         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
-         igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
+         igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1
+         igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1
+         igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1
+         igrdj1(1,2,iobs) = prodatqc%mj(jobs,1)
+         igrdi1(2,1,iobs) = prodatqc%mi(jobs,1)
+         igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1
+         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1)
+         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1)
+         igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1
+         igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1
+         igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1
+         igrdj2(1,2,iobs) = prodatqc%mj(jobs,2)
+         igrdi2(2,1,iobs) = prodatqc%mi(jobs,2)
+         igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1
+         igrdi2(2,2,iobs) = prodatqc%mi(jobs,2)
+         igrdj2(2,2,iobs) = prodatqc%mj(jobs,2)
       END DO
 
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints )
+      ! Initialise depth arrays
+      zgdept(:,:,:,:) = 0.0
+      zgdepw(:,:,:,:) = 0.0
+
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 )
+      
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 )
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 )
+
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 
 
       ! At the end of the day also get interpolated means
-      IF ( idayend == 0 ) THEN
+      IF ( ld_dailyav .AND. idayend == 0 ) THEN
 
          ALLOCATE( &
-            & zinmt(2,2,kpk,ipro),  &
-            & zinms(2,2,kpk,ipro)   &
+            & zinm1(2,2,kpk,ipro),  &
+            & zinm2(2,2,kpk,ipro)   &
             & )
 
-         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
-            &                  prodatqc%vdmean(:,:,:,1), zinmt )
-         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
-            &                  prodatqc%vdmean(:,:,:,2), zinms )
+         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, &
+            &                  prodatqc%vdmean(:,:,:,1), zinm1 )
+         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
+            &                  prodatqc%vdmean(:,:,:,2), zinm2 )
 
       ENDIF
 
+      ! Return if no observations to process 
+      ! Has to be done after comm commands to ensure processors 
+      ! stay in sync 
+      IF ( ipro == 0 ) RETURN 
+
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
 
@@ -283,5 +349,5 @@
 
          IF ( kt /= prodatqc%mstp(jobs) ) THEN
-            
+
             IF(lwp) THEN
                WRITE(numout,*)
@@ -298,162 +364,321 @@
             CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
          ENDIF
-         
+
          zlam = prodatqc%rlam(jobs)
          zphi = prodatqc%rphi(jobs)
+
+         ! Horizontal weights 
+         ! Masked values are calculated later.  
+         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
+
+            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
+               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), &
+               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 )
+
+         ENDIF
+
+         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
+
+            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
+               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), &
+               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 )
+ 
+         ENDIF
+
+         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
+
+            zobsk(:) = obfillflt
+
+            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
+
+               IF ( idayend == 0 )  THEN
+                  ! Daily averaged data
+
+                  ! vertically interpolate all 4 corners 
+                  ista = prodatqc%npvsta(jobs,1) 
+                  iend = prodatqc%npvend(jobs,1) 
+                  inum_obs = iend - ista + 1 
+                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
+
+                  DO iin=1,2 
+                     DO ijn=1,2 
+
+                        IF ( k1dint == 1 ) THEN 
+                           CALL obs_int_z1d_spl( kpk, & 
+                              &     zinm1(iin,ijn,:,iobs), & 
+                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
+                              &     zmask1(iin,ijn,:,iobs)) 
+                        ENDIF 
+       
+                        CALL obs_level_search(kpk, & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                           &    iv_indic) 
+
+                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
+                           &    prodatqc%var(1)%vdep(ista:iend), & 
+                           &    zinm1(iin,ijn,:,iobs), & 
+                           &    zobs2k, interp_corner(iin,ijn,:), & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    zmask1(iin,ijn,:,iobs)) 
+       
+                     ENDDO 
+                  ENDDO 
+
+               ENDIF !idayend
+
+            ELSE   
+
+               ! Point data 
+     
+               ! vertically interpolate all 4 corners 
+               ista = prodatqc%npvsta(jobs,1) 
+               iend = prodatqc%npvend(jobs,1) 
+               inum_obs = iend - ista + 1 
+               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
+               DO iin=1,2  
+                  DO ijn=1,2 
+                    
+                     IF ( k1dint == 1 ) THEN 
+                        CALL obs_int_z1d_spl( kpk, & 
+                           &    zint1(iin,ijn,:,iobs),& 
+                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
+                           &    zmask1(iin,ijn,:,iobs)) 
+  
+                     ENDIF 
+       
+                     CALL obs_level_search(kpk, & 
+                         &        zgdept(iin,ijn,:,iobs),& 
+                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                         &        iv_indic) 
+
+                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
+                         &          prodatqc%var(1)%vdep(ista:iend),     & 
+                         &          zint1(iin,ijn,:,iobs),            & 
+                         &          zobs2k,interp_corner(iin,ijn,:), & 
+                         &          zgdept(iin,ijn,:,iobs),         & 
+                         &          zmask1(iin,ijn,:,iobs) )      
          
-         ! Horizontal weights and vertical mask
-
-         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
-            & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
-
-            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
-               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
-               &                   zmask(:,:,:,iobs), zweig, zobsmask )
-
-         ENDIF
-
-         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
+                  ENDDO 
+               ENDDO 
+             
+            ENDIF 
+
+            !------------------------------------------------------------- 
+            ! Compute the horizontal interpolation for every profile level 
+            !------------------------------------------------------------- 
+             
+            DO ikn=1,inum_obs 
+               iend=ista+ikn-1
+                  
+               zweig(:,:,1) = 0._wp 
+   
+               ! This code forces the horizontal weights to be  
+               ! zero IF the observation is below the bottom of the  
+               ! corners of the interpolation nodes, Or if it is in  
+               ! the mask. This is important for observations near  
+               ! steep bathymetry 
+               DO iin=1,2 
+                  DO ijn=1,2 
+     
+                     depth_loop1: DO ik=kpk,2,-1 
+                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
+                            
+                           zweig(iin,ijn,1) = &  
+                              & zweig1(iin,ijn,1) * & 
+                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
+                              &  - prodatqc%var(1)%vdep(iend)),0._wp) 
+                            
+                           EXIT depth_loop1 
+
+                        ENDIF 
+
+                     ENDDO depth_loop1 
+     
+                  ENDDO 
+               ENDDO 
+   
+               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
+                  &              prodatqc%var(1)%vmod(iend:iend) ) 
+
+                  ! Set QC flag for any observations found below the bottom
+                  ! needed as the check here is more strict than that in obs_prep
+               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
+ 
+            ENDDO 
+ 
+            DEALLOCATE(interp_corner,iv_indic) 
+          
+         ENDIF 
+
+         ! For the second variable
+         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
 
             zobsk(:) = obfillflt
 
-	    IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
+            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
 
                IF ( idayend == 0 )  THEN
+                  ! Daily averaged data
+
+                  ! vertically interpolate all 4 corners 
+                  ista = prodatqc%npvsta(jobs,2) 
+                  iend = prodatqc%npvend(jobs,2) 
+                  inum_obs = iend - ista + 1 
+                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
+
+                  DO iin=1,2 
+                     DO ijn=1,2 
+
+                        IF ( k1dint == 1 ) THEN 
+                           CALL obs_int_z1d_spl( kpk, & 
+                              &     zinm2(iin,ijn,:,iobs), & 
+                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
+                              &     zmask2(iin,ijn,:,iobs)) 
+                        ENDIF 
+       
+                        CALL obs_level_search(kpk, & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
+                           &    iv_indic) 
+
+                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
+                           &    prodatqc%var(2)%vdep(ista:iend), & 
+                           &    zinm2(iin,ijn,:,iobs), & 
+                           &    zobs2k, interp_corner(iin,ijn,:), & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    zmask2(iin,ijn,:,iobs)) 
+       
+                     ENDDO 
+                  ENDDO 
+
+               ENDIF !idayend
+
+            ELSE   
+
+               ! Point data 
+     
+               ! vertically interpolate all 4 corners 
+               ista = prodatqc%npvsta(jobs,2) 
+               iend = prodatqc%npvend(jobs,2) 
+               inum_obs = iend - ista + 1 
+               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
+               DO iin=1,2  
+                  DO ijn=1,2 
+                    
+                     IF ( k1dint == 1 ) THEN 
+                        CALL obs_int_z1d_spl( kpk, & 
+                           &    zint2(iin,ijn,:,iobs),& 
+                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
+                           &    zmask2(iin,ijn,:,iobs)) 
+  
+                     ENDIF 
+       
+                     CALL obs_level_search(kpk, & 
+                         &        zgdept(iin,ijn,:,iobs),& 
+                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
+                         &        iv_indic) 
+
+                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
+                         &          prodatqc%var(2)%vdep(ista:iend),     & 
+                         &          zint2(iin,ijn,:,iobs),            & 
+                         &          zobs2k,interp_corner(iin,ijn,:), & 
+                         &          zgdept(iin,ijn,:,iobs),         & 
+                         &          zmask2(iin,ijn,:,iobs) )      
+         
+                  ENDDO 
+               ENDDO 
+             
+            ENDIF 
+
+            !------------------------------------------------------------- 
+            ! Compute the horizontal interpolation for every profile level 
+            !------------------------------------------------------------- 
+             
+            DO ikn=1,inum_obs 
+               iend=ista+ikn-1
                   
-                  ! Daily averaged moored buoy (MRB) data
-                  
-                  CALL obs_int_h2d( kpk, kpk,      &
-                     &              zweig, zinmt(:,:,:,iobs), zobsk )
-                  
-                  
-               ELSE
-               
-                  CALL ctl_stop( ' A nonzero' //     &
-                     &           ' number of profile T BUOY data should' // &
-                     &           ' only occur at the end of a given day' )
-
-               ENDIF
-	       
-            ELSE 
-               
-               ! Point data
-
-               CALL obs_int_h2d( kpk, kpk,      &
-                  &              zweig, zintt(:,:,:,iobs), zobsk )
-
-            ENDIF
-
-            !-------------------------------------------------------------
-            ! Compute vertical second-derivative of the interpolating 
-            ! polynomial at obs points
-            !-------------------------------------------------------------
-            
-            IF ( k1dint == 1 ) THEN
-               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
-                  &                  pgdept, zobsmask )
-            ENDIF
-            
-            !-----------------------------------------------------------------
-            !  Vertical interpolation to the observation point
-            !-----------------------------------------------------------------
-            ista = prodatqc%npvsta(jobs,1)
-            iend = prodatqc%npvend(jobs,1)
-            CALL obs_int_z1d( kpk,                &
-               & prodatqc%var(1)%mvk(ista:iend),  &
-               & k1dint, iend - ista + 1,         &
-               & prodatqc%var(1)%vdep(ista:iend), &
-               & zobsk, zobs2k,                   &
-               & prodatqc%var(1)%vmod(ista:iend), &
-               & pgdept, zobsmask )
-
-         ENDIF
-
-         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
-
-            zobsk(:) = obfillflt
-
-            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
-
-               IF ( idayend == 0 )  THEN
-
-                  ! Daily averaged moored buoy (MRB) data
-                  
-                  CALL obs_int_h2d( kpk, kpk,      &
-                     &              zweig, zinms(:,:,:,iobs), zobsk )
-                  
-               ELSE
-
-                  CALL ctl_stop( ' A nonzero' //     &
-                     &           ' number of profile S BUOY data should' // &
-                     &           ' only occur at the end of a given day' )
-
-               ENDIF
-
-            ELSE
-               
-               ! Point data
-
-               CALL obs_int_h2d( kpk, kpk,      &
-                  &              zweig, zints(:,:,:,iobs), zobsk )
-
-            ENDIF
-
-
-            !-------------------------------------------------------------
-            ! Compute vertical second-derivative of the interpolating 
-            ! polynomial at obs points
-            !-------------------------------------------------------------
-            
-            IF ( k1dint == 1 ) THEN
-               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
-                  &                  pgdept, zobsmask )
-            ENDIF
-            
-            !----------------------------------------------------------------
-            !  Vertical interpolation to the observation point
-            !----------------------------------------------------------------
-            ista = prodatqc%npvsta(jobs,2)
-            iend = prodatqc%npvend(jobs,2)
-            CALL obs_int_z1d( kpk, &
-               & prodatqc%var(2)%mvk(ista:iend),&
-               & k1dint, iend - ista + 1, &
-               & prodatqc%var(2)%vdep(ista:iend),&
-               & zobsk, zobs2k, &
-               & prodatqc%var(2)%vmod(ista:iend),&
-               & pgdept, zobsmask )
-
-         ENDIF
-
-      END DO
+               zweig(:,:,1) = 0._wp 
+   
+               ! This code forces the horizontal weights to be  
+               ! zero IF the observation is below the bottom of the  
+               ! corners of the interpolation nodes, Or if it is in  
+               ! the mask. This is important for observations near  
+               ! steep bathymetry 
+               DO iin=1,2 
+                  DO ijn=1,2 
+     
+                     depth_loop2: DO ik=kpk,2,-1 
+                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
+                            
+                           zweig(iin,ijn,1) = &  
+                              & zweig2(iin,ijn,1) * & 
+                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
+                              &  - prodatqc%var(2)%vdep(iend)),0._wp) 
+                            
+                           EXIT depth_loop2 
+
+                        ENDIF 
+
+                     ENDDO depth_loop2 
+     
+                  ENDDO 
+               ENDDO 
+   
+               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
+                  &              prodatqc%var(2)%vmod(iend:iend) ) 
+
+                  ! Set QC flag for any observations found below the bottom
+                  ! needed as the check here is more strict than that in obs_prep
+               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
  
+            ENDDO 
+ 
+            DEALLOCATE(interp_corner,iv_indic) 
+          
+         ENDIF 
+
+      ENDDO
+
       ! Deallocate the data for interpolation
       DEALLOCATE( &
-         & igrdi, &
-         & igrdj, &
-         & zglam, &
-         & zgphi, &
-         & zmask, &
-         & zintt, &
-         & zints  &
+         & igrdi1, &
+         & igrdi2, &
+         & igrdj1, &
+         & igrdj2, &
+         & zglam1, &
+         & zglam2, &
+         & zgphi1, &
+         & zgphi2, &
+         & zmask1, &
+         & zmask2, &
+         & zint1,  &
+         & zint2,  &
+         & zgdept, &
+         & zgdepw  &
          & )
+
       ! At the end of the day also get interpolated means
-      IF ( idayend == 0 ) THEN
+      IF ( ld_dailyav .AND. idayend == 0 ) THEN
          DEALLOCATE( &
-            & zinmt,  &
-            & zinms   &
+            & zinm1,  &
+            & zinm2   &
             & )
       ENDIF
 
       prodatqc%nprofup = prodatqc%nprofup + ipro 
-      
-   END SUBROUTINE obs_pro_opt
-
-   SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
-      &                    psshn, psshmask, k2dint )
+
+   END SUBROUTINE obs_prof_opt
+
+   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
+      &                     kit000, kdaystp, psurf, psurfmask,   &
+      &                     k2dint, ldnightav, plamscl, pphiscl, &
+      &                     lindegrees )
+
       !!-----------------------------------------------------------------------
       !!
-      !!                     ***  ROUTINE obs_sla_opt  ***
-      !!
-      !! ** Purpose : Compute the model counterpart of sea level anomaly
+      !!                     ***  ROUTINE obs_surf_opt  ***
+      !!
+      !! ** Purpose : Compute the model counterpart of surface
       !!              data by interpolating from the model grid to the 
       !!              observation point.
@@ -462,5 +687,5 @@
       !!              the model values at the corners of the surrounding grid box.
       !!
-      !!    The now model SSH is first computed at the obs (lon, lat) point.
+      !!    The new model value is first computed at the obs (lon, lat) point.
       !!
       !!    Several horizontal interpolation schemes are available:
@@ -470,7 +695,9 @@
       !!        - bilinear (quadrilateral grid)    (k2dint = 3)
       !!        - polynomial (quadrilateral grid)  (k2dint = 4)
-      !!  
-      !!    The sea level anomaly at the observation points is then computed 
-      !!    by removing a mean dynamic topography (defined at the obs. point).
+      !!
+      !!    Two horizontal averaging schemes are also available:
+      !!        - weighted radial footprint        (k2dint = 5)
+      !!        - weighted rectangular footprint   (k2dint = 6)
+      !!
       !!
       !! ** Action  :
@@ -478,6 +705,8 @@
       !! History :
       !!      ! 07-03 (A. Weaver)
+      !!      ! 15-02 (M. Martin) Combined routine for surface types
+      !!      ! 17-03 (M. Martin) Added horizontal averaging options
       !!-----------------------------------------------------------------------
-  
+
       !! * Modules used
       USE obs_surf_def  ! Definition of storage space for surface observations
@@ -486,15 +715,23 @@
 
       !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening
-      INTEGER, INTENT(IN) :: kt      ! Time step
-      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
+      TYPE(obs_surf), INTENT(INOUT) :: &
+         & surfdataqc                  ! Subset of surface data passing QC
+      INTEGER, INTENT(IN) :: kt        ! Time step
+      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
       INTEGER, INTENT(IN) :: kpj
-      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step 
-                                      !   (kit000-1 = restart time)
-      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & psshn,  &    ! Model SSH field
-         & psshmask     ! Land-sea mask
-         
+      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
+                                       !   (kit000-1 = restart time)
+      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
+      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
+         & psurf,  &                   ! Model surface field
+         & psurfmask                   ! Land-sea mask
+      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions
+         & pphiscl                     ! This is the full width (rather than half-width)
+      LOGICAL, INTENT(IN) :: &
+         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
+
       !! * Local declarations
       INTEGER :: ji
@@ -502,67 +739,189 @@
       INTEGER :: jobs
       INTEGER :: inrc
-      INTEGER :: isla
+      INTEGER :: isurf
       INTEGER :: iobs
-      REAL(KIND=wp) :: zlam
-      REAL(KIND=wp) :: zphi
-      REAL(KIND=wp) :: zext(1), zobsmask(1)
-      REAL(kind=wp), DIMENSION(2,2,1) :: &
-         & zweig
+      INTEGER :: imaxifp, imaxjfp
+      INTEGER :: imodi, imodj
+      INTEGER :: idayend
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
+         & igrdi,   &
+         & igrdj,   &
+         & igrdip1, &
+         & igrdjp1
+      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
+         & icount_night,      &
+         & imask_night
+      REAL(wp) :: zlam
+      REAL(wp) :: zphi
+      REAL(wp), DIMENSION(1) :: zext, zobsmask
+      REAL(wp) :: zdaystp
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zmask, &
-         & zsshl, &
-         & zglam, &
-         & zgphi
-      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi, &
-         & igrdj
+         & zweig,  &
+         & zmask,  &
+         & zsurf,  &
+         & zsurfm, &
+         & zsurftmp, &
+         & zglam,  &
+         & zgphi,  &
+         & zglamf, &
+         & zgphif
+
+      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
+         & zintmp,  &
+         & zouttmp, &
+         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
 
       !------------------------------------------------------------------------
       ! Local initialization 
       !------------------------------------------------------------------------
-      ! ... Record and data counters
+      ! Record and data counters
       inrc = kt - kit000 + 2
-      isla = sladatqc%nsstp(inrc)
+      isurf = surfdataqc%nsstp(inrc)
+
+      ! Work out the maximum footprint size for the 
+      ! interpolation/averaging in model grid-points - has to be even.
+
+      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
+
+
+      IF ( ldnightav ) THEN
+
+      ! Initialize array for night mean
+         IF ( kt == 0 ) THEN
+            ALLOCATE ( icount_night(kpi,kpj) )
+            ALLOCATE ( imask_night(kpi,kpj) )
+            ALLOCATE ( zintmp(kpi,kpj) )
+            ALLOCATE ( zouttmp(kpi,kpj) )
+            ALLOCATE ( zmeanday(kpi,kpj) )
+            nday_qsr = -1   ! initialisation flag for nbc_dcy
+         ENDIF
+
+         ! Night-time means are calculated for night-time values over timesteps:
+         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
+         idayend = MOD( kt - kit000 + 1, kdaystp )
+
+         ! Initialize night-time mean for first timestep of the day
+         IF ( idayend == 1 .OR. kt == 0 ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  surfdataqc%vdmean(ji,jj) = 0.0
+                  zmeanday(ji,jj) = 0.0
+                  icount_night(ji,jj) = 0
+               END DO
+            END DO
+         ENDIF
+
+         zintmp(:,:) = 0.0
+         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
+         imask_night(:,:) = INT( zouttmp(:,:) )
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! Increment the temperature field for computing night mean and counter
+               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
+                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
+               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
+               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
+            END DO
+         END DO
+
+         ! Compute the night-time mean at the end of the day
+         zdaystp = 1.0 / REAL( kdaystp )
+         IF ( idayend == 0 ) THEN
+            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! Test if "no night" point
+                  IF ( icount_night(ji,jj) > 0 ) THEN
+                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
+                       &                        / REAL( icount_night(ji,jj) )
+                  ELSE
+                     !At locations where there is no night (e.g. poles),
+                     ! calculate daily mean instead of night-time mean.
+                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
+                  ENDIF
+               END DO
+            END DO
+         ENDIF
+
+      ENDIF
 
       ! Get the data for interpolation
 
       ALLOCATE( &
-         & igrdi(2,2,isla), &
-         & igrdj(2,2,isla), &
-         & zglam(2,2,isla), &
-         & zgphi(2,2,isla), &
-         & zmask(2,2,isla), &
-         & zsshl(2,2,isla)  &
+         & zweig(imaxifp,imaxjfp,1),      &
+         & igrdi(imaxifp,imaxjfp,isurf), &
+         & igrdj(imaxifp,imaxjfp,isurf), &
+         & zglam(imaxifp,imaxjfp,isurf), &
+         & zgphi(imaxifp,imaxjfp,isurf), &
+         & zmask(imaxifp,imaxjfp,isurf), &
+         & zsurf(imaxifp,imaxjfp,isurf), &
+         & zsurftmp(imaxifp,imaxjfp,isurf),  &
+         & zglamf(imaxifp+1,imaxjfp+1,isurf), &
+         & zgphif(imaxifp+1,imaxjfp+1,isurf), &
+         & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
+         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
          & )
-      
-      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
-         iobs = jobs - sladatqc%nsurfup
-         igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
-         igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
-         igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
-         igrdj(1,2,iobs) = sladatqc%mj(jobs)
-         igrdi(2,1,iobs) = sladatqc%mi(jobs)
-         igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
-         igrdi(2,2,iobs) = sladatqc%mi(jobs)
-         igrdj(2,2,iobs) = sladatqc%mj(jobs)
+
+      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
+         iobs = jobs - surfdataqc%nsurfup
+         DO ji = 0, imaxifp
+            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
+            
+            !Deal with wrap around in longitude
+            IF ( imodi < 1      ) imodi = imodi + jpiglo
+            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
+            
+            DO jj = 0, imaxjfp
+               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
+               !If model values are out of the domain to the north/south then
+               !set them to be the edge of the domain
+               IF ( imodj < 1      ) imodj = 1
+               IF ( imodj > jpjglo ) imodj = jpjglo
+
+               igrdip1(ji+1,jj+1,iobs) = imodi
+               igrdjp1(ji+1,jj+1,iobs) = imodj
+               
+               IF ( ji >= 1 .AND. jj >= 1 ) THEN
+                  igrdi(ji,jj,iobs) = imodi
+                  igrdj(ji,jj,iobs) = imodj
+               ENDIF
+               
+            END DO
+         END DO
       END DO
 
-      CALL obs_int_comm_2d( 2, 2, isla, &
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, glamt, zglam )
-      CALL obs_int_comm_2d( 2, 2, isla, &
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, gphit, zgphi )
-      CALL obs_int_comm_2d( 2, 2, isla, &
-         &                  igrdi, igrdj, psshmask, zmask )
-      CALL obs_int_comm_2d( 2, 2, isla, &
-         &                  igrdi, igrdj, psshn, zsshl )
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
+         &                  igrdi, igrdj, psurfmask, zmask )
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
+         &                  igrdi, igrdj, psurf, zsurf )
+      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
+         &                  igrdip1, igrdjp1, glamf, zglamf )
+      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
+         &                  igrdip1, igrdjp1, gphif, zgphif )
+
+      ! At the end of the day get interpolated means
+      IF ( idayend == 0 .AND. ldnightav ) THEN
+
+         ALLOCATE( &
+            & zsurfm(imaxifp,imaxjfp,isurf)  &
+            & )
+
+         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
+            &               surfdataqc%vdmean(:,:), zsurfm )
+
+      ENDIF
 
       ! Loop over observations
-
-      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
-
-         iobs = jobs - sladatqc%nsurfup
-
-         IF ( kt /= sladatqc%mstp(jobs) ) THEN
-            
+      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
+
+         iobs = jobs - surfdataqc%nsurfup
+
+         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
+
             IF(lwp) THEN
                WRITE(numout,*)
@@ -574,27 +933,58 @@
                WRITE(numout,*) ' Record  = ', jobs,                &
                   &            ' kt      = ', kt,                  &
-                  &            ' mstp    = ', sladatqc%mstp(jobs), &
-                  &            ' ntyp    = ', sladatqc%ntyp(jobs)
+                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
+                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
             ENDIF
-            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
-            
+            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
+
+         ENDIF
+
+         zlam = surfdataqc%rlam(jobs)
+         zphi = surfdataqc%rphi(jobs)
+
+         IF ( ldnightav .AND. idayend == 0 ) THEN
+            ! Night-time averaged data
+            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
+         ELSE
+            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
+         ENDIF
+
+         IF ( k2dint <= 4 ) THEN
+
+            ! Get weights to interpolate the model value to the observation point
+            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
+               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
+               &                   zmask(:,:,iobs), zweig, zobsmask )
+
+            ! Interpolate the model value to the observation point 
+            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
+
+         ELSE
+
+            ! Get weights to average the model SLA to the observation footprint
+            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
+               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
+               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
+               &                   zmask(:,:,iobs), plamscl, pphiscl, &
+               &                   lindegrees, zweig, zobsmask )
+
+            ! Average the model SST to the observation footprint
+            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
+               &              zweig, zsurftmp(:,:,iobs),  zext )
+
+         ENDIF
+
+         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
+            ! ... Remove the MDT from the SSH at the observation point to get the SLA
+            surfdataqc%rext(jobs,1) = zext(1)
+            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
+         ELSE
+            surfdataqc%rmod(jobs,1) = zext(1)
          ENDIF
          
-         zlam = sladatqc%rlam(jobs)
-         zphi = sladatqc%rphi(jobs)
-
-         ! Get weights to interpolate the model SSH to the observation point
-         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
-            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
-            &                   zmask(:,:,iobs), zweig, zobsmask )
-         
-
-         ! Interpolate the model SSH to the observation point
-         CALL obs_int_h2d( 1, 1,      &
-            &              zweig, zsshl(:,:,iobs),  zext )
-         
-         sladatqc%rext(jobs,1) = zext(1)
-         ! ... Remove the MDT at the observation point
-         sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
+         IF ( zext(1) == obfillflt ) THEN
+            ! If the observation value is a fill value, set QC flag to bad
+            surfdataqc%nqc(jobs) = 4
+         ENDIF
 
       END DO
@@ -602,4 +992,5 @@
       ! Deallocate the data for interpolation
       DEALLOCATE( &
+         & zweig, &
          & igrdi, &
          & igrdj, &
@@ -607,836 +998,22 @@
          & zgphi, &
          & zmask, &
-         & zsshl  &
+         & zsurf, &
+         & zsurftmp, &
+         & zglamf, &
+         & zgphif, &
+         & igrdip1,&
+         & igrdjp1 &
          & )
 
-      sladatqc%nsurfup = sladatqc%nsurfup + isla
-
-   END SUBROUTINE obs_sla_opt
-
-   SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
-      &                    psstn, psstmask, k2dint, ld_nightav )
-      !!-----------------------------------------------------------------------
-      !!
-      !!                     ***  ROUTINE obs_sst_opt  ***
-      !!
-      !! ** Purpose : Compute the model counterpart of surface temperature
-      !!              data by interpolating from the model grid to the 
-      !!              observation point.
-      !!
-      !! ** Method  : Linearly interpolate to each observation point using 
-      !!              the model values at the corners of the surrounding grid box.
-      !!
-      !!    The now model SST is first computed at the obs (lon, lat) point.
-      !!
-      !!    Several horizontal interpolation schemes are available:
-      !!        - distance-weighted (great circle) (k2dint = 0)
-      !!        - distance-weighted (small angle)  (k2dint = 1)
-      !!        - bilinear (geographical grid)     (k2dint = 2)
-      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
-      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
-      !!
-      !!
-      !! ** Action  :
-      !!
-      !! History :
-      !!        !  07-07  (S. Ricci ) : Original
-      !!      
-      !!-----------------------------------------------------------------------
-
-      !! * Modules used
-      USE obs_surf_def  ! Definition of storage space for surface observations
-      USE sbcdcy
-
-      IMPLICIT NONE
-
-      !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: &
-         & sstdatqc     ! Subset of surface data not failing screening
-      INTEGER, INTENT(IN) :: kt        ! Time step
-      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
-      INTEGER, INTENT(IN) :: kpj
-      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
-                                       !   (kit000-1 = restart time)
-      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
-      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day  
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & psstn,  &    ! Model SST field
-         & psstmask     ! Land-sea mask
-
-      !! * Local declarations
-      INTEGER :: ji
-      INTEGER :: jj
-      INTEGER :: jobs
-      INTEGER :: inrc
-      INTEGER :: isst
-      INTEGER :: iobs
-      INTEGER :: idayend
-      REAL(KIND=wp) :: zlam
-      REAL(KIND=wp) :: zphi
-      REAL(KIND=wp) :: zext(1), zobsmask(1)
-      REAL(KIND=wp) :: zdaystp
-      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
-         & icount_sstnight,      &
-         & imask_night
-      REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
-         & zintmp, &
-         & zouttmp, & 
-         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
-      REAL(kind=wp), DIMENSION(2,2,1) :: &
-         & zweig
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zmask, &
-         & zsstl, &
-         & zsstm, &
-         & zglam, &
-         & zgphi
-      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi, &
-         & igrdj
-      LOGICAL, INTENT(IN) :: ld_nightav
-
-      !-----------------------------------------------------------------------
-      ! Local initialization 
-      !-----------------------------------------------------------------------
-      ! ... Record and data counters
-      inrc = kt - kit000 + 2
-      isst = sstdatqc%nsstp(inrc)
-
-      IF ( ld_nightav ) THEN
-
-      ! Initialize array for night mean
-
-      IF ( kt .EQ. 0 ) THEN
-         ALLOCATE ( icount_sstnight(kpi,kpj) )
-         ALLOCATE ( imask_night(kpi,kpj) )
-         ALLOCATE ( zintmp(kpi,kpj) )
-         ALLOCATE ( zouttmp(kpi,kpj) )
-         ALLOCATE ( zmeanday(kpi,kpj) )
-         nday_qsr = -1   ! initialisation flag for nbc_dcy
-      ENDIF
-
-      ! Initialize daily mean for first timestep
-      idayend = MOD( kt - kit000 + 1, kdaystp )
-
-      ! Added kt == 0 test to catch restart case 
-      IF ( idayend == 1 .OR. kt == 0) THEN
-         IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               sstdatqc%vdmean(ji,jj) = 0.0
-               zmeanday(ji,jj) = 0.0
-               icount_sstnight(ji,jj) = 0
-            END DO
-         END DO
-      ENDIF
-
-      zintmp(:,:) = 0.0
-      zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
-      imask_night(:,:) = INT( zouttmp(:,:) )
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! Increment the temperature field for computing night mean and counter
-            sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  &
-                   &                        + psstn(ji,jj)*imask_night(ji,jj)
-            zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj)
-            icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
-         END DO
-      END DO
-   
-      ! Compute the daily mean at the end of day
-
-      zdaystp = 1.0 / REAL( kdaystp )
-
-      IF ( idayend == 0 ) THEN 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! Test if "no night" point
-               IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
-                  sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
-                    &                        / icount_sstnight(ji,jj) 
-               ELSE
-                  sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
-               ENDIF
-            END DO
-         END DO
-      ENDIF
-
-      ENDIF
-
-      ! Get the data for interpolation
-      
-      ALLOCATE( &
-         & igrdi(2,2,isst), &
-         & igrdj(2,2,isst), &
-         & zglam(2,2,isst), &
-         & zgphi(2,2,isst), &
-         & zmask(2,2,isst), &
-         & zsstl(2,2,isst)  &
-         & )
-      
-      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
-         iobs = jobs - sstdatqc%nsurfup
-         igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
-         igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
-         igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
-         igrdj(1,2,iobs) = sstdatqc%mj(jobs)
-         igrdi(2,1,iobs) = sstdatqc%mi(jobs)
-         igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
-         igrdi(2,2,iobs) = sstdatqc%mi(jobs)
-         igrdj(2,2,iobs) = sstdatqc%mj(jobs)
-      END DO
-      
-      CALL obs_int_comm_2d( 2, 2, isst, &
-         &                  igrdi, igrdj, glamt, zglam )
-      CALL obs_int_comm_2d( 2, 2, isst, &
-         &                  igrdi, igrdj, gphit, zgphi )
-      CALL obs_int_comm_2d( 2, 2, isst, &
-         &                  igrdi, igrdj, psstmask, zmask )
-      CALL obs_int_comm_2d( 2, 2, isst, &
-         &                  igrdi, igrdj, psstn, zsstl )
-
-      ! At the end of the day get interpolated means
-      IF ( idayend == 0 .AND. ld_nightav ) THEN
-
-         ALLOCATE( &
-            & zsstm(2,2,isst)  &
-            & )
-
-         CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
-            &               sstdatqc%vdmean(:,:), zsstm )
-
-      ENDIF
-
-      ! Loop over observations
-
-      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
-         
-         iobs = jobs - sstdatqc%nsurfup
-         
-         IF ( kt /= sstdatqc%mstp(jobs) ) THEN
-            
-            IF(lwp) THEN
-               WRITE(numout,*)
-               WRITE(numout,*) ' E R R O R : Observation',              &
-                  &            ' time step is not consistent with the', &
-                  &            ' model time step'
-               WRITE(numout,*) ' ========='
-               WRITE(numout,*)
-               WRITE(numout,*) ' Record  = ', jobs,                &
-                  &            ' kt      = ', kt,                  &
-                  &            ' mstp    = ', sstdatqc%mstp(jobs), &
-                  &            ' ntyp    = ', sstdatqc%ntyp(jobs)
-            ENDIF
-            CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
-            
-         ENDIF
-         
-         zlam = sstdatqc%rlam(jobs)
-         zphi = sstdatqc%rphi(jobs)
-         
-         ! Get weights to interpolate the model SST to the observation point
-         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
-            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
-            &                   zmask(:,:,iobs), zweig, zobsmask )
-            
-         ! Interpolate the model SST to the observation point 
-
-         IF ( ld_nightav ) THEN
-
-           IF ( idayend == 0 )  THEN
-               ! Daily averaged/diurnal cycle of SST  data
-               CALL obs_int_h2d( 1, 1,      & 
-                     &              zweig, zsstm(:,:,iobs), zext )
-            ELSE 
-               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     &
-                     &           ' number of night SST data should' // &
-                     &           ' only occur at the end of a given day' )
-            ENDIF
-
-         ELSE
-
-            CALL obs_int_h2d( 1, 1,      &
-            &              zweig, zsstl(:,:,iobs),  zext )
-
-         ENDIF
-         sstdatqc%rmod(jobs,1) = zext(1)
-         
-      END DO
-      
-      ! Deallocate the data for interpolation
-      DEALLOCATE( &
-         & igrdi, &
-         & igrdj, &
-         & zglam, &
-         & zgphi, &
-         & zmask, &
-         & zsstl  &
-         & )
-
-      ! At the end of the day also get interpolated means
-      IF ( idayend == 0 .AND. ld_nightav ) THEN
+      ! At the end of the day also deallocate night-time mean array
+      IF ( idayend == 0 .AND. ldnightav ) THEN
          DEALLOCATE( &
-            & zsstm  &
+            & zsurfm  &
             & )
       ENDIF
-      
-      sstdatqc%nsurfup = sstdatqc%nsurfup + isst
-
-   END SUBROUTINE obs_sst_opt
-
-   SUBROUTINE obs_sss_opt
-      !!-----------------------------------------------------------------------
-      !!
-      !!                     ***  ROUTINE obs_sss_opt  ***
-      !!
-      !! ** Purpose : Compute the model counterpart of sea surface salinity
-      !!              data by interpolating from the model grid to the 
-      !!              observation point.
-      !!
-      !! ** Method  : 
-      !!
-      !! ** Action  :
-      !!
-      !! History :
-      !!      ! ??-?? 
-      !!-----------------------------------------------------------------------
-
-      IMPLICIT NONE
-
-   END SUBROUTINE obs_sss_opt
-
-   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
-      &                    pseaicen, pseaicemask, k2dint )
-
-      !!-----------------------------------------------------------------------
-      !!
-      !!                     ***  ROUTINE obs_seaice_opt  ***
-      !!
-      !! ** Purpose : Compute the model counterpart of surface temperature
-      !!              data by interpolating from the model grid to the 
-      !!              observation point.
-      !!
-      !! ** Method  : Linearly interpolate to each observation point using 
-      !!              the model values at the corners of the surrounding grid box.
-      !!
-      !!    The now model sea ice is first computed at the obs (lon, lat) point.
-      !!
-      !!    Several horizontal interpolation schemes are available:
-      !!        - distance-weighted (great circle) (k2dint = 0)
-      !!        - distance-weighted (small angle)  (k2dint = 1)
-      !!        - bilinear (geographical grid)     (k2dint = 2)
-      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
-      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
-      !!
-      !!
-      !! ** Action  :
-      !!
-      !! History :
-      !!        !  07-07  (S. Ricci ) : Original
-      !!      
-      !!-----------------------------------------------------------------------
-
-      !! * Modules used
-      USE obs_surf_def  ! Definition of storage space for surface observations
-
-      IMPLICIT NONE
-
-      !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
-      INTEGER, INTENT(IN) :: kt       ! Time step
-      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
-      INTEGER, INTENT(IN) :: kpj
-      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step 
-                                      !   (kit000-1 = restart time)
-      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & pseaicen,  &    ! Model sea ice field
-         & pseaicemask     ! Land-sea mask
-         
-      !! * Local declarations
-      INTEGER :: ji
-      INTEGER :: jj
-      INTEGER :: jobs
-      INTEGER :: inrc
-      INTEGER :: iseaice
-      INTEGER :: iobs
-       
-      REAL(KIND=wp) :: zlam
-      REAL(KIND=wp) :: zphi
-      REAL(KIND=wp) :: zext(1), zobsmask(1)
-      REAL(kind=wp), DIMENSION(2,2,1) :: &
-         & zweig
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zmask, &
-         & zseaicel, &
-         & zglam, &
-         & zgphi
-      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi, &
-         & igrdj
-
-      !------------------------------------------------------------------------
-      ! Local initialization 
-      !------------------------------------------------------------------------
-      ! ... Record and data counters
-      inrc = kt - kit000 + 2
-      iseaice = seaicedatqc%nsstp(inrc)
-
-      ! Get the data for interpolation
-      
-      ALLOCATE( &
-         & igrdi(2,2,iseaice), &
-         & igrdj(2,2,iseaice), &
-         & zglam(2,2,iseaice), &
-         & zgphi(2,2,iseaice), &
-         & zmask(2,2,iseaice), &
-         & zseaicel(2,2,iseaice)  &
-         & )
-      
-      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
-         iobs = jobs - seaicedatqc%nsurfup
-         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
-         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
-         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
-         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
-         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
-         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
-         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
-         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
-      END DO
-      
-      CALL obs_int_comm_2d( 2, 2, iseaice, &
-         &                  igrdi, igrdj, glamt, zglam )
-      CALL obs_int_comm_2d( 2, 2, iseaice, &
-         &                  igrdi, igrdj, gphit, zgphi )
-      CALL obs_int_comm_2d( 2, 2, iseaice, &
-         &                  igrdi, igrdj, pseaicemask, zmask )
-      CALL obs_int_comm_2d( 2, 2, iseaice, &
-         &                  igrdi, igrdj, pseaicen, zseaicel )
-      
-      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
-         
-         iobs = jobs - seaicedatqc%nsurfup
-         
-         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
-            
-            IF(lwp) THEN
-               WRITE(numout,*)
-               WRITE(numout,*) ' E R R O R : Observation',              &
-                  &            ' time step is not consistent with the', &
-                  &            ' model time step'
-               WRITE(numout,*) ' ========='
-               WRITE(numout,*)
-               WRITE(numout,*) ' Record  = ', jobs,                &
-                  &            ' kt      = ', kt,                  &
-                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
-                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
-            ENDIF
-            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
-            
-         ENDIF
-         
-         zlam = seaicedatqc%rlam(jobs)
-         zphi = seaicedatqc%rphi(jobs)
-         
-         ! Get weights to interpolate the model sea ice to the observation point
-         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
-            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
-            &                   zmask(:,:,iobs), zweig, zobsmask )
-         
-         ! ... Interpolate the model sea ice to the observation point
-         CALL obs_int_h2d( 1, 1,      &
-            &              zweig, zseaicel(:,:,iobs),  zext )
-         
-         seaicedatqc%rmod(jobs,1) = zext(1)
-         
-      END DO
-      
-      ! Deallocate the data for interpolation
-      DEALLOCATE( &
-         & igrdi,    &
-         & igrdj,    &
-         & zglam,    &
-         & zgphi,    &
-         & zmask,    &
-         & zseaicel  &
-         & )
-      
-      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
-
-   END SUBROUTINE obs_seaice_opt
-
-   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
-      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
-      &                    ld_dailyav )
-      !!-----------------------------------------------------------------------
-      !!
-      !!                     ***  ROUTINE obs_vel_opt  ***
-      !!
-      !! ** Purpose : Compute the model counterpart of velocity profile
-      !!              data by interpolating from the model grid to the 
-      !!              observation point.
-      !!
-      !! ** Method  : Linearly interpolate zonal and meridional components of velocity 
-      !!              to each observation point using the model values at the corners of 
-      !!              the surrounding grid box. The model velocity components are on a 
-      !!              staggered C- grid.
-      !!
-      !!    For velocity data from the TAO array, the model equivalent is
-      !!    a daily mean velocity field. So, we first compute
-      !!    the mean, then interpolate only at the end of the day.
-      !!
-      !! ** Action  :
-      !!
-      !! History :
-      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
-      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
-      !!-----------------------------------------------------------------------
-    
-      !! * Modules used
-      USE obs_profiles_def ! Definition of storage space for profile obs.
-
-      IMPLICIT NONE
-
-      !! * Arguments
-      TYPE(obs_prof), INTENT(INOUT) :: &
-         & prodatqc        ! Subset of profile data not failing screening
-      INTEGER, INTENT(IN) :: kt        ! Time step
-      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
-      INTEGER, INTENT(IN) :: kpj
-      INTEGER, INTENT(IN) :: kpk 
-      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
-                                       !   (kit000-1 = restart time)
-      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
-      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
-      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                    
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
-         & pun,    &    ! Model zonal component of velocity
-         & pvn,    &    ! Model meridional component of velocity
-         & pumask, &    ! Land-sea mask
-         & pvmask       ! Land-sea mask
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
-         & pgdept       ! Model array of depth levels
-      LOGICAL, INTENT(IN) :: ld_dailyav
-         
-      !! * Local declarations
-      INTEGER :: ji
-      INTEGER :: jj
-      INTEGER :: jk
-      INTEGER :: jobs
-      INTEGER :: inrc
-      INTEGER :: ipro
-      INTEGER :: idayend
-      INTEGER :: ista
-      INTEGER :: iend
-      INTEGER :: iobs
-      INTEGER, DIMENSION(imaxavtypes) :: &
-         & idailyavtypes
-      REAL(KIND=wp) :: zlam
-      REAL(KIND=wp) :: zphi
-      REAL(KIND=wp) :: zdaystp
-      REAL(KIND=wp), DIMENSION(kpk) :: &
-         & zobsmasku, &
-         & zobsmaskv, &
-         & zobsmask,  &
-         & zobsk,     &
-         & zobs2k
-      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
-         & zweigu,zweigv
-      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
-         & zumask, zvmask, &
-         & zintu, &
-         & zintv, &
-         & zinmu, &
-         & zinmv
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zglamu, zglamv, &
-         & zgphiu, zgphiv
-      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdiu, &
-         & igrdju, &
-         & igrdiv, &
-         & igrdjv
-
-      !------------------------------------------------------------------------
-      ! Local initialization 
-      !------------------------------------------------------------------------
-      ! ... Record and data counters
-      inrc = kt - kit000 + 2
-      ipro = prodatqc%npstp(inrc)
-
-      ! Initialize daily mean for first timestep
-      idayend = MOD( kt - kit000 + 1, kdaystp )
-
-      ! Added kt == 0 test to catch restart case 
-      IF ( idayend == 1 .OR. kt == 0) THEN
-         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
-         prodatqc%vdmean(:,:,:,1) = 0.0
-         prodatqc%vdmean(:,:,:,2) = 0.0
-      ENDIF
-
-      ! Increment the zonal velocity field for computing daily mean
-      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
-      ! Increment the meridional velocity field for computing daily mean
-      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
-   
-      ! Compute the daily mean at the end of day
-      zdaystp = 1.0 / REAL( kdaystp )
-      IF ( idayend == 0 ) THEN
-         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
-         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
-      ENDIF
-
-      ! Get the data for interpolation
-      ALLOCATE( &
-         & igrdiu(2,2,ipro),      &
-         & igrdju(2,2,ipro),      &
-         & igrdiv(2,2,ipro),      &
-         & igrdjv(2,2,ipro),      &
-         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
-         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
-         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
-         & zintu(2,2,kpk,ipro),  &
-         & zintv(2,2,kpk,ipro)   &
-         & )
-
-      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
-         iobs = jobs - prodatqc%nprofup
-         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
-         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
-         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
-         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
-         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
-         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
-         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
-         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
-         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
-         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
-         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
-         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
-         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
-         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
-      END DO
-
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
-
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
-
-      ! At the end of the day also get interpolated means
-      IF ( idayend == 0 ) THEN
-
-         ALLOCATE( &
-            & zinmu(2,2,kpk,ipro),  &
-            & zinmv(2,2,kpk,ipro)   &
-            & )
-
-         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
-            &                  prodatqc%vdmean(:,:,:,1), zinmu )
-         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
-            &                  prodatqc%vdmean(:,:,:,2), zinmv )
-
-      ENDIF
-
-! loop over observations
-
-      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
-
-         iobs = jobs - prodatqc%nprofup
-
-         IF ( kt /= prodatqc%mstp(jobs) ) THEN
-            
-            IF(lwp) THEN
-               WRITE(numout,*)
-               WRITE(numout,*) ' E R R O R : Observation',              &
-                  &            ' time step is not consistent with the', &
-                  &            ' model time step'
-               WRITE(numout,*) ' ========='
-               WRITE(numout,*)
-               WRITE(numout,*) ' Record  = ', jobs,                    &
-                  &            ' kt      = ', kt,                      &
-                  &            ' mstp    = ', prodatqc%mstp(jobs), &
-                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
-            ENDIF
-            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
-         ENDIF
-         
-         zlam = prodatqc%rlam(jobs)
-         zphi = prodatqc%rphi(jobs)
-
-         ! Initialize observation masks
-
-         zobsmasku(:) = 0.0
-         zobsmaskv(:) = 0.0
-         
-         ! Horizontal weights and vertical mask
-
-         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN
-
-            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
-               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
-               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )
-
-         ENDIF
-
-         
-         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
-
-            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
-               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
-               &                   zvmask(:,:,:,iobs), zweigv, zobsmasku )
-
-         ENDIF
-
-         ! Ensure that the vertical mask on u and v are consistent.
-
-         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
-
-         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
-
-            zobsk(:) = obfillflt
-
-	    IF ( ld_dailyav ) THEN
-
-               IF ( idayend == 0 )  THEN
-                  
-                  ! Daily averaged data
-                  
-                  CALL obs_int_h2d( kpk, kpk,      &
-                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
-                  
-                  
-               ELSE
-               
-                  CALL ctl_stop( ' A nonzero' //     &
-                     &           ' number of U profile data should' // &
-                     &           ' only occur at the end of a given day' )
-
-               ENDIF
-	       
-            ELSE 
-               
-               ! Point data
-
-               CALL obs_int_h2d( kpk, kpk,      &
-                  &              zweigu, zintu(:,:,:,iobs), zobsk )
-
-            ENDIF
-
-            !-------------------------------------------------------------
-            ! Compute vertical second-derivative of the interpolating 
-            ! polynomial at obs points
-            !-------------------------------------------------------------
-            
-            IF ( k1dint == 1 ) THEN
-               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
-                  &                  pgdept, zobsmask )
-            ENDIF
-            
-            !-----------------------------------------------------------------
-            !  Vertical interpolation to the observation point
-            !-----------------------------------------------------------------
-            ista = prodatqc%npvsta(jobs,1)
-            iend = prodatqc%npvend(jobs,1)
-            CALL obs_int_z1d( kpk,                &
-               & prodatqc%var(1)%mvk(ista:iend),  &
-               & k1dint, iend - ista + 1,         &
-               & prodatqc%var(1)%vdep(ista:iend), &
-               & zobsk, zobs2k,                   &
-               & prodatqc%var(1)%vmod(ista:iend), &
-               & pgdept, zobsmask )
-
-         ENDIF
-
-         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
-
-            zobsk(:) = obfillflt
-
-            IF ( ld_dailyav ) THEN
-
-               IF ( idayend == 0 )  THEN
-
-                  ! Daily averaged data
-                  
-                  CALL obs_int_h2d( kpk, kpk,      &
-                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
-                  
-               ELSE
-
-                  CALL ctl_stop( ' A nonzero' //     &
-                     &           ' number of V profile data should' // &
-                     &           ' only occur at the end of a given day' )
-
-               ENDIF
-
-            ELSE
-               
-               ! Point data
-
-               CALL obs_int_h2d( kpk, kpk,      &
-                  &              zweigv, zintv(:,:,:,iobs), zobsk )
-
-            ENDIF
-
-
-            !-------------------------------------------------------------
-            ! Compute vertical second-derivative of the interpolating 
-            ! polynomial at obs points
-            !-------------------------------------------------------------
-            
-            IF ( k1dint == 1 ) THEN
-               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
-                  &                  pgdept, zobsmask )
-            ENDIF
-            
-            !----------------------------------------------------------------
-            !  Vertical interpolation to the observation point
-            !----------------------------------------------------------------
-            ista = prodatqc%npvsta(jobs,2)
-            iend = prodatqc%npvend(jobs,2)
-            CALL obs_int_z1d( kpk, &
-               & prodatqc%var(2)%mvk(ista:iend),&
-               & k1dint, iend - ista + 1, &
-               & prodatqc%var(2)%vdep(ista:iend),&
-               & zobsk, zobs2k, &
-               & prodatqc%var(2)%vmod(ista:iend),&
-               & pgdept, zobsmask )
-
-         ENDIF
-
-      END DO
- 
-      ! Deallocate the data for interpolation
-      DEALLOCATE( &
-         & igrdiu, &
-         & igrdju, &
-         & igrdiv, &
-         & igrdjv, &
-         & zglamu, zglamv, &
-         & zgphiu, zgphiv, &
-         & zumask, zvmask, &
-         & zintu, &
-         & zintv  &
-         & )
-      ! At the end of the day also get interpolated means
-      IF ( idayend == 0 ) THEN
-         DEALLOCATE( &
-            & zinmu,  &
-            & zinmv   &
-            & )
-      ENDIF
-
-      prodatqc%nprofup = prodatqc%nprofup + ipro 
-      
-   END SUBROUTINE obs_vel_opt
+
+      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
+
+   END SUBROUTINE obs_surf_opt
 
 END MODULE obs_oper
-
