Ignore:
Timestamp:
2017-04-18T10:24:44+02:00 (3 years ago)
Author:
jwhile
Message:

Added code for 'Do not assimilate flags'

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/obs_oper_do_not_assim_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r7837 r7915  
    5252CONTAINS 
    5353 
    54    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject ) 
     54   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
     55                            kqc_cutoff ) 
    5556      !!---------------------------------------------------------------------- 
    5657      !!                    ***  ROUTINE obs_pre_sla  *** 
     
    8283      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land 
    8384      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
     85      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    8486      !! * Local declarations 
     87      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    8588      INTEGER :: iyea0        ! Initial date 
    8689      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    130133      ibdysobs = 0  
    131134 
     135      ! Set QC cutoff to optional value if provided 
     136      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
     137 
    132138      ! ----------------------------------------------------------------------- 
    133139      ! Find time coordinate for surface data 
     
    138144         &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
    139145         &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
    140          &              surfdata%nqc,     surfdata%mstp, iotdobs        ) 
     146         &              surfdata%nqc,     surfdata%mstp, iotdobs,       & 
     147         &              kqc_cutoff = iqc_cutoff  ) 
    141148 
    142149      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
     
    179186      ALLOCATE( llvalid(surfdata%nsurf) ) 
    180187       
    181       ! We want all data which has qc flags <= 10 
    182  
    183       llvalid(:)  = ( surfdata%nqc(:)  <= 10 ) 
     188      ! We want all data which has qc flags <= iqc_cutoff 
     189 
     190      llvalid(:)  = ( surfdata%nqc(:)  <= iqc_cutoff ) 
    184191 
    185192      ! The actual copying 
     
    251258      &                     kpi, kpj, kpk, & 
    252259      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
    253       &                     ld_nea, ld_bound_reject, kdailyavtypes ) 
     260      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff ) 
    254261 
    255262!!---------------------------------------------------------------------- 
     
    292299         & pgphi1, & 
    293300         & pgphi2 
     301      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    294302 
    295303      !! * Local declarations 
     304      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    296305      INTEGER :: iyea0        ! Initial date 
    297306      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    361370      iuvchkv   = 0 
    362371 
     372 
     373      ! Set QC cutoff to optional value if provided 
     374      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
     375 
    363376      ! ----------------------------------------------------------------------- 
    364377      ! Find time coordinate for profiles 
     
    371384            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    372385            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    373             &              iotdobs, kdailyavtypes = kdailyavtypes ) 
     386            &              iotdobs, kdailyavtypes = kdailyavtypes,         & 
     387            &              kqc_cutoff = iqc_cutoff ) 
    374388      ELSE 
    375389         CALL obs_coo_tim_prof( icycle, & 
     
    378392            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    379393            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    380             &              iotdobs ) 
     394            &              iotdobs,          kqc_cutoff = iqc_cutoff ) 
    381395      ENDIF 
    382396 
     
    395409 
    396410      ! ----------------------------------------------------------------------- 
    397       ! Reject all observations for profiles with nqc > 10 
    398       ! ----------------------------------------------------------------------- 
    399  
    400       CALL obs_pro_rej( profdata ) 
     411      ! Reject all observations for profiles with nqc > iqc_cutoff 
     412      ! ----------------------------------------------------------------------- 
     413 
     414      CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 
    401415 
    402416      ! ----------------------------------------------------------------------- 
     
    419433         &                 iosdv1obs,             ilanv1obs,            & 
    420434         &                 inlav1obs,             ld_nea,               & 
    421          &                 ibdyv1obs,             ld_bound_reject       ) 
     435         &                 ibdyv1obs,             ld_bound_reject,      & 
     436         &                 iqc_cutoff       ) 
    422437 
    423438      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
     
    440455         &                 iosdv2obs,             ilanv2obs,            & 
    441456         &                 inlav2obs,             ld_nea,               & 
    442          &                 ibdyv2obs,             ld_bound_reject       ) 
     457         &                 ibdyv2obs,             ld_bound_reject,      & 
     458         &                 iqc_cutoff       ) 
    443459 
    444460      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
     
    452468 
    453469      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    454          CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 
     470         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
    455471         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 
    456472         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 
     
    469485      END DO 
    470486 
    471       ! We want all data which has qc flags = 0 
    472  
    473       llvalid%luse(:) = ( profdata%nqc(:)  <= 10 ) 
     487      ! We want all data which has qc flags <= iqc_cutoff 
     488 
     489      llvalid%luse(:) = ( profdata%nqc(:)  <= iqc_cutoff ) 
    474490      DO jvar = 1,profdata%nvar 
    475          llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 
     491         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 
    476492      END DO 
    477493 
     
    783799      &                    kobsno,                                        & 
    784800      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   & 
    785       &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes ) 
     801      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, & 
     802      &                    kqc_cutoff ) 
    786803      !!---------------------------------------------------------------------- 
    787804      !!                    ***  ROUTINE obs_coo_tim *** 
     
    827844      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    828845         & kdailyavtypes    ! Types for daily averages 
     846      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     847 
    829848      !! * Local declarations 
    830849      INTEGER :: jobs 
     850      INTEGER :: iqc_cutoff=255 
    831851 
    832852      !----------------------------------------------------------------------- 
     
    847867         DO jobs = 1, kobsno 
    848868             
    849             IF ( kobsqc(jobs) <= 10 ) THEN 
     869            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 
    850870                
    851871               IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 
    852872                  & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 
    853                   kobsqc(jobs) = kobsqc(jobs) + 14 
     873                  kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    854874                  kotdobs      = kotdobs + 1 
    855875                  CYCLE 
     
    894914      DO jobs = 1, kobsno 
    895915         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 
    896             kobsqc(jobs) = kobsqc(jobs) + 18 
     916            kobsqc(jobs) = IBSET(kobsqc(jobs),12) 
    897917            kgrdobs = kgrdobs + 1 
    898918         ENDIF 
     
    906926      &                       kobsqc, kosdobs, klanobs,          & 
    907927      &                       knlaobs,ld_nea,                    & 
    908       &                       kbdyobs,ld_bound_reject            ) 
     928      &                       kbdyobs,ld_bound_reject,           & 
     929      &                       kqc_cutoff                         ) 
    909930      !!---------------------------------------------------------------------- 
    910931      !!                    ***  ROUTINE obs_coo_spc_2d  *** 
     
    945966      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
    946967      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
     968      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value 
     969 
    947970      !! * Local declarations 
    948971      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
     
    969992         ! For invalid points use 2,2 
    970993 
    971          IF ( kobsqc(jobs) >= 10 ) THEN 
     994         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 
    972995 
    973996            igrdi(1,1,jobs) = 1 
     
    10161039 
    10171040         ! Skip bad observations 
    1018          IF ( kobsqc(jobs) >= 10 ) CYCLE 
     1041         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 
    10191042 
    10201043         ! Flag if the observation falls outside the model spatial domain 
     
    10231046            &  .OR. ( pobsphi(jobs) <  -90. ) & 
    10241047            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN 
    1025             kobsqc(jobs) = kobsqc(jobs) + 11 
     1048            kobsqc(jobs) = IBSET(kobsqc(jobs),11) 
    10261049            kosdobs = kosdobs + 1 
    10271050            CYCLE 
     
    10301053         ! Flag if the observation falls with a model land cell 
    10311054         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
    1032             kobsqc(jobs) = kobsqc(jobs)  + 12 
     1055            kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    10331056            klanobs = klanobs + 1 
    10341057            CYCLE 
     
    10551078         IF (lgridobs) THEN 
    10561079            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
    1057                kobsqc(jobs) = kobsqc(jobs) + 12 
     1080               kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    10581081               klanobs = klanobs + 1 
    10591082               CYCLE 
     
    10661089            knlaobs = knlaobs + 1 
    10671090            IF (ld_nea) THEN 
    1068                kobsqc(jobs) = kobsqc(jobs) + 14 
     1091               kobsqc(jobs) = IBSET(kobsqc(jobs),9) 
    10691092               CYCLE 
    10701093            ENDIF 
     
    11011124      &                       kpobsqc, kobsqc,  kosdobs,        & 
    11021125      &                       klanobs, knlaobs, ld_nea,         & 
    1103       &                       kbdyobs, ld_bound_reject          ) 
     1126      &                       kbdyobs, ld_bound_reject,         & 
     1127      &                       kqc_cutoff                        ) 
    11041128      !!---------------------------------------------------------------------- 
    11051129      !!                    ***  ROUTINE obs_coo_spc_3d  *** 
     
    11691193      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land 
    11701194      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary 
     1195      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value 
     1196 
    11711197      !! * Local declarations 
    11721198      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
     
    11961222         ! For invalid points use 2,2 
    11971223 
    1198          IF ( kpobsqc(jobs) >= 10 ) THEN 
     1224         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 
    11991225             
    12001226            igrdi(1,1,jobs) = 1 
     
    12461272 
    12471273         ! Skip bad profiles 
    1248          IF ( kpobsqc(jobs) >= 10 ) CYCLE 
     1274         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 
    12491275 
    12501276         ! Check if this observation is on a grid point 
     
    12841310               &  .OR. ( pobsdep(jobsp) < 0.0          )       & 
    12851311               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 
    1286                kobsqc(jobsp) = kobsqc(jobsp) + 11 
     1312               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 
    12871313               kosdobs = kosdobs + 1 
    12881314               CYCLE 
     
    12991325               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &  
    13001326                  &  == 0.0_wp ) THEN  
    1301                   kobsqc(jobsp) = kobsqc(jobsp) + 12  
     1327                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
    13021328                  klanobs = klanobs + 1  
    13031329                  CYCLE  
     
    13091335                  knlaobs = knlaobs + 1  
    13101336                  IF (ld_nea) THEN    
    1311                      kobsqc(jobsp) = kobsqc(jobsp) + 14   
     1337                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
    13121338                  ENDIF   
    13131339               ENDIF  
     
    13201346                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
    13211347                  &  == 0.0_wp) ) THEN 
    1322                   kobsqc(jobsp) = kobsqc(jobsp) + 12  
     1348                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
    13231349                  klanobs = klanobs + 1  
    13241350                  CYCLE  
     
    13291355                  knlaobs = knlaobs + 1  
    13301356                  IF (ld_nea) THEN    
    1331                      kobsqc(jobsp) = kobsqc(jobsp) + 14   
     1357                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
    13321358                  ENDIF   
    13331359               ENDIF  
     
    13751401   END SUBROUTINE obs_coo_spc_3d 
    13761402 
    1377    SUBROUTINE obs_pro_rej( profdata ) 
     1403   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 
    13781404      !!---------------------------------------------------------------------- 
    13791405      !!                    ***  ROUTINE obs_pro_rej *** 
     
    13931419      !! * Arguments 
    13941420      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data 
     1421      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value 
     1422 
    13951423      !! * Local declarations 
    13961424      INTEGER :: jprof 
     
    14021430      DO jprof = 1, profdata%nprof 
    14031431 
    1404          IF ( profdata%nqc(jprof) > 10 ) THEN 
     1432         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 
    14051433             
    14061434            DO jvar = 1, profdata%nvar 
     
    14101438                   
    14111439                  profdata%var(jvar)%nvqc(jobs) = & 
    1412                      & profdata%var(jvar)%nvqc(jobs) + 26 
     1440                     & IBSET(profdata%var(jvar)%nvqc(jobs),14) 
    14131441 
    14141442               END DO 
     
    14221450   END SUBROUTINE obs_pro_rej 
    14231451 
    1424    SUBROUTINE obs_uv_rej( profdata, knumu, knumv ) 
     1452   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutof ) 
    14251453      !!---------------------------------------------------------------------- 
    14261454      !!                    ***  ROUTINE obs_uv_rej *** 
     
    14421470      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected 
    14431471      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected 
     1472      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     1473 
    14441474      !! * Local declarations 
    14451475      INTEGER :: jprof 
     
    14611491         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 
    14621492             
    1463             IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. & 
    1464                & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN 
    1465                profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42 
     1493            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1494               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1495               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    14661496               knumv = knumv + 1 
    14671497            ENDIF 
    1468             IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. & 
    1469                & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN 
    1470                profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42 
     1498            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1499               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1500               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    14711501               knumu = knumu + 1 
    14721502            ENDIF 
Note: See TracChangeset for help on using the changeset viewer.