Changeset 8667

Timestamp:

2017-10-30T10:28:45+01:00 (6 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:

• diaobs.F90 (modified) (18 diffs)
• obs_mpp.F90 (modified) (2 diffs)
• obs_oper.F90 (modified) (24 diffs)
• obs_prep.F90 (modified) (57 diffs)
• obs_profiles_def.F90 (modified) (1 diff)
• obs_read_prof.F90 (modified) (19 diffs)
• obs_read_surf.F90 (modified) (4 diffs)
• obs_readmdt.F90 (modified) (1 diff)
• obs_sstbias.F90 (modified) (3 diffs)
• obs_surf_def.F90 (modified) (13 diffs)
• obs_write.F90 (modified) (9 diffs)
• obsinter_h2d.h90 (modified) (1 diff)

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

@@ -155 +155 @@
 
       NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
-         &            ln_sst, ln_sic, ln_vel3d,                       &
-         &            ln_altbias, ln_nea, ln_grid_global,             &
-         &            ln_grid_search_lookup,                          &
-         &            ln_ignmis, ln_s_at_t, ln_sstnight,              &
+         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
+         &            ln_altbias, ln_sstbias, ln_nea,                 &
+         &            ln_grid_global, ln_grid_search_lookup,          &
+         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
+         &            ln_sstnight,                                    &
+         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
+         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
          &            cn_profbfiles, cn_slafbfiles,                   &
          &            cn_sstfbfiles, cn_sicfbfiles,                   &
-         &            cn_velfbfiles, cn_altbiasfile,                  &
+         &            cn_velfbfiles, cn_sssfbfiles,                   &
+         &            cn_sstbiasfiles, cn_altbiasfile,                &
          &            cn_gridsearchfile, rn_gridsearchres,            &
-         &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     &
+         &            rn_dobsini, rn_dobsend,                         &
+         &            rn_sla_avglamscl, rn_sla_avgphiscl,             &
+         &            rn_sst_avglamscl, rn_sst_avgphiscl,             &
+         &            rn_sss_avglamscl, rn_sss_avgphiscl,             &
+         &            rn_sic_avglamscl, rn_sic_avgphiscl,             &
+         &            nn_1dint, nn_2dint,                             &
+         &            nn_2dint_sla, nn_2dint_sst,                     &
+         &            nn_2dint_sss, nn_2dint_sic,                     &
          &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
-         &            nn_profdavtypes, ln_sstbias, cn_sstbias_files
+         &            nn_profdavtypes
 
       INTEGER :: jnumsstbias
@@ -187 +198 @@
       cn_sicfbfiles(:) = ''
       cn_velfbfiles(:) = ''
+      cn_sssfbfiles(:)    = ''
       cn_sstbias_files(:) = ''
       nn_profdavtypes(:) = -1
@@ -208 +220 @@
          RETURN
       ENDIF
-      
-      !-----------------------------------------------------------------------
-      ! Set up list of observation types to be used
-      ! and the files associated with each type
-      !-----------------------------------------------------------------------
-
-      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
-      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) )
-
-      IF (ln_sstbias) THEN 
-         lmask(:) = .FALSE. 
-         WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE. 
-         jnumsstbias = COUNT(lmask) 
-         lmask(:) = .FALSE. 
-      ENDIF      
-
-      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
-         IF(lwp) WRITE(numout,cform_war)
-         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
-            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', &
-            &                    ' are set to .FALSE. so turning off calls to dia_obs'
-         nwarn = nwarn + 1
-         ln_diaobs = .FALSE.
-         RETURN
-      ENDIF
-
-      IF ( nproftypes > 0 ) THEN
-
-         ALLOCATE( cobstypesprof(nproftypes) )
-         ALLOCATE( ifilesprof(nproftypes) )
-         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
-
-         jtype = 0
-         IF (ln_t3d .OR. ln_s3d) THEN
-            jtype = jtype + 1
-            clproffiles(jtype,:) = cn_profbfiles(:)
-            cobstypesprof(jtype) = 'prof  '
-            ifilesprof(jtype) = 0
-            DO jfile = 1, jpmaxnfiles
-               IF ( trim(clproffiles(jtype,jfile)) /= '' ) &
-                  ifilesprof(jtype) = ifilesprof(jtype) + 1
-            END DO
-         ENDIF
-         IF (ln_vel3d) THEN
-            jtype = jtype + 1
-            clproffiles(jtype,:) = cn_velfbfiles(:)
-            cobstypesprof(jtype) = 'vel   '
-            ifilesprof(jtype) = 0
-            DO jfile = 1, jpmaxnfiles
-               IF ( trim(clproffiles(jtype,jfile)) /= '' ) &
-                  ifilesprof(jtype) = ifilesprof(jtype) + 1
-            END DO
-         ENDIF
-
-      ENDIF
-
-      IF ( nsurftypes > 0 ) THEN
-
-         ALLOCATE( cobstypessurf(nsurftypes) )
-         ALLOCATE( ifilessurf(nsurftypes) )
-         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
-
-         jtype = 0
-         IF (ln_sla) THEN
-            jtype = jtype + 1
-            clsurffiles(jtype,:) = cn_slafbfiles(:)
-            cobstypessurf(jtype) = 'sla   '
-            ifilessurf(jtype) = 0
-            DO jfile = 1, jpmaxnfiles
-               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
-                  ifilessurf(jtype) = ifilessurf(jtype) + 1
-            END DO
-         ENDIF
-         IF (ln_sst) THEN
-            jtype = jtype + 1
-            clsurffiles(jtype,:) = cn_sstfbfiles(:)
-            cobstypessurf(jtype) = 'sst   '
-            ifilessurf(jtype) = 0
-            DO jfile = 1, jpmaxnfiles
-               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
-                  ifilessurf(jtype) = ifilessurf(jtype) + 1
-            END DO
-         ENDIF
-#if defined key_lim2 || defined key_lim3
-         IF (ln_sic) THEN
-            jtype = jtype + 1
-            clsurffiles(jtype,:) = cn_sicfbfiles(:)
-            cobstypessurf(jtype) = 'sic   '
-            ifilessurf(jtype) = 0
-            DO jfile = 1, jpmaxnfiles
-               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
-                  ifilessurf(jtype) = ifilessurf(jtype) + 1
-            END DO
-         ENDIF
-#endif
-
-      ENDIF
-
-      !Write namelist settings to stdout
+
       IF(lwp) THEN
          WRITE(numout,*)
@@ -318 +232 @@
          WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
          WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
-         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global
-         WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias 
-         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup
+         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss
+         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global
+         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
          IF (ln_grid_search_lookup) &
             WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
@@ -328 +242 @@
          WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint
          WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
+         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject
          WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
          WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
          WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
          WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
+         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias
          WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
          WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
          WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
-         WRITE(numout,*) '          Number of profile obs types: ',nproftypes
-
-         IF ( nproftypes > 0 ) THEN
-            DO jtype = 1, nproftypes
-               DO jfile = 1, ifilesprof(jtype)
-                  WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', &
-                     TRIM(clproffiles(jtype,jfile))
-               END DO
-            END DO
+      ENDIF
+      !-----------------------------------------------------------------------
+      ! Set up list of observation types to be used
+      ! and the files associated with each type
+      !-----------------------------------------------------------------------
+
+      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
+      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) )
+
+      IF (ln_sstbias) THEN 
+         lmask(:) = .FALSE. 
+         WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE. 
+         jnumsstbias = COUNT(lmask) 
+         lmask(:) = .FALSE. 
+      ENDIF      
+
+      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
+         IF(lwp) WRITE(numout,cform_war)
+         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
+            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', &
+            &                    ' are set to .FALSE. so turning off calls to dia_obs'
+         nwarn = nwarn + 1
+         ln_diaobs = .FALSE.
+         RETURN
+      ENDIF
+
+      IF ( nproftypes > 0 ) THEN
+
+         ALLOCATE( cobstypesprof(nproftypes) )
+         ALLOCATE( ifilesprof(nproftypes) )
+         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
+
+         jtype = 0
+         IF (ln_t3d .OR. ln_s3d) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', &
+               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
          ENDIF
-
-         WRITE(numout,*)'          Number of surface obs types: ',nsurftypes
-         IF ( nsurftypes > 0 ) THEN
-            DO jtype = 1, nsurftypes
-               DO jfile = 1, ifilessurf(jtype)
-                  WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', &
-                     TRIM(clsurffiles(jtype,jfile))
-               END DO
-            END DO
+         IF (ln_vel3d) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', &
+               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
          ENDIF
-         WRITE(numout,*) '~~~~~~~~~~~~'
-
-      ENDIF
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         ALLOCATE( cobstypessurf(nsurftypes) )
+         ALLOCATE( ifilessurf(nsurftypes) )
+         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
+         ALLOCATE(n2dintsurf(nsurftypes))
+         ALLOCATE(ravglamscl(nsurftypes))
+         ALLOCATE(ravgphiscl(nsurftypes))
+         ALLOCATE(lfpindegs(nsurftypes))
+         ALLOCATE(llnightav(nsurftypes))
+
+         jtype = 0
+         IF (ln_sla) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', &
+               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      &
+               &                  nn_2dint, nn_2dint_sla,             &
+               &                  rn_sla_avglamscl, rn_sla_avgphiscl, &
+               &                  ln_sla_fp_indegs, .FALSE.,          &
+               &                  n2dintsurf, ravglamscl, ravgphiscl, &
+               &                  lfpindegs, llnightav )
+         ENDIF
+         IF (ln_sst) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', &
+               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      &
+               &                  nn_2dint, nn_2dint_sst,             &
+               &                  rn_sst_avglamscl, rn_sst_avgphiscl, &
+               &                  ln_sst_fp_indegs, ln_sstnight,      &
+               &                  n2dintsurf, ravglamscl, ravgphiscl, &
+               &                  lfpindegs, llnightav )
+         ENDIF
+#if defined key_lim2 || defined key_lim3 || defined key_cice
+         IF (ln_sic) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', &
+               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      &
+               &                  nn_2dint, nn_2dint_sic,             &
+               &                  rn_sic_avglamscl, rn_sic_avgphiscl, &
+               &                  ln_sic_fp_indegs, .FALSE.,          &
+               &                  n2dintsurf, ravglamscl, ravgphiscl, &
+               &                  lfpindegs, llnightav )
+         ENDIF
+#endif
+         IF (ln_sss) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', &
+               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      &
+               &                  nn_2dint, nn_2dint_sss,             &
+               &                  rn_sss_avglamscl, rn_sss_avgphiscl, &
+               &                  ln_sss_fp_indegs, .FALSE.,          &
+               &                  n2dintsurf, ravglamscl, ravgphiscl, &
+               &                  lfpindegs, llnightav )
+         ENDIF
+
+      ENDIF
+
+
 
       !-----------------------------------------------------------------------
@@ -377 +378 @@
       ENDIF
 
-      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN
+      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN
          CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
             &                    ' is not available')
@@ -442 +443 @@
                &               jpi, jpj, jpk, &
                &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
-               &               ln_nea, kdailyavtypes = nn_profdavtypes )
+               &               ln_nea, ln_bound_reject, &
+               &               kdailyavtypes = nn_profdavtypes )
 
          END DO
@@ -469 +471 @@
                &               clsurffiles(jtype,1:ifilessurf(jtype)), &
                &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
-               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav )
-         
-         
-            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea )
+               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
+
+            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
 
             IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
-               CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint )
-               IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile )
+               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
+               IF ( ln_altbias ) &
+                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
             ENDIF
+
             IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
