Index: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 7837)
+++ branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 7858)
@@ -48,7 +48,25 @@
    LOGICAL :: ln_diaobs           !: Logical switch for the obs operator
    LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs
+   LOGICAL :: ln_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres
+   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres
+   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres
+   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
+
+   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres)
+   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres)
+   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres)
+   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres)
+   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres)
+   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres)
+   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres)
+   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres)
 
    INTEGER :: nn_1dint       !: Vertical interpolation method
-   INTEGER :: nn_2dint       !: Horizontal interpolation method
+   INTEGER :: nn_2dint       !: Default horizontal interpolation method
+   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method 
+   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method 
+   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method 
+   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method 
+ 
    INTEGER, DIMENSION(imaxavtypes) :: &
       & nn_profdavtypes      !: Profile data types representing a daily average
@@ -61,4 +79,12 @@
       & nextrprof, &         !: Number of profile extra variables
       & nextrsurf            !: Number of surface extra variables
+   INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      & n2dintsurf           !: Interpolation option for surface variables
+   REAL(wp), DIMENSION(:), ALLOCATABLE :: &
+      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables
+      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables
+   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
+      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres
+      & llnightav            !: Logical for calculating night-time averages
 
    TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
@@ -132,7 +158,4 @@
          & cn_altbiasfile        ! Altimeter bias input filename
 
-      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
-         & clproffiles, &        ! Profile filenames
-         & clsurffiles           ! Surface filenames
 
       LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
@@ -153,10 +176,15 @@
       LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
       LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary
+
+      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
+      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
+
+      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
+         & clproffiles, &        ! Profile filenames
+         & clsurffiles           ! Surface filenames
+
       LOGICAL :: llvar1          ! Logical for profile variable 1
       LOGICAL :: llvar2          ! Logical for profile variable 1
-      LOGICAL :: llnightav       ! Logical for calculating night-time averages
-
-      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
-      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
+
       REAL(wp), POINTER, DIMENSION(:,:) :: &
          & zglam1, &             ! Model longitudes for profile variable 1
@@ -169,4 +197,5 @@
          & zmask2                ! Model land/sea mask associated with variable 2
 
+
       NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
          &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
@@ -176,4 +205,6 @@
          &            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,                   &
@@ -183,5 +214,12 @@
          &            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
@@ -236,100 +274,4 @@
       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, &
-         &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) )
-
-      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', &
-            &                    ' are set to .FALSE. so turning off calls to dia_obs'
-         nwarn = nwarn + 1
-         lk_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
-         IF (ln_vel3d) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', &
-               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
-         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
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', &
-               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
-         ENDIF
-         IF (ln_sst) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', &
-               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-         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 )
-         ENDIF
-#endif
-         IF (ln_sss) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', &
-               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-         ENDIF
-
-         IF (ln_logchl) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', &
-               &                   cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-         ENDIF
-
-         IF (ln_spm) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm   ', &
-               &                   cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-         ENDIF
-
-         IF (ln_fco2) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2  ', &
-               &                   cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
-         ENDIF
-
-         IF (ln_pco2) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2  ', &
-               &                   cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
-         ENDIF
-
-      ENDIF
-
-      !Write namelist settings to stdout
       IF(lwp) THEN
          WRITE(numout,*)
@@ -348,6 +290,6 @@
          WRITE(numout,*) '             Logical switch for FCO2 observations                    ln_fco2 = ', ln_fco2
          WRITE(numout,*) '             Logical switch for PCO2 observations                    ln_pco2 = ', ln_pco2
-         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
+         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
@@ -366,27 +308,148 @@
          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, &
+         &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) )
+
+      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', &
+            &                    ' are set to .FALSE. so turning off calls to dia_obs'
+         nwarn = nwarn + 1
+         lk_diaobs = .FALSE.
+         RETURN
+      ENDIF
+
+      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes
+      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(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes
+      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
+
+         IF (ln_logchl) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', &
+               &                   cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'logchl',         &
+               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
+               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
+               &                    lfpindegs, llnightav )
+         ENDIF
+
+         IF (ln_spm) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm   ', &
+               &                   cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'spm   ',         &
+               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
+               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
+               &                    lfpindegs, llnightav )
+         ENDIF
+
+         IF (ln_fco2) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2  ', &
+               &                   cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'fco2  ',         &
+               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
+               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
+               &                    lfpindegs, llnightav )
+         ENDIF
+
+         IF (ln_pco2) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2  ', &
+               &                   cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
+            CALL obs_setinterpopts( nsurftypes, jtype, 'pco2  ',         &
+               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
+               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
+               &                    lfpindegs, llnightav )
+         ENDIF
+
+      ENDIF
+
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+
 
       !-----------------------------------------------------------------------
@@ -404,5 +467,5 @@
       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')
@@ -488,7 +551,5 @@
             nvarssurf(jtype) = 1
             nextrsurf(jtype) = 0
-            llnightav = .FALSE.
             IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
-            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight
 
             !Read in surface obs types
@@ -496,11 +557,12 @@
                &               clsurffiles(jtype,1:ifilessurf(jtype)), &
                &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
-               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav )
+               &               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
 
@@ -515,5 +577,5 @@
                ENDIF
 
-               CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 
+               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
                   &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
 
@@ -623,5 +685,4 @@
          & zgphi1,    &            ! Model latitudes for prof variable 1
          & zgphi2                  ! Model latitudes for prof variable 2
-      LOGICAL :: llnightav        ! Logical for calculating night-time average
 
 
@@ -697,10 +758,8 @@
             !Defaults which might be changed
             zsurfmask(:,:) = tmask(:,:,1)
-            llnightav = .FALSE.
 
             SELECT CASE ( TRIM(cobstypessurf(jtype)) )
             CASE('sst')
                zsurfvar(:,:) = tsn(:,:,1,jp_tem)
-               llnightav = ln_sstnight
             CASE('sla')
                zsurfvar(:,:) = sshn(:,:)
