Changeset 12506 for branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

Timestamp:

2020-03-04T13:11:00+01:00 (4 years ago)

Author:

kingr

Message:

Bug-fix to avoid near-bottom profile obs being associated with a level below the bathymetry and leading to reading beyond the end of arrays.

File:

• branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 (modified) (7 diffs)

branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

@@ -16 +16 @@
    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 +1164 @@
       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 +1173 @@
       INTEGER :: iig, ijg           ! i,j of observation on model grid point.
       INTEGER :: jobs, jobsp, jk, ji, jj
+      REAL(KIND=wp) :: maxdept, maxdepw
 
       ! Get grid point indices
@@ -1221 +1224 @@
       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 +1265 @@
          DO jobsp = kpstart(jobs), kpend(jobs)
 
+            ! Calculate max T and W depths of 2x2 grid
+            maxdept=zgdept(1,1,NINT(zbathy(1,1,jobs)),jobs)
+            maxdepw=zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs)
+            DO jj = 1, 2
+               DO ji = 1, 2
+                  IF ( zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs) > maxdept ) THEN
+                     maxdept = zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs)
+                  END IF
+                  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 +1285 @@
                &  .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
@@ -1334 +1352 @@
                pobsdep(jobsp) = MINVAL(zgdept(1:2,1:2,1,jobs))
             ENDIF
+
+            ! Set observation depth equal to that of the last wet T-point
+            IF ( ( pobsdep(jobsp) > maxdept ) .AND. &
+               & ( pobsdep(jobsp) < maxdepw ) ) THEN
+               pobsdep(jobsp) = maxdept
+            END IF
             
 #if defined key_bdy