-               !Read in bias field and correct SST.
-               IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// &
-                                                     "  but no bias"// &
-                                                     " files to read in")   
-                  CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, &
-                                        jnumsstbias, cn_sstbias_files(1:jnumsstbias) )
+               jnumsstbias = 0
+               DO jfile = 1, jpmaxnfiles
+                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
+                     &  jnumsstbias = jnumsstbias + 1
+               END DO
+               IF ( jnumsstbias == 0 ) THEN
+                  CALL ctl_stop("ln_sstbias set but no bias files to read in")    
+               ENDIF
+
+               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
+                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
+
             ENDIF
          END DO
@@ -545 +554 @@
          & frld
 #endif
+#if defined key_cice
+      USE sbc_oce, ONLY : fr_i     ! ice fraction
+#endif
+
       IMPLICIT NONE
 
@@ -561 +574 @@
          & zprofmask2              ! Mask associated with zprofvar2
       REAL(wp), POINTER, DIMENSION(:,:) :: &
-         & zsurfvar                ! Model values equivalent to surface ob.
+         & zsurfvar, &             ! Model values equivalent to surface ob.
+         & zsurfmask               ! Mask associated with surface variable
       REAL(wp), POINTER, DIMENSION(:,:) :: &
          & zglam1,    &            ! Model longitudes for prof variable 1
@@ -567 +581 @@
          & zgphi1,    &            ! Model latitudes for prof variable 1
          & zgphi2                  ! Model latitudes for prof variable 2
-#if ! defined key_lim2 && ! defined key_lim3
-      REAL(wp), POINTER, DIMENSION(:,:) :: frld
-#endif
       LOGICAL :: llnightav        ! Logical for calculating night-time average
 
@@ -578 +589 @@
       CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
       CALL wrk_alloc( jpi, jpj, zsurfvar )
+      CALL wrk_alloc( jpi, jpj, zsurfmask )
       CALL wrk_alloc( jpi, jpj, zglam1 )
       CALL wrk_alloc( jpi, jpj, zglam2 )
       CALL wrk_alloc( jpi, jpj, zgphi1 )
       CALL wrk_alloc( jpi, jpj, zgphi2 )
-#if ! defined key_lim2 && ! defined key_lim3
-      CALL wrk_alloc(jpi,jpj,frld) 
-#endif
 
       IF(lwp) THEN
@@ -594 +603 @@
       idaystp = NINT( rday / rdt )
 
-      !-----------------------------------------------------------------------
-      ! No LIM => frld == 0.0_wp
-      !-----------------------------------------------------------------------
-#if ! defined key_lim2 && ! defined key_lim3
-      frld(:,:) = 0.0_wp
-#endif
       !-----------------------------------------------------------------------
       ! Call the profile and surface observation operators
@@ -627 +630 @@
                zgphi1(:,:) = gphiu(:,:)
                zgphi2(:,:) = gphiv(:,:)
+            CASE DEFAULT
+               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
             END SELECT
 
-            IF( ln_zco .OR. ln_zps ) THEN 
-               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
-                  &               nit000, idaystp,                         &
-                  &               zprofvar1, zprofvar2,                    &
-                  &               gdept_1d, zprofmask1, zprofmask2,        &
-                  &               zglam1, zglam2, zgphi1, zgphi2,          &
-                  &               nn_1dint, nn_2dint,                      &
-                  &               kdailyavtypes = nn_profdavtypes )
-            ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN
-               !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION
-               CALL obs_pro_sco_opt( profdataqc(jtype),                    & 
-                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
-                  &              zprofvar1, zprofvar2,                   & 
-                  &              gdept_n(:,:,:), gdepw_n(:,:,:),           &
-                  &              tmask, nn_1dint, nn_2dint,              & 
-                  &              kdailyavtypes = nn_profdavtypes ) 
-            ELSE
-               CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// &
-                         'yet working for velocity data (turn off velocity observations')
-            ENDIF
+            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
+               &               nit000, idaystp,                         &
+               &               zprofvar1, zprofvar2,                    &
+               &               gdept_n(:,:,:), gdepw_n(:,:,:),            & 
+               &               zprofmask1, zprofmask2,                  &
+               &               zglam1, zglam2, zgphi1, zgphi2,          &
+               &               nn_1dint, nn_2dint,                      &
+               &               kdailyavtypes = nn_profdavtypes )
 
          END DO
@@ -657 +650 @@
 
          DO jtype = 1, nsurftypes
+
+            !Defaults which might be changed
+            zsurfmask(:,:) = tmask(:,:,1)
 
             SELECT CASE ( TRIM(cobstypessurf(jtype)) )
             CASE('sst')
                zsurfvar(:,:) = tsn(:,:,1,jp_tem)
-               llnightav = ln_sstnight
             CASE('sla')
                zsurfvar(:,:) = sshn(:,:)
-               llnightav = .FALSE.
-#if defined key_lim2 || defined key_lim3
+            CASE('sss')
+               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
             CASE('sic')
                IF ( kstp == 0 ) THEN
@@ -678 +673 @@
                   CYCLE
                ELSE
+#if defined key_cice
+                  zsurfvar(:,:) = fr_i(:,:)
+#elif defined key_lim2 || defined key_lim3
                   zsurfvar(:,:) = 1._wp - frld(:,:)
+#else
+               CALL ctl_stop( ' Trying to run sea-ice observation operator', &
+                  &           ' but no sea-ice model appears to have been defined' )
+#endif
                ENDIF
 
-               llnightav = .FALSE.
-#endif
             END SELECT
 
             CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
-               &               nit000, idaystp, zsurfvar, tmask(:,:,1), &
-               &               nn_2dint, llnightav )
+               &               nit000, idaystp, zsurfvar, zsurfmask,    &
+               &               n2dintsurf(jtype), llnightav(jtype),     &
+               &               ravglamscl(jtype), ravgphiscl(jtype),     &
+               &               lfpindegs(jtype) )
 
          END DO
@@ -698 +700 @@
       CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
       CALL wrk_dealloc( jpi, jpj, zsurfvar )
+      CALL wrk_dealloc( jpi, jpj, zsurfmask )
       CALL wrk_dealloc( jpi, jpj, zglam1 )
       CALL wrk_dealloc( jpi, jpj, zglam2 )
       CALL wrk_dealloc( jpi, jpj, zgphi1 )
       CALL wrk_dealloc( jpi, jpj, zgphi2 )
-#if ! defined key_lim2 && ! defined key_lim3
-      CALL wrk_dealloc(jpi,jpj,frld)
-#endif
 
    END SUBROUTINE dia_obs
@@ -821 +821 @@
 
       IF ( nsurftypes > 0 ) &
-         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf )
+         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
+         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
 
    END SUBROUTINE dia_obs_dealloc

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90

@@ -142 +142 @@
       !!
 
-      iobsp=kobsp
+      iobsp(:)=kobsp(:)
 
       WHERE( iobsp(:) == -1 )
@@ -148 +148 @@
       END WHERE
 
-      iobsp=-1*iobsp
+      iobsp(:)=-1*iobsp(:)
 
       CALL obs_mpp_max_integer( iobsp, kno )
 
-      kobsp=-1*iobsp
+      kobsp(:)=-1*iobsp(:)
 
       isum=0

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

@@ -9 +9 @@
    !!   obs_prof_opt :    Compute the model counterpart of profile data
    !!   obs_surf_opt :    Compute the model counterpart of surface data
-   !!   obs_pro_sco_opt: Compute the model counterpart of temperature and 
-   !!                    salinity observations from profiles in generalised 
-   !!                    vertical coordinates 
    !!----------------------------------------------------------------------
 
@@ -22 +19 @@
       & obs_int_h2d, &
       & obs_int_h2d_init
+   USE obs_averg_h2d, ONLY : &    ! Horizontal averaging to the obs footprint
+      & obs_avg_h2d, &
+      & obs_avg_h2d_init, &
+      & obs_max_fpsize
    USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt
       & obs_int_z1d,    &
@@ -28 +29 @@
       & obfillflt      ! Fillvalue   
    USE dom_oce,       ONLY : &
-      & glamt, glamu, glamv, &
-      & gphit, gphiu, gphiv, & 
+      & glamt, glamf, &
+      & gphit, gphif, & 
       & gdept_n, gdept_0 
    USE lib_mpp,       ONLY : &
@@ -44 +45 @@
 
    PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs
-      &   obs_pro_sco_opt, &  ! Compute the model counterpart of profile observations 
       &   obs_surf_opt     ! Compute the model counterpart of surface obs
 
@@ -290 +290 @@
       END DO
 
+      ! Initialise depth arrays
+      zgdept(:,:,:,:) = 0.0
+      zgdepw(:,:,:,:) = 0.0
+
       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
@@ -300 +304 @@
       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 )
 
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 
+
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
@@ -314 +321 @@
 
       ENDIF
+
+      ! Return if no observations to process 
+      ! Has to be done after comm commands to ensure processors 
+      ! stay in sync 
+      IF ( ipro == 0 ) RETURN 
 
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
@@ -339 +351 @@
          zphi = prodatqc%rphi(jobs)
 
-         ! Horizontal weights and vertical mask
-
+         ! Horizontal weights 
+         ! Masked values are calculated later.  
          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
 
-            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
+            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
                &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), &
-               &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 )
+               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 )
 
          ENDIF
@@ -351 +363 @@
          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
 
-            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
+            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
                &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), &
-               &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 )
+               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 )
  
          ENDIF
@@ -365 +377 @@
                IF ( idayend == 0 )  THEN
                   ! Daily averaged data
-                  CALL obs_int_h2d( kpk, kpk,      &
-                     &              zweig1, zinm1(:,:,:,iobs), zobsk )
-
-               ENDIF
+
+                  ! vertically interpolate all 4 corners 
+                  ista = prodatqc%npvsta(jobs,1) 
+                  iend = prodatqc%npvend(jobs,1) 
+                  inum_obs = iend - ista + 1 
+                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
+
+                  DO iin=1,2 
+                     DO ijn=1,2 
+
+                        IF ( k1dint == 1 ) THEN 
+                           CALL obs_int_z1d_spl( kpk, & 
+                              &     zinm1(iin,ijn,:,iobs), & 
+                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
+                              &     zmask1(iin,ijn,:,iobs)) 
+                        ENDIF 
+       
+                        CALL obs_level_search(kpk, & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                           &    iv_indic) 
+
+                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
+                           &    prodatqc%var(1)%vdep(ista:iend), & 
+                           &    zinm1(iin,ijn,:,iobs), & 
+                           &    zobs2k, interp_corner(iin,ijn,:), & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    zmask1(iin,ijn,:,iobs)) 
+       
+                     ENDDO 
+                  ENDDO 
+
+               ENDIF !idayend
 
             ELSE 
 