@@ -741,5 +800,4 @@
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               llnightav = .FALSE.
                zsurfmask(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
                ! Take the log10 where we can, otherwise exclude
@@ -846,5 +904,7 @@
             CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
                &               nit000, idaystp, zsurfvar, zsurfmask,    &
-               &               nn_2dint, llnightav )
+               &               n2dintsurf(jtype), llnightav(jtype),     &
+               &               ravglamscl(jtype), ravgphiscl(jtype),     &
+               &               lfpindegs(jtype) )
 
          END DO
@@ -978,5 +1038,6 @@
 
       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
@@ -1177,5 +1238,85 @@
     ENDIF
 
+    IF(lwp) THEN    
+       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
+       DO jfile = 1, ifiles(jtype)
+          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
+       END DO
+    ENDIF
+
     END SUBROUTINE obs_settypefiles
 
+    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
+               &                  n2dint_default, n2dint_type,        &
+               &                  ravglamscl_type, ravgphiscl_type,   &
+               &                  lfp_indegs_type, lavnight_type,     &
+               &                  n2dint, ravglamscl, ravgphiscl,     &
+               &                  lfpindegs, lavnight )
+
+    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
+    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
+    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
+    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
+    REAL(wp), INTENT(IN) :: &
+       &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
+       &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
+    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
+    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
+    CHARACTER(len=6), INTENT(IN) :: ctypein 
+
+    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
+       &                    n2dint 
+    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
+       &                    ravglamscl, ravgphiscl
+    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
+       &                    lfpindegs, lavnight
+
+    lavnight(jtype) = lavnight_type
+
+    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
+       n2dint(jtype) = n2dint_type
+    ELSE
+       n2dint(jtype) = n2dint_default
+    ENDIF
+
+    ! For averaging observation footprints set options for size of footprint 
+    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
+       IF ( ravglamscl_type > 0._wp ) THEN
+          ravglamscl(jtype) = ravglamscl_type
+       ELSE
+          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
+                         'scale (ravglamscl) for observation type '//TRIM(ctypein) )      
+       ENDIF
+
+       IF ( ravgphiscl_type > 0._wp ) THEN
+          ravgphiscl(jtype) = ravgphiscl_type
+       ELSE
+          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
+                         'scale (ravgphiscl) for observation type '//TRIM(ctypein) )      
+       ENDIF
+
+       lfpindegs(jtype) = lfp_indegs_type 
+
+    ENDIF
+
+    ! Write out info 
+    IF(lwp) THEN
+       IF ( n2dint(jtype) <= 4 ) THEN
+          WRITE(numout,*) '             '//TRIM(ctypein)// &
+             &            ' model counterparts will be interpolated horizontally'
+       ELSE IF ( n2dint(jtype) <= 6 ) THEN
+          WRITE(numout,*) TRIM(ctypein)// &
+             &            ' model counterparts will be averaged horizontally'
+          WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype)
+          WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype)
+          IF ( lfpindegs(jtype) ) THEN
+              WRITE(numout,*) '             '//'    (in degrees)'
+          ELSE
+              WRITE(numout,*) '             '//'    (in metres)'
+          ENDIF
+       ENDIF
+    ENDIF
+
+    END SUBROUTINE obs_setinterpopts
+
 END MODULE diaobs
Index: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_averg_h2d.F90
===================================================================
--- branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_averg_h2d.F90	(revision 7858)
+++ branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_averg_h2d.F90	(revision 7858)
@@ -0,0 +1,54 @@
+MODULE obs_averg_h2d
+   !!======================================================================
+   !!                       ***  MODULE obs_averg_h2d   ***
+   !! Observation diagnostics: Perform the horizontal averaging
+   !!                          from model grid to observation footprint
+   !!=====================================================================
+
+   !!----------------------------------------------------------------------
+   !!   obs_averg_h2d     : Horizontal averaging to the observation footprint
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_kind, ONLY : &  ! Precision variables
+      & wp
+   USE phycst,   ONLY : &  ! Physical constants
+      & rad,  &
+      & ra,   &
+      & rpi
+   USE dom_oce,   ONLY : &
+      & jpi, jpj, &
+      & e1t, e2t, &
+      & e1f, e2f, &
+      & glamt, gphit, &
+      & nproc
+   USE in_out_manager
+   USE obs_const, ONLY : &
+      & obfillflt		! Fillvalue
+   USE obs_utils           ! Utility functions
+   USE lib_mpp,   ONLY : &
+      & ctl_warn, ctl_stop, &
+      & mpp_min, lk_mpp
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE obs_avg_h2d_rad, & ! Horizontal averaging using a radial footprint
+      &    obs_avg_h2d_rec, & ! Horizontal averaging using a rectangular footprint
+      &    obs_deg2dist,    & ! Conversion of distance in degrees to distance in metres
+      &    obs_dist2corners   ! Distance from the centre of obs footprint to the corners of a grid box
+   
+   PUBLIC obs_avg_h2d,      & ! Horizontal averaging to the observation footprint
+      &   obs_avg_h2d_init, & ! Set up weights for the averaging
+      &   obs_max_fpsize      ! Works out the maximum number of grid points required for the averaging
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obs_averg_h2d.F90 2715 2011-03-30 15:58:35Z rblod $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+ 
+#include "obsaverg_h2d.h90"
+
+END MODULE obs_averg_h2d
Index: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 7837)
+++ branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 7858)
@@ -19,4 +19,8 @@
       & 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,    &
@@ -24,7 +28,7 @@
    USE obs_const,  ONLY :    &    ! Obs fill value
       & obfillflt
-   USE dom_oce,       ONLY : &    ! Model lats/lons
-      & glamt, &
-      & gphit
+   USE dom_oce,       ONLY : &
+      & glamt, glamf, &
+      & gphit, gphif
    USE lib_mpp,       ONLY : &    ! Warning and stopping routines
       & ctl_warn, ctl_stop
