Changeset 8667


Ignore:
Timestamp:
2017-10-30T10:28:45+01:00 (3 years ago)
Author:
timgraham
Message:

Update of OBS code from local v3.6 branch to head of trunk

Location:
branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r6140 r8667  
    155155 
    156156      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    157          &            ln_sst, ln_sic, ln_vel3d,                       & 
    158          &            ln_altbias, ln_nea, ln_grid_global,             & 
    159          &            ln_grid_search_lookup,                          & 
    160          &            ln_ignmis, ln_s_at_t, ln_sstnight,              & 
     157         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     158         &            ln_altbias, ln_sstbias, ln_nea,                 & 
     159         &            ln_grid_global, ln_grid_search_lookup,          & 
     160         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
     161         &            ln_sstnight,                                    & 
     162         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
     163         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
    161164         &            cn_profbfiles, cn_slafbfiles,                   & 
    162165         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    163          &            cn_velfbfiles, cn_altbiasfile,                  & 
     166         &            cn_velfbfiles, cn_sssfbfiles,                   & 
     167         &            cn_sstbiasfiles, cn_altbiasfile,                & 
    164168         &            cn_gridsearchfile, rn_gridsearchres,            & 
    165          &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     & 
     169         &            rn_dobsini, rn_dobsend,                         & 
     170         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
     171         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
     172         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     173         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     174         &            nn_1dint, nn_2dint,                             & 
     175         &            nn_2dint_sla, nn_2dint_sst,                     & 
     176         &            nn_2dint_sss, nn_2dint_sic,                     & 
    166177         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    167          &            nn_profdavtypes, ln_sstbias, cn_sstbias_files 
     178         &            nn_profdavtypes 
    168179 
    169180      INTEGER :: jnumsstbias 
     
    187198      cn_sicfbfiles(:) = '' 
    188199      cn_velfbfiles(:) = '' 
     200      cn_sssfbfiles(:)    = '' 
    189201      cn_sstbias_files(:) = '' 
    190202      nn_profdavtypes(:) = -1 
     
    208220         RETURN 
    209221      ENDIF 
    210        
    211       !----------------------------------------------------------------------- 
    212       ! Set up list of observation types to be used 
    213       ! and the files associated with each type 
    214       !----------------------------------------------------------------------- 
    215  
    216       nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
    217       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 
    218  
    219       IF (ln_sstbias) THEN  
    220          lmask(:) = .FALSE.  
    221          WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE.  
    222          jnumsstbias = COUNT(lmask)  
    223          lmask(:) = .FALSE.  
    224       ENDIF       
    225  
    226       IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
    227          IF(lwp) WRITE(numout,cform_war) 
    228          IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
    229             &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
    230             &                    ' are set to .FALSE. so turning off calls to dia_obs' 
    231          nwarn = nwarn + 1 
    232          ln_diaobs = .FALSE. 
    233          RETURN 
    234       ENDIF 
    235  
    236       IF ( nproftypes > 0 ) THEN 
    237  
    238          ALLOCATE( cobstypesprof(nproftypes) ) 
    239          ALLOCATE( ifilesprof(nproftypes) ) 
    240          ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
    241  
    242          jtype = 0 
    243          IF (ln_t3d .OR. ln_s3d) THEN 
    244             jtype = jtype + 1 
    245             clproffiles(jtype,:) = cn_profbfiles(:) 
    246             cobstypesprof(jtype) = 'prof  ' 
    247             ifilesprof(jtype) = 0 
    248             DO jfile = 1, jpmaxnfiles 
    249                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    250                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    251             END DO 
    252          ENDIF 
    253          IF (ln_vel3d) THEN 
    254             jtype = jtype + 1 
    255             clproffiles(jtype,:) = cn_velfbfiles(:) 
    256             cobstypesprof(jtype) = 'vel   ' 
    257             ifilesprof(jtype) = 0 
    258             DO jfile = 1, jpmaxnfiles 
    259                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    260                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    261             END DO 
    262          ENDIF 
    263  
    264       ENDIF 
    265  
    266       IF ( nsurftypes > 0 ) THEN 
    267  
    268          ALLOCATE( cobstypessurf(nsurftypes) ) 
    269          ALLOCATE( ifilessurf(nsurftypes) ) 
    270          ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
    271  
    272          jtype = 0 
    273          IF (ln_sla) THEN 
    274             jtype = jtype + 1 
    275             clsurffiles(jtype,:) = cn_slafbfiles(:) 
    276             cobstypessurf(jtype) = 'sla   ' 
    277             ifilessurf(jtype) = 0 
    278             DO jfile = 1, jpmaxnfiles 
    279                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    280                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    281             END DO 
    282          ENDIF 
    283          IF (ln_sst) THEN 
    284             jtype = jtype + 1 
    285             clsurffiles(jtype,:) = cn_sstfbfiles(:) 
    286             cobstypessurf(jtype) = 'sst   ' 
    287             ifilessurf(jtype) = 0 
    288             DO jfile = 1, jpmaxnfiles 
    289                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    290                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    291             END DO 
    292          ENDIF 
    293 #if defined key_lim2 || defined key_lim3 
    294          IF (ln_sic) THEN 
    295             jtype = jtype + 1 
    296             clsurffiles(jtype,:) = cn_sicfbfiles(:) 
    297             cobstypessurf(jtype) = 'sic   ' 
    298             ifilessurf(jtype) = 0 
    299             DO jfile = 1, jpmaxnfiles 
    300                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    301                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    302             END DO 
    303          ENDIF 
    304 #endif 
    305  
    306       ENDIF 
    307  
    308       !Write namelist settings to stdout 
     222 
    309223      IF(lwp) THEN 
    310224         WRITE(numout,*) 
     
    318232         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
    319233         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    320          WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global 
    321          WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias  
    322          WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 
     234         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     235         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     236         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    323237         IF (ln_grid_search_lookup) & 
    324238            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
     
    328242         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
    329243         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     244         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
    330245         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
    331246         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
    332247         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
    333248         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     249         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
    334250         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
    335251         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    336252         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
    337          WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
    338  
    339          IF ( nproftypes > 0 ) THEN 
    340             DO jtype = 1, nproftypes 
    341                DO jfile = 1, ifilesprof(jtype) 
    342                   WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', & 
    343                      TRIM(clproffiles(jtype,jfile)) 
    344                END DO 
    345             END DO 
     253      ENDIF 
     254      !----------------------------------------------------------------------- 
     255      ! Set up list of observation types to be used 
     256      ! and the files associated with each type 
     257      !----------------------------------------------------------------------- 
     258 
     259      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     260      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) ) 
     261 
     262      IF (ln_sstbias) THEN  
     263         lmask(:) = .FALSE.  
     264         WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE.  
     265         jnumsstbias = COUNT(lmask)  
     266         lmask(:) = .FALSE.  
     267      ENDIF       
     268 
     269      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     270         IF(lwp) WRITE(numout,cform_war) 
     271         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     272            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
     273            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     274         nwarn = nwarn + 1 
     275         ln_diaobs = .FALSE. 
     276         RETURN 
     277      ENDIF 
     278 
     279      IF ( nproftypes > 0 ) THEN 
     280 
     281         ALLOCATE( cobstypesprof(nproftypes) ) 
     282         ALLOCATE( ifilesprof(nproftypes) ) 
     283         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     284 
     285         jtype = 0 
     286         IF (ln_t3d .OR. ln_s3d) THEN 
     287            jtype = jtype + 1 
     288            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
     289               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    346290         ENDIF 
    347  
    348          WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
    349          IF ( nsurftypes > 0 ) THEN 
    350             DO jtype = 1, nsurftypes 
    351                DO jfile = 1, ifilessurf(jtype) 
    352                   WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', & 
    353                      TRIM(clsurffiles(jtype,jfile)) 
    354                END DO 
    355             END DO 
     291         IF (ln_vel3d) THEN 
     292            jtype = jtype + 1 
     293            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
     294               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    356295         ENDIF 
    357          WRITE(numout,*) '~~~~~~~~~~~~' 
    358  
    359       ENDIF 
     296 
     297      ENDIF 
     298 
     299      IF ( nsurftypes > 0 ) THEN 
     300 
     301         ALLOCATE( cobstypessurf(nsurftypes) ) 
     302         ALLOCATE( ifilessurf(nsurftypes) ) 
     303         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     304         ALLOCATE(n2dintsurf(nsurftypes)) 
     305         ALLOCATE(ravglamscl(nsurftypes)) 
     306         ALLOCATE(ravgphiscl(nsurftypes)) 
     307         ALLOCATE(lfpindegs(nsurftypes)) 
     308         ALLOCATE(llnightav(nsurftypes)) 
     309 
     310         jtype = 0 
     311         IF (ln_sla) THEN 
     312            jtype = jtype + 1 
     313            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', & 
     314               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     315            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      & 
     316               &                  nn_2dint, nn_2dint_sla,             & 
     317               &                  rn_sla_avglamscl, rn_sla_avgphiscl, & 
     318               &                  ln_sla_fp_indegs, .FALSE.,          & 
     319               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     320               &                  lfpindegs, llnightav ) 
     321         ENDIF 
     322         IF (ln_sst) THEN 
     323            jtype = jtype + 1 
     324            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', & 
     325               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     326            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      & 
     327               &                  nn_2dint, nn_2dint_sst,             & 
     328               &                  rn_sst_avglamscl, rn_sst_avgphiscl, & 
     329               &                  ln_sst_fp_indegs, ln_sstnight,      & 
     330               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     331               &                  lfpindegs, llnightav ) 
     332         ENDIF 
     333#if defined key_lim2 || defined key_lim3 || defined key_cice 
     334         IF (ln_sic) THEN 
     335            jtype = jtype + 1 
     336            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', & 
     337               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     338            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      & 
     339               &                  nn_2dint, nn_2dint_sic,             & 
     340               &                  rn_sic_avglamscl, rn_sic_avgphiscl, & 
     341               &                  ln_sic_fp_indegs, .FALSE.,          & 
     342               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     343               &                  lfpindegs, llnightav ) 
     344         ENDIF 
     345#endif 
     346         IF (ln_sss) THEN 
     347            jtype = jtype + 1 
     348            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', & 
     349               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     350            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      & 
     351               &                  nn_2dint, nn_2dint_sss,             & 
     352               &                  rn_sss_avglamscl, rn_sss_avgphiscl, & 
     353               &                  ln_sss_fp_indegs, .FALSE.,          & 
     354               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     355               &                  lfpindegs, llnightav ) 
     356         ENDIF 
     357 
     358      ENDIF 
     359 
     360 
    360361 
    361362      !----------------------------------------------------------------------- 
     
    377378      ENDIF 
    378379 
    379       IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 
     380      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 
    380381         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
    381382            &                    ' is not available') 
     
    442443               &               jpi, jpj, jpk, & 
    443444               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  & 
    444                &               ln_nea, kdailyavtypes = nn_profdavtypes ) 
     445               &               ln_nea, ln_bound_reject, & 
     446               &               kdailyavtypes = nn_profdavtypes ) 
    445447 
    446448         END DO 
     
    469471               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    470472               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    471                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 
    472           
    473           
    474             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 
     473               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     474 
     475            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
    475476 
    476477            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    477                CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 
    478                IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 
     478               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
     479               IF ( ln_altbias ) & 
     480                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
    479481            ENDIF 
     482 
    480483            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
    481                !Read in bias field and correct SST. 
    482                IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
    483                                                      "  but no bias"// & 
    484                                                      " files to read in")    
    485                   CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 
    486                                         jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 
     484               jnumsstbias = 0 
     485               DO jfile = 1, jpmaxnfiles 
     486                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
     487                     &  jnumsstbias = jnumsstbias + 1 
     488               END DO 
     489               IF ( jnumsstbias == 0 ) THEN 
     490                  CALL ctl_stop("ln_sstbias set but no bias files to read in")     
     491               ENDIF 
     492 
     493               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
     494                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
     495 
    487496            ENDIF 
    488497         END DO 
     
    545554         & frld 
    546555#endif 
     556#if defined key_cice 
     557      USE sbc_oce, ONLY : fr_i     ! ice fraction 
     558#endif 
     559 
    547560      IMPLICIT NONE 
    548561 
     
    561574         & zprofmask2              ! Mask associated with zprofvar2 
    562575      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    563          & zsurfvar                ! Model values equivalent to surface ob. 
     576         & zsurfvar, &             ! Model values equivalent to surface ob. 
     577         & zsurfmask               ! Mask associated with surface variable 
    564578      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    565579         & zglam1,    &            ! Model longitudes for prof variable 1 
     
    567581         & zgphi1,    &            ! Model latitudes for prof variable 1 
    568582         & zgphi2                  ! Model latitudes for prof variable 2 
    569 #if ! defined key_lim2 && ! defined key_lim3 
    570       REAL(wp), POINTER, DIMENSION(:,:) :: frld 
    571 #endif 
    572583      LOGICAL :: llnightav        ! Logical for calculating night-time average 
    573584 
     
    578589      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 
    579590      CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     591      CALL wrk_alloc( jpi, jpj, zsurfmask ) 
    580592      CALL wrk_alloc( jpi, jpj, zglam1 ) 
    581593      CALL wrk_alloc( jpi, jpj, zglam2 ) 
    582594      CALL wrk_alloc( jpi, jpj, zgphi1 ) 
    583595      CALL wrk_alloc( jpi, jpj, zgphi2 ) 
    584 #if ! defined key_lim2 && ! defined key_lim3 
    585       CALL wrk_alloc(jpi,jpj,frld)  
    586 #endif 
    587596 
    588597      IF(lwp) THEN 
     
    594603      idaystp = NINT( rday / rdt ) 
    595604 
    596       !----------------------------------------------------------------------- 
    597       ! No LIM => frld == 0.0_wp 
    598       !----------------------------------------------------------------------- 
    599 #if ! defined key_lim2 && ! defined key_lim3 
    600       frld(:,:) = 0.0_wp 
    601 #endif 
    602605      !----------------------------------------------------------------------- 
    603606      ! Call the profile and surface observation operators 
     
    627630               zgphi1(:,:) = gphiu(:,:) 
    628631               zgphi2(:,:) = gphiv(:,:) 
     632            CASE DEFAULT 
     633               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
    629634            END SELECT 
    630635 
    631             IF( ln_zco .OR. ln_zps ) THEN  
    632                CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
    633                   &               nit000, idaystp,                         & 
    634                   &               zprofvar1, zprofvar2,                    & 
    635                   &               gdept_1d, zprofmask1, zprofmask2,        & 
    636                   &               zglam1, zglam2, zgphi1, zgphi2,          & 
    637                   &               nn_1dint, nn_2dint,                      & 
    638                   &               kdailyavtypes = nn_profdavtypes ) 
    639             ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 
    640                !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 
    641                CALL obs_pro_sco_opt( profdataqc(jtype),                    &  
    642                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
    643                   &              zprofvar1, zprofvar2,                   &  
    644                   &              gdept_n(:,:,:), gdepw_n(:,:,:),           & 
    645                   &              tmask, nn_1dint, nn_2dint,              &  
    646                   &              kdailyavtypes = nn_profdavtypes )  
    647             ELSE 
    648                CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 
    649                          'yet working for velocity data (turn off velocity observations') 
    650             ENDIF 
     636            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     637               &               nit000, idaystp,                         & 
     638               &               zprofvar1, zprofvar2,                    & 
     639               &               gdept_n(:,:,:), gdepw_n(:,:,:),            &  
     640               &               zprofmask1, zprofmask2,                  & 
     641               &               zglam1, zglam2, zgphi1, zgphi2,          & 
     642               &               nn_1dint, nn_2dint,                      & 
     643               &               kdailyavtypes = nn_profdavtypes ) 
    651644 
    652645         END DO 
     
    657650 
    658651         DO jtype = 1, nsurftypes 
     652 
     653            !Defaults which might be changed 
     654            zsurfmask(:,:) = tmask(:,:,1) 
    659655 
    660656            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    661657            CASE('sst') 
    662658               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
    663                llnightav = ln_sstnight 
    664659            CASE('sla') 
    665660               zsurfvar(:,:) = sshn(:,:) 
    666                llnightav = .FALSE. 
    667 #if defined key_lim2 || defined key_lim3 
     661            CASE('sss') 
     662               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
    668663            CASE('sic') 
    669664               IF ( kstp == 0 ) THEN 
     
    678673                  CYCLE 
    679674               ELSE 
     675#if defined key_cice 
     676                  zsurfvar(:,:) = fr_i(:,:) 
     677#elif defined key_lim2 || defined key_lim3 
    680678                  zsurfvar(:,:) = 1._wp - frld(:,:) 
     679#else 
     680               CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     681                  &           ' but no sea-ice model appears to have been defined' ) 
     682#endif 
    681683               ENDIF 
    682684 
    683                llnightav = .FALSE. 
    684 #endif 
    685685            END SELECT 
    686686 
    687687            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    688                &               nit000, idaystp, zsurfvar, tmask(:,:,1), & 
    689                &               nn_2dint, llnightav ) 
     688               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
     689               &               n2dintsurf(jtype), llnightav(jtype),     & 
     690               &               ravglamscl(jtype), ravgphiscl(jtype),     & 
     691               &               lfpindegs(jtype) ) 
    690692 
    691693         END DO 
     
    698700      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 
    699701      CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     702      CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
    700703      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
    701704      CALL wrk_dealloc( jpi, jpj, zglam2 ) 
    702705      CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
    703706      CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
    704 #if ! defined key_lim2 && ! defined key_lim3 
    705       CALL wrk_dealloc(jpi,jpj,frld) 
    706 #endif 
    707707 
    708708   END SUBROUTINE dia_obs 
     
    821821 
    822822      IF ( nsurftypes > 0 ) & 
    823          &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 
     823         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 
     824         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 
    824825 
    825826   END SUBROUTINE dia_obs_dealloc 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90

    r6140 r8667  
    142142      !! 
    143143 
    144       iobsp=kobsp 
     144      iobsp(:)=kobsp(:) 
    145145 
    146146      WHERE( iobsp(:) == -1 ) 
     
    148148      END WHERE 
    149149 
    150       iobsp=-1*iobsp 
     150      iobsp(:)=-1*iobsp(:) 
    151151 
    152152      CALL obs_mpp_max_integer( iobsp, kno ) 
    153153 
    154       kobsp=-1*iobsp 
     154      kobsp(:)=-1*iobsp(:) 
    155155 
    156156      isum=0 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r7646 r8667  
    99   !!   obs_prof_opt :    Compute the model counterpart of profile data 
    1010   !!   obs_surf_opt :    Compute the model counterpart of surface data 
    11    !!   obs_pro_sco_opt: Compute the model counterpart of temperature and  
    12    !!                    salinity observations from profiles in generalised  
    13    !!                    vertical coordinates  
    1411   !!---------------------------------------------------------------------- 
    1512 
     
    2219      & obs_int_h2d, & 
    2320      & obs_int_h2d_init 
     21   USE obs_averg_h2d, ONLY : &    ! Horizontal averaging to the obs footprint 
     22      & obs_avg_h2d, & 
     23      & obs_avg_h2d_init, & 
     24      & obs_max_fpsize 
    2425   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt 
    2526      & obs_int_z1d,    & 
     
    2829      & obfillflt      ! Fillvalue    
    2930   USE dom_oce,       ONLY : & 
    30       & glamt, glamu, glamv, & 
    31       & gphit, gphiu, gphiv, &  
     31      & glamt, glamf, & 
     32      & gphit, gphif, &  
    3233      & gdept_n, gdept_0  
    3334   USE lib_mpp,       ONLY : & 
     
    4445 
    4546   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs 
    46       &   obs_pro_sco_opt, &  ! Compute the model counterpart of profile observations  
    4747      &   obs_surf_opt     ! Compute the model counterpart of surface obs 
    4848 
     
    290290      END DO 
    291291 
     292      ! Initialise depth arrays 
     293      zgdept(:,:,:,:) = 0.0 
     294      zgdepw(:,:,:,:) = 0.0 
     295 
    292296      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 
    293297      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 
     
    300304      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    301305 
     306      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
     307      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )  
     308 
    302309      ! At the end of the day also get interpolated means 
    303310      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     
    314321 
    315322      ENDIF 
     323 
     324      ! Return if no observations to process  
     325      ! Has to be done after comm commands to ensure processors  
     326      ! stay in sync  
     327      IF ( ipro == 0 ) RETURN  
    316328 
    317329      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
     
    339351         zphi = prodatqc%rphi(jobs) 
    340352 
    341          ! Horizontal weights and vertical mask 
    342  
     353         ! Horizontal weights  
     354         ! Masked values are calculated later.   
    343355         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    344356 
    345             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     357            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    346358               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
    347                &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
     359               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 
    348360 
    349361         ENDIF 
     
    351363         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    352364 
    353             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     365            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    354366               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    355                &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
     367               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
    356368  
    357369         ENDIF 
     
    365377               IF ( idayend == 0 )  THEN 
    366378                  ! Daily averaged data 
    367                   CALL obs_int_h2d( kpk, kpk,      & 
    368                      &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    369  
    370                ENDIF 
     379 
     380                  ! vertically interpolate all 4 corners  
     381                  ista = prodatqc%npvsta(jobs,1)  
     382                  iend = prodatqc%npvend(jobs,1)  
     383                  inum_obs = iend - ista + 1  
     384                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     385 
     386                  DO iin=1,2  
     387                     DO ijn=1,2  
     388 
     389                        IF ( k1dint == 1 ) THEN  
     390                           CALL obs_int_z1d_spl( kpk, &  
     391                              &     zinm1(iin,ijn,:,iobs), &  
     392                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     393                              &     zmask1(iin,ijn,:,iobs))  
     394                        ENDIF  
     395        
     396                        CALL obs_level_search(kpk, &  
     397                           &    zgdept(iin,ijn,:,iobs), &  
     398                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     399                           &    iv_indic)  
     400 
     401                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     402                           &    prodatqc%var(1)%vdep(ista:iend), &  
     403                           &    zinm1(iin,ijn,:,iobs), &  
     404                           &    zobs2k, interp_corner(iin,ijn,:), &  
     405                           &    zgdept(iin,ijn,:,iobs), &  
     406                           &    zmask1(iin,ijn,:,iobs))  
     407        
     408                     ENDDO  
     409                  ENDDO  
     410 
     411               ENDIF !idayend 
    371412 
    372413            ELSE  
    373414 
    374                ! Point data 
    375                CALL obs_int_h2d( kpk, kpk,      & 
    376                   &              zweig1, zint1(:,:,:,iobs), zobsk ) 
    377  
    378             ENDIF 
    379  
    380             !------------------------------------------------------------- 
    381             ! Compute vertical second-derivative of the interpolating  
    382             ! polynomial at obs points 
    383             !------------------------------------------------------------- 
    384  
    385             IF ( k1dint == 1 ) THEN 
    386                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    387                   &                  pgdept, zobsmask1 ) 
    388             ENDIF 
    389  
    390             !----------------------------------------------------------------- 
    391             !  Vertical interpolation to the observation point 
    392             !----------------------------------------------------------------- 
    393             ista = prodatqc%npvsta(jobs,1) 
    394             iend = prodatqc%npvend(jobs,1) 
    395             CALL obs_int_z1d( kpk,                & 
    396                & prodatqc%var(1)%mvk(ista:iend),  & 
    397                & k1dint, iend - ista + 1,         & 
    398                & prodatqc%var(1)%vdep(ista:iend), & 
    399                & zobsk, zobs2k,                   & 
    400                & prodatqc%var(1)%vmod(ista:iend), & 
    401                & pgdept, zobsmask1 ) 
    402  
    403          ENDIF 
    404  
     415               ! Point data  
     416      
     417               ! vertically interpolate all 4 corners  
     418               ista = prodatqc%npvsta(jobs,1)  
     419               iend = prodatqc%npvend(jobs,1)  
     420               inum_obs = iend - ista + 1  
     421               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     422               DO iin=1,2   
     423                  DO ijn=1,2  
     424                     
     425                     IF ( k1dint == 1 ) THEN  
     426                        CALL obs_int_z1d_spl( kpk, &  
     427                           &    zint1(iin,ijn,:,iobs),&  
     428                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     429                           &    zmask1(iin,ijn,:,iobs))  
     430   
     431                     ENDIF  
     432        
     433                     CALL obs_level_search(kpk, &  
     434                         &        zgdept(iin,ijn,:,iobs),&  
     435                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     436                         &        iv_indic)  
     437 
     438                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     439                         &          prodatqc%var(1)%vdep(ista:iend),     &  
     440                         &          zint1(iin,ijn,:,iobs),            &  
     441                         &          zobs2k,interp_corner(iin,ijn,:), &  
     442                         &          zgdept(iin,ijn,:,iobs),         &  
     443                         &          zmask1(iin,ijn,:,iobs) )       
     444          
     445                  ENDDO  
     446               ENDDO  
     447              
     448            ENDIF  
     449 
     450            !-------------------------------------------------------------  
     451            ! Compute the horizontal interpolation for every profile level  
     452            !-------------------------------------------------------------  
     453              
     454            DO ikn=1,inum_obs  
     455               iend=ista+ikn-1 
     456                   
     457               zweig(:,:,1) = 0._wp  
     458    
     459               ! This code forces the horizontal weights to be   
     460               ! zero IF the observation is below the bottom of the   
     461               ! corners of the interpolation nodes, Or if it is in   
     462               ! the mask. This is important for observations near   
     463               ! steep bathymetry  
     464               DO iin=1,2  
     465                  DO ijn=1,2  
     466      
     467                     depth_loop1: DO ik=kpk,2,-1  
     468                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     469                             
     470                           zweig(iin,ijn,1) = &   
     471                              & zweig1(iin,ijn,1) * &  
     472                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     473                              &  - prodatqc%var(1)%vdep(iend)),0._wp)  
     474                             
     475                           EXIT depth_loop1  
     476 
     477                        ENDIF  
     478                     ENDDO depth_loop1  
     479      
     480                  ENDDO  
     481               ENDDO  
     482    
     483               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     484                  &              prodatqc%var(1)%vmod(iend:iend) )  
     485 
     486                  ! Set QC flag for any observations found below the bottom 
     487                  ! needed as the check here is more strict than that in obs_prep 
     488               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 
     489  
     490            ENDDO  
     491  
     492            DEALLOCATE(interp_corner,iv_indic)  
     493           
     494         ENDIF  
     495 
     496         ! For the second variable 
    405497         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    406498 
     
    410502 
    411503               IF ( idayend == 0 )  THEN 
    412  
    413504                  ! Daily averaged data 
    414                   CALL obs_int_h2d( kpk, kpk,      & 
    415                      &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    416  
    417                ENDIF 
    418  
    419             ELSE 
    420  
    421                ! Point data 
    422                CALL obs_int_h2d( kpk, kpk,      & 
    423                   &              zweig2, zint2(:,:,:,iobs), zobsk ) 
    424  
    425             ENDIF 
    426  
    427  
    428             !------------------------------------------------------------- 
    429             ! Compute vertical second-derivative of the interpolating  
    430             ! polynomial at obs points 
    431             !------------------------------------------------------------- 
    432  
    433             IF ( k1dint == 1 ) THEN 
    434                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    435                   &                  pgdept, zobsmask2 ) 
    436             ENDIF 
    437  
    438             !---------------------------------------------------------------- 
    439             !  Vertical interpolation to the observation point 
    440             !---------------------------------------------------------------- 
    441             ista = prodatqc%npvsta(jobs,2) 
    442             iend = prodatqc%npvend(jobs,2) 
    443             CALL obs_int_z1d( kpk, & 
    444                & prodatqc%var(2)%mvk(ista:iend),& 
    445                & k1dint, iend - ista + 1, & 
    446                & prodatqc%var(2)%vdep(ista:iend),& 
    447                & zobsk, zobs2k, & 
    448                & prodatqc%var(2)%vmod(ista:iend),& 
    449                & pgdept, zobsmask2 ) 
    450  
    451          ENDIF 
    452  
    453       END DO 
     505 
     506                  ! vertically interpolate all 4 corners  
     507                  ista = prodatqc%npvsta(jobs,2)  
     508                  iend = prodatqc%npvend(jobs,2)  
     509                  inum_obs = iend - ista + 1  
     510                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     511 
     512                  DO iin=1,2  
     513                     DO ijn=1,2  
     514 
     515                        IF ( k1dint == 1 ) THEN  
     516                           CALL obs_int_z1d_spl( kpk, &  
     517                              &     zinm2(iin,ijn,:,iobs), &  
     518                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     519                              &     zmask2(iin,ijn,:,iobs))  
     520                        ENDIF  
     521        
     522                        CALL obs_level_search(kpk, &  
     523                           &    zgdept(iin,ijn,:,iobs), &  
     524                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     525                           &    iv_indic)  
     526 
     527                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     528                           &    prodatqc%var(2)%vdep(ista:iend), &  
     529                           &    zinm2(iin,ijn,:,iobs), &  
     530                           &    zobs2k, interp_corner(iin,ijn,:), &  
     531                           &    zgdept(iin,ijn,:,iobs), &  
     532                           &    zmask2(iin,ijn,:,iobs))  
     533        
     534                     ENDDO  
     535                  ENDDO  
     536 
     537               ENDIF !idayend 
     538 
     539            ELSE    
     540 
     541               ! Point data  
     542      
     543               ! vertically interpolate all 4 corners  
     544               ista = prodatqc%npvsta(jobs,2)  
     545               iend = prodatqc%npvend(jobs,2)  
     546               inum_obs = iend - ista + 1  
     547               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     548               DO iin=1,2   
     549                  DO ijn=1,2  
     550                     
     551                     IF ( k1dint == 1 ) THEN  
     552                        CALL obs_int_z1d_spl( kpk, &  
     553                           &    zint2(iin,ijn,:,iobs),&  
     554                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     555                           &    zmask2(iin,ijn,:,iobs))  
     556   
     557                     ENDIF  
     558        
     559                     CALL obs_level_search(kpk, &  
     560                         &        zgdept(iin,ijn,:,iobs),&  
     561                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     562                         &        iv_indic)  
     563 
     564                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     565                         &          prodatqc%var(2)%vdep(ista:iend),     &  
     566                         &          zint2(iin,ijn,:,iobs),            &  
     567                         &          zobs2k,interp_corner(iin,ijn,:), &  
     568                         &          zgdept(iin,ijn,:,iobs),         &  
     569                         &          zmask2(iin,ijn,:,iobs) )       
     570          
     571                  ENDDO  
     572               ENDDO  
     573              
     574            ENDIF  
     575 
     576            !-------------------------------------------------------------  
     577            ! Compute the horizontal interpolation for every profile level  
     578            !-------------------------------------------------------------  
     579              
     580            DO ikn=1,inum_obs  
     581               iend=ista+ikn-1 
     582                   
     583               zweig(:,:,1) = 0._wp  
     584    
     585               ! This code forces the horizontal weights to be   
     586               ! zero IF the observation is below the bottom of the   
     587               ! corners of the interpolation nodes, Or if it is in   
     588               ! the mask. This is important for observations near   
     589               ! steep bathymetry  
     590               DO iin=1,2  
     591                  DO ijn=1,2  
     592      
     593                     depth_loop2: DO ik=kpk,2,-1  
     594                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     595                             
     596                           zweig(iin,ijn,1) = &   
     597                              & zweig2(iin,ijn,1) * &  
     598                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     599                              &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     600                             
     601                           EXIT depth_loop2  
     602 
     603                        ENDIF  
     604 
     605                     ENDDO depth_loop2  
     606      
     607                  ENDDO  
     608               ENDDO  
     609    
     610               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     611                  &              prodatqc%var(2)%vmod(iend:iend) )  
     612 
     613                  ! Set QC flag for any observations found below the bottom 
     614                  ! needed as the check here is more strict than that in obs_prep 
     615               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
     616  
     617            ENDDO  
     618  
     619            DEALLOCATE(interp_corner,iv_indic)  
     620           
     621         ENDIF  
     622 
     623      ENDDO 
    454624 
    455625      ! Deallocate the data for interpolation 
     
    466636         & zmask2, & 
    467637         & zint1,  & 
    468          & zint2   & 
     638         & zint2,  & 
     639         & zgdept, & 
     640         & zgdepw  & 
    469641         & ) 
    470642 
     
    481653   END SUBROUTINE obs_prof_opt 
    482654 
    483    SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &  
    484       &                    ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, &  
    485       &                    kdailyavtypes )  
    486       !!-----------------------------------------------------------------------  
    487       !!  
    488       !!                     ***  ROUTINE obs_pro_opt  ***  
    489       !!  
    490       !! ** Purpose : Compute the model counterpart of profiles  
    491       !!              data by interpolating from the model grid to the   
    492       !!              observation point. Generalised vertical coordinate version  
    493       !!  
    494       !! ** Method  : Linearly interpolate to each observation point using   
    495       !!              the model values at the corners of the surrounding grid box.  
    496       !!  
    497       !!          First, model values on the model grid are interpolated vertically to the  
    498       !!          Depths of the profile observations.  Two vertical interpolation schemes are  
    499       !!          available:  
    500       !!          - linear       (k1dint = 0)  
    501       !!          - Cubic spline (k1dint = 1)     
    502       !!  
    503       !!  
    504       !!         Secondly the interpolated values are interpolated horizontally to the   
    505       !!         obs (lon, lat) point.  
    506       !!         Several horizontal interpolation schemes are available:  
    507       !!        - distance-weighted (great circle) (k2dint = 0)  
    508       !!        - distance-weighted (small angle)  (k2dint = 1)  
    509       !!        - bilinear (geographical grid)     (k2dint = 2)  
    510       !!        - bilinear (quadrilateral grid)    (k2dint = 3)  
    511       !!        - polynomial (quadrilateral grid)  (k2dint = 4)  
    512       !!  
    513       !!    For the cubic spline the 2nd derivative of the interpolating   
    514       !!    polynomial is computed before entering the vertical interpolation   
    515       !!    routine.  
    516       !!  
    517       !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is  
    518       !!    a daily mean model temperature field. So, we first compute  
    519       !!    the mean, then interpolate only at the end of the day.  
    520       !!  
    521       !!    This is the procedure to be used with generalised vertical model   
    522       !!    coordinates (ie s-coordinates. It is ~4x slower than the equivalent  
    523       !!    horizontal then vertical interpolation algorithm, but can deal with situations  
    524       !!    where the model levels are not flat.  
    525       !!    ONLY PERFORMED if ln_sco=.TRUE.   
    526       !!        
    527       !!    Note: the in situ temperature observations must be converted  
    528       !!    to potential temperature (the model variable) prior to  
    529       !!    assimilation.   
    530       !!??????????????????????????????????????????????????????????????  
    531       !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???  
    532       !!??????????????????????????????????????????????????????????????  
    533       !!  
    534       !! ** Action  :  
    535       !!  
    536       !! History :  
    537       !!      ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised  
    538       !!                           vertical coordinates 
    539       !!-----------------------------------------------------------------------  
    540     
    541       !! * Modules used  
    542       USE obs_profiles_def   ! Definition of storage space for profile obs.  
    543         
    544       IMPLICIT NONE  
    545   
    546       !! * Arguments  
    547       TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening  
    548       INTEGER, INTENT(IN) :: kt        ! Time step  
    549       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters  
    550       INTEGER, INTENT(IN) :: kpj  
    551       INTEGER, INTENT(IN) :: kpk  
    552       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step   
    553                                        !   (kit000-1 = restart time)  
    554       INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)  
    555       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)  
    556       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                      
    557       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
    558          & ptn,    &    ! Model temperature field  
    559          & psn,    &    ! Model salinity field  
    560          & ptmask       ! Land-sea mask  
    561       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
    562          & pgdept,  &    ! Model array of depth T levels     
    563          & pgdepw       ! Model array of depth W levels  
    564       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &  
    565          & kdailyavtypes   ! Types for daily averages  
    566        
    567       !! * Local declarations  
    568       INTEGER ::   ji  
    569       INTEGER ::   jj  
    570       INTEGER ::   jk  
    571       INTEGER ::   iico, ijco  
    572       INTEGER ::   jobs  
    573       INTEGER ::   inrc  
    574       INTEGER ::   ipro  
    575       INTEGER ::   idayend  
    576       INTEGER ::   ista  
    577       INTEGER ::   iend  
    578       INTEGER ::   iobs  
    579       INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
    580       INTEGER, DIMENSION(imaxavtypes) :: &  
    581          & idailyavtypes  
    582       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &  
    583          & igrdi, &  
    584          & igrdj  
    585       INTEGER :: &  
    586          & inum_obs 
    587       INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic     
    588       REAL(KIND=wp) :: zlam  
    589       REAL(KIND=wp) :: zphi  
    590       REAL(KIND=wp) :: zdaystp  
    591       REAL(KIND=wp), DIMENSION(kpk) :: &  
    592          & zobsmask, &  
    593          & zobsk,    &  
    594          & zobs2k  
    595       REAL(KIND=wp), DIMENSION(2,2,1) :: &  
    596          & zweig, &  
    597          & l_zweig  
    598       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &  
    599          & zmask, &  
    600          & zintt, &  
    601          & zints, &  
    602          & zinmt, &  
    603          & zgdept,&  
    604          & zgdepw,&  
    605          & zinms  
    606       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &  
    607          & zglam, &  
    608          & zgphi     
    609       REAL(KIND=wp), DIMENSION(1) :: zmsk_1        
    610       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner        
    611   
    612       !------------------------------------------------------------------------  
    613       ! Local initialization   
    614       !------------------------------------------------------------------------  
    615       ! ... Record and data counters  
    616       inrc = kt - kit000 + 2  
    617       ipro = prodatqc%npstp(inrc)  
    618    
    619       ! Daily average types  
    620       IF ( PRESENT(kdailyavtypes) ) THEN  
    621          idailyavtypes(:) = kdailyavtypes(:)  
    622       ELSE  
    623          idailyavtypes(:) = -1  
    624       ENDIF  
    625   
    626       ! Initialize daily mean for first time-step  
    627       idayend = MOD( kt - kit000 + 1, kdaystp )  
    628   
    629       ! Added kt == 0 test to catch restart case   
    630       IF ( idayend == 1 .OR. kt == 0) THEN  
    631            
    632          IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt  
    633          DO jk = 1, jpk  
    634             DO jj = 1, jpj  
    635                DO ji = 1, jpi  
    636                   prodatqc%vdmean(ji,jj,jk,1) = 0.0  
    637                   prodatqc%vdmean(ji,jj,jk,2) = 0.0  
    638                END DO  
    639             END DO  
    640          END DO  
    641         
    642       ENDIF  
    643         
    644       DO jk = 1, jpk  
    645          DO jj = 1, jpj  
    646             DO ji = 1, jpi  
    647                ! Increment the temperature field for computing daily mean  
    648                prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
    649                &                        + ptn(ji,jj,jk)  
    650                ! Increment the salinity field for computing daily mean  
    651                prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
    652                &                        + psn(ji,jj,jk)  
    653             END DO  
    654          END DO  
    655       END DO  
    656      
    657       ! Compute the daily mean at the end of day  
    658       zdaystp = 1.0 / REAL( kdaystp )  
    659       IF ( idayend == 0 ) THEN  
    660          DO jk = 1, jpk  
    661             DO jj = 1, jpj  
    662                DO ji = 1, jpi  
    663                   prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
    664                   &                        * zdaystp  
    665                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
    666                   &                           * zdaystp  
    667                END DO  
    668             END DO  
    669          END DO  
    670       ENDIF  
    671   
    672       ! Get the data for interpolation  
    673       ALLOCATE( &  
    674          & igrdi(2,2,ipro),      &  
    675          & igrdj(2,2,ipro),      &  
    676          & zglam(2,2,ipro),      &  
    677          & zgphi(2,2,ipro),      &  
    678          & zmask(2,2,kpk,ipro),  &  
    679          & zintt(2,2,kpk,ipro),  &  
    680          & zints(2,2,kpk,ipro),  &  
    681          & zgdept(2,2,kpk,ipro), &  
    682          & zgdepw(2,2,kpk,ipro)  &  
    683          & )  
    684   
    685       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
    686          iobs = jobs - prodatqc%nprofup  
    687          igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1  
    688          igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1  
    689          igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1  
    690          igrdj(1,2,iobs) = prodatqc%mj(jobs,1)  
    691          igrdi(2,1,iobs) = prodatqc%mi(jobs,1)  
    692          igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1  
    693          igrdi(2,2,iobs) = prodatqc%mi(jobs,1)  
    694          igrdj(2,2,iobs) = prodatqc%mj(jobs,1)  
    695       END DO  
    696       
    697       ! Initialise depth arrays 
    698       zgdept = 0.0 
    699       zgdepw = 0.0 
    700   
    701       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam )  
    702       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi )  
    703       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask )  
    704       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn,   zintt )  
    705       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn,   zints )  
    706       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), &  
    707         &                     zgdept )  
    708       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), &  
    709         &                     zgdepw )  
    710   
    711       ! At the end of the day also get interpolated means  
    712       IF ( idayend == 0 ) THEN  
    713   
    714          ALLOCATE( &  
    715             & zinmt(2,2,kpk,ipro),  &  
    716             & zinms(2,2,kpk,ipro)   &  
    717             & )  
    718   
    719          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &  
    720             &                  prodatqc%vdmean(:,:,:,1), zinmt )  
    721          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &  
    722             &                  prodatqc%vdmean(:,:,:,2), zinms )  
    723   
    724       ENDIF  
    725         
    726       ! Return if no observations to process  
    727       ! Has to be done after comm commands to ensure processors  
    728       ! stay in sync  
    729       IF ( ipro == 0 ) RETURN  
    730   
    731       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
    732      
    733          iobs = jobs - prodatqc%nprofup  
    734      
    735          IF ( kt /= prodatqc%mstp(jobs) ) THEN  
    736               
    737             IF(lwp) THEN  
    738                WRITE(numout,*)  
    739                WRITE(numout,*) ' E R R O R : Observation',              &  
    740                   &            ' time step is not consistent with the', &  
    741                   &            ' model time step'  
    742                WRITE(numout,*) ' ========='  
    743                WRITE(numout,*)  
    744                WRITE(numout,*) ' Record  = ', jobs,                    &  
    745                   &            ' kt      = ', kt,                      &  
    746                   &            ' mstp    = ', prodatqc%mstp(jobs), &  
    747                   &            ' ntyp    = ', prodatqc%ntyp(jobs)  
    748             ENDIF  
    749             CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )  
    750          ENDIF  
    751            
    752          zlam = prodatqc%rlam(jobs)  
    753          zphi = prodatqc%rphi(jobs)  
    754            
    755          ! Horizontal weights  
    756          ! Only calculated once, for both T and S.  
    757          ! Masked values are calculated later.   
    758   
    759          IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &  
    760             & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN  
    761   
    762             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &  
    763                &                   zglam(:,:,iobs), zgphi(:,:,iobs), &  
    764                &                   zmask(:,:,1,iobs), zweig, zmsk_1 )  
    765   
    766          ENDIF  
    767           
    768          ! IF zmsk_1 = 0; then ob is on land  
    769          IF (zmsk_1(1) < 0.1) THEN  
    770             WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'  
    771     
    772          ELSE   
    773               
    774             ! Temperature  
    775               
    776             IF ( prodatqc%npvend(jobs,1) > 0 ) THEN   
    777      
    778                zobsk(:) = obfillflt  
    779      
    780                IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
    781      
    782                   IF ( idayend == 0 )  THEN  
    783                     
    784                      ! Daily averaged moored buoy (MRB) data  
    785                     
    786                      ! vertically interpolate all 4 corners  
    787                      ista = prodatqc%npvsta(jobs,1)  
    788                      iend = prodatqc%npvend(jobs,1)  
    789                      inum_obs = iend - ista + 1  
    790                      ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    791        
    792                      DO iin=1,2  
    793                         DO ijn=1,2  
    794                                         
    795                                         
    796             
    797                            IF ( k1dint == 1 ) THEN  
    798                               CALL obs_int_z1d_spl( kpk, &  
    799                                  &     zinmt(iin,ijn,:,iobs), &  
    800                                  &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    801                                  &     zmask(iin,ijn,:,iobs))  
    802                            ENDIF  
    803         
    804                            CALL obs_level_search(kpk, &  
    805                               &    zgdept(iin,ijn,:,iobs), &  
    806                               &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
    807                               &    iv_indic)  
    808                            CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    809                               &    prodatqc%var(1)%vdep(ista:iend), &  
    810                               &    zinmt(iin,ijn,:,iobs), &  
    811                               &    zobs2k, interp_corner(iin,ijn,:), &  
    812                               &    zgdept(iin,ijn,:,iobs), &  
    813                               &    zmask(iin,ijn,:,iobs))  
    814         
    815                         ENDDO  
    816                      ENDDO  
    817                     
    818                     
    819                   ELSE  
    820                  
    821                      CALL ctl_stop( ' A nonzero' //     &  
    822                         &           ' number of profile T BUOY data should' // &  
    823                         &           ' only occur at the end of a given day' )  
    824      
    825                   ENDIF  
    826           
    827                ELSE   
    828                  
    829                   ! Point data  
    830       
    831                   ! vertically interpolate all 4 corners  
    832                   ista = prodatqc%npvsta(jobs,1)  
    833                   iend = prodatqc%npvend(jobs,1)  
    834                   inum_obs = iend - ista + 1  
    835                   ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    836                   DO iin=1,2   
    837                      DO ijn=1,2  
    838                                      
    839                                      
    840                         IF ( k1dint == 1 ) THEN  
    841                            CALL obs_int_z1d_spl( kpk, &  
    842                               &    zintt(iin,ijn,:,iobs),&  
    843                               &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    844                               &    zmask(iin,ijn,:,iobs))  
    845    
    846                         ENDIF  
    847         
    848                         CALL obs_level_search(kpk, &  
    849                             &        zgdept(iin,ijn,:,iobs),&  
    850                             &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
    851                             &         iv_indic)  
    852                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
    853                             &          prodatqc%var(1)%vdep(ista:iend),     &  
    854                             &          zintt(iin,ijn,:,iobs),            &  
    855                             &          zobs2k,interp_corner(iin,ijn,:), &  
    856                             &          zgdept(iin,ijn,:,iobs),         &  
    857                             &          zmask(iin,ijn,:,iobs) )       
    858           
    859                      ENDDO  
    860                   ENDDO  
    861               
    862                ENDIF  
    863         
    864                !-------------------------------------------------------------  
    865                ! Compute the horizontal interpolation for every profile level  
    866                !-------------------------------------------------------------  
    867               
    868                DO ikn=1,inum_obs  
    869                   iend=ista+ikn-1  
    870  
    871                   l_zweig(:,:,1) = 0._wp  
    872  
    873                   ! This code forces the horizontal weights to be   
    874                   ! zero IF the observation is below the bottom of the   
    875                   ! corners of the interpolation nodes, Or if it is in   
    876                   ! the mask. This is important for observations are near   
    877                   ! steep bathymetry  
    878                   DO iin=1,2  
    879                      DO ijn=1,2  
    880       
    881                         depth_loop1: DO ik=kpk,2,-1  
    882                            IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    883                              
    884                               l_zweig(iin,ijn,1) = &   
    885                                  & zweig(iin,ijn,1) * &  
    886                                  & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    887                                  &  - prodatqc%var(1)%vdep(iend)),0._wp)  
    888                              
    889                               EXIT depth_loop1  
    890                            ENDIF  
    891                         ENDDO depth_loop1  
    892       
    893                      ENDDO  
    894                   ENDDO  
    895     
    896                   CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
    897                   &          prodatqc%var(1)%vmod(iend:iend) )  
    898   
    899                ENDDO  
    900   
    901   
    902                DEALLOCATE(interp_corner,iv_indic)  
    903            
    904             ENDIF  
    905         
    906   
    907             ! Salinity   
    908            
    909             IF ( prodatqc%npvend(jobs,2) > 0 ) THEN   
    910      
    911                zobsk(:) = obfillflt  
    912      
    913                IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
    914      
    915                   IF ( idayend == 0 )  THEN  
    916                     
    917                      ! Daily averaged moored buoy (MRB) data  
    918                     
    919                      ! vertically interpolate all 4 corners  
    920                      ista = prodatqc%npvsta(iobs,2)  
    921                      iend = prodatqc%npvend(iobs,2)  
    922                      inum_obs = iend - ista + 1  
    923                      ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    924        
    925                      DO iin=1,2  
    926                         DO ijn=1,2  
    927                                         
    928                                         
    929             
    930                            IF ( k1dint == 1 ) THEN  
    931                               CALL obs_int_z1d_spl( kpk, &  
    932                                  &     zinms(iin,ijn,:,iobs), &  
    933                                  &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    934                                  &     zmask(iin,ijn,:,iobs))  
    935                            ENDIF  
    936         
    937                            CALL obs_level_search(kpk, &  
    938                               &    zgdept(iin,ijn,:,iobs), &  
    939                               &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    940                               &    iv_indic)  
    941                            CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    942                               &    prodatqc%var(2)%vdep(ista:iend), &  
    943                               &    zinms(iin,ijn,:,iobs), &  
    944                               &    zobs2k, interp_corner(iin,ijn,:), &  
    945                               &    zgdept(iin,ijn,:,iobs), &  
    946                               &    zmask(iin,ijn,:,iobs))  
    947         
    948                         ENDDO  
    949                      ENDDO  
    950                     
    951                     
    952                   ELSE  
    953                  
    954                      CALL ctl_stop( ' A nonzero' //     &  
    955                         &           ' number of profile T BUOY data should' // &  
    956                         &           ' only occur at the end of a given day' )  
    957      
    958                   ENDIF  
    959           
    960                ELSE   
    961                  
    962                   ! Point data  
    963       
    964                   ! vertically interpolate all 4 corners  
    965                   ista = prodatqc%npvsta(jobs,2)  
    966                   iend = prodatqc%npvend(jobs,2)  
    967                   inum_obs = iend - ista + 1  
    968                   ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    969                     
    970                   DO iin=1,2      
    971                      DO ijn=1,2   
    972                                   
    973                                   
    974                         IF ( k1dint == 1 ) THEN  
    975                            CALL obs_int_z1d_spl( kpk, &  
    976                               &    zints(iin,ijn,:,iobs),&  
    977                               &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    978                               &    zmask(iin,ijn,:,iobs))  
    979    
    980                         ENDIF  
    981         
    982                         CALL obs_level_search(kpk, &  
    983                            &        zgdept(iin,ijn,:,iobs),&  
    984                            &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    985                            &         iv_indic)  
    986                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,  &  
    987                            &          prodatqc%var(2)%vdep(ista:iend),     &  
    988                            &          zints(iin,ijn,:,iobs),               &  
    989                            &          zobs2k,interp_corner(iin,ijn,:),     &  
    990                            &          zgdept(iin,ijn,:,iobs),              &  
    991                            &          zmask(iin,ijn,:,iobs) )       
    992           
    993                      ENDDO  
    994                   ENDDO  
    995               
    996                ENDIF  
    997         
    998                !-------------------------------------------------------------  
    999                ! Compute the horizontal interpolation for every profile level  
    1000                !-------------------------------------------------------------  
    1001               
    1002                DO ikn=1,inum_obs  
    1003                   iend=ista+ikn-1  
    1004  
    1005                   l_zweig(:,:,1) = 0._wp 
    1006     
    1007                   ! This code forces the horizontal weights to be   
    1008                   ! zero IF the observation is below the bottom of the   
    1009                   ! corners of the interpolation nodes, Or if it is in   
    1010                   ! the mask. This is important for observations are near   
    1011                   ! steep bathymetry  
    1012                   DO iin=1,2  
    1013                      DO ijn=1,2  
    1014       
    1015                         depth_loop2: DO ik=kpk,2,-1  
    1016                            IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    1017                              
    1018                               l_zweig(iin,ijn,1) = &   
    1019                                  &  zweig(iin,ijn,1) * &  
    1020                                  &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    1021                                  &  - prodatqc%var(2)%vdep(iend)),0._wp)  
    1022                              
    1023                               EXIT depth_loop2  
    1024                            ENDIF  
    1025                         ENDDO depth_loop2  
    1026       
    1027                      ENDDO  
    1028                   ENDDO  
    1029     
    1030                   CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
    1031                   &          prodatqc%var(2)%vmod(iend:iend) )  
    1032   
    1033                ENDDO  
    1034   
    1035   
    1036                DEALLOCATE(interp_corner,iv_indic)  
    1037            
    1038             ENDIF  
    1039            
    1040          ENDIF  
    1041         
    1042       END DO  
    1043       
    1044       ! Deallocate the data for interpolation  
    1045       DEALLOCATE( &  
    1046          & igrdi, &  
    1047          & igrdj, &  
    1048          & zglam, &  
    1049          & zgphi, &  
    1050          & zmask, &  
    1051          & zintt, &  
    1052          & zints, &  
    1053          & zgdept,& 
    1054          & zgdepw & 
    1055          & )  
    1056       ! At the end of the day also get interpolated means  
    1057       IF ( idayend == 0 ) THEN  
    1058          DEALLOCATE( &  
    1059             & zinmt,  &  
    1060             & zinms   &  
    1061             & )  
    1062       ENDIF  
    1063      
    1064       prodatqc%nprofup = prodatqc%nprofup + ipro   
    1065         
    1066    END SUBROUTINE obs_pro_sco_opt  
    1067   
    1068    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         & 
    1069       &                    kit000, kdaystp, psurf, psurfmask, & 
    1070       &                    k2dint, ldnightav ) 
     655   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
     656      &                     kit000, kdaystp, psurf, psurfmask,   & 
     657      &                     k2dint, ldnightav, plamscl, pphiscl, & 
     658      &                     lindegrees ) 
    1071659 
    1072660      !!----------------------------------------------------------------------- 
     
    1090678      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    1091679      !! 
     680      !!    Two horizontal averaging schemes are also available: 
     681      !!        - weighted radial footprint        (k2dint = 5) 
     682      !!        - weighted rectangular footprint   (k2dint = 6) 
     683      !! 
    1092684      !! 
    1093685      !! ** Action  : 
     
    1096688      !!      ! 07-03 (A. Weaver) 
    1097689      !!      ! 15-02 (M. Martin) Combined routine for surface types 
     690      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
    1098691      !!----------------------------------------------------------------------- 
    1099692 
     
    1117710         & psurfmask                   ! Land-sea mask 
    1118711      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     712      REAL(KIND=wp), INTENT(IN) :: & 
     713         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions 
     714         & pphiscl                     ! This is the full width (rather than half-width) 
     715      LOGICAL, INTENT(IN) :: & 
     716         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 
    1119717 
    1120718      !! * Local declarations 
     
    1125723      INTEGER :: isurf 
    1126724      INTEGER :: iobs 
     725      INTEGER :: imaxifp, imaxjfp 
     726      INTEGER :: imodi, imodj 
    1127727      INTEGER :: idayend 
    1128728      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    1129          & igrdi, & 
    1130          & igrdj 
     729         & igrdi,   & 
     730         & igrdj,   & 
     731         & igrdip1, & 
     732         & igrdjp1 
    1131733      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1132734         & icount_night,      & 
     
    1136738      REAL(wp), DIMENSION(1) :: zext, zobsmask 
    1137739      REAL(wp) :: zdaystp 
    1138       REAL(wp), DIMENSION(2,2,1) :: & 
    1139          & zweig 
    1140740      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     741         & zweig,  & 
    1141742         & zmask,  & 
    1142743         & zsurf,  & 
    1143744         & zsurfm, & 
     745         & zsurftmp, & 
    1144746         & zglam,  & 
    1145          & zgphi 
     747         & zgphi,  & 
     748         & zglamf, & 
     749         & zgphif 
     750 
    1146751      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1147752         & zintmp,  & 
     
    1155760      inrc = kt - kit000 + 2 
    1156761      isurf = surfdataqc%nsstp(inrc) 
     762 
     763      ! Work out the maximum footprint size for the  
     764      ! interpolation/averaging in model grid-points - has to be even. 
     765 
     766      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     767 
    1157768 
    1158769      IF ( ldnightav ) THEN 
     
    1221832 
    1222833      ALLOCATE( & 
    1223          & igrdi(2,2,isurf), & 
    1224          & igrdj(2,2,isurf), & 
    1225          & zglam(2,2,isurf), & 
    1226          & zgphi(2,2,isurf), & 
    1227          & zmask(2,2,isurf), & 
    1228          & zsurf(2,2,isurf)  & 
     834         & zweig(imaxifp,imaxjfp,1),      & 
     835         & igrdi(imaxifp,imaxjfp,isurf), & 
     836         & igrdj(imaxifp,imaxjfp,isurf), & 
     837         & zglam(imaxifp,imaxjfp,isurf), & 
     838         & zgphi(imaxifp,imaxjfp,isurf), & 
     839         & zmask(imaxifp,imaxjfp,isurf), & 
     840         & zsurf(imaxifp,imaxjfp,isurf), & 
     841         & zsurftmp(imaxifp,imaxjfp,isurf),  & 
     842         & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
     843         & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
     844         & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
     845         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
    1229846         & ) 
    1230847 
    1231848      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    1232849         iobs = jobs - surfdataqc%nsurfup 
    1233          igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 
    1234          igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 
    1235          igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 
    1236          igrdj(1,2,iobs) = surfdataqc%mj(jobs) 
    1237          igrdi(2,1,iobs) = surfdataqc%mi(jobs) 
    1238          igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 
    1239          igrdi(2,2,iobs) = surfdataqc%mi(jobs) 
    1240          igrdj(2,2,iobs) = surfdataqc%mj(jobs) 
     850         DO ji = 0, imaxifp 
     851            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
     852             
     853            !Deal with wrap around in longitude 
     854            IF ( imodi < 1      ) imodi = imodi + jpiglo 
     855            IF ( imodi > jpiglo ) imodi = imodi - jpiglo 
     856             
     857            DO jj = 0, imaxjfp 
     858               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     859               !If model values are out of the domain to the north/south then 
     860               !set them to be the edge of the domain 
     861               IF ( imodj < 1      ) imodj = 1 
     862               IF ( imodj > jpjglo ) imodj = jpjglo 
     863 
     864               igrdip1(ji+1,jj+1,iobs) = imodi 
     865               igrdjp1(ji+1,jj+1,iobs) = imodj 
     866                
     867               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
     868                  igrdi(ji,jj,iobs) = imodi 
     869                  igrdj(ji,jj,iobs) = imodj 
     870               ENDIF 
     871                
     872            END DO 
     873         END DO 
    1241874      END DO 
    1242875 
    1243       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     876      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1244877         &                  igrdi, igrdj, glamt, zglam ) 
    1245       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     878      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1246879         &                  igrdi, igrdj, gphit, zgphi ) 
    1247       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     880      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1248881         &                  igrdi, igrdj, psurfmask, zmask ) 
    1249       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     882      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1250883         &                  igrdi, igrdj, psurf, zsurf ) 
     884      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     885         &                  igrdip1, igrdjp1, glamf, zglamf ) 
     886      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     887         &                  igrdip1, igrdjp1, gphif, zgphif ) 
    1251888 
    1252889      ! At the end of the day get interpolated means 
    1253       IF (ldnightav ) THEN 
    1254          IF ( idayend == 0 ) THEN 
    1255  
    1256             ALLOCATE( & 
    1257                & zsurfm(2,2,isurf)  & 
    1258                & ) 
    1259  
    1260             CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 
    1261                &               surfdataqc%vdmean(:,:), zsurfm ) 
    1262  
    1263          ENDIF 
     890      IF ( idayend == 0 .AND. ldnightav ) THEN 
     891 
     892         ALLOCATE( & 
     893            & zsurfm(imaxifp,imaxjfp,isurf)  & 
     894            & ) 
     895 
     896         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 
     897            &               surfdataqc%vdmean(:,:), zsurfm ) 
     898 
    1264899      ENDIF 
    1265900 
     
    1290925         zphi = surfdataqc%rphi(jobs) 
    1291926 
    1292          ! Get weights to interpolate the model value to the observation point 
    1293          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    1294             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    1295             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    1296  
    1297          ! Interpolate the model field to the observation point 
    1298927         IF ( ldnightav .AND. idayend == 0 ) THEN 
    1299928            ! Night-time averaged data 
    1300             CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext ) 
     929            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    1301930         ELSE 
    1302             CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext ) 
     931            zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
     932         ENDIF 
     933 
     934         IF ( k2dint <= 4 ) THEN 
     935 
     936            ! Get weights to interpolate the model value to the observation point 
     937            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     938               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     939               &                   zmask(:,:,iobs), zweig, zobsmask ) 
     940 
     941            ! Interpolate the model value to the observation point  
     942            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     943 
     944         ELSE 
     945 
     946            ! Get weights to average the model SLA to the observation footprint 
     947            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     948               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     949               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     950               &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     951               &                   lindegrees, zweig, zobsmask ) 
     952 
     953            ! Average the model SST to the observation footprint 
     954            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     955               &              zweig, zsurftmp(:,:,iobs),  zext ) 
     956 
    1303957         ENDIF 
    1304958 
     
    1310964            surfdataqc%rmod(jobs,1) = zext(1) 
    1311965         ENDIF 
     966          
     967         IF ( zext(1) == obfillflt ) THEN 
     968            ! If the observation value is a fill value, set QC flag to bad 
     969            surfdataqc%nqc(jobs) = 4 
     970         ENDIF 
    1312971 
    1313972      END DO 
     
    1315974      ! Deallocate the data for interpolation 
    1316975      DEALLOCATE( & 
     976         & zweig, & 
    1317977         & igrdi, & 
    1318978         & igrdj, & 
     
    1320980         & zgphi, & 
    1321981         & zmask, & 
    1322          & zsurf  & 
     982         & zsurf, & 
     983         & zsurftmp, & 
     984         & zglamf, & 
     985         & zgphif, & 
     986         & igrdip1,& 
     987         & igrdjp1 & 
    1323988         & ) 
    1324989 
    1325990      ! At the end of the day also deallocate night-time mean array 
    1326       IF ( ldnightav ) THEN 
    1327          IF ( idayend == 0 ) THEN 
    1328             DEALLOCATE( & 
    1329                & zsurfm  & 
    1330                & ) 
    1331          ENDIF 
     991      IF ( idayend == 0 .AND. ldnightav ) THEN 
     992         DEALLOCATE( & 
     993            & zsurfm  & 
     994            & ) 
    1332995      ENDIF 
    1333996 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r7646 r8667  
    2323   USE obs_oper           ! Observation operators 
    2424   USE lib_mpp, ONLY :   ctl_warn, ctl_stop 
     25   USE bdy_oce, ONLY : &        ! Boundary information 
     26      idx_bdy, nb_bdy, ln_bdy 
    2527 
    2628   IMPLICIT NONE 
     
    4042CONTAINS 
    4143 
    42    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 
     44   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
     45                            kqc_cutoff ) 
    4346      !!---------------------------------------------------------------------- 
    4447      !!                    ***  ROUTINE obs_pre_sla  *** 
     
    5760      !!        !  2015-02  (M. Martin) Combined routine for surface types. 
    5861      !!---------------------------------------------------------------------- 
     62      !! * Modules used 
     63      USE domstp              ! Domain: set the time-step 
    5964      USE par_oce             ! Ocean parameters 
    6065      USE dom_oce, ONLY       :   glamt, gphit, tmask, nproc   ! Geographical information 
     
    6368      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening 
    6469      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
    65       ! 
     70      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
     71      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
     72      !! * Local declarations 
     73      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    6674      INTEGER :: iyea0        ! Initial date 
    6775      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    7684      INTEGER :: inlasobs     !  - close to land 
    7785      INTEGER :: igrdobs      !  - fail the grid search 
     86      INTEGER :: ibdysobs     !  - close to open boundary 
    7887                              ! Global counters for observations that 
    7988      INTEGER :: iotdobsmpp     !  - outside time domain 
     
    8291      INTEGER :: inlasobsmpp    !  - close to land 
    8392      INTEGER :: igrdobsmpp     !  - fail the grid search 
    84       LOGICAL, DIMENSION(:), ALLOCATABLE ::   llvalid            ! SLA data selection 
     93      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     94      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     95         & llvalid            ! SLA data selection 
    8596      INTEGER :: jobs         ! Obs. loop variable 
    8697      INTEGER :: jstp         ! Time loop variable 
     
    107118      ilansobs = 0 
    108119      inlasobs = 0 
     120      ibdysobs = 0  
     121 
     122      ! Set QC cutoff to optional value if provided 
     123      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    109124 
    110125      ! ----------------------------------------------------------------------- 
     
    140155         &                 tmask(:,:,1), surfdata%nqc,  & 
    141156         &                 iosdsobs,     ilansobs,     & 
    142          &                 inlasobs,     ld_nea        ) 
     157         &                 inlasobs,     ld_nea,       & 
     158         &                 ibdysobs,     ld_bound_reject, & 
     159         &                 iqc_cutoff                     ) 
    143160 
    144161      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    145162      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    146163      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     164      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
    147165 
    148166      ! ----------------------------------------------------------------------- 
     
    155173      ALLOCATE( llvalid(surfdata%nsurf) ) 
    156174       
    157       ! We want all data which has qc flags <= 10 
    158  
    159       llvalid(:)  = ( surfdata%nqc(:)  <= 10 ) 
     175      ! We want all data which has qc flags <= iqc_cutoff 
     176 
     177      llvalid(:)  = ( surfdata%nqc(:)  <= iqc_cutoff ) 
    160178 
    161179      ! The actual copying 
     
    190208               &            inlasobsmpp 
    191209         ENDIF 
     210         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
     211            &            ibdysobsmpp   
    192212         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    193213            &            surfdataqc%nsurfmpp 
     
    225245      &                     kpi, kpj, kpk, & 
    226246      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
    227       &                     ld_nea, kdailyavtypes ) 
     247      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff ) 
    228248 
    229249!!---------------------------------------------------------------------- 
     
    241261      !! 
    242262      !!---------------------------------------------------------------------- 
    243       USE par_oce           ! Ocean parameters 
    244       USE dom_oce, ONLY :   gdept_1d, nproc   ! Geographical information 
     263      !! * Modules used 
     264      USE domstp              ! Domain: set the time-step 
     265      USE par_oce             ! Ocean parameters 
     266      USE dom_oce, ONLY : &   ! Geographical information 
     267         & gdept_1d,             & 
     268         & nproc 
    245269 
    246270      !! * Arguments 
     
    250274      LOGICAL, INTENT(IN) :: ld_var2 
    251275      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land 
     276      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary 
    252277      INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes 
    253278      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
     
    261286         & pgphi1, & 
    262287         & pgphi2 
     288      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    263289 
    264290      !! * Local declarations 
     291      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    265292      INTEGER :: iyea0        ! Initial date 
    266293      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    277304      INTEGER :: inlav1obs    !  - close to land (variable 1) 
    278305      INTEGER :: inlav2obs    !  - close to land (variable 2) 
     306      INTEGER :: ibdyv1obs    !  - boundary (variable 1)  
     307      INTEGER :: ibdyv2obs    !  - boundary (variable 2)       
    279308      INTEGER :: igrdobs      !  - fail the grid search 
    280309      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa 
     
    288317      INTEGER :: inlav1obsmpp !  - close to land (variable 1) 
    289318      INTEGER :: inlav2obsmpp !  - close to land (variable 2) 
     319      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1)  
     320      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)       
    290321      INTEGER :: igrdobsmpp   !  - fail the grid search 
    291322      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa 
     
    322353      inlav1obs = 0 
    323354      inlav2obs = 0 
    324       iuvchku  = 0 
    325       iuvchkv = 0 
     355      ibdyv1obs = 0 
     356      ibdyv2obs = 0 
     357      iuvchku   = 0 
     358      iuvchkv   = 0 
     359 
     360 
     361      ! Set QC cutoff to optional value if provided 
     362      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    326363 
    327364      ! ----------------------------------------------------------------------- 
     
    335372            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    336373            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    337             &              iotdobs, kdailyavtypes = kdailyavtypes ) 
     374            &              iotdobs, kdailyavtypes = kdailyavtypes,         & 
     375            &              kqc_cutoff = iqc_cutoff ) 
    338376      ELSE 
    339377         CALL obs_coo_tim_prof( icycle, & 
     
    342380            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    343381            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    344             &              iotdobs ) 
     382            &              iotdobs,          kqc_cutoff = iqc_cutoff ) 
    345383      ENDIF 
    346384 
     
    359397 
    360398      ! ----------------------------------------------------------------------- 
    361       ! Reject all observations for profiles with nqc > 10 
    362       ! ----------------------------------------------------------------------- 
    363  
    364       CALL obs_pro_rej( profdata ) 
     399      ! Reject all observations for profiles with nqc > iqc_cutoff 
     400      ! ----------------------------------------------------------------------- 
     401 
     402      CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 
    365403 
    366404      ! ----------------------------------------------------------------------- 
     
    381419         &                 gdept_1d,              zmask1,               & 
    382420         &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    383          &                 iosdv1obs,              ilanv1obs,           & 
    384          &                 inlav1obs,              ld_nea                ) 
     421         &                 iosdv1obs,             ilanv1obs,            & 
     422         &                 inlav1obs,             ld_nea,               & 
     423         &                 ibdyv1obs,             ld_bound_reject,      & 
     424         &                 iqc_cutoff       ) 
    385425 
    386426      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
    387427      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 
    388428      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 
     429      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 
    389430 
    390431      ! Variable 2 
     
    400441         &                 gdept_1d,              zmask2,               & 
    401442         &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    402          &                 iosdv2obs,              ilanv2obs,           & 
    403          &                 inlav2obs,              ld_nea                ) 
     443         &                 iosdv2obs,             ilanv2obs,            & 
     444         &                 inlav2obs,             ld_nea,               & 
     445         &                 ibdyv2obs,             ld_bound_reject,      & 
     446         &                 iqc_cutoff       ) 
    404447 
    405448      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    406449      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    407450      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
     451      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
    408452 
    409453      ! ----------------------------------------------------------------------- 
     
    412456 
    413457      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    414          CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 
     458         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
    415459         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 
    416460         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 
     
    429473      END DO 
    430474 
    431       ! We want all data which has qc flags = 0 
    432  
    433       llvalid%luse(:) = ( profdata%nqc(:)  <= 10 ) 
     475      ! We want all data which has qc flags <= iqc_cutoff 
     476 
     477      llvalid%luse(:) = ( profdata%nqc(:)  <= iqc_cutoff ) 
    434478      DO jvar = 1,profdata%nvar 
    435          llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 
     479         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 
    436480      END DO 
    437481 
     
    475519               &            iuvchku 
    476520         ENDIF 
     521         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 
     522               &            ibdyv1obsmpp 
    477523         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    478524            &            prodatqc%nvprotmpp(1) 
     
    492538               &            iuvchkv 
    493539         ENDIF 
     540         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
     541               &            ibdyv2obsmpp 
    494542         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    495543            &            prodatqc%nvprotmpp(2) 
     
    644692            &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 
    645693            kobsstp(jobs) = -1 
    646             kobsqc(jobs)  = kobsqc(jobs) + 11 
     694            kobsqc(jobs)  = IBSET(kobsqc(jobs),13) 
    647695            kotdobs       = kotdobs + 1 
    648696            CYCLE 
     
    695743         IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 
    696744            & .OR.( kobsstp(jobs) > nitend ) ) THEN 
    697             kobsqc(jobs) = kobsqc(jobs) + 12 
     745            kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    698746            kotdobs = kotdobs + 1 
    699747            CYCLE 
     
    739787      &                    kobsno,                                        & 
    740788      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   & 
    741       &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes ) 
     789      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, & 
     790      &                    kqc_cutoff ) 
    742791      !!---------------------------------------------------------------------- 
    743792      !!                    ***  ROUTINE obs_coo_tim *** 
     
    783832      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    784833         & kdailyavtypes    ! Types for daily averages 
     834      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     835 
    785836      !! * Local declarations 
    786837      INTEGER :: jobs 
     838      INTEGER :: iqc_cutoff=255 
    787839 
    788840      !----------------------------------------------------------------------- 
     
    803855         DO jobs = 1, kobsno 
    804856             
    805             IF ( kobsqc(jobs) <= 10 ) THEN 
     857            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 
    806858                
    807859               IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 
    808860                  & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 
    809                   kobsqc(jobs) = kobsqc(jobs) + 14 
     861                  kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    810862                  kotdobs      = kotdobs + 1 
    811863                  CYCLE 
     
    850902      DO jobs = 1, kobsno 
    851903         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 
    852             kobsqc(jobs) = kobsqc(jobs) + 18 
     904            kobsqc(jobs) = IBSET(kobsqc(jobs),12) 
    853905            kgrdobs = kgrdobs + 1 
    854906         ENDIF 
     
    861913      &                       plam,   pphi,    pmask,            & 
    862914      &                       kobsqc, kosdobs, klanobs,          & 
    863       &                       knlaobs,ld_nea                     ) 
     915      &                       knlaobs,ld_nea,                    & 
     916      &                       kbdyobs,ld_bound_reject,           & 
     917      &                       kqc_cutoff                         ) 
    864918      !!---------------------------------------------------------------------- 
    865919      !!                    ***  ROUTINE obs_coo_spc_2d  *** 
     
    894948      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 
    895949         & kobsqc             ! Observation quality control 
    896       INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain 
    897       INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell 
    898       INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land 
    899       LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land 
     950      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain 
     951      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell 
     952      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land 
     953      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary 
     954      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
     955      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
     956      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value 
     957 
    900958      !! * Local declarations 
    901959      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    902960         & zgmsk              ! Grid mask 
     961 
     962      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
     963         & zbmsk              ! Boundary mask 
     964      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
    903965      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    904966         & zglam, &           ! Model longitude at grid points 
     
    917979         ! For invalid points use 2,2 
    918980 
    919          IF ( kobsqc(jobs) >= 10 ) THEN 
     981         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 
    920982 
    921983            igrdi(1,1,jobs) = 1 
     
    9421004 
    9431005      END DO 
     1006 
     1007      IF (ln_bdy) THEN 
     1008        ! Create a mask grid points in boundary rim 
     1009        IF (ld_bound_reject) THEN 
     1010           zbdymask(:,:) = 1.0_wp 
     1011           DO ji = 1, nb_bdy 
     1012              DO jj = 1, idx_bdy(ji)%nblen(1) 
     1013                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1014              ENDDO 
     1015           ENDDO 
     1016 
     1017           CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     1018        ENDIF 
     1019      ENDIF 
     1020 
    9441021       
    9451022      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 
     
    9501027 
    9511028         ! Skip bad observations 
    952          IF ( kobsqc(jobs) >= 10 ) CYCLE 
     1029         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 
    9531030 
    9541031         ! Flag if the observation falls outside the model spatial domain 
     
    9571034            &  .OR. ( pobsphi(jobs) <  -90. ) & 
    9581035            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN 
    959             kobsqc(jobs) = kobsqc(jobs) + 11 
     1036            kobsqc(jobs) = IBSET(kobsqc(jobs),11) 
    9601037            kosdobs = kosdobs + 1 
    9611038            CYCLE 
     
    9641041         ! Flag if the observation falls with a model land cell 
    9651042         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
    966             kobsqc(jobs) = kobsqc(jobs)  + 12 
     1043            kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9671044            klanobs = klanobs + 1 
    9681045            CYCLE 
     
    9781055               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 
    9791056                  & .AND. & 
    980                   & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) & 
    981                   & ) THEN 
     1057                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 
     1058                  & < 1.0e-6_wp ) ) THEN 
    9821059                  lgridobs = .TRUE. 
    9831060                  iig = ji 
     
    9921069         IF (lgridobs) THEN 
    9931070            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
    994                kobsqc(jobs) = kobsqc(jobs) + 12 
     1071               kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9951072               klanobs = klanobs + 1 
    9961073               CYCLE 
     
    10001077         ! Flag if the observation falls is close to land 
    10011078         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 
    1002             IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14 
    10031079            knlaobs = knlaobs + 1 
    1004             CYCLE 
     1080            IF (ld_nea) THEN 
     1081               kobsqc(jobs) = IBSET(kobsqc(jobs),9) 
     1082               CYCLE 
     1083            ENDIF 
     1084         ENDIF 
     1085 
     1086         IF (ln_bdy) THEN 
     1087         ! Flag if the observation falls close to the boundary rim 
     1088           IF (ld_bound_reject) THEN 
     1089              IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1090                 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1091                 kbdyobs = kbdyobs + 1 
     1092                 CYCLE 
     1093              ENDIF 
     1094              ! for observations on the grid... 
     1095              IF (lgridobs) THEN 
     1096                 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1097                    kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1098                    kbdyobs = kbdyobs + 1 
     1099                    CYCLE 
     1100                 ENDIF 
     1101              ENDIF 
     1102            ENDIF 
    10051103         ENDIF 
    10061104             
     
    10151113      &                       plam,    pphi,    pdep,    pmask, & 
    10161114      &                       kpobsqc, kobsqc,  kosdobs,        & 
    1017       &                       klanobs, knlaobs, ld_nea          ) 
     1115      &                       klanobs, knlaobs, ld_nea,         & 
     1116      &                       kbdyobs, ld_bound_reject,         & 
     1117      &                       kqc_cutoff                        ) 
    10181118      !!---------------------------------------------------------------------- 
    10191119      !!                    ***  ROUTINE obs_coo_spc_3d  *** 
     
    10771177      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell 
    10781178      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land 
     1179      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary 
    10791180      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land 
     1181      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary 
     1182      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value 
     1183 
    10801184      !! * Local declarations 
    10811185      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10821186         & zgmsk              ! Grid mask 
     1187      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
     1188         & zbmsk              ! Boundary mask 
     1189      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
    10831190      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10841191         & zgdepw 
     
    11001207         ! For invalid points use 2,2 
    11011208 
    1102          IF ( kpobsqc(jobs) >= 10 ) THEN 
     1209         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 
    11031210             
    11041211            igrdi(1,1,jobs) = 1 
     
    11251232          
    11261233      END DO 
     1234 
     1235      IF (ln_bdy) THEN 
     1236        ! Create a mask grid points in boundary rim 
     1237        IF (ld_bound_reject) THEN            
     1238           zbdymask(:,:) = 1.0_wp 
     1239           DO ji = 1, nb_bdy 
     1240              DO jj = 1, idx_bdy(ji)%nblen(1) 
     1241                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1242              ENDDO 
     1243           ENDDO 
     1244        ENDIF 
     1245    
     1246        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     1247      ENDIF 
    11271248       
    11281249      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 
    11291250      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 
    11301251      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
    1131       IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 
    1132         ! Need to know the bathy depth for each observation for sco 
    1133         CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 
     1252      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 
    11341253        &                     zgdepw ) 
    1135       ENDIF 
    11361254 
    11371255      DO jobs = 1, kprofno 
    11381256 
    11391257         ! Skip bad profiles 
    1140          IF ( kpobsqc(jobs) >= 10 ) CYCLE 
     1258         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 
    11411259 
    11421260         ! Check if this observation is on a grid point 
     
    11491267               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 
    11501268                  & .AND. & 
    1151                   & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) & 
     1269                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) & 
    11521270                  & ) THEN 
    11531271                  lgridobs = .TRUE. 
     
    11761294               &  .OR. ( pobsdep(jobsp) < 0.0          )       & 
    11771295               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 
    1178                kobsqc(jobsp) = kobsqc(jobsp) + 11 
     1296               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 
    11791297               kosdobs = kosdobs + 1 
    11801298               CYCLE 
    11811299            ENDIF 
    11821300 
    1183             ! To check if an observations falls within land there are two cases: 
    1184             ! 1: z-coordibnates, where the check uses the mask 
    1185             ! 2: terrain following (eg s-coordinates),  
    1186             !    where we use the depth of the bottom cell to mask observations 
     1301            ! To check if an observations falls within land: 
    11871302              
    1188             IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 
     1303            ! Flag if the observation is deeper than the bathymetry 
     1304            ! Or if it is within the mask 
     1305            IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 
     1306               &     .OR. & 
     1307               &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
     1308               &  == 0.0_wp) ) THEN 
     1309               kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1310               klanobs = klanobs + 1 
     1311               CYCLE 
     1312            ENDIF 
    11891313                
    1190                ! Flag if the observation falls with a model land cell 
    1191                IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
    1192                   &  == 0.0_wp ) THEN 
    1193                   kobsqc(jobsp) = kobsqc(jobsp) + 12 
    1194                   klanobs = klanobs + 1 
    1195                   CYCLE 
    1196                ENDIF 
    1197               
    1198                ! Flag if the observation is close to land 
    1199               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
    1200                   &  0.0_wp) THEN 
    1201                   knlaobs = knlaobs + 1 
    1202                   IF (ld_nea) THEN    
    1203                      kobsqc(jobsp) = kobsqc(jobsp) + 14  
    1204                   ENDIF  
    1205                ENDIF 
    1206               
    1207             ELSE ! Case 2 
    1208   
    1209                ! Flag if the observation is deeper than the bathymetry 
    1210                ! Or if it is within the mask 
    1211                IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 
    1212                   &     .OR. & 
    1213                   &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
    1214                   &  == 0.0_wp) ) THEN 
    1215                   kobsqc(jobsp) = kobsqc(jobsp) + 12 
    1216                   klanobs = klanobs + 1 
    1217                   CYCLE 
    1218                ENDIF 
    1219                 
    1220                ! Flag if the observation is close to land 
    1221                IF ( ll_next_to_land ) THEN 
    1222                   knlaobs = knlaobs + 1 
    1223                   IF (ld_nea) THEN    
    1224                      kobsqc(jobsp) = kobsqc(jobsp) + 14  
    1225                   ENDIF  
    1226                ENDIF 
     1314            ! Flag if the observation is close to land 
     1315            IF ( ll_next_to_land ) THEN 
     1316               knlaobs = knlaobs + 1 
     1317               IF (ld_nea) THEN    
     1318                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1319               ENDIF  
    12271320            ENDIF 
    12281321             
     
    12321325            IF (lgridobs) THEN 
    12331326               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 
    1234                   kobsqc(jobsp) = kobsqc(jobsp) + 12 
     1327                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 
    12351328                  klanobs = klanobs + 1 
    12361329                  CYCLE 
     
    12501343            ENDIF 
    12511344             
     1345            IF (ln_bdy) THEN 
     1346               ! Flag if the observation falls close to the boundary rim 
     1347               IF (ld_bound_reject) THEN 
     1348                  IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1349                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 
     1350                     kbdyobs = kbdyobs + 1 
     1351                     CYCLE 
     1352                  ENDIF 
     1353                  ! for observations on the grid... 
     1354                  IF (lgridobs) THEN 
     1355                     IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1356                        kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 
     1357                        kbdyobs = kbdyobs + 1 
     1358                        CYCLE 
     1359                     ENDIF 
     1360                  ENDIF 
     1361               ENDIF 
     1362            ENDIF 
     1363             
    12521364         END DO 
    12531365      END DO 
     
    12551367   END SUBROUTINE obs_coo_spc_3d 
    12561368 
    1257    SUBROUTINE obs_pro_rej( profdata ) 
     1369   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 
    12581370      !!---------------------------------------------------------------------- 
    12591371      !!                    ***  ROUTINE obs_pro_rej *** 
     
    12731385      !! * Arguments 
    12741386      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data 
     1387      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value 
     1388 
    12751389      !! * Local declarations 
    12761390      INTEGER :: jprof 
     
    12821396      DO jprof = 1, profdata%nprof 
    12831397 
    1284          IF ( profdata%nqc(jprof) > 10 ) THEN 
     1398         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 
    12851399             
    12861400            DO jvar = 1, profdata%nvar 
     
    12901404                   
    12911405                  profdata%var(jvar)%nvqc(jobs) = & 
    1292                      & profdata%var(jvar)%nvqc(jobs) + 26 
     1406                     & IBSET(profdata%var(jvar)%nvqc(jobs),14) 
    12931407 
    12941408               END DO 
     
    13021416   END SUBROUTINE obs_pro_rej 
    13031417 
    1304    SUBROUTINE obs_uv_rej( profdata, knumu, knumv ) 
     1418   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 
    13051419      !!---------------------------------------------------------------------- 
    13061420      !!                    ***  ROUTINE obs_uv_rej *** 
     
    13221436      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected 
    13231437      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected 
     1438      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     1439 
    13241440      !! * Local declarations 
    13251441      INTEGER :: jprof 
     
    13411457         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 
    13421458             
    1343             IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. & 
    1344                & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN 
    1345                profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42 
     1459            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1460               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1461               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    13461462               knumv = knumv + 1 
    13471463            ENDIF 
    1348             IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. & 
    1349                & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN 
    1350                profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42 
     1464            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1465               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1466               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    13511467               knumu = knumu + 1 
    13521468            ENDIF 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90

    r6140 r8667  
    104104      ! Bookkeeping arrays with sizes equal to number of variables 
    105105 
    106       CHARACTER(len=6), POINTER, DIMENSION(:) :: & 
     106      CHARACTER(len=8), POINTER, DIMENSION(:) :: & 
    107107         & cvars          !: Variable names 
    108108 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

    r6140 r8667  
    8787      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 
    8888      CHARACTER(len=8) :: clrefdate 
    89       CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars 
     89      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
    9090      INTEGER :: jvar 
    9191      INTEGER :: ji 
     
    307307            inowin = 0 
    308308            DO ji = 1, inpfiles(jj)%nobs 
    309                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    310                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    311                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     309               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     310               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     311                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    312312               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    313313                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    325325            inowin = 0 
    326326            DO ji = 1, inpfiles(jj)%nobs 
    327                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    328                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    329                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     327               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     328               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     329                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    330330               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    331331                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    351351            inowin = 0 
    352352            DO ji = 1, inpfiles(jj)%nobs 
    353                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    354                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    355                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     353               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     354               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     355                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    356356               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    357357                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    373373 
    374374            DO ji = 1, inpfiles(jj)%nobs 
    375                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    376                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    377                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     375               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     376               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     377                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    378378               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    379379                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    388388                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    389389                           & CYCLE 
    390                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    391                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     390                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     391                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    392392                           ivar1t0 = ivar1t0 + 1 
    393393                        ENDIF 
     
    398398                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    399399                           & CYCLE 
    400                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    401                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     400                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     401                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    402402                           ivar2t0 = ivar2t0 + 1 
    403403                        ENDIF 
     
    407407                     IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    408408                        & CYCLE 
    409                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    410                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    411                         &     ldvar1 ) .OR. & 
    412                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    413                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
     409                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     410                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     411                        &    ldvar1 ) .OR. & 
     412                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     413                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
    414414                        &     ldvar2 ) ) THEN 
    415415                        ip3dt = ip3dt + 1 
     
    437437      DO jj = 1, inobf 
    438438         DO ji = 1, inpfiles(jj)%nobs 
    439             IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    440             IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    441                & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     439            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     440            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     441               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    442442            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    443443               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    452452      DO jj = 1, inobf 
    453453         DO ji = 1, inpfiles(jj)%nobs 
    454             IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    455             IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    456                & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     454            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     455            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     456               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    457457            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    458458               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    501501         ji = iprofidx(iindx(jk)) 
    502502 
    503          IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    504          IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    505             & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     503            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     504            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     505               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    506506 
    507507         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  & 
     
    518518            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    519519 
    520             IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    521                & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE 
     520            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     521            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     522               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    522523 
    523524            loop_prof : DO ij = 1, inpfiles(jj)%nlev 
     
    526527                  & CYCLE 
    527528 
    528                IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    529                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     529               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     530                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    530531 
    531532                  llvalprof = .TRUE.  
     
    534535               ENDIF 
    535536 
    536                IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    537                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     537               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     538                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    538539 
    539540                  llvalprof = .TRUE.  
     
    615616                  IF (ldsatt) THEN 
    616617 
    617                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    618                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    619                         &     ldvar1 ) .OR. & 
    620                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    621                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    622                         &     ldvar2 ) ) THEN 
     618                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     619                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     620                        &    ldvar1 ) .OR. & 
     621                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     622                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     623                        &   ldvar2 ) ) THEN 
    623624                        ip3dt = ip3dt + 1 
    624625                     ELSE 
     
    628629                  ENDIF 
    629630 
    630                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    631                      &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    632                      &       ldvar1 ) .OR. ldsatt ) THEN 
     631                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     632                    &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     633                    &    ldvar1 ) .OR. ldsatt ) THEN 
    633634 
    634635                     IF (ldsatt) THEN 
     
    661662 
    662663                     ! Profile var1 value 
    663                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    664                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     664                     IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     665                        & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    665666                        profdata%var(1)%vobs(ivar1t) = & 
    666667                           &                inpfiles(jj)%pob(ij,ji,1) 
     
    692693                  ENDIF 
    693694 
    694                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    695                      &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
     695                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     696                     &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    .AND. & 
    696697                     &   ldvar2 ) .OR. ldsatt ) THEN 
    697698 
     
    725726 
    726727                     ! Profile var2 value 
    727                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    728                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     728                  IF (  ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. & 
     729                    &   ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    ) ) THEN 
    729730                        profdata%var(2)%vobs(ivar2t) = & 
    730731                           &                inpfiles(jj)%pob(ij,ji,2) 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

    r6140 r8667  
    7777      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 
    7878      CHARACTER(len=8) :: clrefdate 
    79       CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars 
     79      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
    8080      INTEGER :: ji 
    8181      INTEGER :: jj 
     
    172172 
    173173            !------------------------------------------------------------------ 
    174             !  Read the profile file into inpfiles 
     174            !  Read the surface file into inpfiles 
    175175            !------------------------------------------------------------------ 
    176176            CALL init_obfbdata( inpfiles(jj) ) 
     
    296296                  ENDIF 
    297297                  llvalprof = .FALSE. 
    298                   IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 
    299                      & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 
     298                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 
    300299                     iobs = iobs + 1 
    301300                  ENDIF 
     
    370369            ! Set observation information 
    371370 
    372             IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 
    373                & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 
     371            IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 
    374372 
    375373               iobs = iobs + 1 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90

    r6140 r8667  
    154154 
    155155! mark any masked data with a QC flag 
    156          IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = 11 
     156         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15) 
    157157 
    158158         END DO 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_sstbias.F90

    r6140 r8667  
    11MODULE obs_sstbias 
    22   !!====================================================================== 
    3    !!                       ***  MODULE obs_readsstbias  *** 
    4    !! Observation diagnostics: Read the bias for SLA data 
     3   !!                       ***  MODULE obs_sstbias  *** 
     4   !! Observation diagnostics: Read the bias for SST data 
    55   !!====================================================================== 
    66   !!---------------------------------------------------------------------- 
    7    !!   obs_rea_sstbias : Driver for reading altimeter bias 
     7   !!   obs_app_sstbias : Driver for reading and applying the SST bias 
    88   !!---------------------------------------------------------------------- 
    99   !! * Modules used    
     
    139139               cl_bias_files(jtype) ) 
    140140            ! Get the SST bias data 
    141             CALL iom_get( numsstbias, jpdom_data, 'sst', z_sstbias_2d(:,:), 1 ) 
     141            CALL iom_get( numsstbias, jpdom_data, 'tn', z_sstbias_2d(:,:), 1 ) 
    142142            z_sstbias(:,:,jtype) = z_sstbias_2d(:,:)        
    143143            ! Close the file 
     
    190190               igrdj_tmp(:,:,jt) = igrdj(:,:,jobs) 
    191191               zglam_tmp(:,:,jt) = zglam(:,:,jobs) 
    192                zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
    193192               zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
    194193               zmask_tmp(:,:,jt) = zmask(:,:,jobs) 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90

    r6140 r8667  
    5050      INTEGER :: npj 
    5151      INTEGER :: nsurfup    !: Observation counter used in obs_oper 
     52      INTEGER :: nrec       !: Number of surface observation records in window 
    5253 
    5354      ! Arrays with size equal to the number of surface observations 
     
    5657         & mi,   &        !: i-th grid coord. for interpolating to surface observation 
    5758         & mj,   &        !: j-th grid coord. for interpolating to surface observation 
     59         & mt,   &        !: time record number for gridded data 
    5860         & nsidx,&        !: Surface observation number 
    5961         & nsfil,&        !: Surface observation number in file 
     
    6769         & ntyp           !: Type of surface observation product 
    6870 
    69       CHARACTER(len=6), POINTER, DIMENSION(:) :: & 
     71      CHARACTER(len=8), POINTER, DIMENSION(:) :: & 
    7072         & cvars          !: Variable names 
    7173 
     
    9395         & nsstpmpp       !: Global number of surface observations per time step 
    9496 
     97      ! Arrays with size equal to the number of observation records in the window 
     98      INTEGER, POINTER, DIMENSION(:) :: & 
     99         & mrecstp   ! Time step of the records 
     100 
    95101      ! Arrays used to store source indices when  
    96102      ! compressing obs_surf derived types 
     
    100106      INTEGER, POINTER, DIMENSION(:) :: & 
    101107         & nsind          !: Source indices of surface data in compressed data 
     108 
     109      ! Is this a gridded product? 
     110      
     111      LOGICAL :: lgrid 
    102112 
    103113   END TYPE obs_surf 
     
    160170         & surf%mi(ksurf),      & 
    161171         & surf%mj(ksurf),      & 
     172         & surf%mt(ksurf),      & 
    162173         & surf%nsidx(ksurf),   & 
    163174         & surf%nsfil(ksurf),   & 
     
    176187         & ) 
    177188 
     189      surf%mt(:) = -1 
     190 
    178191 
    179192      ! Allocate arrays of number of surface data size * number of variables 
     
    190203         & ) 
    191204 
     205      surf%rext(:,:) = 0.0_wp  
     206 
    192207      ! Allocate arrays of number of time step size 
    193208 
     
    217232 
    218233      surf%nsurfup     = 0 
     234       
     235      ! Not gridded by default 
     236           
     237      surf%lgrid       = .FALSE. 
    219238               
    220239   END SUBROUTINE obs_surf_alloc 
     
    242261         & surf%mi,      & 
    243262         & surf%mj,      & 
     263         & surf%mt,      & 
    244264         & surf%nsidx,   & 
    245265         & surf%nsfil,   & 
     
    370390            newsurf%mi(insurf)    = surf%mi(ji) 
    371391            newsurf%mj(insurf)    = surf%mj(ji) 
     392            newsurf%mt(insurf)    = surf%mt(ji) 
    372393            newsurf%nsidx(insurf) = surf%nsidx(ji) 
    373394            newsurf%nsfil(insurf) = surf%nsfil(ji) 
     
    414435      newsurf%nstp     = surf%nstp 
    415436      newsurf%cvars(:) = surf%cvars(:) 
     437       
     438      ! Set gridded stuff 
     439       
     440      newsurf%mt(insurf)    = surf%mt(ji) 
    416441  
    417442      ! Deallocate temporary data 
     
    454479         oldsurf%mi(jj)    = surf%mi(ji) 
    455480         oldsurf%mj(jj)    = surf%mj(ji) 
     481         oldsurf%mt(jj)    = surf%mt(ji) 
    456482         oldsurf%nsidx(jj) = surf%nsidx(ji) 
    457483         oldsurf%nsfil(jj) = surf%nsfil(ji) 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r6140 r8667  
    88   !!   obs_wri_prof   : Write profile observations in feedback format 
    99   !!   obs_wri_surf   : Write surface observations in feedback format 
    10    !!   obs_wri_stats : Print basic statistics on the data being written out 
     10   !!   obs_wri_stats  : Print basic statistics on the data being written out 
    1111   !!---------------------------------------------------------------------- 
    1212 
     
    8383      TYPE(obfbdata) :: fbdata 
    8484      CHARACTER(LEN=40) :: clfname 
    85       CHARACTER(LEN=6) :: clfiletype 
     85      CHARACTER(LEN=10) :: clfiletype 
    8686      INTEGER :: ilevel 
    8787      INTEGER :: jvar 
     
    196196         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:) 
    197197         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:) 
    198          IF ( profdata%nqc(jo) > 10 ) THEN 
    199             fbdata%ioqc(jo)    = 4 
     198         IF ( profdata%nqc(jo) > 255 ) THEN 
     199            fbdata%ioqc(jo)    = IBSET(profdata%nqc(jo),2) 
    200200            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo) 
    201             fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10 
     201            fbdata%ioqcf(2,jo) = profdata%nqc(jo) 
    202202         ELSE 
    203203            fbdata%ioqc(jo)    = profdata%nqc(jo) 
     
    236236               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk) 
    237237               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk) 
    238                IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN 
    239                   fbdata%ivlqc(ik,jo,jvar) = 4 
     238               IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN 
     239                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2) 
    240240                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 
    241                   fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10 
     241                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111') 
    242242               ELSE 
    243243                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) 
     
    320320      TYPE(obfbdata) :: fbdata 
    321321      CHARACTER(LEN=40) :: clfname         ! netCDF filename 
    322       CHARACTER(LEN=6) :: clfiletype 
     322      CHARACTER(LEN=10) :: clfiletype 
    323323      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 
    324324      INTEGER :: jo 
     
    395395         END DO 
    396396 
    397       CASE('ICECON') 
     397      CASE('ICECONC') 
    398398 
    399399         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     
    418418         END DO 
    419419 
     420      CASE('SSS') 
     421 
     422         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     423            &                 1 + iadd, iext, .TRUE. ) 
     424 
     425         clfiletype = 'sssfb' 
     426         fbdata%cname(1)      = surfdata%cvars(1) 
     427         fbdata%coblong(1)    = 'Sea surface salinity' 
     428         fbdata%cobunit(1)    = 'psu' 
     429         DO je = 1, iext 
     430            fbdata%cextname(je) = pext%cdname(je) 
     431            fbdata%cextlong(je) = pext%cdlong(je,1) 
     432            fbdata%cextunit(je) = pext%cdunit(je,1) 
     433         END DO 
     434         fbdata%caddlong(1,1) = 'Model interpolated SSS' 
     435         fbdata%caddunit(1,1) = 'psu' 
     436         fbdata%cgrid(1)      = 'T' 
     437         DO ja = 1, iadd 
     438            fbdata%caddname(1+ja) = padd%cdname(ja) 
     439            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     440            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     441         END DO 
     442 
     443      CASE DEFAULT 
     444 
     445         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 
     446 
    420447      END SELECT 
    421448 
     
    439466         fbdata%ivqc(jo,:)    = 0 
    440467         fbdata%ivqcf(:,jo,:) = 0 
    441          IF ( surfdata%nqc(jo) > 10 ) THEN 
     468         IF ( surfdata%nqc(jo) > 255 ) THEN 
    442469            fbdata%ioqc(jo)    = 4 
    443470            fbdata%ioqcf(1,jo) = 0 
    444             fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10 
     471            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 
    445472         ELSE 
    446473            fbdata%ioqc(jo)    = surfdata%nqc(jo) 
     
    474501         fbdata%idqc(1,jo)     = 0 
    475502         fbdata%idqcf(:,1,jo)  = 0 
    476          IF ( surfdata%nqc(jo) > 10 ) THEN 
     503         IF ( surfdata%nqc(jo) > 255 ) THEN 
    477504            fbdata%ivqc(jo,1)       = 4 
    478505            fbdata%ivlqc(1,jo,1)    = 4 
    479506            fbdata%ivlqcf(1,1,jo,1) = 0 
    480             fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10 
     507            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 
    481508         ELSE 
    482509            fbdata%ivqc(jo,1)       = surfdata%nqc(jo) 
  • branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90

    r2474 r8667  
    12401240         & zdum,  & 
    12411241         & zaamax 
    1242         
     1242       
     1243      imax = -1  
    12431244      ! Main computation 
    12441245      pflt = 1.0_wp 
Note: See TracChangeset for help on using the changeset viewer.