-               ! Point data
-               CALL obs_int_h2d( kpk, kpk,      &
-                  &              zweig1, zint1(:,:,:,iobs), zobsk )
-
-            ENDIF
-
-            !-------------------------------------------------------------
-            ! Compute vertical second-derivative of the interpolating 
-            ! polynomial at obs points
-            !-------------------------------------------------------------
-
-            IF ( k1dint == 1 ) THEN
-               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
-                  &                  pgdept, zobsmask1 )
-            ENDIF
-
-            !-----------------------------------------------------------------
-            !  Vertical interpolation to the observation point
-            !-----------------------------------------------------------------
-            ista = prodatqc%npvsta(jobs,1)
-            iend = prodatqc%npvend(jobs,1)
-            CALL obs_int_z1d( kpk,                &
-               & prodatqc%var(1)%mvk(ista:iend),  &
-               & k1dint, iend - ista + 1,         &
-               & prodatqc%var(1)%vdep(ista:iend), &
-               & zobsk, zobs2k,                   &
-               & prodatqc%var(1)%vmod(ista:iend), &
-               & pgdept, zobsmask1 )
-
-         ENDIF
-
+               ! Point data 
+     
+               ! vertically interpolate all 4 corners 
+               ista = prodatqc%npvsta(jobs,1) 
+               iend = prodatqc%npvend(jobs,1) 
+               inum_obs = iend - ista + 1 
+               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
+               DO iin=1,2  
+                  DO ijn=1,2 
+                    
+                     IF ( k1dint == 1 ) THEN 
+                        CALL obs_int_z1d_spl( kpk, & 
+                           &    zint1(iin,ijn,:,iobs),& 
+                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
+                           &    zmask1(iin,ijn,:,iobs)) 
+  
+                     ENDIF 
+       
+                     CALL obs_level_search(kpk, & 
+                         &        zgdept(iin,ijn,:,iobs),& 
+                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                         &        iv_indic) 
+
+                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
+                         &          prodatqc%var(1)%vdep(ista:iend),     & 
+                         &          zint1(iin,ijn,:,iobs),            & 
+                         &          zobs2k,interp_corner(iin,ijn,:), & 
+                         &          zgdept(iin,ijn,:,iobs),         & 
+                         &          zmask1(iin,ijn,:,iobs) )      
+         
+                  ENDDO 
+               ENDDO 
+             
+            ENDIF 
+
+            !------------------------------------------------------------- 
+            ! Compute the horizontal interpolation for every profile level 
+            !------------------------------------------------------------- 
+             
+            DO ikn=1,inum_obs 
+               iend=ista+ikn-1
+                  
+               zweig(:,:,1) = 0._wp 
+   
+               ! This code forces the horizontal weights to be  
+               ! zero IF the observation is below the bottom of the  
+               ! corners of the interpolation nodes, Or if it is in  
+               ! the mask. This is important for observations near  
+               ! steep bathymetry 
+               DO iin=1,2 
+                  DO ijn=1,2 
+     
+                     depth_loop1: DO ik=kpk,2,-1 
+                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
+                            
+                           zweig(iin,ijn,1) = &  
+                              & zweig1(iin,ijn,1) * & 
+                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
+                              &  - prodatqc%var(1)%vdep(iend)),0._wp) 
+                            
+                           EXIT depth_loop1 
+
+                        ENDIF 
+                     ENDDO depth_loop1 
+     
+                  ENDDO 
+               ENDDO 
+   
+               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
+                  &              prodatqc%var(1)%vmod(iend:iend) ) 
+
+                  ! Set QC flag for any observations found below the bottom
+                  ! needed as the check here is more strict than that in obs_prep
+               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
+ 
+            ENDDO 
+ 
+            DEALLOCATE(interp_corner,iv_indic) 
+          
+         ENDIF 
+
+         ! For the second variable
          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
 
@@ -410 +502 @@
 
                IF ( idayend == 0 )  THEN
-
                   ! Daily averaged data
-                  CALL obs_int_h2d( kpk, kpk,      &
-                     &              zweig2, zinm2(:,:,:,iobs), zobsk )
-
-               ENDIF
-
-            ELSE
-
-               ! Point data
-               CALL obs_int_h2d( kpk, kpk,      &
-                  &              zweig2, zint2(:,:,:,iobs), zobsk )
-
-            ENDIF
-
-
-            !-------------------------------------------------------------
-            ! Compute vertical second-derivative of the interpolating 
-            ! polynomial at obs points
-            !-------------------------------------------------------------
-
-            IF ( k1dint == 1 ) THEN
-               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
-                  &                  pgdept, zobsmask2 )
-            ENDIF
-
-            !----------------------------------------------------------------
-            !  Vertical interpolation to the observation point
-            !----------------------------------------------------------------
-            ista = prodatqc%npvsta(jobs,2)
-            iend = prodatqc%npvend(jobs,2)
-            CALL obs_int_z1d( kpk, &
-               & prodatqc%var(2)%mvk(ista:iend),&
-               & k1dint, iend - ista + 1, &
-               & prodatqc%var(2)%vdep(ista:iend),&
-               & zobsk, zobs2k, &
-               & prodatqc%var(2)%vmod(ista:iend),&
-               & pgdept, zobsmask2 )
-
-         ENDIF
-
-      END DO
+
+                  ! vertically interpolate all 4 corners 
+                  ista = prodatqc%npvsta(jobs,2) 
+                  iend = prodatqc%npvend(jobs,2) 
+                  inum_obs = iend - ista + 1 
+                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
+
+                  DO iin=1,2 
+                     DO ijn=1,2 
+
+                        IF ( k1dint == 1 ) THEN 
+                           CALL obs_int_z1d_spl( kpk, & 
+                              &     zinm2(iin,ijn,:,iobs), & 
+                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
+                              &     zmask2(iin,ijn,:,iobs)) 
+                        ENDIF 
+       
+                        CALL obs_level_search(kpk, & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
+                           &    iv_indic) 
+
+                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
+                           &    prodatqc%var(2)%vdep(ista:iend), & 
+                           &    zinm2(iin,ijn,:,iobs), & 
+                           &    zobs2k, interp_corner(iin,ijn,:), & 
+                           &    zgdept(iin,ijn,:,iobs), & 
+                           &    zmask2(iin,ijn,:,iobs)) 
+       
+                     ENDDO 
+                  ENDDO 
+
+               ENDIF !idayend
+
+            ELSE   
+
+               ! Point data 
+     
+               ! vertically interpolate all 4 corners 
+               ista = prodatqc%npvsta(jobs,2) 
+               iend = prodatqc%npvend(jobs,2) 
+               inum_obs = iend - ista + 1 
+               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
+               DO iin=1,2  
+                  DO ijn=1,2 
+                    
+                     IF ( k1dint == 1 ) THEN 
+                        CALL obs_int_z1d_spl( kpk, & 
+                           &    zint2(iin,ijn,:,iobs),& 
+                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
+                           &    zmask2(iin,ijn,:,iobs)) 
+  
+                     ENDIF 
+       
+                     CALL obs_level_search(kpk, & 
+                         &        zgdept(iin,ijn,:,iobs),& 
+                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
+                         &        iv_indic) 
+
+                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
+                         &          prodatqc%var(2)%vdep(ista:iend),     & 
+                         &          zint2(iin,ijn,:,iobs),            & 
+                         &          zobs2k,interp_corner(iin,ijn,:), & 
+                         &          zgdept(iin,ijn,:,iobs),         & 
+                         &          zmask2(iin,ijn,:,iobs) )      
+         
+                  ENDDO 
+               ENDDO 
+             
+            ENDIF 
+
+            !------------------------------------------------------------- 
+            ! Compute the horizontal interpolation for every profile level 
+            !------------------------------------------------------------- 
+             
+            DO ikn=1,inum_obs 
+               iend=ista+ikn-1
+                  
+               zweig(:,:,1) = 0._wp 
+   
+               ! This code forces the horizontal weights to be  
+               ! zero IF the observation is below the bottom of the  
+               ! corners of the interpolation nodes, Or if it is in  
+               ! the mask. This is important for observations near  
+               ! steep bathymetry 
+               DO iin=1,2 
+                  DO ijn=1,2 
+     
+                     depth_loop2: DO ik=kpk,2,-1 
+                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
+                            
+                           zweig(iin,ijn,1) = &  
+                              & zweig2(iin,ijn,1) * & 
+                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
+                              &  - prodatqc%var(2)%vdep(iend)),0._wp) 
+                            
+                           EXIT depth_loop2 
+
+                        ENDIF 
+
+                     ENDDO depth_loop2 
+     
+                  ENDDO 
+               ENDDO 
+   
+               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
+                  &              prodatqc%var(2)%vmod(iend:iend) ) 
+
+                  ! Set QC flag for any observations found below the bottom
+                  ! needed as the check here is more strict than that in obs_prep
+               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
+ 
+            ENDDO 
+ 
+            DEALLOCATE(interp_corner,iv_indic) 
+          
+         ENDIF 
+
+      ENDDO
 
       ! Deallocate the data for interpolation
@@ -466 +636 @@
          & zmask2, &
          & zint1,  &
-         & zint2   &
+         & zint2,  &
+         & zgdept, &
+         & zgdepw  &
          & )
 
@@ -481 +653 @@
    END SUBROUTINE obs_prof_opt
 
