Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 4745)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 4746)
@@ -118,4 +118,6 @@
    !!----------------------------------------------------------------------
 
+   !! * Substitutions 
+#  include "domzgr_substitute.h90"
 CONTAINS
 
@@ -207,7 +209,30 @@
       ! Read namelist parameters
       !-----------------------------------------------------------------------
-
+      
+      !Initalise all values in namelist arrays
+      enactfiles(:) = ''
+      coriofiles(:) = ''
+      profbfiles(:) = ''
+      slafilesact(:) = ''
+      slafilespas(:) = ''
+      slafbfiles(:) = ''
+      sstfiles(:)   = ''
+      sstfbfiles(:) = ''
+      seaicefiles(:) = ''
       velcurfiles(:) = ''
       veladcpfiles(:) = ''
+      velavcurfiles(:) = ''
+      velhrcurfiles(:) = ''
+      velavadcpfiles(:) = ''
+      velhradcpfiles(:) = ''
+      velfbfiles(:) = ''
+      velcurfiles(:) = ''
+      veladcpfiles(:) = ''
+      endailyavtypes(:) = -1
+      endailyavtypes(1) = 820
+      ln_profb_ena(:) = .FALSE.
+      ln_profb_enatim(:) = .TRUE.
+      ln_velfb_av(:) = .FALSE.
+      ln_ignmis = .FALSE.      
       CALL ini_date( dobsini )
       CALL fin_date( dobsend )
@@ -972,4 +997,6 @@
       !!        !  07-04  (G. Smith) Generalized surface operators
       !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
+      !!        !  14-08  (J. While) observation operator for profiles in 
+      !!                             generalised vertical coordinates
       !!----------------------------------------------------------------------
       !! * Modules used
@@ -977,4 +1004,9 @@
          & rdt,           &                       
          & gdept_1d,       &             
+#if defined key_vvl 
+         & gdept_n,       & 
+#else 
+         & gdept_1d,         & 
+#endif                                        
          & tmask, umask, vmask                            
       USE phycst, ONLY : &              ! Physical constants
@@ -1034,15 +1066,30 @@
       IF ( ln_t3d .OR. ln_s3d ) THEN
          DO jprofset = 1, nprofsets
-            IF ( ld_enact(jprofset) ) THEN
-               CALL obs_pro_opt( prodatqc(jprofset),                     &
-                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   &
-                  &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &
-                  &              gdept_1d, tmask, n1dint, n2dint,        &
-                  &              kdailyavtypes = endailyavtypes )
+            IF( ln_zco .OR. ln_zps ) THEN 
+               IF ( ld_enact(jprofset) ) THEN 
+                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
+                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
+                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
+                     &              gdept_1d, tmask, n1dint, n2dint,         & 
+                     &              kdailyavtypes = endailyavtypes ) 
+               ELSE 
+                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
+                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
+                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
+                     &              gdept_1d, tmask, n1dint, n2dint               ) 
+               ENDIF 
             ELSE
-               CALL obs_pro_opt( prodatqc(jprofset),                     &
-                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   &
-                  &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &
-                  &              gdept_1d, tmask, n1dint, n2dint              )
+               IF ( ld_enact(jprofset) ) THEN 
+                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
+                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
+                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
+                     &              fsdept(:,:,:), tmask, n1dint, n2dint,   & 
+                     &              kdailyavtypes = endailyavtypes ) 
+               ELSE 
+                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
+                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
+                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
+                     &              fsdept(:,:,:), tmask, n1dint, n2dint ) 
+               ENDIF 
             ENDIF
          END DO
Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 4745)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 4746)
@@ -9,4 +9,7 @@
    !!   obs_pro_opt :    Compute the model counterpart of temperature and
    !!                    salinity observations from profiles
+   !!   obs_pro_sco_opt: Compute the model counterpart of temperature and 
+   !!                    salinity observations from profiles in generalised 
+   !!                    vertical coordinates 
    !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly
    !!                    observations
@@ -37,8 +40,15 @@
    USE dom_oce,       ONLY : &
       & glamt, glamu, glamv, &
-      & gphit, gphiu, gphiv
+      & gphit, gphiu, gphiv, & 
+#if defined key_vvl 
+      & gdept_n 
+#else 
+      & gdept_0 
+#endif  
    USE lib_mpp,       ONLY : &
       & ctl_warn, ctl_stop
-
+   USE obs_grid,      ONLY : & 
+      & obs_level_search     
+      
    IMPLICIT NONE
 
@@ -47,4 +57,6 @@
 
    PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations
+      &   obs_pro_sco_opt, &  ! Compute the model counterpart of profile observations 
+                              ! in generalised vertical coordinates 
       &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations
       &   obs_sst_opt, &  ! Compute the model counterpart of SST observations
