Index: branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90	(revision 13645)
+++ branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90	(revision 15415)
@@ -13,5 +13,8 @@
       !! ** Method  : Straightforward search
       !!
-      !! ** Action  : 
+      !! ** Action  : Will return level associated with T-point below the obs
+      !!              depth, except when observation is in the top box will 
+      !!              return level 2. Also, if obs depth greater than depth 
+      !!              of last wet T-point (kpk-1) will return level kpk.
       !!
       !! History :
@@ -43,6 +46,6 @@
       DO ji = 1, kobs 
          kobsk(ji) = 1
-         depk: DO jk = 2, kgrd
-            IF ( pgrddep(jk) >= pobsdep(ji) ) EXIT depk
+         depk: DO jk = 2, kgrd-1
+            IF ( pgrddep(jk) > pobsdep(ji) ) EXIT depk
          END DO depk
          kobsk(ji) = jk
Index: branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 13645)
+++ branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 15415)
@@ -16,4 +16,5 @@
    USE par_kind, ONLY : & ! Precision variables
       & wp   
+   USE dom_oce            ! ocean space and time domain
    USE in_out_manager     ! I/O manager
    USE obs_profiles_def   ! Definitions for storage arrays for profiles
@@ -1163,5 +1164,6 @@
       REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
          & zglam, &           ! Model longitude at grid points
-         & zgphi              ! Model latitude at grid points
+         & zgphi, &           ! Model latitude at grid points
+         & zbathy             ! Index of deepest wet level at grid points
       INTEGER, DIMENSION(2,2,kprofno) :: &
          & igrdi, &           ! Grid i,j
@@ -1171,4 +1173,5 @@
       INTEGER :: iig, ijg           ! i,j of observation on model grid point.
       INTEGER :: jobs, jobsp, jk, ji, jj
+      REAL(KIND=wp) :: maxdepw
 
       ! Get grid point indices
@@ -1221,4 +1224,5 @@
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
+      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, REAL(mbathy), zbathy )
       ! Need to know the bathy depth for each observation for sco
       CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), &
@@ -1261,4 +1265,14 @@
          DO jobsp = kpstart(jobs), kpend(jobs)
 
+            ! Calculate max T and W depths of 2x2 grid
+            maxdepw=zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs)
+            DO jj = 1, 2
+               DO ji = 1, 2
+                  IF ( zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) > maxdepw ) THEN
+                     maxdepw = zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs)
+                  END IF
+               END DO
+            END DO
+
             ! Flag if the observation falls outside the model spatial domain
             IF (       ( pobslam(jobs) < -180.         )       &
@@ -1267,5 +1281,5 @@
                &  .OR. ( pobsphi(jobs) >   90.         )       &
                &  .OR. ( pobsdep(jobsp) < 0.0          )       &
-               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
+               &  .OR. ( pobsdep(jobsp) >= maxdepw ) ) THEN
                kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
                kosdobs = kosdobs + 1
@@ -1329,10 +1343,5 @@
                ENDIF
             ENDIF
-            
-            ! Set observation depth equal to that of the first model depth
-            IF ( pobsdep(jobsp) < MINVAL(zgdept(1:2,1:2,1,jobs) ) ) THEN
-               pobsdep(jobsp) = MINVAL(zgdept(1:2,1:2,1,jobs))
-            ENDIF
-            
+           
 #if defined key_bdy
             ! Flag if the observation falls close to the boundary rim
Index: branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90	(revision 13645)
+++ branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90	(revision 15415)
@@ -62,11 +62,31 @@
          z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
          z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
+
          
-         ! If kkco(jdep) is masked then set pobs(jdep) to the lowest value located above bathymetry
-         IF ( pobsmask(kkco(jdep)) == 0.0_wp ) THEN
-            pobs(jdep) = pobsk(kkco(jdep)-1)
+         ! Where both levels are masked, return a fill value
+         IF ( ( pobsmask(kkco(jdep)-1) == 0.0_wp ) .AND. (pobsmask(kkco(jdep)) == 0.0_wp) ) THEN
+            pobs(jdep)  = 99999.
          ELSE
+         
+            ! Where upper level is masked (e.g., under ice cavity), only use deeper level
+            ! otherwise where ob is at or above upper level model T-point, 
+            ! use upper model level rather than extrapolate
+            IF ( pobsmask(kkco(jdep)-1) == 0.0_wp ) THEN
+               z1dm = 0.0_wp
+            ELSE IF ( pobsdep(jdep) <= pdep(kkco(jdep)-1) ) THEN
+               z1dp = 0.0_wp   
+            END IF   
+
+            ! Where deeper level is masked (e.g., near sea bed), only use upper level
+            ! otherwise where ob is at or below deeper level model T-point, 
+            ! use deeper model level rather than extrapolate
+            IF ( pobsmask(kkco(jdep)) == 0.0_wp ) THEN
+               z1dp = 0.0_wp
+            ELSE IF ( pobsdep(jdep) >= pdep(kkco(jdep)) ) THEN
+               z1dm = 0.0_wp   
+            END IF   
+
             zsum = z1dm + z1dp
-
+         
             IF ( k1dint == 0 ) THEN
 