-   SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
-      &                    ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, & 
-      &                    kdailyavtypes ) 
-      !!----------------------------------------------------------------------- 
-      !! 
-      !!                     ***  ROUTINE obs_pro_opt  *** 
-      !! 
-      !! ** Purpose : Compute the model counterpart of profiles 
-      !!              data by interpolating from the model grid to the  
-      !!              observation point. Generalised vertical coordinate version 
-      !! 
-      !! ** Method  : Linearly interpolate to each observation point using  
-      !!              the model values at the corners of the surrounding grid box. 
-      !! 
-      !!          First, model values on the model grid are interpolated vertically to the 
-      !!          Depths of the profile observations.  Two vertical interpolation schemes are 
-      !!          available: 
-      !!          - linear       (k1dint = 0) 
-      !!          - Cubic spline (k1dint = 1)    
-      !! 
-      !! 
-      !!         Secondly the interpolated values are interpolated horizontally to the  
-      !!         obs (lon, lat) point. 
-      !!         Several horizontal interpolation schemes are available: 
-      !!        - distance-weighted (great circle) (k2dint = 0) 
-      !!        - distance-weighted (small angle)  (k2dint = 1) 
-      !!        - bilinear (geographical grid)     (k2dint = 2) 
-      !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
-      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
-      !! 
-      !!    For the cubic spline the 2nd derivative of the interpolating  
-      !!    polynomial is computed before entering the vertical interpolation  
-      !!    routine. 
-      !! 
-      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is 
-      !!    a daily mean model temperature field. So, we first compute 
-      !!    the mean, then interpolate only at the end of the day. 
-      !! 
-      !!    This is the procedure to be used with generalised vertical model  
-      !!    coordinates (ie s-coordinates. It is ~4x slower than the equivalent 
-      !!    horizontal then vertical interpolation algorithm, but can deal with situations 
-      !!    where the model levels are not flat. 
-      !!    ONLY PERFORMED if ln_sco=.TRUE.  
-      !!       
-      !!    Note: the in situ temperature observations must be converted 
-      !!    to potential temperature (the model variable) prior to 
-      !!    assimilation.  
-      !!?????????????????????????????????????????????????????????????? 
-      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 
-      !!?????????????????????????????????????????????????????????????? 
-      !! 
-      !! ** Action  : 
-      !! 
-      !! History : 
-      !!      ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised 
-      !!                           vertical coordinates
-      !!----------------------------------------------------------------------- 
-   
-      !! * Modules used 
-      USE obs_profiles_def   ! Definition of storage space for profile obs. 
-       
-      IMPLICIT NONE 
- 
-      !! * Arguments 
-      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening 
-      INTEGER, INTENT(IN) :: kt        ! Time step 
-      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
-      INTEGER, INTENT(IN) :: kpj 
-      INTEGER, INTENT(IN) :: kpk 
-      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
-                                       !   (kit000-1 = restart time) 
-      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header) 
-      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
-      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
-         & ptn,    &    ! Model temperature field 
-         & psn,    &    ! Model salinity field 
-         & ptmask       ! Land-sea mask 
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
-         & pgdept,  &    ! Model array of depth T levels    
-         & pgdepw       ! Model array of depth W levels 
-      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
-         & kdailyavtypes   ! Types for daily averages 
-      
-      !! * Local declarations 
-      INTEGER ::   ji 
-      INTEGER ::   jj 
-      INTEGER ::   jk 
-      INTEGER ::   iico, ijco 
-      INTEGER ::   jobs 
-      INTEGER ::   inrc 
-      INTEGER ::   ipro 
-      INTEGER ::   idayend 
-      INTEGER ::   ista 
-      INTEGER ::   iend 
-      INTEGER ::   iobs 
-      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes 
-      INTEGER, DIMENSION(imaxavtypes) :: & 
-         & idailyavtypes 
-      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
-         & igrdi, & 
-         & igrdj 
-      INTEGER :: & 
-         & inum_obs
-      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic    
-      REAL(KIND=wp) :: zlam 
-      REAL(KIND=wp) :: zphi 
-      REAL(KIND=wp) :: zdaystp 
-      REAL(KIND=wp), DIMENSION(kpk) :: & 
-         & zobsmask, & 
-         & zobsk,    & 
-         & zobs2k 
-      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
-         & zweig, & 
-         & l_zweig 
-      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
-         & zmask, & 
-         & zintt, & 
-         & zints, & 
-         & zinmt, & 
-         & zgdept,& 
-         & zgdepw,& 
-         & zinms 
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
-         & zglam, & 
-         & zgphi    
-      REAL(KIND=wp), DIMENSION(1) :: zmsk_1       
-      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner       
- 
-      !------------------------------------------------------------------------ 
-      ! Local initialization  
-      !------------------------------------------------------------------------ 
-      ! ... Record and data counters 
-      inrc = kt - kit000 + 2 
-      ipro = prodatqc%npstp(inrc) 
-  
-      ! Daily average types 
-      IF ( PRESENT(kdailyavtypes) ) THEN 
-         idailyavtypes(:) = kdailyavtypes(:) 
-      ELSE 
-         idailyavtypes(:) = -1 
-      ENDIF 
- 
-      ! Initialize daily mean for first time-step 
-      idayend = MOD( kt - kit000 + 1, kdaystp ) 
- 
-      ! Added kt == 0 test to catch restart case  
-      IF ( idayend == 1 .OR. kt == 0) THEN 
-          
-         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 
-         DO jk = 1, jpk 
-            DO jj = 1, jpj 
-               DO ji = 1, jpi 
-                  prodatqc%vdmean(ji,jj,jk,1) = 0.0 
-                  prodatqc%vdmean(ji,jj,jk,2) = 0.0 
-               END DO 
-            END DO 
-         END DO 
-       
-      ENDIF 
-       
-      DO jk = 1, jpk 
-         DO jj = 1, jpj 
-            DO ji = 1, jpi 
-               ! Increment the temperature field for computing daily mean 
-               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
-               &                        + ptn(ji,jj,jk) 
-               ! Increment the salinity field for computing daily mean 
-               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
-               &                        + psn(ji,jj,jk) 
-            END DO 
-         END DO 
-      END DO 
-    
-      ! Compute the daily mean at the end of day 
-      zdaystp = 1.0 / REAL( kdaystp ) 
-      IF ( idayend == 0 ) THEN 
-         DO jk = 1, jpk 
-            DO jj = 1, jpj 
-               DO ji = 1, jpi 
-                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
-                  &                        * zdaystp 
-                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
-                  &                           * zdaystp 
-               END DO 
-            END DO 
-         END DO 
-      ENDIF 
- 
-      ! Get the data for interpolation 
-      ALLOCATE( & 
-         & igrdi(2,2,ipro),      & 
-         & igrdj(2,2,ipro),      & 
-         & zglam(2,2,ipro),      & 
-         & zgphi(2,2,ipro),      & 
-         & zmask(2,2,kpk,ipro),  & 
-         & zintt(2,2,kpk,ipro),  & 
-         & zints(2,2,kpk,ipro),  & 
-         & zgdept(2,2,kpk,ipro), & 
-         & zgdepw(2,2,kpk,ipro)  & 
-         & ) 
- 
-      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
-         iobs = jobs - prodatqc%nprofup 
-         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 
-         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 
-         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 
-         igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 
-         igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 
-         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 
-         igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 
-         igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 
-      END DO 
-     
-      ! Initialise depth arrays
-      zgdept = 0.0
-      zgdepw = 0.0
- 
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam ) 
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi ) 
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask ) 
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn,   zintt ) 
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn,   zints ) 
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), & 
-        &                     zgdept ) 
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), & 
-        &                     zgdepw ) 
- 
-      ! At the end of the day also get interpolated means 
-      IF ( idayend == 0 ) THEN 
- 
-         ALLOCATE( & 
-            & zinmt(2,2,kpk,ipro),  & 
-            & zinms(2,2,kpk,ipro)   & 
-            & ) 
- 
-         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 
-            &                  prodatqc%vdmean(:,:,:,1), zinmt ) 
-         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 
-            &                  prodatqc%vdmean(:,:,:,2), zinms ) 
- 
-      ENDIF 
-       
-      ! Return if no observations to process 
-      ! Has to be done after comm commands to ensure processors 
-      ! stay in sync 
-      IF ( ipro == 0 ) RETURN 
- 
-      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
-    
-         iobs = jobs - prodatqc%nprofup 
-    
-         IF ( kt /= prodatqc%mstp(jobs) ) THEN 
-             
-            IF(lwp) THEN 
-               WRITE(numout,*) 
-               WRITE(numout,*) ' E R R O R : Observation',              & 
-                  &            ' time step is not consistent with the', & 
-                  &            ' model time step' 
-               WRITE(numout,*) ' =========' 
-               WRITE(numout,*) 
-               WRITE(numout,*) ' Record  = ', jobs,                    & 
-                  &            ' kt      = ', kt,                      & 
-                  &            ' mstp    = ', prodatqc%mstp(jobs), & 
-                  &            ' ntyp    = ', prodatqc%ntyp(jobs) 
-            ENDIF 
-            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 
-         ENDIF 
-          
-         zlam = prodatqc%rlam(jobs) 
-         zphi = prodatqc%rphi(jobs) 
-          
-         ! Horizontal weights 
-         ! Only calculated once, for both T and S. 
-         ! Masked values are calculated later.  
- 
-         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 
-            & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 
- 
-            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
-               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
-               &                   zmask(:,:,1,iobs), zweig, zmsk_1 ) 
- 
-         ENDIF 
-         
-         ! IF zmsk_1 = 0; then ob is on land 
-         IF (zmsk_1(1) < 0.1) THEN 
-            WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 
-   
-         ELSE  
-             
-            ! Temperature 
-             
-            IF ( prodatqc%npvend(jobs,1) > 0 ) THEN  
-    
-               zobsk(:) = obfillflt 
-    
-               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
-    
-                  IF ( idayend == 0 )  THEN 
-                   
-                     ! Daily averaged moored buoy (MRB) data 
-                   
-                     ! vertically interpolate all 4 corners 
-                     ista = prodatqc%npvsta(jobs,1) 
-                     iend = prodatqc%npvend(jobs,1) 
-                     inum_obs = iend - ista + 1 
-                     ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
-      
-                     DO iin=1,2 
-                        DO ijn=1,2 
-                                       
-                                       
-           
-                           IF ( k1dint == 1 ) THEN 
-                              CALL obs_int_z1d_spl( kpk, & 
-                                 &     zinmt(iin,ijn,:,iobs), & 
-                                 &     zobs2k, zgdept(iin,ijn,:,iobs), & 
-                                 &     zmask(iin,ijn,:,iobs)) 
-                           ENDIF 
-       
-                           CALL obs_level_search(kpk, & 
-                              &    zgdept(iin,ijn,:,iobs), & 
-                              &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
-                              &    iv_indic) 
-                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
-                              &    prodatqc%var(1)%vdep(ista:iend), & 
-                              &    zinmt(iin,ijn,:,iobs), & 
-                              &    zobs2k, interp_corner(iin,ijn,:), & 
-                              &    zgdept(iin,ijn,:,iobs), & 
-                              &    zmask(iin,ijn,:,iobs)) 
-       
-                        ENDDO 
-                     ENDDO 
-                   
-                   
-                  ELSE 
-                
-                     CALL ctl_stop( ' A nonzero' //     & 
-                        &           ' number of profile T BUOY data should' // & 
-                        &           ' only occur at the end of a given day' ) 
-    
-                  ENDIF 
-         
-               ELSE  
-                
-                  ! Point data 
-     
-                  ! vertically interpolate all 4 corners 
-                  ista = prodatqc%npvsta(jobs,1) 
-                  iend = prodatqc%npvend(jobs,1) 
-                  inum_obs = iend - ista + 1 
-                  ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
-                  DO iin=1,2  
-                     DO ijn=1,2 
-                                    
-                                    
-                        IF ( k1dint == 1 ) THEN 
-                           CALL obs_int_z1d_spl( kpk, & 
-                              &    zintt(iin,ijn,:,iobs),& 
-                              &    zobs2k, zgdept(iin,ijn,:,iobs), & 
-                              &    zmask(iin,ijn,:,iobs)) 
-  
-                        ENDIF 
-       
-                        CALL obs_level_search(kpk, & 
-                            &        zgdept(iin,ijn,:,iobs),& 
-                            &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
-                            &         iv_indic) 
-                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
-                            &          prodatqc%var(1)%vdep(ista:iend),     & 
-                            &          zintt(iin,ijn,:,iobs),            & 
-                            &          zobs2k,interp_corner(iin,ijn,:), & 
-                            &          zgdept(iin,ijn,:,iobs),         & 
-                            &          zmask(iin,ijn,:,iobs) )      
-         
-                     ENDDO 
-                  ENDDO 
-             
-               ENDIF 
-       
-               !------------------------------------------------------------- 
-               ! Compute the horizontal interpolation for every profile level 
-               !------------------------------------------------------------- 
-             
-               DO ikn=1,inum_obs 
-                  iend=ista+ikn-1 
-
-                  l_zweig(:,:,1) = 0._wp 
-
-                  ! This code forces the horizontal weights to be  
-                  ! zero IF the observation is below the bottom of the  
-                  ! corners of the interpolation nodes, Or if it is in  
-                  ! the mask. This is important for observations are near  
-                  ! steep bathymetry 
-                  DO iin=1,2 
-                     DO ijn=1,2 
-     
-                        depth_loop1: DO ik=kpk,2,-1 
-                           IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
-                            
-                              l_zweig(iin,ijn,1) = &  
-                                 & zweig(iin,ijn,1) * & 
-                                 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
-                                 &  - prodatqc%var(1)%vdep(iend)),0._wp) 
-                            
-                              EXIT depth_loop1 
-                           ENDIF 
-                        ENDDO depth_loop1 
-     
-                     ENDDO 
-                  ENDDO 
-   
-                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 
-                  &          prodatqc%var(1)%vmod(iend:iend) ) 
- 
-               ENDDO 
- 
- 
-               DEALLOCATE(interp_corner,iv_indic) 
-          
-            ENDIF 
-       
- 
-            ! Salinity  
-          
-            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN  
-    
-               zobsk(:) = obfillflt 
-    
-               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
-    
-                  IF ( idayend == 0 )  THEN 
-                   
-                     ! Daily averaged moored buoy (MRB) data 
-                   
-                     ! vertically interpolate all 4 corners 
-                     ista = prodatqc%npvsta(iobs,2) 
-                     iend = prodatqc%npvend(iobs,2) 
-                     inum_obs = iend - ista + 1 
-                     ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
-      
-                     DO iin=1,2 
-                        DO ijn=1,2 
-                                       
-                                       
-           
-                           IF ( k1dint == 1 ) THEN 
-                              CALL obs_int_z1d_spl( kpk, & 
-                                 &     zinms(iin,ijn,:,iobs), & 
-                                 &     zobs2k, zgdept(iin,ijn,:,iobs), & 
-                                 &     zmask(iin,ijn,:,iobs)) 
-                           ENDIF 
-       
-                           CALL obs_level_search(kpk, & 
-                              &    zgdept(iin,ijn,:,iobs), & 
-                              &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
-                              &    iv_indic) 
-                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
-                              &    prodatqc%var(2)%vdep(ista:iend), & 
-                              &    zinms(iin,ijn,:,iobs), & 
-                              &    zobs2k, interp_corner(iin,ijn,:), & 
-                              &    zgdept(iin,ijn,:,iobs), & 
-                              &    zmask(iin,ijn,:,iobs)) 
-       
-                        ENDDO 
-                     ENDDO 
-                   
-                   
-                  ELSE 
-                
-                     CALL ctl_stop( ' A nonzero' //     & 
-                        &           ' number of profile T BUOY data should' // & 
-                        &           ' only occur at the end of a given day' ) 
-    
-                  ENDIF 
-         
-               ELSE  
-                
-                  ! Point data 
-     
-                  ! vertically interpolate all 4 corners 
-                  ista = prodatqc%npvsta(jobs,2) 
-                  iend = prodatqc%npvend(jobs,2) 
-                  inum_obs = iend - ista + 1 
-                  ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
-                   
-                  DO iin=1,2     
-                     DO ijn=1,2  
-                                 
-                                 
-                        IF ( k1dint == 1 ) THEN 
-                           CALL obs_int_z1d_spl( kpk, & 
-                              &    zints(iin,ijn,:,iobs),& 
-                              &    zobs2k, zgdept(iin,ijn,:,iobs), & 
-                              &    zmask(iin,ijn,:,iobs)) 
-  
-                        ENDIF 
-       
-                        CALL obs_level_search(kpk, & 
-                           &        zgdept(iin,ijn,:,iobs),& 
-                           &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
-                           &         iv_indic) 
-                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,  & 
-                           &          prodatqc%var(2)%vdep(ista:iend),     & 
-                           &          zints(iin,ijn,:,iobs),               & 
-                           &          zobs2k,interp_corner(iin,ijn,:),     & 
-                           &          zgdept(iin,ijn,:,iobs),              & 
-                           &          zmask(iin,ijn,:,iobs) )      
-         
-                     ENDDO 
-                  ENDDO 
-             
-               ENDIF 
-       
-               !------------------------------------------------------------- 
-               ! Compute the horizontal interpolation for every profile level 
-               !------------------------------------------------------------- 
-             
-               DO ikn=1,inum_obs 
-                  iend=ista+ikn-1 
-
-                  l_zweig(:,:,1) = 0._wp
-   
-                  ! This code forces the horizontal weights to be  
-                  ! zero IF the observation is below the bottom of the  
-                  ! corners of the interpolation nodes, Or if it is in  
-                  ! the mask. This is important for observations are near  
-                  ! steep bathymetry 
-                  DO iin=1,2 
-                     DO ijn=1,2 
-     
-                        depth_loop2: DO ik=kpk,2,-1 
-                           IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
-                            
-                              l_zweig(iin,ijn,1) = &  
-                                 &  zweig(iin,ijn,1) * & 
-                                 &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
-                                 &  - prodatqc%var(2)%vdep(iend)),0._wp) 
-                            
-                              EXIT depth_loop2 
-                           ENDIF 
-                        ENDDO depth_loop2 
-     
-                     ENDDO 
-                  ENDDO 
-   
-                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 
-                  &          prodatqc%var(2)%vmod(iend:iend) ) 
- 
-               ENDDO 
- 
- 
-               DEALLOCATE(interp_corner,iv_indic) 
-          
-            ENDIF 
-          
-         ENDIF 
-       
-      END DO 
-     
-      ! Deallocate the data for interpolation 
-      DEALLOCATE( & 
-         & igrdi, & 
-         & igrdj, & 
-         & zglam, & 
-         & zgphi, & 
-         & zmask, & 
-         & zintt, & 
-         & zints, & 
-         & zgdept,&
-         & zgdepw &
-         & ) 
-      ! At the end of the day also get interpolated means 
-      IF ( idayend == 0 ) THEN 
-         DEALLOCATE( & 
-            & zinmt,  & 
-            & zinms   & 
-            & ) 
-      ENDIF 
-    
-      prodatqc%nprofup = prodatqc%nprofup + ipro  
-       
-   END SUBROUTINE obs_pro_sco_opt 
- 
-   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         &
-      &                    kit000, kdaystp, psurf, psurfmask, &
-      &                    k2dint, ldnightav )
+   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
+      &                     kit000, kdaystp, psurf, psurfmask,   &
+      &                     k2dint, ldnightav, plamscl, pphiscl, &
+      &                     lindegrees )
 
       !!-----------------------------------------------------------------------
