Changeset 11256 for branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis3/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90

Timestamp:

2019-07-11T16:06:20+02:00 (5 years ago)

Author:

rrenshaw

Message:

bug fix for obs below model bathymetry

File:

• branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis3/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90 (modified) (5 diffs)

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

@@ -26 +26 @@
       !!      ! 06-10 (A. Weaver) Cleanup
       !!      ! 07-01 (K. Mogensen) Use profile rather than single level
+      !!      ! 19-07 (R. Renshaw) Checks on kkco and zsum to avoid occasional NaNs
       !!---------------------------------------------------------------------
 
@@ -50 +51 @@
       REAL(KIND=wp) :: zsum2
       INTEGER :: jdep             ! Observation depths loop variable
+      INTEGER :: kkco_fix        ! kkco(jdep) limited to number of model levels
     
       !------------------------------------------------------------------------
@@ -60 +62 @@
          ! Initialization
          !---------------------------------------------------------------------
-         z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
-         z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
-         IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp
+!        kkco (from obs_level_search) sometimes has model level kpk+1 which doesn't exist
+         kkco_fix = MIN( kkco(jdep), kpk )
+         z1dm = ( pdep(kkco_fix) - pobsdep(jdep)      )
+         z1dp = ( pobsdep(jdep)    - pdep(kkco_fix-1) )
+!        Where ob is below model bottom, use model bottom and don't extrapolate
+         IF ( pdep(kkco_fix) <= pobsdep(jdep) ) z1dm = 0.0_wp
+!        Where lower level is missing, use higher level
+         IF ( pobsmask(kkco_fix) == 0.0_wp ) z1dp = 0.0_wp
 
          zsum = z1dm + z1dp
+         ! if pobsmask==0 and model level depth==observed depth, we get zsum=0
+         IF ( zsum > 0.0_wp ) THEN
          
          IF ( k1dint == 0 ) THEN
@@ -71 +80 @@
             !  Linear interpolation
             !-----------------------------------------------------------------
-            pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &
-               &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum
+            pobs(jdep) = (   z1dm * pobsk(kkco_fix-1) &
+               &           + z1dp * pobsk(kkco_fix  ) ) / zsum
 
          ELSEIF ( k1dint == 1 ) THEN
@@ -80 +89 @@
             !-----------------------------------------------------------------
             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)  ) &
+            pobs(jdep)  = (  z1dm                             * pobsk (kkco_fix-1) &
+               &           + z1dp                             * pobsk (kkco_fix  ) &
+               &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco_fix-1) &
+               &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco_fix  ) &
                &             ) / 6.0_wp                                              &
                &          ) / zsum
 
          ENDIF
+
+         ELSE  ! take value directly from the higher model level
+            pobs(jdep)  = pobsk(kkco_fix-1)
+         ENDIF
+
       END DO