@@ -667,7 +671,8 @@
    END SUBROUTINE obs_prof_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 )
 
       !!-----------------------------------------------------------------------
@@ -691,4 +696,8 @@
       !!        - polynomial (quadrilateral grid)  (k2dint = 4)
       !!
+      !!    Two horizontal averaging schemes are also available:
+      !!        - weighted radial footprint        (k2dint = 5)
+      !!        - weighted rectangular footprint   (k2dint = 6)
+      !!
       !!
       !! ** Action  :
@@ -697,4 +706,5 @@
       !!      ! 07-03 (A. Weaver)
       !!      ! 15-02 (M. Martin) Combined routine for surface types
+      !!      ! 17-03 (M. Martin) Added horizontal averaging options
       !!-----------------------------------------------------------------------
 
@@ -718,4 +728,9 @@
          & 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
@@ -726,8 +741,12 @@
       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,      &
@@ -737,12 +756,15 @@
       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,  &
@@ -756,4 +778,10 @@
       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
@@ -822,32 +850,58 @@
 
       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
@@ -855,8 +909,8 @@
 
          ALLOCATE( &
-            & zsurfm(2,2,isurf)  &
+            & zsurfm(imaxifp,imaxjfp,isurf)  &
             & )
 
-         CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, &
+         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
             &               surfdataqc%vdmean(:,:), zsurfm )
 
@@ -889,15 +943,34 @@
          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
 
@@ -914,4 +987,5 @@
       ! Deallocate the data for interpolation
       DEALLOCATE( &
+         & zweig, &
          & igrdi, &
          & igrdj, &
@@ -919,5 +993,10 @@
          & zgphi, &
          & zmask, &
-         & zsurf  &
+         & zsurf, &
+         & zsurftmp, &
+         & zglamf, &
+         & zgphif, &
+         & igrdip1,&
+         & igrdjp1 &
          & )
 