@@ -1090 +678 @@
       !!        - polynomial (quadrilateral grid)  (k2dint = 4)
       !!
+      !!    Two horizontal averaging schemes are also available:
+      !!        - weighted radial footprint        (k2dint = 5)
+      !!        - weighted rectangular footprint   (k2dint = 6)
+      !!
       !!
       !! ** Action  :
@@ -1096 +688 @@
       !!      ! 07-03 (A. Weaver)
       !!      ! 15-02 (M. Martin) Combined routine for surface types
+      !!      ! 17-03 (M. Martin) Added horizontal averaging options
       !!-----------------------------------------------------------------------
 
@@ -1117 +710 @@
          & psurfmask                   ! Land-sea mask
       LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions
+         & pphiscl                     ! This is the full width (rather than half-width)
+      LOGICAL, INTENT(IN) :: &
+         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
 
       !! * Local declarations
@@ -1125 +723 @@
       INTEGER :: isurf
       INTEGER :: iobs
+      INTEGER :: imaxifp, imaxjfp
+      INTEGER :: imodi, imodj
       INTEGER :: idayend
       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi, &
-         & igrdj
+         & igrdi,   &
+         & igrdj,   &
+         & igrdip1, &
+         & igrdjp1
       INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
          & icount_night,      &
@@ -1136 +738 @@
       REAL(wp), DIMENSION(1) :: zext, zobsmask
       REAL(wp) :: zdaystp
-      REAL(wp), DIMENSION(2,2,1) :: &
-         & zweig
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
+         & zweig,  &
          & zmask,  &
          & zsurf,  &
          & zsurfm, &
+         & zsurftmp, &
          & zglam,  &
-         & zgphi
+         & zgphi,  &
+         & zglamf, &
+         & zgphif
+
       REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
          & zintmp,  &
@@ -1155 +760 @@
       inrc = kt - kit000 + 2
       isurf = surfdataqc%nsstp(inrc)
+
+      ! Work out the maximum footprint size for the 
+      ! interpolation/averaging in model grid-points - has to be even.
+
+      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
+
 
       IF ( ldnightav ) THEN
@@ -1221 +832 @@
 
       ALLOCATE( &
-         & igrdi(2,2,isurf), &
-         & igrdj(2,2,isurf), &
-         & zglam(2,2,isurf), &
-         & zgphi(2,2,isurf), &
-         & zmask(2,2,isurf), &
-         & zsurf(2,2,isurf)  &
+         & zweig(imaxifp,imaxjfp,1),      &
+         & igrdi(imaxifp,imaxjfp,isurf), &
+         & igrdj(imaxifp,imaxjfp,isurf), &
+         & zglam(imaxifp,imaxjfp,isurf), &
+         & zgphi(imaxifp,imaxjfp,isurf), &
+         & zmask(imaxifp,imaxjfp,isurf), &
+         & zsurf(imaxifp,imaxjfp,isurf), &
+         & zsurftmp(imaxifp,imaxjfp,isurf),  &
+         & zglamf(imaxifp+1,imaxjfp+1,isurf), &
+         & zgphif(imaxifp+1,imaxjfp+1,isurf), &
+         & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
+         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
          & )
 
       DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
          iobs = jobs - surfdataqc%nsurfup
-         igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1
-         igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1
-         igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1
-         igrdj(1,2,iobs) = surfdataqc%mj(jobs)
-         igrdi(2,1,iobs) = surfdataqc%mi(jobs)
-         igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1
-         igrdi(2,2,iobs) = surfdataqc%mi(jobs)
-         igrdj(2,2,iobs) = surfdataqc%mj(jobs)
+         DO ji = 0, imaxifp
+            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
+            
+            !Deal with wrap around in longitude
+            IF ( imodi < 1      ) imodi = imodi + jpiglo
+            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
+            
+            DO jj = 0, imaxjfp
+               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
+               !If model values are out of the domain to the north/south then
+               !set them to be the edge of the domain
+               IF ( imodj < 1      ) imodj = 1
+               IF ( imodj > jpjglo ) imodj = jpjglo
+
+               igrdip1(ji+1,jj+1,iobs) = imodi
+               igrdjp1(ji+1,jj+1,iobs) = imodj
+               
+               IF ( ji >= 1 .AND. jj >= 1 ) THEN
+                  igrdi(ji,jj,iobs) = imodi
+                  igrdj(ji,jj,iobs) = imodj
+               ENDIF
+               
+            END DO
+         END DO
       END DO
 
-      CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, glamt, zglam )
-      CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, gphit, zgphi )
-      CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, psurfmask, zmask )
-      CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
+      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, psurf, zsurf )
+      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
+         &                  igrdip1, igrdjp1, glamf, zglamf )
+      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
+         &                  igrdip1, igrdjp1, gphif, zgphif )
 
       ! At the end of the day get interpolated means
-      IF (ldnightav ) THEN
-         IF ( idayend == 0 ) THEN
-
-            ALLOCATE( &
-               & zsurfm(2,2,isurf)  &
-               & )
-
-            CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, &
-               &               surfdataqc%vdmean(:,:), zsurfm )
-
-         ENDIF
+      IF ( idayend == 0 .AND. ldnightav ) THEN
+
+         ALLOCATE( &
+            & zsurfm(imaxifp,imaxjfp,isurf)  &
+            & )
+
+         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
+            &               surfdataqc%vdmean(:,:), zsurfm )
+
       ENDIF
 
@@ -1290 +925 @@
          zphi = surfdataqc%rphi(jobs)
 
-         ! Get weights to interpolate the model value to the observation point
-         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
-            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
-            &                   zmask(:,:,iobs), zweig, zobsmask )
-
-         ! Interpolate the model field to the observation point
          IF ( ldnightav .AND. idayend == 0 ) THEN
             ! Night-time averaged data
-            CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext )
+            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
          ELSE
-            CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext )
+            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
+         ENDIF
+
+         IF ( k2dint <= 4 ) THEN
+
+            ! Get weights to interpolate the model value to the observation point
+            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
+               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
+               &                   zmask(:,:,iobs), zweig, zobsmask )
+
+            ! Interpolate the model value to the observation point 
+            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
+
+         ELSE
+
+            ! Get weights to average the model SLA to the observation footprint
+            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
+               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
+               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
+               &                   zmask(:,:,iobs), plamscl, pphiscl, &
+               &                   lindegrees, zweig, zobsmask )
+
+            ! Average the model SST to the observation footprint
+            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
+               &              zweig, zsurftmp(:,:,iobs),  zext )
+
          ENDIF
 
@@ -1310 +964 @@
             surfdataqc%rmod(jobs,1) = zext(1)
          ENDIF
