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 7773 for branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 – NEMO

Ignore:
Timestamp:
2017-03-09T13:52:43+01:00 (7 years ago)
Author:
mattmartin
Message:

Committing updates after doing the following:

  • merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
  • updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
    • generalised vertical coordinate for profile obs. This was done so that is now the default option.
    • sst bias correction implemented with the new simplified obs code.
    • included the biogeochemical obs types int he new simplified obs code.
    • included the changes to exclude obs in the boundary for limited area models
    • included other changes for the efficiency of the obs operator to remove global arrays.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    • Property svn:keywords deleted
    r5785 r7773  
    2424   USE obs_inter_sup      ! Interpolation support 
    2525   USE obs_oper           ! Observation operators 
     26#if defined key_bdy 
     27   USE bdy_oce, ONLY : &        ! Boundary information 
     28      idx_bdy, nb_bdy 
     29#endif 
    2630   USE lib_mpp, ONLY : & 
    2731      & ctl_warn, ctl_stop 
    … …  
    4549CONTAINS 
    4650 
    47    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 
     51   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject ) 
    4852      !!---------------------------------------------------------------------- 
    4953      !!                    ***  ROUTINE obs_pre_sla  *** 
    … …  
    7276      !! * Arguments 
    7377      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data 
    74       TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening 
    75       LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
     78      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc  ! Subset of surface data not failing screening 
     79      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land 
     80      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
    7681      !! * Local declarations 
    7782      INTEGER :: iyea0        ! Initial date 
    … …  
    8792      INTEGER :: inlasobs     !  - close to land 
    8893      INTEGER :: igrdobs      !  - fail the grid search 
     94      INTEGER :: ibdysobs     !  - close to open boundary 
    8995                              ! Global counters for observations that 
    9096      INTEGER :: iotdobsmpp     !  - outside time domain 
    … …  
    9399      INTEGER :: inlasobsmpp    !  - close to land 
    94100      INTEGER :: igrdobsmpp     !  - fail the grid search 
     101      INTEGER :: ibdysobsmpp  !  - close to open boundary 
    95102      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
    96103         & llvalid            ! SLA data selection 
    … …  
    118125      ilansobs = 0 
    119126      inlasobs = 0 
     127      ibdysobs = 0  
    120128 
    121129      ! ----------------------------------------------------------------------- 
    … …  
    151159         &                 tmask(:,:,1), surfdata%nqc,  & 
    152160         &                 iosdsobs,     ilansobs,     & 
    153          &                 inlasobs,     ld_nea        ) 
     161         &                 inlasobs,     ld_nea,       & 
     162         &                 ibdysobs,     ld_bound_reject        ) 
    154163 
    155164      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    156165      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    157166      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     167      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
    158168 
    159169      ! ----------------------------------------------------------------------- 
    … …  
    201211               &            inlasobsmpp 
    202212         ENDIF 
     213         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
     214            &            ibdysobsmpp   
    203215         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    204216            &            surfdataqc%nsurfmpp 
    … …  
    236248      &                     kpi, kpj, kpk, & 
    237249      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
    238       &                     ld_nea, kdailyavtypes ) 
     250      &                     ld_nea, ld_bound_reject, kdailyavtypes ) 
    239251 
    240252!!---------------------------------------------------------------------- 
    … …  
    265277      LOGICAL, INTENT(IN) :: ld_var2 
    266278      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land 
     279      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary 
    267280      INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes 
    268281      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    … …  
    292305      INTEGER :: inlav1obs    !  - close to land (variable 1) 
    293306      INTEGER :: inlav2obs    !  - close to land (variable 2) 
     307      INTEGER :: ibdyv1obs    !  - boundary (variable 1)  
     308      INTEGER :: ibdyv2obs    !  - boundary (variable 2)       
    294309      INTEGER :: igrdobs      !  - fail the grid search 
    295310      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa 
    … …  
    303318      INTEGER :: inlav1obsmpp !  - close to land (variable 1) 
    304319      INTEGER :: inlav2obsmpp !  - close to land (variable 2) 
     320      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1)  
     321      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)       
    305322      INTEGER :: igrdobsmpp   !  - fail the grid search 
    306323      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa 
    … …  
    328345      ! Diagnotics counters for various failures. 
    329346 
    330       iotdobs  = 0 
    331       igrdobs  = 0 
     347      iotdobs   = 0 
     348      igrdobs   = 0 
    332349      iosdv1obs = 0 
    333350      iosdv2obs = 0 
    … …  
    336353      inlav1obs = 0 
    337354      inlav2obs = 0 
    338       iuvchku  = 0 
    339       iuvchkv = 0 
     355      ibdyv1obs = 0 
     356      ibdyv2obs = 0 
     357      iuvchku   = 0 
     358      iuvchkv   = 0 
    340359 
    341360      ! ----------------------------------------------------------------------- 
    … …  
    395414         &                 gdept_1d,              zmask1,               & 
    396415         &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    397          &                 iosdv1obs,              ilanv1obs,           & 
    398          &                 inlav1obs,              ld_nea                ) 
     416         &                 iosdv1obs,             ilanv1obs,            & 
     417         &                 inlav1obs,             ld_nea,               & 
     418         &                 ibdyv1obs,             ld_bound_reject       ) 
    399419 
    400420      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
    401421      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 
    402422      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 
     423      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 
    403424 
    404425      ! Variable 2 
    … …  
    414435         &                 gdept_1d,              zmask2,               & 
    415436         &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    416          &                 iosdv2obs,              ilanv2obs,           & 
    417          &                 inlav2obs,              ld_nea                ) 
     437         &                 iosdv2obs,             ilanv2obs,            & 
     438         &                 inlav2obs,             ld_nea,               & 
     439         &                 ibdyv2obs,             ld_bound_reject       ) 
    418440 
    419441      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    420442      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    421443      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
     444      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
    422445 
    423446      ! ----------------------------------------------------------------------- 
    … …  
    489512               &            iuvchku 
    490513         ENDIF 
     514         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 
     515               &            ibdyv1obsmpp 
    491516         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    492517            &            prodatqc%nvprotmpp(1) 
    … …  
    506531               &            iuvchkv 
    507532         ENDIF 
     533         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
     534               &            ibdyv2obsmpp 
    508535         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    509536            &            prodatqc%nvprotmpp(2) 
    … …  
    875902      &                       plam,   pphi,    pmask,            & 
    876903      &                       kobsqc, kosdobs, klanobs,          & 
    877       &                       knlaobs,ld_nea                     ) 
     904      &                       knlaobs,ld_nea,                    & 
     905      &                       kbdyobs,ld_bound_reject            ) 
    878906      !!---------------------------------------------------------------------- 
    879907      !!                    ***  ROUTINE obs_coo_spc_2d  *** 
    … …  
    908936      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 
    909937         & kobsqc             ! Observation quality control 
    910       INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain 
    911       INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell 
    912       INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land 
    913       LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land 
     938      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain 
     939      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell 
     940      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land 
     941      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary 
     942      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
     943      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
    914944      !! * Local declarations 
    915945      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    916946         & zgmsk              ! Grid mask 
     947#if defined key_bdy  
     948      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
     949         & zbmsk              ! Boundary mask 
     950      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
     951#endif  
    917952      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    918953         & zglam, &           ! Model longitude at grid points 
    … …  
    956991 
    957992      END DO 
     993 
     994#if defined key_bdy              
     995      ! Create a mask grid points in boundary rim 
     996      IF (ld_bound_reject) THEN 
     997         zbdymask(:,:) = 1.0_wp 
     998         DO ji = 1, nb_bdy 
     999            DO jj = 1, idx_bdy(ji)%nblen(1) 
     1000               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1001            ENDDO 
     1002         ENDDO 
     1003  
     1004         CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, zbdymask, zbmsk )        
     1005      ENDIF 
     1006#endif        
    9581007       
    9591008      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 
    … …  
    10001049            END DO 
    10011050         END DO 
    1002    
    1003          ! For observations on the grid reject them if their are at 
    1004          ! a masked point 
    1005           
     1051  
    10061052         IF (lgridobs) THEN 
    10071053            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
    … …  
    10111057            ENDIF 
    10121058         ENDIF 
    1013                        
     1059 
     1060  
    10141061         ! Flag if the observation falls is close to land 
    10151062         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 
    1016             IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14 
    10171063            knlaobs = knlaobs + 1 
    1018             CYCLE 
    1019          ENDIF 
     1064            IF (ld_nea) THEN 
     1065               kobsqc(jobs) = kobsqc(jobs) + 14 
     1066               CYCLE 
     1067            ENDIF 
     1068         ENDIF 
     1069 
     1070#if defined key_bdy 
     1071         ! Flag if the observation falls close to the boundary rim 
     1072         IF (ld_bound_reject) THEN 
     1073            IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1074               kobsqc(jobs) = kobsqc(jobs) + 15 
     1075               kbdyobs = kbdyobs + 1 
     1076               CYCLE 
     1077            ENDIF 
     1078            ! for observations on the grid... 
     1079            IF (lgridobs) THEN 
     1080               IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1081                  kobsqc(jobs) = kobsqc(jobs) + 15 
     1082                  kbdyobs = kbdyobs + 1 
     1083                  CYCLE 
     1084               ENDIF 
     1085            ENDIF 
     1086         ENDIF 
     1087#endif  
    10201088             
    10211089      END DO 
    … …  
    10291097      &                       plam,    pphi,    pdep,    pmask, & 
    10301098      &                       kpobsqc, kobsqc,  kosdobs,        & 
    1031       &                       klanobs, knlaobs, ld_nea          ) 
     1099      &                       klanobs, knlaobs, ld_nea,         & 
     1100      &                       kbdyobs, ld_bound_reject          ) 
    10321101      !!---------------------------------------------------------------------- 
    10331102      !!                    ***  ROUTINE obs_coo_spc_3d  *** 
    … …  
    10521121      !! * Modules used 
    10531122      USE dom_oce, ONLY : &       ! Geographical information 
    1054          & gdepw_1d                         
     1123         & gdepw_1d,      & 
     1124         & gdepw_0,       &                        
     1125#if defined key_vvl 
     1126         & gdepw_n,       &  
     1127         & gdept_n,       & 
     1128#endif 
     1129         & ln_zco,        & 
     1130         & ln_zps,        & 
     1131         & lk_vvl                         
    10551132 
    10561133      !! * Arguments 
    … …  
    10861163      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell 
    10871164      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land 
     1165      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary 
    10881166      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land 
     1167      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary 
    10891168      !! * Local declarations 
    10901169      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10911170         & zgmsk              ! Grid mask 
     1171#if defined key_bdy  
     1172      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
     1173         & zbmsk              ! Boundary mask 
     1174      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
     1175#endif  
     1176      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
     1177         & zgdepw          
    10921178      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
    10931179         & zglam, &           ! Model longitude at grid points 
    … …  
    10971183         & igrdj 
    10981184      LOGICAL :: lgridobs           ! Is observation on a model grid point. 
     1185      LOGICAL :: ll_next_to_land    ! Is a profile next to land  
    10991186      INTEGER :: iig, ijg           ! i,j of observation on model grid point. 
    11001187      INTEGER :: jobs, jobsp, jk, ji, jj 
    … …  
    11311218          
    11321219      END DO 
     1220 
     1221#if defined key_bdy  
     1222      ! Create a mask grid points in boundary rim 
     1223      IF (ld_bound_reject) THEN            
     1224         zbdymask(:,:) = 1.0_wp 
     1225         DO ji = 1, nb_bdy 
     1226            DO jj = 1, idx_bdy(ji)%nblen(1) 
     1227               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1228            ENDDO 
     1229         ENDDO 
     1230      ENDIF 
     1231  
     1232      CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, zbdymask, zbmsk ) 
     1233#endif  
    11331234       
    11341235      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 
    … …  
    11591260         END DO 
    11601261 
     1262         ! Check if next to land  
     1263         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN  
     1264            ll_next_to_land=.TRUE.  
     1265         ELSE  
     1266            ll_next_to_land=.FALSE.  
     1267         ENDIF  
     1268          
    11611269         ! Reject observations 
    11621270 
    … …  
    11751283            ENDIF 
    11761284 
    1177             ! Flag if the observation falls with a model land cell 
    1178             IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
    1179                &  == 0.0_wp ) THEN 
    1180                kobsqc(jobsp) = kobsqc(jobsp) + 12 
    1181                klanobs = klanobs + 1 
    1182                CYCLE 
     1285            ! To check if an observations falls within land there are two cases:  
     1286            ! 1: z-coordibnates, where the check uses the mask  
     1287            ! 2: terrain following (eg s-coordinates),   
     1288            !    where we use the depth of the bottom cell to mask observations  
     1289              
     1290            IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco )  ) THEN !(CASE 1)  
     1291                 
     1292               ! Flag if the observation falls with a model land cell  
     1293               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &  
     1294                  &  == 0.0_wp ) THEN  
     1295                  kobsqc(jobsp) = kobsqc(jobsp) + 12  
     1296                  klanobs = klanobs + 1  
     1297                  CYCLE  
     1298               ENDIF  
     1299              
     1300               ! Flag if the observation is close to land  
     1301               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &  
     1302                  &  0.0_wp) THEN  
     1303                  knlaobs = knlaobs + 1  
     1304                  IF (ld_nea) THEN    
     1305                     kobsqc(jobsp) = kobsqc(jobsp) + 14   
     1306                  ENDIF   
     1307               ENDIF  
     1308              
     1309            ELSE ! Case 2  
     1310               ! Flag if the observation is deeper than the bathymetry  
     1311               ! Or if it is within the mask  
     1312               IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 
     1313                  &     .OR. &  
     1314                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
     1315                  &  == 0.0_wp) ) THEN 
     1316                  kobsqc(jobsp) = kobsqc(jobsp) + 12  
     1317                  klanobs = klanobs + 1  
     1318                  CYCLE  
     1319               ENDIF  
     1320                 
     1321               ! Flag if the observation is close to land  
     1322               IF ( ll_next_to_land ) THEN  
     1323                  knlaobs = knlaobs + 1  
     1324                  IF (ld_nea) THEN    
     1325                     kobsqc(jobsp) = kobsqc(jobsp) + 14   
     1326                  ENDIF   
     1327               ENDIF  
     1328              
    11831329            ENDIF 
    11841330 
    … …  
    11941340            ENDIF 
    11951341             
    1196             ! Flag if the observation falls is close to land 
    1197             IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
    1198                &  0.0_wp) THEN 
    1199                IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14 
    1200                knlaobs = knlaobs + 1 
    1201             ENDIF 
    1202  
    12031342            ! Set observation depth equal to that of the first model depth 
    12041343            IF ( pobsdep(jobsp) <= pdep(1) ) THEN 
    12051344               pobsdep(jobsp) = pdep(1)   
    12061345            ENDIF 
     1346             
     1347#if defined key_bdy 
     1348            ! Flag if the observation falls close to the boundary rim 
     1349            IF (ld_bound_reject) THEN 
     1350               IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1351                  kobsqc(jobsp) = kobsqc(jobsp) + 15 
     1352                  kbdyobs = kbdyobs + 1 
     1353                  CYCLE 
     1354               ENDIF 
     1355               ! for observations on the grid... 
     1356               IF (lgridobs) THEN 
     1357                  IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1358                     kobsqc(jobsp) = kobsqc(jobsp) + 15 
     1359                     kbdyobs = kbdyobs + 1 
     1360                     CYCLE 
     1361                  ENDIF 
     1362               ENDIF 
     1363            ENDIF 
     1364#endif  
    12071365             
    12081366         END DO 
Note: See TracChangeset for help on using the changeset viewer.