Index: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsaverg_h2d.h90
===================================================================
--- branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsaverg_h2d.h90	(revision 7858)
+++ branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsaverg_h2d.h90	(revision 7858)
@@ -0,0 +1,777 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obs_averg_h2d.h90 2474 2010-12-16 15:32:33Z djlea $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE obs_avg_h2d_init( kpk, kpk2, kmaxifp, kmaxjfp, k2dint, plam,  pphi, &
+      &                         pglam, pgphi, pglamf, pgphif, pmask, plamscl, pphiscl, lindegrees, &
+      &                         pweig, pobsmask, iminpoints )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_avg_h2d_init  ***
+      !!
+      !! ** Purpose : Computes weights for horizontal averaging to the 
+      !!              observation footprint.
+      !!
+      !! ** Method  : Horizontal averaging to the observation footprint using
+      !!              model values at a defined area.
+      !!
+      !!    Averaging schemes :
+      !!
+      !!    Two horizontal averaging schemes are available:
+      !!        - weighted radial footprint        (k2dint = 5)
+      !!        - weighted rectangular footprint   (k2dint = 6)
+      !!
+      !! History :
+      !!        ! 13-10 (M. Martin)
+      !!-----------------------------------------------------------------------
+      !! * Modules used
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kpk,   &             ! Parameter values for automatic arrays
+         & kpk2,  &
+         & kmaxifp,  &          ! Max size of model points in i-direction for obs footprint
+         & kmaxjfp,  &          ! Max size of model points in j-direction for obs footprint
+         & k2dint               ! Averaging scheme options - see header
+      REAL(KIND=wp), INTENT(INOUT) :: &
+         & plam, &              ! Geographical (lat,lon) coordinates of
+         & pphi                 ! observation
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
+         & pglam, &             ! Model variable lon
+         & pgphi                ! Model variable lat
+      REAL(KIND=wp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
+         & pglamf, &            ! Model variable lon at corners of grid-boxes
+         & pgphif               ! Model variable lat at corners of grid-boxes
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
+         & pmask                ! Model variable mask
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plamscl, &           ! Diameter (lat,lon) of obs footprint in metres
+         & pphiscl              ! This is the full width (rather than half-width)
+      LOGICAL, INTENT(IN) :: &
+         & lindegrees           ! T=> obs footprint specified in degrees, F=> in metres
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) ::  &
+         & pweig                ! Weights for averaging
+      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) ::  &
+         & pobsmask             ! Vertical mask for observations
+      INTEGER, INTENT(IN), OPTIONAL :: &
+         & iminpoints           ! Reject point which is not surrounded
+                                ! by at least iminpoints sea points
+
+      !! * Local declarations
+      INTEGER :: &
+         & jk
+      INTEGER :: &
+         & ikmax
+
+
+      !------------------------------------------------------------------------
+      !
+      !------------------------------------------------------------------------
+
+      !------------------------------------------------------------------------
+      ! Initialize number of levels
+      !------------------------------------------------------------------------
+      IF ( kpk2 == 1 ) THEN
+         ikmax = 1
+      ELSEIF ( kpk2 == kpk) THEN
+         ikmax = kpk-1
+      ENDIF
+
+
+         SELECT CASE (k2dint)
+         CASE(5)
+            CALL obs_avg_h2d_rad( kpk2, ikmax, kmaxifp, kmaxjfp, plam, pphi, &
+      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
+         CASE(6)
+            CALL obs_avg_h2d_rec( kpk2, ikmax, kmaxifp, kmaxjfp, plam, pphi, &
+      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
+         END SELECT
+
+
+      END SUBROUTINE obs_avg_h2d_init
+
+
+      SUBROUTINE obs_avg_h2d_rad( kpk2, kmax, kmaxifp, kmaxjfp, plam, pphi, &
+      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_avg_h2d_rad  ***
+      !!
+      !! ** Purpose : Computes weights for horizontal averaging to the 
+      !!              observation using a radial footprint.
+      !!
+      !! ** Method  : Calculate whether each grid box is completely or 
+      !!              partially within the observation footprint.
+      !!              If it is partially in the footprint then calculate
+      !!              the ratio of the area inside the footprint to the total 
+      !!              area of the grid box.
+      !!
+      !! History :
+      !!        ! 14-01 (M. Martin)
+      !!-----------------------------------------------------------------------
+      !! * Modules used
+      USE phycst,   ONLY : &  ! Physical constants
+         & ra,  &
+         & rpi
+
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kpk2, &             ! Parameter values for automatic arrays
+         & kmax
+
+      INTEGER, INTENT(IN) :: &
+         & kmaxifp,   &         ! Max size of model points in i-direction for obs footprint
+         & kmaxjfp              ! Max size of model points in j-direction for obs footprint
+
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plam, &
+         & pphi                 ! Geographical (lat,lon) coordinates of
+                                ! observation
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plamscl, &           ! Diameter (lat,lon) of obs footprint in metres or degrees (see below)
+         & pphiscl              ! This is the full width (rather than half-width)
+      LOGICAL, INTENT(IN) :: &
+         & lindegrees           ! T=>scales specified in degrees, F=> in metres
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
+         & pmask                ! Model variable mask
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
+         & pglam, &             ! Model variable lon
+         & pgphi                ! Model variable lat
+      REAL(KIND=wp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
+         & pglamf, &             ! Model variable lon at corners of grid boxes
+         & pgphif                ! Model variable lat at corners of grid boxes
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) ::  &
+         & pweig                ! Weights for interpolation
+
+      !! Local declarations
+      INTEGER :: ji, jj, jk
+      INTEGER :: jvert, jis, jjs
+      INTEGER :: jnumvert, jnumvertbig
+      INTEGER, PARAMETER :: &
+         & jnumsubgrid = 20     ! The number of sub grid-boxes (in x and y directions) used to approximate area of obs fp
+
+      REAL(KIND=wp), DIMENSION(4) :: &
+         & zxvert, zyvert, &    ! The lon/lat of the vertices(corners) of the grid box in m relative to centre of obs fp
+         & zdist                ! Distance of each vertex to the centre of the obs footprint
+      REAL(KIND=wp), DIMENSION(4) :: &
+         & zxgrid, zygrid, &      ! Distance of each vertex of grid box to the centre of the grid box in x/y directions
+         & zdgrid
+      REAL(KIND=wp) :: &
+         & zdx, zdy, &          ! The sub grid-box sizes (in metres)
+         & zarea_subbox, &      ! The area of each sub grid-box (in metres squared)
+         & zxpos, zypos, &      ! The x,y position (relative to centre of obs footprint) of the centre of each sub grid-box
+         & zsubdist, &          ! The distance of the centre of each sub grid-box from the centre of the obs footprint
+         & zarea_fp, &          ! Total area of obs footprint within the grid box
+         & zareabox             ! Total area of the grid box
+      REAL(KIND=wp) :: &
+         & zphiscl_m, &         ! Diameter of obs footprint in metres
+         & zlamscl_m
+      !---------------------------------------------------------------------------------------------------
+      !Initialise weights to zero.
+      pweig(:,:,:) = 0.0_wp
+      
+      !Two footprint sizes can be specified in the namelist but this routine assumes a circular footprint. 
+      !If the two sizes are different then write out a warning.
+      IF ( pphiscl /= plamscl ) THEN
+               CALL ctl_warn( 'obs_avg_h2d_rad:',   &
+                  &           'The two components of the obs footprint size are not equal', &
+                  &           'yet the radial option has been selected - using pphiscl here' )
+      ENDIF
+      
+      DO jk = 1, kmax
+         DO ji = 1, kmaxifp
+            DO jj = 1, kmaxjfp
+            
+               IF ( pmask(ji,jj,jk) == 1.0_wp ) THEN
+
+                  IF ( lindegrees ) THEN
+                     !If the scales are specified in degrees, work out the
+                     !scales (metres) in x/y directions
+                     CALL obs_deg2dist( 1, 1, pglam(ji,jj), pgphi(ji,jj), &
+                        &               plamscl, pphiscl, zlamscl_m, zphiscl_m )
+                  ELSE
+                     zphiscl_m = pphiscl
+                  ENDIF
+
+
+                  ! Work out the area of the grid box using distance of corners relative to centre of grid box
+                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
+                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
+                     &                  pglam(ji,jj), pgphi(ji,jj), zxgrid, zygrid, zdgrid)
+                  zareabox = ABS( zxgrid(1) - zxgrid(2) ) * ABS( zygrid(1) - zygrid(4) )
+
+                  !1. Determine how many of the vertices of the grid box lie within the circle
+                  
+                  !For each vertex, calculate its location and distance relative 
+                  !to the centre of the observation footprint
+                  
+                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
+                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
+                     &                  plam, pphi, zxvert, zyvert, zdist)
+
+                  jnumvert = 0
+                  jnumvertbig = 0
+                  DO jvert = 1, 4
+
+                     !If the distance to the center to the observation footprint is less
+                     !than the radius of the footprint (half the diameter) then this
+                     !vertex is within the observation footprint
+                     IF ( zdist(jvert) <= ( zphiscl_m / 2.0_wp ) ) jnumvert = jnumvert + 1
+                     
+                     !For expediency, check if the vertices are "nearly" within the obs
+                     !footprint as if none of them are close to the edge of the footprint
+                     !then the footprint is unlikely to be intersecting the grid box
+                     IF ( zdist(jvert) - ( 0.5_wp * zareabox ) <= ( zphiscl_m / 2.0 ) ) &
+                        & jnumvertbig = jnumvertbig + 1
+                     
+                  END DO
+                  
+                  !2. If none of the vertices are even close to the edge of the obs
+                  !footprint then leave weight as zero and cycle to next grid box.
+                  IF ( jnumvertbig == 0 ) CYCLE
+
+                  !3. If all the vertices of the box are within the observation footprint then the 
+                  !      whole grid box is within the footprint so set the weight to one and 
+                  !      move to the next grid box.
+                  IF ( jnumvert == 4 ) THEN
+                     pweig(ji,jj,jk) = 1.0_wp
+                     CYCLE
+                  ENDIF
+
+
+                  !4. Use a brute force technique for calculating the area within 
+                  !   the grid box covered by the obs footprint.
+                  !  (alternative could be to use formulae on 
+                  !         http://mathworld.wolfram.com/Circle-LineIntersection.html)
+                  !   For now split the grid box into a specified number of smaller 
+                  !   boxes and add up the area of those whose centre is within the obs footprint.
+                  !   Order of vertices is 1=topleft, 2=topright, 3=bottomright, 4=bottomleft
+                  zdx = ABS( zxvert(3) - zxvert(4) ) / REAL(jnumsubgrid, wp)
+                  zdy = ABS( zyvert(1) - zyvert(4) ) / REAL(jnumsubgrid, wp)
+                  zarea_subbox = zdx * zdy
+
+                  zarea_fp = 0.0_wp
+                  DO jis = 1, jnumsubgrid
+                     zxpos = zxvert(4) + ( REAL(jis, wp) * zdx ) - (0.5_wp * zdx )
+                     DO jjs = 1, jnumsubgrid
+                        !Find the distance of the centre of this sub grid box to the
+                        !centre of the obs footprint
+                        zypos = zyvert(4) + ( REAL(jjs, wp) * zdy ) - ( 0.5_wp * zdy )
+                        zsubdist = SQRT( (zxpos * zxpos) + (zypos * zypos) )
+                        IF ( zsubdist < ( zphiscl_m / 2.0_wp ) ) &
+                           &  zarea_fp = zarea_fp + zarea_subbox
+                     END DO
+                  END DO
+
+                  !6. Calculate the ratio of the area of the footprint within the box 
+                  !   to the total area of the grid box and use this fraction to weight 
+                  !   the model value in this grid box.
+                  pweig(ji,jj,jk) = MIN( zarea_fp / zareabox, 1.0_wp )
+
+               END IF  !pmask
+            END DO
+         END DO
+      END DO
+      
+      END SUBROUTINE obs_avg_h2d_rad
+
+
+      SUBROUTINE obs_avg_h2d_rec( kpk2, kmax, kmaxifp, kmaxjfp, plam, pphi, &
+      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_avg_h2d_rec  ***
+      !!
+      !! ** Purpose : Computes weights for horizontal averaging to the
+      !!              observation using a rectangular footprint which
+      !!              is aligned with lines of lat/lon.
+      !!
+      !! ** Method  : Horizontal averaging to the observation footprint using
+      !!              model values at a defined area.
+      !!
+      !! History :
+      !!        ! 14-01 (M. Martin)
+      !!-----------------------------------------------------------------------
+      !! * Modules used
+      USE phycst,   ONLY : &    ! Physical constants
+         & ra,  &
+         & rpi
+
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kpk2, &              ! Parameter values for automatic arrays
+         & kmax
+
+      INTEGER, INTENT(IN) :: &
+         & kmaxifp,   &         ! Max size of model points in i-direction for obs footprint
+         & kmaxjfp              ! Max size of model points in j-direction for obs footprint
+
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plam, &
+         & pphi                 ! Geographical (lat,lon) coordinates of
+                                ! observation
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plamscl, &
+         & pphiscl              ! Width in x/y directions of obs footprint in metres
+                                ! This is the full width (rather than half-width)
+      LOGICAL, INTENT(IN) :: &
+         & lindegrees           !T=> scales specified in degrees, F=> in metres
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
+         & pmask                ! Model variable mask
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
+         & pglam, &             ! Model variable lat at centre of grid boxes
+         & pgphi                ! Model variable lon at centre of grid boxes
+      REAL(KIND=wp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
+         & pglamf, &             ! Model variable lat at corners of grid boxes
+         & pgphif                ! Model variable lon at corners of grid boxes
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) ::  &
+         & pweig                ! Weights for interpolation
+
+      !! Local declarations
+      INTEGER :: ji, jj, jk
+      INTEGER :: jvert
+      INTEGER, DIMENSION(4) :: &
+         & jnumvert
+      REAL(KIND=wp), DIMENSION(4) :: &
+         & zxvert, zyvert         ! The lon/lat of the vertices(corners) of the grid box in m relative to centre of obs fp
+      REAL(KIND=wp), DIMENSION(4) :: &
+         & zdist                  ! Distance of each vertex to the centre of the obs footprint
+      REAL(KIND=wp), DIMENSION(4) :: &
+         & zxgrid, zygrid, &      ! Distance of each vertex of grid box to the centre of the grid box in x/y directions
+         & zdgrid
+      REAL(KIND=wp) :: &
+         & zareabox, &            ! Total area of grid box
+         & zarea_fp, &            ! Total area of obs footprint
+         & zarea_intersect        ! Area of the intersection between the grid box and the obs footprint
+      REAL(KIND=wp) :: &
+         & zlamscl_m, &
+         & zphiscl_m              ! Total width (lat,lon) of obs footprint in metres
+      REAL(KIND=wp) :: &
+         & z_awidth, z_aheight, & ! Width and height of model grid box
+         & z_cwidth, z_cheight    ! Width and height of union of model grid box and obs footprint
+      REAL(KIND=wp) :: &
+         & zleft, zright, &       ! Distance (metres) of corners area of intersection
+         & ztop, zbottom          ! between grid box and obs footprint
+
+      !-----------------------------------------------------------------------
+
+      !Initialise weights to zero
+      pweig(:,:,:) = 0.0_wp
+      
+      !Loop over the grid boxes which have been identified as potentially being within the
+      !observation footprint
+      DO jk = 1, kmax
+         DO ji = 1, kmaxifp
+            DO jj = 1, kmaxjfp
+
+               IF ( pmask(ji,jj,jk) == 1.0_wp ) THEN
+
+
+                  IF ( lindegrees ) THEN
+                     !If the scales are specified in degrees, work out the
+                     !scales (metres) in x/y directions
+                     CALL obs_deg2dist( 1, 1, pglam(ji,jj), pgphi(ji,jj), &
+                        &               plamscl, pphiscl, zlamscl_m, zphiscl_m )
+                  ELSE
+                     zlamscl_m = plamscl
+                     zphiscl_m = pphiscl
+                  ENDIF
+
+                  ! Work out the area of the grid box using distance of corners relative to centre of grid box
+                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
+                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
+                     &                  pglam(ji,jj), pgphi(ji,jj), zxgrid, zygrid, zdgrid)
+
+                  !Calculate width and height of model grid box
+                  z_awidth  = ABS( zxgrid(1) - zxgrid(2) )
+                  z_aheight = ABS( zygrid(1) - zygrid(4) )
+                  zareabox = z_awidth * z_aheight
+
+                  ! Work out area of the observation footprint
+                  zarea_fp = zlamscl_m * zphiscl_m
+
+                  ! For each corner of the grid box, calculate its location and distance relative
+                  ! to the centre of the observation footprint
+                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
+                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
+                     &                  plam, pphi, zxvert, zyvert, zdist)
+
+                  !Work out maximum width and height of rectangle covered by corners of obs fp and grid box
+                  z_cwidth  = MAX( zxvert(1), zxvert(2), -zlamscl_m/2.0_wp, zlamscl_m/2.0_wp ) - &
+                     &       MIN( zxvert(1), zxvert(2), -zlamscl_m/2.0_wp, zlamscl_m/2.0_wp )
+                     
+                  z_cheight = MAX( zyvert(1), zyvert(4), zphiscl_m/2.0_wp, -zphiscl_m/2.0_wp ) - &
+                     &       MIN( zyvert(1), zyvert(4), zphiscl_m/2.0_wp, -zphiscl_m/2.0_wp )
+                  
+                  IF ( ( z_cwidth  >= z_awidth  + zlamscl_m  ) .OR. &
+                     & ( z_cheight >= z_aheight + zphiscl_m ) ) THEN
+                     !The obs footprint and the model grid box don't overlap so set weight to zero
+                     pweig(ji,jj,jk) = 0.0_wp
+                  ELSE IF ( ( z_cwidth  == zlamscl_m ) .AND. &
+                     &      ( z_cheight == zphiscl_m ) ) THEN
+                     !The grid box is totally contained within the obs footprint so set weight to one
+                     pweig(ji,jj,jk) = 1.0_wp
+                  ELSE IF ( ( z_cwidth  == z_awidth  ) .AND. &
+                     &      ( z_cheight == z_aheight ) ) THEN
+                     !The obs footprint is totally contained within the grid box so set weight as ratio of the two
+                     pweig(ji,jj,jk) = zarea_fp / zareabox
+                  ELSE 
+                     !The obs footprint and the grid box overlap so calculate the area of the intersection of the two
+                     zleft   = max(zxvert(1), -zlamscl_m/2.0_wp)
+                     zright  = min(zxvert(2),  zlamscl_m/2.0_wp)
+                     zbottom = max(zyvert(4), -zphiscl_m/2.0_wp)
+                     ztop    = min(zyvert(1),  zphiscl_m/2.0_wp)
+                     
+                     IF (  ( zleft < zright ) .AND. ( zbottom < ztop ) ) THEN
+                        zarea_intersect = ( zright - zleft ) * ( ztop - zbottom )
+                        pweig(ji,jj,jk) = zarea_intersect / zareabox
+                     ENDIF
+                  ENDIF
+
+               END IF !pmask
+            END DO
+         END DO
+      END DO
+
+   END SUBROUTINE obs_avg_h2d_rec
+
+   SUBROUTINE obs_avg_h2d( kpk, kpk2, kmaxifp, kmaxjfp, pweig, pmod, pobsk )
+
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_int_h2d  ***
+      !!
+      !! ** Purpose : Horizontal averaging to the observation footprint.
+      !!
+      !! ** Method  : Average the model points based on the weights already calculated.
+      !!
+      !! ** Action  :
+      !!                   
+      !! References : 
+      !!
+      !! History :
+      !!        ! 13/10. M. Martin.
+      !!-----------------------------------------------------------------------
+      !! * Modules used
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kpk,   &             ! Parameter values for automatic arrays
+         & kpk2
+      INTEGER, INTENT(IN) :: &
+         & kmaxifp,   &         ! Max size of model points in i-direction for obs footprint
+         & kmaxjfp              ! Max size of model points in j-direction for obs footprint
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
+         & pweig                ! Interpolation weights
+      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
+         & pmod                 ! Model variable to interpolate
+      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) ::  &
+         & pobsk                ! Model profile interpolated to obs (i,j) pt
+
+      INTEGER :: &
+         & jk
+      INTEGER :: &
+         & ikmax
+      REAL(KIND=wp) :: &
+         & zsum
+
+      !------------------------------------------------------------------------
+      ! Initialize number of levels
+      !------------------------------------------------------------------------
+      IF ( kpk2 == 1 ) THEN
+         ikmax = 1
+      ELSEIF ( kpk2 == kpk) THEN
+         ikmax = kpk-1
+      ENDIF
+
+      !------------------------------------------------------------------------
+      ! Average model values to the observation footprint
+      !------------------------------------------------------------------------
+      pobsk = obfillflt
+
+      DO jk = 1, ikmax
+
+         zsum = SUM( pweig(:,:,jk) )
+
+         IF ( zsum /= 0.0_wp ) THEN
+            pobsk(jk) = SUM ( pweig(:,:,jk) * pmod(:,:,jk), Mask=pweig(:,:,jk) > 0.0_wp )
+            pobsk(jk) = pobsk(jk) / zsum
+         END IF
+
+      END DO
+
+   END SUBROUTINE obs_avg_h2d
+
+   SUBROUTINE obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, pmask, kmaxifp, kmaxjfp )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_max_fpsize  ***
+      !!
+      !! ** Purpose : Calculate maximum number of grid points which may
+      !!              need to be used in the averaging in the global domain.
+      !!
+      !!
+      !! ** Method  : Work out the minimum grid size and work out
+      !!              how many of the smallest grid points would be needed
+      !!              to cover the scale of the observation footprint.
+      !!              This needs to be done using the max/min of the global domain
+      !!              as the obs can be distributed from other parts of the grid.
+      !!
+      !! ** Action  :
+      !!                   
+      !! References : 
+      !!
+      !! History :
+      !!        ! 14/01. M. Martin. 
+      !!-----------------------------------------------------------------------
+      !! * Modules used
+      !! * Arguments
+      INTEGER , INTENT(IN) :: &
+         & k2dint                   !Type of interpolation/averaging used
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plamscl,   &             !Total width/radius in metres of the observation footprint
+         & pphiscl                  !
+      LOGICAL, INTENT(IN) :: &
+         & lindegrees               !T=> plamscl and pphiscl are specified in degrees
+      REAL(KIND=wp), DIMENSION(jpi,jpj), INTENT(IN) :: &
+         & pmask                    !Land/sea mask
+                                    !F=> plamscl and pphiscl are specified in metres
+      INTEGER, INTENT(OUT)  :: &
+         & kmaxifp,   &             !Max number of grid points in i,j directions to use in averaging
+         & kmaxjfp                  !these have to be even so that the footprint is centred
+
+      !! * Local variables
+      REAL(KIND=wp) :: &
+         & ze1min,     &            !Minimum global grid-size in i,j directions
+         & ze2min
+      REAL(KIND=wp) :: &
+         & zphiscl_m, &
+         & zlamscl_m
+      !------------------------------------------------------------------------
+
+      IF ( k2dint <= 4 ) THEN
+         !If interpolation is being used then only need to use a 2x2 footprint
+         kmaxifp = 2
+         kmaxjfp = 2
+
+      ELSE
+
+         IF ( lindegrees ) THEN
+            !If the scales are specified in degrees, work out the max 
+            !distance (metres) in x/y directions
+            CALL obs_deg2dist( jpi, jpj, glamt, gphit, &
+               &               plamscl, pphiscl, zlamscl_m, zphiscl_m )
+         ELSE
+            zlamscl_m = plamscl
+            zphiscl_m = pphiscl
+         ENDIF
+
+         ze1min = MINVAL( e1t(:,:), mask = pmask(:,:) == 1._wp )
+         ze2min = MINVAL( e2t(:,:), mask = pmask(:,:) == 1._wp )
+         
+         IF(lk_mpp) THEN
+            CALL mpp_min( ze1min )
+            CALL mpp_min( ze2min )
+         ENDIF
+
+         kmaxifp = ceiling(zlamscl_m/ze1min) + 1
+         kmaxjfp = ceiling(zphiscl_m/ze2min) + 1
+         
+         !Ensure that these numbers are even
+         kmaxifp = kmaxifp + MOD(kmaxifp,2)
+         kmaxjfp = kmaxjfp + MOD(kmaxjfp,2)
+         
+
+      ENDIF
+
+   END SUBROUTINE obs_max_fpsize
+
+   SUBROUTINE obs_deg2dist( ki, kj, pglam, pgphi, plamscl_deg, pphiscl_deg, &
+      &                     plamscl_max, pphiscl_max )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_deg2dist  ***
+      !!
+      !! ** Purpose : Calculate the maximum distance in m of the length scale
+      !!              in degrees.
+      !!
+      !! ** Method  : At each lon/lat point, work out the distances in the 
+      !!              zonal and meridional directions.
+      !!
+      !! ** Action  :
+      !!                   
+      !! References : 
+      !!
+      !! History :
+      !!        ! 14/01. M. Martin. 
+      !!-----------------------------------------------------------------------
+      !! * Modules used
+      !! * Arguments
+      INTEGER , INTENT(IN) :: &
+         & ki, kj                   !x/y dimensions of input lat/lon variables
+      REAL(KIND=wp), INTENT(IN), DIMENSION(ki,kj) :: &
+         & pglam, pgphi             !Longitude and latitudes of grid points
+      REAL(KIND=wp), INTENT(IN) :: &
+         & plamscl_deg,   &         !Size in degrees of the observation footprint
+         & pphiscl_deg              !
+      REAL(KIND=wp), INTENT(OUT) :: &
+         & plamscl_max, &           !Maximum size of obs footprint in metres
+         & pphiscl_max
+      
+      !! * Local declarations
+      INTEGER :: &
+         & ji, jj                   !Counters
+      REAL(KIND=wp) :: &
+         & zlon1, zlon2, &          !Lon values surrounding centre of grid point
+         & zlat1, zlat2, &          !Lat values surrounding centre of grid point
+         & zdlat, zdlon             !Distance in radians in lat/lon directions
+      REAL(KIND=wp) :: &
+         & za1, za2, za, zc, zd
+      
+      plamscl_max = -1.0_wp
+      pphiscl_max = -1.0_wp
+
+      DO ji = 1, ki
+         DO jj = 1, kj
+
+            !Calculate distance in metres in zonal(x) direction
+
+            zlon1 = rad * ( pglam(ji,jj) + ( 0.5_wp * plamscl_deg ) )
+            zlon2 = rad * ( pglam(ji,jj) - ( 0.5_wp * plamscl_deg ) )
+            zlat1 = rad * pgphi(ji,jj)
+            zlat2 = rad * pgphi(ji,jj)
+            zdlon = zlon2 - zlon1
+            zdlat = zlat2 - zlat1
+
+            za1 = sin( zdlat/2.0_wp )
+            za2 = sin( zdlon/2.0_wp )
+            za = ( za1 * za1 ) + ( COS( zlat1 ) * COS( zlat2 ) * ( za2 * za2 ) )
+            zc = 2.0_wp * atan2( SQRT( za ), SQRT( 1.0_wp-za ) )
+            zd = ra * zc
+
+            IF ( zd > plamscl_max ) plamscl_max = zd
+
+            !Calculate distance in metres in meridional(y) direction
+
+            zlon1 = rad * pglam(ji,jj)
+            zlon2 = rad * pglam(ji,jj)
+            zlat1 = rad * ( pgphi(ji,jj) + ( 0.5_wp * pphiscl_deg ) )
+            zlat2 = rad * ( pgphi(ji,jj) - ( 0.5_wp * pphiscl_deg ) )
+            zdlon = zlon2 - zlon1
+            zdlat = zlat2 - zlat1
+
+            za1 = sin( zdlat/2.0_wp )
+            za2 = sin( zdlon/2.0_wp )
+            za = ( za1 * za1 ) + ( COS( zlat1 ) * COS( zlat2 ) * ( za2 * za2 ) )
+            zc = 2.0_wp * atan2( SQRT( za ), SQRT( 1.0_wp-za ) )
+            zd = ra * zc
+
+            IF ( zd > pphiscl_max ) pphiscl_max = zd
+
+         END DO
+      END DO
+         
+   END SUBROUTINE obs_deg2dist
+
+   SUBROUTINE obs_dist2corners(pglam_bl, pglam_br, pglam_tl, pglam_tr, &
+      &                        pgphi_bl, pgphi_br, pgphi_tl, pgphi_tr, &
+      &                        plam, pphi, pxvert, pyvert, pdist)
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_dist2corners  ***
+      !!
+      !! ** Purpose : Calculate distance from centre of obs footprint to the corners of a grid box
+      !!
+      !! ** Method  : Use great circle distance formulae.
+      !!              Order of corners is 1=topleft, 2=topright, 3=bottomright, 4=bottomleft
+      !!
+      !! ** Action  :
+      !!                   
+      !! References : 
+      !!
+      !! History :
+      !!        ! 14/01. M. Martin.
+      !!-----------------------------------------------------------------------
+      !! * Modules used
+      !! * Arguments
+         REAL(KIND=wp), INTENT(IN) :: &
+            &  pglam_bl, pglam_br, &       !lon at corners of grid box
+            &  pglam_tl, pglam_tr
+         REAL(KIND=wp), INTENT(IN) :: &
+            &  pgphi_bl, pgphi_br, &       !lat at corners of grid box
+            &  pgphi_tl, pgphi_tr
+         REAL(KIND=wp), INTENT(IN) :: &
+            &  pphi, plam                  !lat/lon of centre of obs footprint
+         REAL(KIND=wp), DIMENSION(4), INTENT(OUT) :: &
+            &  pxvert, pyvert              !x/y location (in metres relative to centre of obs footprint) of corners
+         REAL(KIND=wp), DIMENSION(4), INTENT(OUT) :: &
+            &  pdist                       !distance (in metres) of each corner relative to centre of obs footprint
+
+      !! * Local variables
+         INTEGER :: &
+            &  jvert                       !Counter for corners
+         REAL(KIND=wp) :: &
+            &  zphi, zlam                  !Local values for lon/lat of corners
+         REAL(KIND=wp) :: & 
+            &  za1, za2,  &                !For great circle distance calculations
+            &  zb1, zb2,  &
+            &  zc1, zc2
+         REAL(KIND=wp) :: &
+            &  zdist_centre_lat, &         !Distances in lat/lon directions (in metres)
+            &  zdist_centre_lon
+
+      !!-----------------------------------------------------------------------
+
+         ! Work out latitudinal and longitudinal distance from centre of 
+         ! obs fp to corners of grid box
+         DO jvert = 1, 4
+            SELECT CASE(jvert)
+            CASE(1)
+               zphi = pgphi_tl
+               zlam = pglam_tl
+            CASE(2)
+               zphi = pgphi_tr
+               zlam = pglam_tr
+            CASE(3)
+               zphi = pgphi_br
+               zlam = pglam_br
+            CASE(4)
+               zphi = pgphi_bl
+               zlam = pglam_bl
+            END SELECT
+            
+            IF (zlam == plam ) THEN
+               pxvert(jvert) = 0.0_wp
+            ELSE
+               za1 = SIN( zphi * rad )
+               za2 = SIN( zphi * rad )
+               zb1 = COS( zphi * rad ) * COS( zlam * rad )
+               zb2 = COS( zphi * rad ) * COS( plam  * rad )
+               zc1 = COS( zphi * rad ) * SIN( zlam * rad )
+               zc2 = COS( zphi * rad ) * SIN( plam  * rad )
+               pxvert(jvert) = grt_cir_dis( za1, za2, zb1, zb2, zc1, zc2 )
+               pxvert(jvert) =  ra * pxvert(jvert)
+               IF ( zlam < plam ) pxvert(jvert) = - pxvert(jvert)
+            ENDIF
+            
+            IF ( zphi == pphi ) THEN
+               pyvert(jvert) = 0.0_wp
+            ELSE
+               za1 = SIN( zphi * rad )
+               za2 = SIN( pphi * rad )
+               zb1 = COS( zphi * rad ) * COS( zlam * rad )
+               zb2 = COS( pphi * rad ) * COS( zlam * rad )
+               zc1 = COS( zphi * rad ) * SIN( zlam * rad )
+               zc2 = COS( pphi * rad ) * SIN( zlam * rad )
+               pyvert(jvert) = grt_cir_dis( za1, za2, zb1, zb2, zc1, zc2 )
+               pyvert(jvert) =  ra * pyvert(jvert)
+               IF ( zphi < pphi ) pyvert(jvert) = - pyvert(jvert)
+            ENDIF
+
+            !Calculate the distance of each vertex relative to centre of obs footprint
+            pdist(jvert) = SQRT( ( pxvert(jvert) * pxvert(jvert) ) + &
+               &                 ( pyvert(jvert) * pyvert(jvert) ) )
+
+         END DO
+
+      END SUBROUTINE obs_dist2corners