+         
+         IF ( zext(1) == obfillflt ) THEN
+            ! If the observation value is a fill value, set QC flag to bad
+            surfdataqc%nqc(jobs) = 4
+         ENDIF
 
       END DO
@@ -1315 +974 @@
       ! Deallocate the data for interpolation
       DEALLOCATE( &
+         & zweig, &
          & igrdi, &
          & igrdj, &
@@ -1320 +980 @@
          & zgphi, &
          & zmask, &
-         & zsurf  &
+         & zsurf, &
+         & zsurftmp, &
+         & zglamf, &
+         & zgphif, &
+         & igrdip1,&
+         & igrdjp1 &
          & )
 
       ! At the end of the day also deallocate night-time mean array
-      IF ( ldnightav ) THEN
-         IF ( idayend == 0 ) THEN
-            DEALLOCATE( &
-               & zsurfm  &
-               & )
-         ENDIF
+      IF ( idayend == 0 .AND. ldnightav ) THEN
+         DEALLOCATE( &
+            & zsurfm  &
+            & )
       ENDIF
 

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

@@ -23 +23 @@
    USE obs_oper           ! Observation operators
    USE lib_mpp, ONLY :   ctl_warn, ctl_stop
+   USE bdy_oce, ONLY : &        ! Boundary information
+      idx_bdy, nb_bdy, ln_bdy
 
    IMPLICIT NONE
@@ -40 +42 @@
 CONTAINS
 
-   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea )
+   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, &
+                            kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_pre_sla  ***
@@ -57 +60 @@
       !!        !  2015-02  (M. Martin) Combined routine for surface types.
       !!----------------------------------------------------------------------
+      !! * Modules used
+      USE domstp              ! Domain: set the time-step
       USE par_oce             ! Ocean parameters
       USE dom_oce, ONLY       :   glamt, gphit, tmask, nproc   ! Geographical information
@@ -63 +68 @@
       TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening
       LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land
-      !
+      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary
+      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
+      !! * Local declarations
+      INTEGER :: iqc_cutoff = 255   ! cut off for QC value
       INTEGER :: iyea0        ! Initial date
       INTEGER :: imon0        !  - (year, month, day, hour, minute)
@@ -76 +84 @@
       INTEGER :: inlasobs     !  - close to land
       INTEGER :: igrdobs      !  - fail the grid search
+      INTEGER :: ibdysobs     !  - close to open boundary
                               ! Global counters for observations that
       INTEGER :: iotdobsmpp     !  - outside time domain
@@ -82 +91 @@
       INTEGER :: inlasobsmpp    !  - close to land
       INTEGER :: igrdobsmpp     !  - fail the grid search
-      LOGICAL, DIMENSION(:), ALLOCATABLE ::   llvalid            ! SLA data selection
+      INTEGER :: ibdysobsmpp  !  - close to open boundary
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
+         & llvalid            ! SLA data selection
       INTEGER :: jobs         ! Obs. loop variable
       INTEGER :: jstp         ! Time loop variable
@@ -107 +118 @@
       ilansobs = 0
       inlasobs = 0
+      ibdysobs = 0 
+
+      ! Set QC cutoff to optional value if provided
+      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
 
       ! -----------------------------------------------------------------------
@@ -140 +155 @@
          &                 tmask(:,:,1), surfdata%nqc,  &
          &                 iosdsobs,     ilansobs,     &
-         &                 inlasobs,     ld_nea        )
+         &                 inlasobs,     ld_nea,       &
+         &                 ibdysobs,     ld_bound_reject, &
+         &                 iqc_cutoff                     )
 
       CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
       CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
       CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
+      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
 
       ! -----------------------------------------------------------------------
@@ -155 +173 @@
       ALLOCATE( llvalid(surfdata%nsurf) )
       
-      ! We want all data which has qc flags <= 10
-
-      llvalid(:)  = ( surfdata%nqc(:)  <= 10 )
+      ! We want all data which has qc flags <= iqc_cutoff
+
+      llvalid(:)  = ( surfdata%nqc(:)  <= iqc_cutoff )
 
       ! The actual copying
@@ -190 +208 @@
                &            inlasobsmpp
          ENDIF
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', &
+            &            ibdysobsmpp  
          WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', &
             &            surfdataqc%nsurfmpp
@@ -225 +245 @@
       &                     kpi, kpj, kpk, &
       &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  &
-      &                     ld_nea, kdailyavtypes )
+      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff )
 
 !!----------------------------------------------------------------------
@@ -241 +261 @@
       !!
       !!----------------------------------------------------------------------
-      USE par_oce           ! Ocean parameters
-      USE dom_oce, ONLY :   gdept_1d, nproc   ! Geographical information
+      !! * Modules used
+      USE domstp              ! Domain: set the time-step
+      USE par_oce             ! Ocean parameters
+      USE dom_oce, ONLY : &   ! Geographical information
+         & gdept_1d,             &
+         & nproc
 
       !! * Arguments
@@ -250 +274 @@
       LOGICAL, INTENT(IN) :: ld_var2
       LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land
+      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary
       INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
@@ -261 +286 @@
          & pgphi1, &
          & pgphi2
+      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
 
       !! * Local declarations
+      INTEGER :: iqc_cutoff = 255   ! cut off for QC value
       INTEGER :: iyea0        ! Initial date
       INTEGER :: imon0        !  - (year, month, day, hour, minute)
@@ -277 +304 @@
       INTEGER :: inlav1obs    !  - close to land (variable 1)
       INTEGER :: inlav2obs    !  - close to land (variable 2)
+      INTEGER :: ibdyv1obs    !  - boundary (variable 1) 
+      INTEGER :: ibdyv2obs    !  - boundary (variable 2)      
       INTEGER :: igrdobs      !  - fail the grid search
       INTEGER :: iuvchku      !  - reject u if v rejected and vice versa
@@ -288 +317 @@
       INTEGER :: inlav1obsmpp !  - close to land (variable 1)
       INTEGER :: inlav2obsmpp !  - close to land (variable 2)
+      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1) 
+      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)      
       INTEGER :: igrdobsmpp   !  - fail the grid search
       INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa
@@ -322 +353 @@
       inlav1obs = 0
       inlav2obs = 0
-      iuvchku  = 0
-      iuvchkv = 0
+      ibdyv1obs = 0
+      ibdyv2obs = 0
+      iuvchku   = 0
+      iuvchkv   = 0
+
+
+      ! Set QC cutoff to optional value if provided
+      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
 
       ! -----------------------------------------------------------------------
@@ -335 +372 @@
             &              profdata%nday,    profdata%nhou, profdata%nmin, &
             &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
-            &              iotdobs, kdailyavtypes = kdailyavtypes )
+            &              iotdobs, kdailyavtypes = kdailyavtypes,         &
+            &              kqc_cutoff = iqc_cutoff )
       ELSE
          CALL obs_coo_tim_prof( icycle, &
@@ -342 +380 @@
             &              profdata%nday,    profdata%nhou, profdata%nmin, &
             &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
-            &              iotdobs )
+            &              iotdobs,          kqc_cutoff = iqc_cutoff )
       ENDIF
 
@@ -359 +397 @@
 
       ! -----------------------------------------------------------------------
-      ! Reject all observations for profiles with nqc > 10
-      ! -----------------------------------------------------------------------
-
-      CALL obs_pro_rej( profdata )
+      ! Reject all observations for profiles with nqc > iqc_cutoff
+      ! -----------------------------------------------------------------------
+
+      CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff )
 
       ! -----------------------------------------------------------------------
@@ -381 +419 @@
          &                 gdept_1d,              zmask1,               &
          &                 profdata%nqc,          profdata%var(1)%nvqc, &
-         &                 iosdv1obs,              ilanv1obs,           &
-         &                 inlav1obs,              ld_nea                )
+         &                 iosdv1obs,             ilanv1obs,            &
+         &                 inlav1obs,             ld_nea,               &
+         &                 ibdyv1obs,             ld_bound_reject,      &
+         &                 iqc_cutoff       )
 
       CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp )
       CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp )
       CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp )
+      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp )
 
       ! Variable 2
@@ -400 +441 @@
          &                 gdept_1d,              zmask2,               &
          &                 profdata%nqc,          profdata%var(2)%nvqc, &
-         &                 iosdv2obs,              ilanv2obs,           &
-         &                 inlav2obs,              ld_nea                )
+         &                 iosdv2obs,             ilanv2obs,            &
+         &                 inlav2obs,             ld_nea,               &
+         &                 ibdyv2obs,             ld_bound_reject,      &
+         &                 iqc_cutoff       )
 
       CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp )
       CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp )
       CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp )
+      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp )
 
       ! -----------------------------------------------------------------------
@@ -412 +456 @@
 
       IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
-         CALL obs_uv_rej( profdata, iuvchku, iuvchkv )
+         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff )
          CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
          CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
@@ -429 +473 @@
       END DO
 
-      ! We want all data which has qc flags = 0
-
-      llvalid%luse(:) = ( profdata%nqc(:)  <= 10 )
+      ! We want all data which has qc flags <= iqc_cutoff
+
+      llvalid%luse(:) = ( profdata%nqc(:)  <= iqc_cutoff )
       DO jvar = 1,profdata%nvar
-         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
+         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff )
       END DO
 
@@ -475 +519 @@
                &            iuvchku
          ENDIF
+         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',&
+               &            ibdyv1obsmpp
          WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', &
             &            prodatqc%nvprotmpp(1)
@@ -492 +538 @@
                &            iuvchkv
          ENDIF
+         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',&
+               &            ibdyv2obsmpp
          WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', &
             &            prodatqc%nvprotmpp(2)
@@ -644 +692 @@
             &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
             kobsstp(jobs) = -1
-            kobsqc(jobs)  = kobsqc(jobs) + 11
+            kobsqc(jobs)  = IBSET(kobsqc(jobs),13)
             kotdobs       = kotdobs + 1
             CYCLE
@@ -695 +743 @@
          IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
             & .OR.( kobsstp(jobs) > nitend ) ) THEN
-            kobsqc(jobs) = kobsqc(jobs) + 12
+            kobsqc(jobs) = IBSET(kobsqc(jobs),13)
             kotdobs = kotdobs + 1
             CYCLE
@@ -739 +787 @@
       &                    kobsno,                                        &
       &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
-      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes )
+      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, &
+      &                    kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_tim ***
@@ -783 +832 @@
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes    ! Types for daily averages
+      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
+
       !! * Local declarations
       INTEGER :: jobs
+      INTEGER :: iqc_cutoff=255
 
       !-----------------------------------------------------------------------
@@ -803 +855 @@
          DO jobs = 1, kobsno
             
-            IF ( kobsqc(jobs) <= 10 ) THEN
+            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
                
                IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
                   & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
-                  kobsqc(jobs) = kobsqc(jobs) + 14
+                  kobsqc(jobs) = IBSET(kobsqc(jobs),13)
                   kotdobs      = kotdobs + 1
                   CYCLE
@@ -850 +902 @@
       DO jobs = 1, kobsno
          IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
-            kobsqc(jobs) = kobsqc(jobs) + 18
+            kobsqc(jobs) = IBSET(kobsqc(jobs),12)
             kgrdobs = kgrdobs + 1
          ENDIF
@@ -861 +913 @@
       &                       plam,   pphi,    pmask,            &
       &                       kobsqc, kosdobs, klanobs,          &
-      &                       knlaobs,ld_nea                     )
+      &                       knlaobs,ld_nea,                    &
+      &                       kbdyobs,ld_bound_reject,           &
+      &                       kqc_cutoff                         )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_spc_2d  ***
@@ -894 +948 @@
       INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
          & kobsqc             ! Observation quality control
