Changeset 15356

Timestamp:

2021-10-12T15:52:10+02:00 (3 years ago)

Author:

kingr

Message:

Chnages to fix test in obs_level_search to ensure index of level below obdepth is returned, and to obs_prep to remove code which altered obdepths, and to vertical interpolation routine to deal with obs in top half of top box or bottom half of bottom box.

Location:

branches/UKMO/dev_r5518_obs_oper_update_DepthBug/NEMOGCM/NEMO/OPA_SRC/OBS

Files:

• obs_level_search.h90 (modified) (1 diff)
• obs_prep.F90 (modified) (3 diffs)
• obsinter_z1d.h90 (modified) (1 diff)

branches/UKMO/dev_r5518_obs_oper_update_DepthBug/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90

@@ -47 +47 @@
          kobsk(ji) = 1
          depk: DO jk = 2, kgrd-1
-            IF ( pgrddep(jk) >= pobsdep(ji) ) EXIT depk
+            IF ( pgrddep(jk) > pobsdep(ji) ) EXIT depk
          END DO depk
          kobsk(ji) = jk

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

@@ -1224 +1224 @@
       INTEGER :: iig, ijg           ! i,j of observation on model grid point.
       INTEGER :: jobs, jobsp, jk, ji, jj
-      REAL(KIND=wp) :: maxdept, maxdepw
+      REAL(KIND=wp) :: maxdepw
 
       ! Get grid point indices
@@ -1317 +1317 @@
 
             ! 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)
@@ -1398 +1394 @@
                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
-
-            ! 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
             ! Flag if the observation falls close to the boundary rim

branches/UKMO/dev_r5518_obs_oper_update_DepthBug/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90

@@ -62 +62 @@
          z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
          z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
+         
+         ! Where ob is at or above top model T-point, use model top rather than extrapolate
+         IF ( pobsdep(jdep) <= pdep(kkco(jdep)-1) ) z1dp = 0.0_wp
+         
+         ! Where ob is at or below bottom model T-point, use model bottom rather than extrapolate
+         IF ( pobsdep(jdep) >= pdep(kkco(jdep)) ) z1dm = 0.0_wp
+
+         ! Where deeper level is masked, only use higher level
          IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp
 
          zsum = z1dm + z1dp
          
-         IF ( k1dint == 0 ) THEN
-
-            !-----------------------------------------------------------------
-            !  Linear interpolation
-            !-----------------------------------------------------------------
-            pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &
-               &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum
-
-         ELSEIF ( k1dint == 1 ) THEN
-
-            !-----------------------------------------------------------------
-            ! Cubic spline interpolation
-            !-----------------------------------------------------------------
-            zsum2 = zsum * zsum
-            pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) &
-               &           + z1dp                             * pobsk (kkco(jdep)  ) &
-               &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) &
-               &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) &
-               &             ) / 6.0_wp                                              &
-               &          ) / zsum
-
+         IF ( zsum > 0.0_wp ) THEN
+            IF ( k1dint == 0 ) THEN
+
+               !-----------------------------------------------------------------
+               !  Linear interpolation
+               !-----------------------------------------------------------------
+               pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &
+                  &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum
+
+            ELSEIF ( k1dint == 1 ) THEN
+
+               !-----------------------------------------------------------------
+               ! Cubic spline interpolation
+               !-----------------------------------------------------------------
+               zsum2 = zsum * zsum
+               pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) &
+                  &           + z1dp                             * pobsk (kkco(jdep)  ) &
+                  &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) &
+                  &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) &
+                  &             ) / 6.0_wp                                              &
+                  &          ) / zsum
+
+            ENDIF
+         ELSE ! take value directly from the higher model level
+            pobs(jdep)  = pobsk(kkco(jdep)-1)
          ENDIF
       END DO