New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15400 for branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90 – NEMO

Ignore:
Timestamp:
2021-10-19T11:34:55+02:00 (11 months ago)
Author:
kingr
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90

    r7960 r15400  
    6262         z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      ) 
    6363         z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) ) 
    64          IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp 
    65  
    66          zsum = z1dm + z1dp 
     64 
    6765          
    68          IF ( k1dint == 0 ) THEN 
    69  
    70             !----------------------------------------------------------------- 
    71             !  Linear interpolation 
    72             !----------------------------------------------------------------- 
    73             pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) & 
    74                &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum 
    75  
    76          ELSEIF ( k1dint == 1 ) THEN 
    77  
    78             !----------------------------------------------------------------- 
    79             ! Cubic spline interpolation 
    80             !----------------------------------------------------------------- 
    81             zsum2 = zsum * zsum 
    82             pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) & 
    83                &           + z1dp                             * pobsk (kkco(jdep)  ) & 
    84                &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) & 
    85                &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) & 
    86                &             ) / 6.0_wp                                              & 
    87                &          ) / zsum 
    88  
    89          ENDIF 
     66         ! Where both levels are masked, return a fill value 
     67         IF ( ( pobsmask(kkco(jdep)-1) == 0.0_wp ) .AND. (pobsmask(kkco(jdep)) == 0.0_wp) ) THEN 
     68            pobs(jdep)  = 99999. 
     69         ELSE 
     70          
     71            ! Where upper level is masked (e.g., under ice cavity), only use deeper level 
     72            ! otherwise where ob is at or above upper level model T-point,  
     73            ! use upper model level rather than extrapolate 
     74            IF ( pobsmask(kkco(jdep)-1) == 0.0_wp ) THEN 
     75               z1dm = 0.0_wp 
     76            ELSE IF ( pobsdep(jdep) <= pdep(kkco(jdep)-1) ) THEN 
     77               z1dp = 0.0_wp    
     78            END IF    
     79 
     80            ! Where deeper level is masked (e.g., near sea bed), only use upper level 
     81            ! otherwise where ob is at or below deeper level model T-point,  
     82            ! use deeper model level rather than extrapolate 
     83            IF ( pobsmask(kkco(jdep)) == 0.0_wp ) THEN 
     84               z1dp = 0.0_wp 
     85            ELSE IF ( pobsdep(jdep) >= pdep(kkco(jdep)) ) THEN 
     86               z1dm = 0.0_wp    
     87            END IF    
     88 
     89            zsum = z1dm + z1dp 
     90          
     91            IF ( k1dint == 0 ) THEN 
     92 
     93               !----------------------------------------------------------------- 
     94               !  Linear interpolation 
     95               !----------------------------------------------------------------- 
     96               pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) & 
     97                  &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum 
     98 
     99            ELSEIF ( k1dint == 1 ) THEN 
     100 
     101               !----------------------------------------------------------------- 
     102               ! Cubic spline interpolation 
     103               !----------------------------------------------------------------- 
     104               zsum2 = zsum * zsum 
     105               pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) & 
     106                  &           + z1dp                             * pobsk (kkco(jdep)  ) & 
     107                  &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) & 
     108                  &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) & 
     109                  &             ) / 6.0_wp                                              & 
     110                  &          ) / zsum 
     111 
     112            ENDIF 
     113         ENDIF    
    90114      END DO 
    91115 
Note: See TracChangeset for help on using the changeset viewer.