-      INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain
-      INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell
-      INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land
-      LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land
+      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain
+      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell
+      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land
+      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary
+      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land
+      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary 
+      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value
+
       !! * Local declarations
       REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
          & zgmsk              ! Grid mask
+
+      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
+         & zbmsk              ! Boundary mask
+      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
       REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
          & zglam, &           ! Model longitude at grid points
@@ -917 +979 @@
          ! For invalid points use 2,2
 
-         IF ( kobsqc(jobs) >= 10 ) THEN
+         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
 
             igrdi(1,1,jobs) = 1
@@ -942 +1004 @@
 
       END DO
+
+      IF (ln_bdy) THEN
+        ! Create a mask grid points in boundary rim
+        IF (ld_bound_reject) THEN
+           zbdymask(:,:) = 1.0_wp
+           DO ji = 1, nb_bdy
+              DO jj = 1, idx_bdy(ji)%nblen(1)
+                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
+              ENDDO
+           ENDDO
+
+           CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
+        ENDIF
+      ENDIF
+
       
       CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk )
@@ -950 +1027 @@
 
          ! Skip bad observations
-         IF ( kobsqc(jobs) >= 10 ) CYCLE
+         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
 
          ! Flag if the observation falls outside the model spatial domain
@@ -957 +1034 @@
             &  .OR. ( pobsphi(jobs) <  -90. ) &
             &  .OR. ( pobsphi(jobs) >   90. ) ) THEN
-            kobsqc(jobs) = kobsqc(jobs) + 11
+            kobsqc(jobs) = IBSET(kobsqc(jobs),11)
             kosdobs = kosdobs + 1
             CYCLE
@@ -964 +1041 @@
          ! Flag if the observation falls with a model land cell
          IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
-            kobsqc(jobs) = kobsqc(jobs)  + 12
+            kobsqc(jobs) = IBSET(kobsqc(jobs),10)
             klanobs = klanobs + 1
             CYCLE
@@ -978 +1055 @@
                IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
                   & .AND. &
-                  & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
-                  & ) THEN
+                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) )  &
+                  & < 1.0e-6_wp ) ) THEN
                   lgridobs = .TRUE.
                   iig = ji
@@ -992 +1069 @@
          IF (lgridobs) THEN
             IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
-               kobsqc(jobs) = kobsqc(jobs) + 12
+               kobsqc(jobs) = IBSET(kobsqc(jobs),10)
                klanobs = klanobs + 1
                CYCLE
@@ -1000 +1077 @@
          ! Flag if the observation falls is close to land
          IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
-            IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14
             knlaobs = knlaobs + 1
-            CYCLE
+            IF (ld_nea) THEN
+               kobsqc(jobs) = IBSET(kobsqc(jobs),9)
+               CYCLE
+            ENDIF
+         ENDIF
+
+         IF (ln_bdy) THEN
+         ! Flag if the observation falls close to the boundary rim
+           IF (ld_bound_reject) THEN
+              IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
+                 kobsqc(jobs) = IBSET(kobsqc(jobs),8)
+                 kbdyobs = kbdyobs + 1
+                 CYCLE
+              ENDIF
+              ! for observations on the grid...
+              IF (lgridobs) THEN
+                 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
+                    kobsqc(jobs) = IBSET(kobsqc(jobs),8)
+                    kbdyobs = kbdyobs + 1
+                    CYCLE
+                 ENDIF
+              ENDIF
+            ENDIF
          ENDIF
             
@@ -1015 +1113 @@
       &                       plam,    pphi,    pdep,    pmask, &
       &                       kpobsqc, kobsqc,  kosdobs,        &
-      &                       klanobs, knlaobs, ld_nea          )
+      &                       klanobs, knlaobs, ld_nea,         &
+      &                       kbdyobs, ld_bound_reject,         &
+      &                       kqc_cutoff                        )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_spc_3d  ***
@@ -1077 +1177 @@
       INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell
       INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land
+      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary
       LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land
+      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary
+      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value
+
       !! * Local declarations
       REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
          & zgmsk              ! Grid mask
+      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
+         & zbmsk              ! Boundary mask
+      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
       REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
          & zgdepw
@@ -1100 +1207 @@
          ! For invalid points use 2,2
 
-         IF ( kpobsqc(jobs) >= 10 ) THEN
+         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
             
             igrdi(1,1,jobs) = 1
@@ -1125 +1232 @@
          
       END DO
+
+      IF (ln_bdy) THEN
+        ! Create a mask grid points in boundary rim
+        IF (ld_bound_reject) THEN           
+           zbdymask(:,:) = 1.0_wp
+           DO ji = 1, nb_bdy
+              DO jj = 1, idx_bdy(ji)%nblen(1)
+                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
+              ENDDO
+           ENDDO
+        ENDIF
+   
+        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
+      ENDIF
       
       CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
-      IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN
-        ! Need to know the bathy depth for each observation for sco
-        CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
+      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
         &                     zgdepw )
-      ENDIF
 
       DO jobs = 1, kprofno
 
          ! Skip bad profiles
-         IF ( kpobsqc(jobs) >= 10 ) CYCLE
+         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
 
          ! Check if this observation is on a grid point
@@ -1149 +1267 @@
                IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
                   & .AND. &
-                  & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
+                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
                   & ) THEN
                   lgridobs = .TRUE.
@@ -1176 +1294 @@
                &  .OR. ( pobsdep(jobsp) < 0.0          )       &
                &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
-               kobsqc(jobsp) = kobsqc(jobsp) + 11
+               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
                kosdobs = kosdobs + 1
                CYCLE
             ENDIF
 
-            ! To check if an observations falls within land there are two cases:
-            ! 1: z-coordibnates, where the check uses the mask
-            ! 2: terrain following (eg s-coordinates), 
-            !    where we use the depth of the bottom cell to mask observations
+            ! To check if an observations falls within land:
              
-            IF( ln_zps .OR. ln_zco ) THEN !(CASE 1)
+            ! Flag if the observation is deeper than the bathymetry
+            ! Or if it is within the mask
+            IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
+               &     .OR. &
+               &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
+               &  == 0.0_wp) ) THEN
+               kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
+               klanobs = klanobs + 1
+               CYCLE
+            ENDIF
                
-               ! Flag if the observation falls with a model land cell
-               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
-                  &  == 0.0_wp ) THEN
-                  kobsqc(jobsp) = kobsqc(jobsp) + 12
-                  klanobs = klanobs + 1
-                  CYCLE
-               ENDIF
-             
-               ! Flag if the observation is close to land
-              IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
-                  &  0.0_wp) THEN
-                  knlaobs = knlaobs + 1
-                  IF (ld_nea) THEN   
-                     kobsqc(jobsp) = kobsqc(jobsp) + 14 
-                  ENDIF 
-               ENDIF
-             
-            ELSE ! Case 2
- 
-               ! Flag if the observation is deeper than the bathymetry
-               ! Or if it is within the mask
-               IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
-                  &     .OR. &
-                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
-                  &  == 0.0_wp) ) THEN
-                  kobsqc(jobsp) = kobsqc(jobsp) + 12
-                  klanobs = klanobs + 1
-                  CYCLE
-               ENDIF
-               
-               ! Flag if the observation is close to land
-               IF ( ll_next_to_land ) THEN
-                  knlaobs = knlaobs + 1
-                  IF (ld_nea) THEN   
-                     kobsqc(jobsp) = kobsqc(jobsp) + 14 
-                  ENDIF 
-               ENDIF
+            ! Flag if the observation is close to land
+            IF ( ll_next_to_land ) THEN
+               knlaobs = knlaobs + 1
+               IF (ld_nea) THEN   
+                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
+               ENDIF 
             ENDIF
             
@@ -1232 +1325 @@
             IF (lgridobs) THEN
                IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
-                  kobsqc(jobsp) = kobsqc(jobsp) + 12
+                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
                   klanobs = klanobs + 1
                   CYCLE
@@ -1250 +1343 @@
             ENDIF
             
+            IF (ln_bdy) THEN
+               ! Flag if the observation falls close to the boundary rim
+               IF (ld_bound_reject) THEN
+                  IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
+                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
+                     kbdyobs = kbdyobs + 1
+                     CYCLE
+                  ENDIF
+                  ! for observations on the grid...
+                  IF (lgridobs) THEN
+                     IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
+                        kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
+                        kbdyobs = kbdyobs + 1
+                        CYCLE
+                     ENDIF
+                  ENDIF
+               ENDIF
+            ENDIF
+            
          END DO
       END DO
@@ -1255 +1367 @@
    END SUBROUTINE obs_coo_spc_3d
 
-   SUBROUTINE obs_pro_rej( profdata )
+   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_pro_rej ***
@@ -1273 +1385 @@
       !! * Arguments
       TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data
+      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value
+
       !! * Local declarations
       INTEGER :: jprof
@@ -1282 +1396 @@
       DO jprof = 1, profdata%nprof
 
-         IF ( profdata%nqc(jprof) > 10 ) THEN
+         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
             
             DO jvar = 1, profdata%nvar
@@ -1290 +1404 @@
                   
                   profdata%var(jvar)%nvqc(jobs) = &
-                     & profdata%var(jvar)%nvqc(jobs) + 26
+                     & IBSET(profdata%var(jvar)%nvqc(jobs),14)
 
                END DO
@@ -1302 +1416 @@
    END SUBROUTINE obs_pro_rej
 
-   SUBROUTINE obs_uv_rej( profdata, knumu, knumv )
+   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_uv_rej ***
@@ -1322 +1436 @@
       INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected
       INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected
+      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
+
       !! * Local declarations
       INTEGER :: jprof
@@ -1341 +1457 @@
          DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
             
-            IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. &
-               & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN
-               profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42
+            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. &
+               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN
+               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
                knumv = knumv + 1
             ENDIF
-            IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. &
-               & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN
-               profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42
+            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. &
+               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN
+               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
                knumu = knumu + 1
             ENDIF

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90

@@ -104 +104 @@
       ! Bookkeeping arrays with sizes equal to number of variables
 
-      CHARACTER(len=6), POINTER, DIMENSION(:) :: &
+      CHARACTER(len=8), POINTER, DIMENSION(:) :: &
          & cvars          !: Variable names
 

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

@@ -87 +87 @@
       CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
       CHARACTER(len=8) :: clrefdate
-      CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars
       INTEGER :: jvar
       INTEGER :: ji
@@ -307 +307 @@
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
-               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
-               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
+               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -325 +325 @@
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
-               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
-               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
+               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -351 +351 @@
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
-               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
-               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
+               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -373 +373 @@
 
             DO ji = 1, inpfiles(jj)%nobs
-               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
-               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
+               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -388 +388 @@
                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                            & CYCLE
-                        IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
-                           & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
+                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
                            ivar1t0 = ivar1t0 + 1
                         ENDIF
@@ -398 +398 @@
                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                            & CYCLE
-                        IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
-                           & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
+                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
+                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
                            ivar2t0 = ivar2t0 + 1
                         ENDIF
@@ -407 +407 @@
                      IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                         & CYCLE
-                     IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
-                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     ldvar1 ) .OR. &
-                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
-                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
+                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                        &    ldvar1 ) .OR. &
+                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
+                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
                         &     ldvar2 ) ) THEN
                         ip3dt = ip3dt + 1
@@ -437 +437 @@
       DO jj = 1, inobf
          DO ji = 1, inpfiles(jj)%nobs
-            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
-            IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-               & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
+            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -452 +452 @@
       DO jj = 1, inobf
          DO ji = 1, inpfiles(jj)%nobs