@@ -61,4 +73,6 @@
    !!----------------------------------------------------------------------
 
+   !! * Substitutions 
+#  include "domzgr_substitute.h90" 
 CONTAINS
 
@@ -449,4 +463,583 @@
    END SUBROUTINE obs_pro_opt
 
+   SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
+      &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, & 
+      &                    kdailyavtypes ) 
+      !!----------------------------------------------------------------------- 
+      !! 
+      !!                     ***  ROUTINE obs_pro_opt  *** 
+      !! 
+      !! ** Purpose : Compute the model counterpart of profiles 
+      !!              data by interpolating from the model grid to the  
+      !!              observation point. Generalised vertical coordinate version 
+      !! 
+      !! ** Method  : Linearly interpolate to each observation point using  
+      !!              the model values at the corners of the surrounding grid box. 
+      !! 
+      !!          First, model values on the model grid are interpolated vertically to the 
+      !!          Depths of the profile observations.  Two vertical interpolation schemes are 
+      !!          available: 
+      !!          - linear       (k1dint = 0) 
+      !!          - Cubic spline (k1dint = 1)    
+      !! 
+      !! 
+      !!         Secondly the interpolated values are interpolated horizontally to 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) 
+      !! 
+      !!    For the cubic spline the 2nd derivative of the interpolating  
+      !!    polynomial is computed before entering the vertical interpolation  
+      !!    routine. 
+      !! 
+      !!    For ENACT moored buoy data (e.g., TAO), 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. 
+      !! 
+      !!    This is the procedure to be used with generalised vertical model  
+      !!    coordinates (ie s-coordinates. It is ~4x slower than the equivalent 
+      !!    horizontal then vertical interpolation algorithm, but can deal with situations 
+      !!    where the model levels are not flat. 
+      !!    ONLY PERFORMED if ln_sco=.TRUE.  
+      !!       
+      !!    Note: the 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  : 
+      !! 
+      !! History : 
+      !!      ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised 
+      !!                           vertical coordinates
+      !!----------------------------------------------------------------------- 
+   
+      !! * Modules used 
+      USE obs_profiles_def   ! Definition of storage space for profile obs. 
+      USE dom_oce,  ONLY : & 
+#if defined key_vvl  
+      gdepw_n 
+#else 
+      gdepw_0 
+#endif 
+       
+      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) :: & 
+         & ptn,    &    ! Model temperature field 
+         & psn,    &    ! Model salinity field 
+         & ptmask       ! Land-sea mask 
+      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,jpj,kpk) :: & 
+         & pgdept       ! Model array of depth levels 
+      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
+         & kdailyavtypes   ! Types for daily averages 
+      !! * Local declarations 
+      INTEGER ::   ji 
+      INTEGER ::   jj 
+      INTEGER ::   jk 
+      INTEGER ::   iico, ijco 
+      INTEGER ::   jobs 
+      INTEGER ::   inrc 
+      INTEGER ::   ipro 
+      INTEGER ::   idayend 
+      INTEGER ::   ista 
+      INTEGER ::   iend 
+      INTEGER ::   iobs 
+      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes 
+      INTEGER, DIMENSION(imaxavtypes) :: & 
+         & idailyavtypes 
+      REAL(KIND=wp) :: zlam 
+      REAL(KIND=wp) :: zphi 
+      REAL(KIND=wp) :: zdaystp 
+      REAL(KIND=wp), DIMENSION(kpk) :: & 
+         & zobsmask, & 
+         & zobsk,    & 
+         & zobs2k 
+      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
+         & zweig, & 
+         & l_zweig 
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
+         & zmask, & 
+         & zintt, & 
+         & zints, & 
+         & zinmt, & 
+         & zgdept,& 
+         & zgdepw,& 
+         & zinms 
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
+         & zglam, & 
+         & zgphi 
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
+         & igrdi, & 
+         & igrdj 
+      INTEGER :: & 
+         & inum_obs    
+      REAL(KIND=wp), DIMENSION(1) :: zmsk_1       
+      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner       
+      INTEGER, ALLOCATABLE, DIMENSION(:) ::           v_indic    
+ 
+      !------------------------------------------------------------------------ 
+      ! Local initialization  
+      !------------------------------------------------------------------------ 
+      ! ... Record and data counters 
+      inrc = kt - kit000 + 2 
+      ipro = prodatqc%npstp(inrc) 
+  
+      ! Daily average types 
+      IF ( PRESENT(kdailyavtypes) ) THEN 
+         idailyavtypes(:) = kdailyavtypes(:) 
+      ELSE 
+         idailyavtypes(:) = -1 
+      ENDIF 
+ 
+      ! Initialize daily mean for first time-step 
+      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 
+         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 
+               ! 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 
+               END DO 
+            END DO 
+         END DO 
+      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),  & 
+      & 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) 
+      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 ) 
+      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdept(:,:,:), & 
+        &                     zgdept ) 
+      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdepw(:,:,:), & 
+        &                     zgdepw ) 
+ 
+      ! At the end of the day also get interpolated means 
+      IF ( idayend == 0 ) THEN 
+ 
+         ALLOCATE( & 
+         & zinmt(2,2,kpk,ipro),  & 
+         & zinms(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 ) 
+ 
+      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 
+    
+         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) 
+          
+         ! Horizontal weights 
+         ! Only calculated once, for both T and S. 
+         ! Masked values are calculated later.  
+ 
+         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 
+         & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 
+ 
+            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
+            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
+            &                   zmask(:,:,1,iobs), zweig, zmsk_1 ) 
+ 
+         ENDIF 
+         
+         ! IF zmsk_1 = 0; then ob is on land 
+         IF (zmsk_1(1) < 0.1) THEN 
+            WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 
+   
+         ELSE  
+             
+            ! Temperature 
+             
+            IF ( prodatqc%npvend(jobs,1) > 0 ) THEN  
+    
+               zobsk(:) = obfillflt 
+    
+               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
+    
+                  IF ( idayend == 0 )  THEN 
+                   
+                     ! Daily averaged moored buoy (MRB) 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),v_indic(inum_obs)) 
+      
+                     DO iin=1,2 
+                        DO ijn=1,2 
+                                       
+                                       
+           
+                           IF ( k1dint == 1 ) THEN 
+                              CALL obs_int_z1d_spl( kpk, & 
+                              &     zinmt(iin,ijn,:,jobs), & 
+                              &     zobs2k, zgdept(iin,ijn,:,jobs), & 
+                              &     zmask(iin,ijn,:,jobs)) 
+                           ENDIF 
+       
+                           CALL obs_level_search(kpk, & 
+                           &    zgdept(iin,ijn,:,jobs), & 
+                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                           &    v_indic) 
+                           CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, & 
+                           &    prodatqc%var(1)%vdep(ista:iend), & 
+                           &    zinmt(iin,ijn,:,jobs), & 
+                           &    zobs2k, interp_corner(iin,ijn,:), & 
+                           &    zgdept(iin,ijn,:,jobs), & 
+                           &    zmask(iin,ijn,:,jobs)) 
+       
+                        ENDDO 
+                     ENDDO 
+                   
+                   
+                  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 
+     
+                  ! 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), v_indic(inum_obs)) 
+                  DO iin=1,2  
+                     DO ijn=1,2 
+                                    
+                                    
+                        IF ( k1dint == 1 ) THEN 
+                           CALL obs_int_z1d_spl( kpk, & 
+                           &    zintt(iin,ijn,:,jobs),& 
+                           &    zobs2k, zgdept(iin,ijn,:,jobs), & 
+                           &    zmask(iin,ijn,:,jobs)) 
+  
+                        ENDIF 
+       
+                        CALL obs_level_search(kpk, & 
+                         &        zgdept(iin,ijn,:,jobs),& 
+                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                         &         v_indic) 
+                        CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs,     & 
+                         &          prodatqc%var(1)%vdep(ista:iend),     & 
+                         &          zintt(iin,ijn,:,jobs),            & 
+                         &          zobs2k,interp_corner(iin,ijn,:), & 
+                         &          zgdept(iin,ijn,:,jobs),         & 
+                         &          zmask(iin,ijn,:,jobs) )      
+         
+                     ENDDO 
+                  ENDDO 
+             
+               ENDIF 
+       
+               !------------------------------------------------------------- 
+               ! Compute the horizontal interpolation for every profile level 
+               !------------------------------------------------------------- 
+             
+               DO ikn=1,inum_obs 
+                  iend=ista+ikn-1 
+   
+                  ! 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 are near  
+                  ! steep bathymetry 
+                  DO iin=1,2 
+                     DO ijn=1,2 
+     
+                        depth_loop1: DO ik=kpk,2,-1 
+                           IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN   
+                            
+                              l_zweig(iin,ijn,1) = &  
+                           &  zweig(iin,ijn,1) *   & 
+                           &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) & 
+                           &  - prodatqc%var(1)%vdep(iend)),0._wp) 
+                            
+                              EXIT depth_loop1 
+                           ENDIF 
+                        ENDDO depth_loop1 
+     
+                     ENDDO 
+                  ENDDO 
+   
+                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 
+                  &          prodatqc%var(1)%vmod(iend:iend) ) 
+ 
+               ENDDO 
+ 
+ 
+               DEALLOCATE(interp_corner,v_indic) 
+          
+            ENDIF 
+       
+ 
+            ! Salinity  
+          
+            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 
+                   
+                     ! 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),v_indic(inum_obs)) 
+      
+                     DO iin=1,2 
+                        DO ijn=1,2 
+                                       
+                                       
+           
+                           IF ( k1dint == 1 ) THEN 
+                              CALL obs_int_z1d_spl( kpk, & 
+                              &     zinms(iin,ijn,:,jobs), & 
+                              &     zobs2k, zgdept(iin,ijn,:,jobs), & 
+                              &     zmask(iin,ijn,:,jobs)) 
+                           ENDIF 
+       
+                           CALL obs_level_search(kpk, & 
+                           &    zgdept(iin,ijn,:,jobs), & 
+                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
+                           &    v_indic) 
+                           CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, & 
+                           &    prodatqc%var(2)%vdep(ista:iend), & 
+                           &    zinms(iin,ijn,:,jobs), & 
+                           &    zobs2k, interp_corner(iin,ijn,:), & 
+                           &    zgdept(iin,ijn,:,jobs), & 
+                           &    zmask(iin,ijn,:,jobs)) 
+       
+                        ENDDO 
+                     ENDDO 
+                   
+                   
+                  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 
+     
+                  ! 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), v_indic(inum_obs)) 
+                   
+                  DO iin=1,2     
+                     DO ijn=1,2  
+                                 
+                                 
+                        IF ( k1dint == 1 ) THEN 
+                           CALL obs_int_z1d_spl( kpk, & 
+                           &    zints(iin,ijn,:,jobs),& 
+                           &    zobs2k, zgdept(iin,ijn,:,jobs), & 
+                           &    zmask(iin,ijn,:,jobs)) 
+  
+                        ENDIF 
+       
+                        CALL obs_level_search(kpk, & 
+                         &        zgdept(iin,ijn,:,jobs),& 
+                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
+                         &         v_indic) 
+                        CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs,     & 
+                         &          prodatqc%var(2)%vdep(ista:iend),     & 
+                         &          zints(iin,ijn,:,jobs),            & 
+                         &          zobs2k,interp_corner(iin,ijn,:), & 
+                         &          zgdept(iin,ijn,:,jobs),         & 
+                         &          zmask(iin,ijn,:,jobs) )      
+         
+                     ENDDO 
+                  ENDDO 
+             
+               ENDIF 
+       
+               !------------------------------------------------------------- 
+               ! Compute the horizontal interpolation for every profile level 
+               !------------------------------------------------------------- 
+             
+               DO ikn=1,inum_obs 
+                  iend=ista+ikn-1 
+   
+                  ! 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 are near  
+                  ! steep bathymetry 
+                  DO iin=1,2 
+                     DO ijn=1,2 
+     
+                        depth_loop2: DO ik=kpk,2,-1 
+                           IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN   
+                            
+                              l_zweig(iin,ijn,1) = &  
+                           &  zweig(iin,ijn,1) *   & 
+                           &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) & 
+                           &  - prodatqc%var(2)%vdep(iend)),0._wp) 
+                            
+                              EXIT depth_loop2 
+                           ENDIF 
+                        ENDDO depth_loop2 
+     
+                     ENDDO 
+                  ENDDO 
+   
+                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 
+                  &          prodatqc%var(2)%vmod(iend:iend) ) 
+ 
+               ENDDO 
+ 
+ 
+               DEALLOCATE(interp_corner,v_indic) 
+          
+            ENDIF 
+          
+         ENDIF 
+       
+      END DO 
+     
+      ! Deallocate the data for interpolation 
+      DEALLOCATE( & 
+         & igrdi, & 
+         & igrdj, & 
+         & zglam, & 
+         & zgphi, & 
+         & zmask, & 
+         & zintt, & 
+         & zints  & 
+         & ) 
+      ! At the end of the day also get interpolated means 
+      IF ( idayend == 0 ) THEN 
+         DEALLOCATE( & 
+            & zinmt,  & 
+            & zinms   & 
+            & ) 
+      ENDIF 
+    
+      prodatqc%nprofup = prodatqc%nprofup + ipro  
+       
+   END SUBROUTINE obs_pro_sco_opt 
+ 
    SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
       &                    psshn, psshmask, k2dint )
@@ -1293,5 +1886,5 @@
             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
                &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
-               &                   zvmask(:,:,:,iobs), zweigv, zobsmasku )
+               &                   zvmask(:,:,:,iobs), zweigv, zobsmaskv )
 
          ENDIF
