Changeset 12506


Ignore:
Timestamp:
2020-03-04T13:11:00+01:00 (8 months 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.

Location:
branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
2 edited

Legend:

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

    r7960 r12506  
    1313      !! ** Method  : Straightforward search 
    1414      !! 
    15       !! ** Action  :  
     15      !! ** Action  : Will return level associated with T-point below the obs 
     16      !!              depth, except when observation is in the top box will  
     17      !!              return level 2. Also, if obs depth greater than depth  
     18      !!              of last wet T-point (kpk-1) will return level kpk. 
    1619      !! 
    1720      !! History : 
     
    4346      DO ji = 1, kobs  
    4447         kobsk(ji) = 1 
    45          depk: DO jk = 2, kgrd 
     48         depk: DO jk = 2, kgrd-1 
    4649            IF ( pgrddep(jk) >= pobsdep(ji) ) EXIT depk 
    4750         END DO depk 
  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r11461 r12506  
    1616   USE par_kind, ONLY : & ! Precision variables 
    1717      & wp    
     18   USE dom_oce            ! ocean space and time domain 
    1819   USE in_out_manager     ! I/O manager 
    1920   USE obs_profiles_def   ! Definitions for storage arrays for profiles 
     
    11631164      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
    11641165         & zglam, &           ! Model longitude at grid points 
    1165          & zgphi              ! Model latitude at grid points 
     1166         & zgphi, &           ! Model latitude at grid points 
     1167         & zbathy             ! Index of deepest wet level at grid points 
    11661168      INTEGER, DIMENSION(2,2,kprofno) :: & 
    11671169         & igrdi, &           ! Grid i,j 
     
    11711173      INTEGER :: iig, ijg           ! i,j of observation on model grid point. 
    11721174      INTEGER :: jobs, jobsp, jk, ji, jj 
     1175      REAL(KIND=wp) :: maxdept, maxdepw 
    11731176 
    11741177      ! Get grid point indices 
     
    12211224      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 
    12221225      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
     1226      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, REAL(mbathy), zbathy ) 
    12231227      ! Need to know the bathy depth for each observation for sco 
    12241228      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), & 
     
    12611265         DO jobsp = kpstart(jobs), kpend(jobs) 
    12621266 
     1267            ! Calculate max T and W depths of 2x2 grid 
     1268            maxdept=zgdept(1,1,NINT(zbathy(1,1,jobs)),jobs) 
     1269            maxdepw=zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs) 
     1270            DO jj = 1, 2 
     1271               DO ji = 1, 2 
     1272                  IF ( zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs) > maxdept ) THEN 
     1273                     maxdept = zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs) 
     1274                  END IF 
     1275                  IF ( zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) > maxdepw ) THEN 
     1276                     maxdepw = zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) 
     1277                  END IF 
     1278               END DO 
     1279            END DO 
     1280 
    12631281            ! Flag if the observation falls outside the model spatial domain 
    12641282            IF (       ( pobslam(jobs) < -180.         )       & 
     
    12671285               &  .OR. ( pobsphi(jobs) >   90.         )       & 
    12681286               &  .OR. ( pobsdep(jobsp) < 0.0          )       & 
    1269                &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 
     1287               &  .OR. ( pobsdep(jobsp) >= maxdepw ) ) THEN 
    12701288               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 
    12711289               kosdobs = kosdobs + 1 
     
    13341352               pobsdep(jobsp) = MINVAL(zgdept(1:2,1:2,1,jobs)) 
    13351353            ENDIF 
     1354 
     1355            ! Set observation depth equal to that of the last wet T-point 
     1356            IF ( ( pobsdep(jobsp) > maxdept ) .AND. & 
     1357               & ( pobsdep(jobsp) < maxdepw ) ) THEN 
     1358               pobsdep(jobsp) = maxdept 
     1359            END IF 
    13361360             
    13371361#if defined key_bdy 
Note: See TracChangeset for help on using the changeset viewer.