-            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
-            IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-               & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
+            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -501 +501 @@
          ji = iprofidx(iindx(jk))
 
-         IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
-         IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-            & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
+            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
 
          IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
@@ -518 +518 @@
             IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
 
-            IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
-               & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE
+            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
 
             loop_prof : DO ij = 1, inpfiles(jj)%nlev
@@ -526 +527 @@
                   & CYCLE
 
-               IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
-                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
+               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
 
                   llvalprof = .TRUE. 
@@ -534 +535 @@
                ENDIF
 
-               IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
-                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
+               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
+                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
 
                   llvalprof = .TRUE. 
@@ -615 +616 @@
                   IF (ldsatt) THEN
 
-                     IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
-                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     ldvar1 ) .OR. &
-                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
-                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     ldvar2 ) ) THEN
+                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                        &    ldvar1 ) .OR. &
+                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
+                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                        &   ldvar2 ) ) THEN
                         ip3dt = ip3dt + 1
                      ELSE
@@ -628 +629 @@
                   ENDIF
 
-                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
-                     &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                     &       ldvar1 ) .OR. ldsatt ) THEN
+                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                    &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                    &    ldvar1 ) .OR. ldsatt ) THEN
 
                      IF (ldsatt) THEN
@@ -661 +662 @@
 
                      ! Profile var1 value
-                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
-                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
+                     IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                        & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
                         profdata%var(1)%vobs(ivar1t) = &
                            &                inpfiles(jj)%pob(ij,ji,1)
@@ -692 +693 @@
                   ENDIF
 
-                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
-                     &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
+                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
+                     &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    .AND. &
                      &   ldvar2 ) .OR. ldsatt ) THEN
 
@@ -725 +726 @@
 
                      ! Profile var2 value
-                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
-                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
+                  IF (  ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. &
+                    &   ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    )  ) THEN
                         profdata%var(2)%vobs(ivar2t) = &
                            &                inpfiles(jj)%pob(ij,ji,2)

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

@@ -77 +77 @@
       CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
       CHARACTER(len=8) :: clrefdate
-      CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars
       INTEGER :: ji
       INTEGER :: jj
@@ -172 +172 @@
 
             !------------------------------------------------------------------
-            !  Read the profile file into inpfiles
+            !  Read the surface file into inpfiles
             !------------------------------------------------------------------
             CALL init_obfbdata( inpfiles(jj) )
@@ -296 +296 @@
                   ENDIF
                   llvalprof = .FALSE.
-                  IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
-                     & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
+                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
                      iobs = iobs + 1
                   ENDIF
@@ -370 +369 @@
             ! Set observation information
 
-            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
-               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
+            IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
 
                iobs = iobs + 1

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90

@@ -154 +154 @@
 
 ! mark any masked data with a QC flag
-         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = 11
+         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15)
 
          END DO

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_sstbias.F90

@@ -1 +1 @@
 MODULE obs_sstbias
    !!======================================================================
-   !!                       ***  MODULE obs_readsstbias  ***
-   !! Observation diagnostics: Read the bias for SLA data
+   !!                       ***  MODULE obs_sstbias  ***
+   !! Observation diagnostics: Read the bias for SST data
    !!======================================================================
    !!----------------------------------------------------------------------
-   !!   obs_rea_sstbias : Driver for reading altimeter bias
+   !!   obs_app_sstbias : Driver for reading and applying the SST bias
    !!----------------------------------------------------------------------
    !! * Modules used   
@@ -139 +139 @@
                cl_bias_files(jtype) )
             ! Get the SST bias data
-            CALL iom_get( numsstbias, jpdom_data, 'sst', z_sstbias_2d(:,:), 1 )
+            CALL iom_get( numsstbias, jpdom_data, 'tn', z_sstbias_2d(:,:), 1 )
             z_sstbias(:,:,jtype) = z_sstbias_2d(:,:)       
             ! Close the file
@@ -190 +190 @@
                igrdj_tmp(:,:,jt) = igrdj(:,:,jobs)
                zglam_tmp(:,:,jt) = zglam(:,:,jobs)
-               zgphi_tmp(:,:,jt) = zgphi(:,:,jobs)
                zgphi_tmp(:,:,jt) = zgphi(:,:,jobs)
                zmask_tmp(:,:,jt) = zmask(:,:,jobs)

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90

@@ -50 +50 @@
       INTEGER :: npj
       INTEGER :: nsurfup    !: Observation counter used in obs_oper
+      INTEGER :: nrec       !: Number of surface observation records in window
 
       ! Arrays with size equal to the number of surface observations
@@ -56 +57 @@
          & mi,   &        !: i-th grid coord. for interpolating to surface observation
          & mj,   &        !: j-th grid coord. for interpolating to surface observation
+         & mt,   &        !: time record number for gridded data
          & nsidx,&        !: Surface observation number
          & nsfil,&        !: Surface observation number in file
@@ -67 +69 @@
          & ntyp           !: Type of surface observation product
 
-      CHARACTER(len=6), POINTER, DIMENSION(:) :: &
+      CHARACTER(len=8), POINTER, DIMENSION(:) :: &
          & cvars          !: Variable names
 
@@ -93 +95 @@
          & nsstpmpp       !: Global number of surface observations per time step
 
+      ! Arrays with size equal to the number of observation records in the window
+      INTEGER, POINTER, DIMENSION(:) :: &
+         & mrecstp   ! Time step of the records
+
       ! Arrays used to store source indices when 
       ! compressing obs_surf derived types
@@ -100 +106 @@
       INTEGER, POINTER, DIMENSION(:) :: &
          & nsind          !: Source indices of surface data in compressed data
+
+      ! Is this a gridded product?
+     
+      LOGICAL :: lgrid
 
    END TYPE obs_surf
@@ -160 +170 @@
          & surf%mi(ksurf),      &
          & surf%mj(ksurf),      &
+         & surf%mt(ksurf),      &
          & surf%nsidx(ksurf),   &
          & surf%nsfil(ksurf),   &
@@ -176 +187 @@
          & )
 
+      surf%mt(:) = -1
+
 
       ! Allocate arrays of number of surface data size * number of variables
@@ -190 +203 @@
          & )
 
+      surf%rext(:,:) = 0.0_wp 
+
       ! Allocate arrays of number of time step size
 
@@ -217 +232 @@
 
       surf%nsurfup     = 0
+      
+      ! Not gridded by default
+          
+      surf%lgrid       = .FALSE.
               
    END SUBROUTINE obs_surf_alloc
@@ -242 +261 @@
          & surf%mi,      &
          & surf%mj,      &
+         & surf%mt,      &
          & surf%nsidx,   &
          & surf%nsfil,   &
@@ -370 +390 @@
             newsurf%mi(insurf)    = surf%mi(ji)
             newsurf%mj(insurf)    = surf%mj(ji)
+            newsurf%mt(insurf)    = surf%mt(ji)
             newsurf%nsidx(insurf) = surf%nsidx(ji)
             newsurf%nsfil(insurf) = surf%nsfil(ji)
@@ -414 +435 @@
       newsurf%nstp     = surf%nstp
       newsurf%cvars(:) = surf%cvars(:)
+      
+      ! Set gridded stuff
+      
+      newsurf%mt(insurf)    = surf%mt(ji)
  
       ! Deallocate temporary data
@@ -454 +479 @@
          oldsurf%mi(jj)    = surf%mi(ji)
          oldsurf%mj(jj)    = surf%mj(ji)
+         oldsurf%mt(jj)    = surf%mt(ji)
          oldsurf%nsidx(jj) = surf%nsidx(ji)
          oldsurf%nsfil(jj) = surf%nsfil(ji)

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

@@ -8 +8 @@
    !!   obs_wri_prof   : Write profile observations in feedback format
    !!   obs_wri_surf   : Write surface observations in feedback format
-   !!   obs_wri_stats : Print basic statistics on the data being written out
+   !!   obs_wri_stats  : Print basic statistics on the data being written out
    !!----------------------------------------------------------------------
 
@@ -83 +83 @@
       TYPE(obfbdata) :: fbdata
       CHARACTER(LEN=40) :: clfname
-      CHARACTER(LEN=6) :: clfiletype
+      CHARACTER(LEN=10) :: clfiletype
       INTEGER :: ilevel
       INTEGER :: jvar
@@ -196 +196 @@
          fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
          fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
-         IF ( profdata%nqc(jo) > 10 ) THEN
-            fbdata%ioqc(jo)    = 4
+         IF ( profdata%nqc(jo) > 255 ) THEN
+            fbdata%ioqc(jo)    = IBSET(profdata%nqc(jo),2)
             fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
-            fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10
+            fbdata%ioqcf(2,jo) = profdata%nqc(jo)
          ELSE
             fbdata%ioqc(jo)    = profdata%nqc(jo)
@@ -236 +236 @@
                fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
                fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
-               IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN
-                  fbdata%ivlqc(ik,jo,jvar) = 4
+               IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN
+                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
                   fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
-                  fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10
+                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111')
                ELSE
                   fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
@@ -320 +320 @@
       TYPE(obfbdata) :: fbdata
       CHARACTER(LEN=40) :: clfname         ! netCDF filename
-      CHARACTER(LEN=6)  :: clfiletype
+      CHARACTER(LEN=10) :: clfiletype
       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
       INTEGER :: jo
@@ -395 +395 @@
          END DO
 
-      CASE('ICECON')
+      CASE('ICECONC')
 
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
@@ -418 +418 @@
          END DO
 
+      CASE('SSS')
+
+         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
+            &                 1 + iadd, iext, .TRUE. )
+
+         clfiletype = 'sssfb'
+         fbdata%cname(1)      = surfdata%cvars(1)
+         fbdata%coblong(1)    = 'Sea surface salinity'
+         fbdata%cobunit(1)    = 'psu'
+         DO je = 1, iext
+            fbdata%cextname(je) = pext%cdname(je)
+            fbdata%cextlong(je) = pext%cdlong(je,1)
+            fbdata%cextunit(je) = pext%cdunit(je,1)
+         END DO
+         fbdata%caddlong(1,1) = 'Model interpolated SSS'
+         fbdata%caddunit(1,1) = 'psu'
+         fbdata%cgrid(1)      = 'T'
+         DO ja = 1, iadd
+            fbdata%caddname(1+ja) = padd%cdname(ja)
+            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
+            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
+         END DO
+
+      CASE DEFAULT
+
+         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
+
       END SELECT
 
@@ -439 +466 @@
          fbdata%ivqc(jo,:)    = 0
          fbdata%ivqcf(:,jo,:) = 0
-         IF ( surfdata%nqc(jo) > 10 ) THEN
+         IF ( surfdata%nqc(jo) > 255 ) THEN
             fbdata%ioqc(jo)    = 4
             fbdata%ioqcf(1,jo) = 0
-            fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10
+            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
          ELSE
             fbdata%ioqc(jo)    = surfdata%nqc(jo)
@@ -474 +501 @@
          fbdata%idqc(1,jo)     = 0
          fbdata%idqcf(:,1,jo)  = 0
-         IF ( surfdata%nqc(jo) > 10 ) THEN
+         IF ( surfdata%nqc(jo) > 255 ) THEN
             fbdata%ivqc(jo,1)       = 4
             fbdata%ivlqc(1,jo,1)    = 4
             fbdata%ivlqcf(1,1,jo,1) = 0
-            fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10
+            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
          ELSE
             fbdata%ivqc(jo,1)       = surfdata%nqc(jo)

branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90

@@ -1240 +1240 @@
          & zdum,  &
          & zaamax
-       
+      
+      imax = -1 
       ! Main computation
       pflt = 1.0_wp