Changeset 15400

Timestamp:

2021-10-19T11:34:55+02:00 (3 years ago)

Author:

kingr

Message:

Bug-fix to obsoper to avoid division by zero in rare cases.

Location:

branches/UKMO/dev_r5518_obs_oper_update/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/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/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/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90

@@ -62 +62 @@
          z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
          z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
-         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
-
-         ENDIF
+         ! 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
+
+               !-----------------------------------------------------------------
+               !  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
+         ENDIF   
       END DO