Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 11138)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 11202)
@@ -6,6 +6,4 @@
    !!======================================================================
 
-   !!----------------------------------------------------------------------
-   !!   'key_diaobs' : Switch on the observation diagnostic computation
    !!----------------------------------------------------------------------
    !!   dia_obs_init : Reading and prepare observations
@@ -15,5 +13,5 @@
    !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
    !!----------------------------------------------------------------------
-   !! * Modules used   
+   !! * Modules used
    USE wrk_nemo                 ! Memory Allocation
    USE par_kind                 ! Precision variables
@@ -21,11 +19,7 @@
    USE par_oce
    USE dom_oce                  ! Ocean space and time domain variables
-   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch
-   USE obs_read_prof            ! Reading and allocation of observations (Coriolis)
-   USE obs_read_sla             ! Reading and allocation of SLA observations  
-   USE obs_read_sst             ! Reading and allocation of SST observations  
+   USE obs_read_prof            ! Reading and allocation of profile obs
+   USE obs_read_surf            ! Reading and allocation of surface obs
    USE obs_readmdt              ! Reading and allocation of MDT for SLA.
-   USE obs_read_seaice          ! Reading and allocation of Sea Ice observations  
-   USE obs_read_vel             ! Reading and allocation of velocity component observations
    USE obs_prep                 ! Preparation of obs. (grid search etc).
    USE obs_oper                 ! Observation operators
@@ -33,10 +27,7 @@
    USE obs_grid                 ! Grid searching
    USE obs_read_altbias         ! Bias treatment for altimeter
+   USE obs_sstbias              ! Bias correction routine for SST
    USE obs_profiles_def         ! Profile data definitions
-   USE obs_profiles             ! Profile data storage
    USE obs_surf_def             ! Surface data definitions
-   USE obs_sla                  ! SLA data storage
-   USE obs_sst                  ! SST data storage
-   USE obs_seaice               ! Sea Ice data storage
    USE obs_types                ! Definitions for observation types
    USE mpp_map                  ! MPP mapping
@@ -52,63 +43,62 @@
       &   dia_obs_dealloc  ! Deallocate dia_obs data
 
-   !! * Shared Module variables
-   LOGICAL, PUBLIC, PARAMETER :: &
-#if defined key_diaobs
-      & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics
-#else
-      & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics
-#endif
-
    !! * Module variables
-   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles
-   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles
-   LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set
-   LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set
-   LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles
-   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies 
-   LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files
-   LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files
-   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature
-   LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature
-   LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data
-   LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files
-   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration
-   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations
-   LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data 
-   LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data 
-   LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data 
-   LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data 
-   LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files
-   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
-   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
-   LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations
-   LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land
-   LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias  
-   LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files
-   LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations
-
-   REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS
-   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS
-  
-   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method
-   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method 
-
+   LOGICAL, PUBLIC :: &
+      &       lk_diaobs = .TRUE.   !: Include this for backwards compatibility at NEMO 3.6.
+   LOGICAL :: ln_diaobs            !: Logical switch for the obs operator
+   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs
+   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres
+   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_default_avglamscl !: Default E/W diameter of observation footprint
+   REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint
+   REAL(wp) :: rn_sla_avglamscl     !: E/W diameter of SLA observation footprint
+   REAL(wp) :: rn_sla_avgphiscl     !: N/S diameter of SLA observation footprint
+   REAL(wp) :: rn_sst_avglamscl     !: E/W diameter of SST observation footprint
+   REAL(wp) :: rn_sst_avgphiscl     !: N/S diameter of SST observation footprint
+   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint
+   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint
+   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint
+   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint
+
+   INTEGER :: nn_1dint         !: Vertical interpolation method
+   INTEGER :: nn_2dint_default !: Default horizontal interpolation method
+   INTEGER :: nn_2dint_sla     !: SLA horizontal interpolation method (-1 = default)
+   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default)
+   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default)
+   INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default)
+ 
    INTEGER, DIMENSION(imaxavtypes) :: &
-      & endailyavtypes !: ENACT data types which are daily average
-
-   INTEGER, PARAMETER :: MaxNumFiles = 1000
-   LOGICAL, DIMENSION(MaxNumFiles) :: &
-      & ln_profb_ena, & !: Is the feedback files from ENACT data ?
-   !                    !: If so use endailyavtypes
-      & ln_profb_enatim !: Change tim for 820 enact data set.
-   
-   LOGICAL, DIMENSION(MaxNumFiles) :: &
-      & ln_velfb_av   !: Is the velocity feedback files daily average?
+      & nn_profdavtypes      !: Profile data types representing a daily average
+   INTEGER :: nproftypes     !: Number of profile obs types
+   INTEGER :: nsurftypes     !: Number of surface obs types
+   INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      & nvarsprof, &         !: Number of profile variables
+      & nvarssurf            !: Number of surface variables
+   INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      & 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 :: &
-      & ld_enact     !: Profile data is ENACT so use endailyavtypes
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_velav     !: Velocity data is daily averaged
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_sstnight  !: SST observation corresponds to night mean
+      & 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(:) :: &
+      & surfdata, &          !: Initial surface data
+      & surfdataqc           !: Surface data after quality control
+   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
+      & profdata, &          !: Initial profile data
+      & profdataqc           !: Profile data after quality control
+
+   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
+      & cobstypesprof, &     !: Profile obs types
+      & cobstypessurf        !: Surface obs types
 
    !!----------------------------------------------------------------------
@@ -118,4 +108,6 @@
    !!----------------------------------------------------------------------
 
+   !! * Substitutions 
+#  include "domzgr_substitute.h90"
 CONTAINS
 
@@ -135,4 +127,6 @@
       !!        !  06-10  (A. Weaver) Cleaning and add controls
       !!        !  07-03  (K. Mogensen) General handling of profiles
+      !!        !  14-08  (J.While) Incorporated SST bias correction
+      !!        !  15-02  (M. Martin) Simplification of namelist and code
       !!----------------------------------------------------------------------
 
@@ -140,67 +134,154 @@
 
       !! * Local declarations
-      CHARACTER(len=128) :: enactfiles(MaxNumFiles)
-      CHARACTER(len=128) :: coriofiles(MaxNumFiles)
-      CHARACTER(len=128) :: profbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: sstfiles(MaxNumFiles)      
-      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 
-      CHARACTER(len=128) :: slafilesact(MaxNumFiles)      
-      CHARACTER(len=128) :: slafilespas(MaxNumFiles)      
-      CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: seaicefiles(MaxNumFiles)           
-      CHARACTER(len=128) :: velcurfiles(MaxNumFiles)  
-      CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)    
-      CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
-      CHARACTER(LEN=128) :: reysstname
-      CHARACTER(LEN=12)  :: reysstfmt
-      CHARACTER(LEN=128) :: bias_file
-      CHARACTER(LEN=20)  :: datestr=" ", timestr=" "
-      NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       &
-         &            ln_sla, ln_sladt, ln_slafb,                     &
-         &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       &
-         &            enactfiles, coriofiles, profbfiles,             &
-         &            slafilesact, slafilespas, slafbfiles,           &
-         &            sstfiles, sstfbfiles,                           &
-         &            ln_seaice, seaicefiles,                         &
-         &            dobsini, dobsend, n1dint, n2dint,               &
-         &            nmsshc, mdtcorr, mdtcutoff,                     &
-         &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    &
-         &            ln_sstnight,                                    &
-         &            ln_grid_search_lookup,                          &
-         &            grid_search_file, grid_search_res,              &
-         &            ln_grid_global, bias_file, ln_altbias,          &
-         &            endailyavtypes, ln_s_at_t, ln_profb_ena,        &
-         &            ln_vel3d, ln_velavcur, velavcurfiles,           &
-         &            ln_velhrcur, velhrcurfiles,                     &
-         &            ln_velavadcp, velavadcpfiles,                   &
-         &            ln_velhradcp, velhradcpfiles,                   &
-         &            ln_velfb, velfbfiles, ln_velfb_av,              &
-         &            ln_profb_enatim, ln_ignmis, ln_cl4
-
-      INTEGER :: jprofset
-      INTEGER :: jveloset
-      INTEGER :: jvar
-      INTEGER :: jnumenact
-      INTEGER :: jnumcorio
-      INTEGER :: jnumprofb
-      INTEGER :: jnumslaact
-      INTEGER :: jnumslapas
-      INTEGER :: jnumslafb
-      INTEGER :: jnumsst
-      INTEGER :: jnumsstfb
-      INTEGER :: jnumseaice
-      INTEGER :: jnumvelavcur
-      INTEGER :: jnumvelhrcur  
-      INTEGER :: jnumvelavadcp
-      INTEGER :: jnumvelhradcp   
-      INTEGER :: jnumvelfb
-      INTEGER :: ji
-      INTEGER :: jset
-      INTEGER :: ios                 ! Local integer output status for namelist read
-      LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
+      INTEGER, PARAMETER :: &
+         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+         & ifilesprof, &         ! Number of profile files
+         & ifilessurf            ! Number of surface files
+      INTEGER :: ios             ! Local integer output status for namelist read
+      INTEGER :: jtype           ! Counter for obs types
+      INTEGER :: jvar            ! Counter for variables
+      INTEGER :: jfile           ! Counter for files
+      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply
+      INTEGER :: n2dint_type     ! Local version of nn_2dint*
+
+      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
+         & cn_profbfiles,      & ! T/S profile input filenames
+         & cn_sstfbfiles,      & ! Sea surface temperature input filenames
+         & cn_slafbfiles,      & ! Sea level anomaly input filenames
+         & cn_sicfbfiles,      & ! Seaice concentration input filenames
+         & cn_velfbfiles,      & ! Velocity profile input filenames
+         & cn_sssfbfiles,      & ! Sea surface salinity input filenames
+         & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames
+         & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames
+         & cn_slchlnonfbfiles, & ! Surface non-diatom         log10(chlorophyll) input filenames
+         & cn_slchldinfbfiles, & ! Surface dinoflagellate     log10(chlorophyll) input filenames
+         & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames
+         & cn_slchlnanfbfiles, & ! Surface nanophytoplankton  log10(chlorophyll) input filenames
+         & cn_slchlpicfbfiles, & ! Surface picophytoplankton  log10(chlorophyll) input filenames
+         & cn_schltotfbfiles,  & ! Surface total              chlorophyll        input filenames
+         & cn_slphytotfbfiles, & ! Surface total      log10(phytoplankton carbon) input filenames
+         & cn_slphydiafbfiles, & ! Surface diatom     log10(phytoplankton carbon) input filenames
+         & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames
+         & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames
+         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames
+         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames
+         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames
+         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames
+         & cn_pno3fbfiles,     & ! Profile nitrate input filenames
+         & cn_psi4fbfiles,     & ! Profile silicate input filenames
+         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames
+         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames
+         & cn_palkfbfiles,     & ! Profile alkalinity input filenames
+         & cn_pphfbfiles,      & ! Profile pH input filenames
+         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames
+         & cn_sstbiasfiles       ! SST bias input filenames
+
+      CHARACTER(LEN=128) :: &
+         & cn_altbiasfile        ! Altimeter bias input filename
+
+
+      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
+      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
+      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies 
+      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
+      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
+      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs
+      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
+      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs
+      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs
+      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs
+      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs
+      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs
+      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs
+      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs
+      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs
+      LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs
+      LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs
+      LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs
+      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs
+      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs
+      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs
+      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs
+      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs
+      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs
+      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs
+      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs
+      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs
+      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs
+      LOGICAL :: ln_pph          ! Logical switch for profile pH obs
+      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs
+      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
+      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
+      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST
+      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
+      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
+
+      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl
+      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl
+
+      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
+         & clproffiles, &        ! Profile filenames
+         & clsurffiles           ! Surface filenames
+
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read
+      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
+      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables)
+
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zglam                 ! Model longitudes for profile variables
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zgphi                 ! Model latitudes for profile variables
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
+         & zmask                 ! Model land/sea mask associated with variables
+
+
+      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
+         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
+         &            ln_slchltot, ln_slchldia, ln_slchlnon,          &
+         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          &
+         &            ln_slchlpic, ln_schltot,                        &
+         &            ln_slphytot, ln_slphydia, ln_slphynon,          &
+         &            ln_sspm,     ln_sfco2,    ln_spco2,             &
+         &            ln_plchltot, ln_pchltot,  ln_pno3,              &
+         &            ln_psi4,     ln_ppo4,     ln_pdic,              &
+         &            ln_palk,     ln_pph,      ln_po2,               &
+         &            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_default_fp_indegs,              &
+         &            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_sssfbfiles,                   &
+         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         &
+         &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         &
+         &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         &
+         &            cn_slchlpicfbfiles, cn_schltotfbfiles,          &
+         &            cn_slphytotfbfiles, cn_slphydiafbfiles,         &
+         &            cn_slphynonfbfiles, cn_sspmfbfiles,             &
+         &            cn_sfco2fbfiles, cn_spco2fbfiles,               &
+         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          &
+         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, &
+         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  &
+         &            cn_po2fbfiles,                                  &
+         &            cn_sstbiasfiles, cn_altbiasfile,                &
+         &            cn_gridsearchfile, rn_gridsearchres,            &
+         &            rn_dobsini, rn_dobsend,                         &
+         &            rn_default_avglamscl, rn_default_avgphiscl,     &
+         &            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_default,                     &
+         &            nn_2dint_sla, nn_2dint_sst,                     &
+         &            nn_2dint_sss, nn_2dint_sic,                     &
+         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
+         &            nn_profdavtypes
 
       !-----------------------------------------------------------------------
@@ -208,116 +289,61 @@
       !-----------------------------------------------------------------------
 
-      enactfiles(:) = ''
-      coriofiles(:) = ''
-      profbfiles(:) = ''
-      slafilesact(:) = ''
-      slafilespas(:) = ''
-      slafbfiles(:) = ''
-      sstfiles(:)   = ''
-      sstfbfiles(:) = ''
-      seaicefiles(:) = ''
-      velcurfiles(:) = ''
-      veladcpfiles(:) = ''
-      velavcurfiles(:) = ''
-      velhrcurfiles(:) = ''
-      velavadcpfiles(:) = ''
-      velhradcpfiles(:) = ''
-      velfbfiles(:) = ''
-      velcurfiles(:) = ''
-      veladcpfiles(:) = ''
-      endailyavtypes(:) = -1
-      endailyavtypes(1) = 820
-      ln_profb_ena(:) = .FALSE.
-      ln_profb_enatim(:) = .TRUE.
-      ln_velfb_av(:) = .FALSE.
-      ln_ignmis = .FALSE.
-      
-      CALL ini_date( dobsini )
-      CALL fin_date( dobsend )
- 
-      ! Read Namelist namobs : control observation diagnostics
-      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation
+      ! Some namelist arrays need initialising
+      cn_profbfiles(:)      = ''
+      cn_slafbfiles(:)      = ''
+      cn_sstfbfiles(:)      = ''
+      cn_sicfbfiles(:)      = ''
+      cn_velfbfiles(:)      = ''
+      cn_sssfbfiles(:)      = ''
+      cn_slchltotfbfiles(:) = ''
+      cn_slchldiafbfiles(:) = ''
+      cn_slchlnonfbfiles(:) = ''
+      cn_slchldinfbfiles(:) = ''
+      cn_slchlmicfbfiles(:) = ''
+      cn_slchlnanfbfiles(:) = ''
+      cn_slchlpicfbfiles(:) = ''
+      cn_schltotfbfiles(:)  = ''
+      cn_slphytotfbfiles(:) = ''
+      cn_slphydiafbfiles(:) = ''
+      cn_slphynonfbfiles(:) = ''
+      cn_sspmfbfiles(:)     = ''
+      cn_sfco2fbfiles(:)    = ''
+      cn_spco2fbfiles(:)    = ''
+      cn_plchltotfbfiles(:) = ''
+      cn_pchltotfbfiles(:)  = ''
+      cn_pno3fbfiles(:)     = ''
+      cn_psi4fbfiles(:)     = ''
+      cn_ppo4fbfiles(:)     = ''
+      cn_pdicfbfiles(:)     = ''
+      cn_palkfbfiles(:)     = ''
+      cn_pphfbfiles(:)      = ''
+      cn_po2fbfiles(:)      = ''
+      cn_sstbiasfiles(:)    = ''
+      nn_profdavtypes(:)    = -1
+
+      CALL ini_date( rn_dobsini )
+      CALL fin_date( rn_dobsend )
+
+      ! Read namelist namobs : control observation diagnostics
+      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
       READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
 
-      REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation
+      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
       READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namobs )
 
-      ! Count number of files for each type
-      IF (ln_ena) THEN
-         lmask(:) = .FALSE.
-         WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
-         jnumenact = COUNT(lmask)
+      lk_diaobs = .FALSE.
+#if defined key_diaobs
+      IF ( ln_diaobs ) lk_diaobs = .TRUE.
+#endif
+
+      IF ( .NOT. lk_diaobs ) THEN
+         IF(lwp) WRITE(numout,cform_war)
+         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs'
+         RETURN
       ENDIF
-      IF (ln_cor) THEN
-         lmask(:) = .FALSE.
-         WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
-         jnumcorio = COUNT(lmask)
-      ENDIF
-      IF (ln_profb) THEN
-         lmask(:) = .FALSE.
-         WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumprofb = COUNT(lmask)
-      ENDIF
-      IF (ln_sladt) THEN
-         lmask(:) = .FALSE.
-         WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
-         jnumslaact = COUNT(lmask)
-         lmask(:) = .FALSE.
-         WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
-         jnumslapas = COUNT(lmask)
-      ENDIF
-      IF (ln_slafb) THEN
-         lmask(:) = .FALSE.
-         WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumslafb = COUNT(lmask)
-         lmask(:) = .FALSE.
-      ENDIF
-      IF (ln_ghrsst) THEN
-         lmask(:) = .FALSE.
-         WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
-         jnumsst = COUNT(lmask)
-      ENDIF      
-      IF (ln_sstfb) THEN
-         lmask(:) = .FALSE.
-         WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumsstfb = COUNT(lmask)
-         lmask(:) = .FALSE.
-      ENDIF
-      IF (ln_seaice) THEN
-         lmask(:) = .FALSE.
-         WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
-         jnumseaice = COUNT(lmask)
-      ENDIF
-      IF (ln_velavcur) THEN
-         lmask(:) = .FALSE.
-         WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelavcur = COUNT(lmask)
-      ENDIF
-      IF (ln_velhrcur) THEN
-         lmask(:) = .FALSE.
-         WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelhrcur = COUNT(lmask)
-      ENDIF
-      IF (ln_velavadcp) THEN
-         lmask(:) = .FALSE.
-         WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelavadcp = COUNT(lmask)
-      ENDIF
-      IF (ln_velhradcp) THEN
-         lmask(:) = .FALSE.
-         WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelhradcp = COUNT(lmask)
-      ENDIF
-      IF (ln_velfb) THEN
-         lmask(:) = .FALSE.
-         WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelfb = COUNT(lmask)
-         lmask(:) = .FALSE.
-      ENDIF
-      
-      ! Control print
+
       IF(lwp) THEN
          WRITE(numout,*)
@@ -325,138 +351,334 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
-         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d
-         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d
-         WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena
-         WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor
-         WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb
-         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla
-         WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt
-         WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb
-         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh
-         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst
-         WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst    
-         WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst
-         WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb
-         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight
-         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss
-         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice
-         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d
-         WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur
-         WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur
-         WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp
-         WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp
-         WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb
-         WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
-         WRITE(numout,*) &
-   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup
+         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
+         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
+         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
+         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
+         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,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss
+         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot
+         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia
+         WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon
+         WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin
+         WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic
+         WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan
+         WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic
+         WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot
+         WRITE(numout,*) '             Logical switch for surface total log(phyC) obs      ln_slphytot = ', ln_slphytot
+         WRITE(numout,*) '             Logical switch for surface diatom log(phyC) obs     ln_slphydia = ', ln_slphydia
+         WRITE(numout,*) '             Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon
+         WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm
+         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2
+         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2
+         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot
+         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot
+         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3
+         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4
+         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4
+         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic
+         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk
+         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph
+         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2
+         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       grid_search_file = ', grid_search_file
-         IF (ln_ena) THEN
-            DO ji = 1, jnumenact
-               WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', &
-                  TRIM(enactfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_cor) THEN
-            DO ji = 1, jnumcorio
-               WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', &
-                  TRIM(coriofiles(ji))
-            END DO
-         ENDIF
-         IF (ln_profb) THEN
-            DO ji = 1, jnumprofb
-               IF (ln_profb_ena(ji)) THEN
-                  WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', &
-                     TRIM(profbfiles(ji))
+            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
+         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
+         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
+         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
+         WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default
+         WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla
+         WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst
+         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss
+         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic
+         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl
+         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl
+         WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs
+         WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl
+         WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl
+         WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs
+         WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl
+         WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl
+         WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs
+         WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl
+         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl
+         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs
+         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
+      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, ln_plchltot,          &
+         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     &
+         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) )
+      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     &
+         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, &
+         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  &
+         &                  ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm,     &
+         &                  ln_sfco2,    ln_spco2 /) )
+
+      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
+            cobstypesprof(jtype) = 'prof'
+            clproffiles(jtype,:) = cn_profbfiles
+         ENDIF
+         IF (ln_vel3d) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) =  'vel'
+            clproffiles(jtype,:) = cn_velfbfiles
+         ENDIF
+         IF (ln_plchltot) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'plchltot'
+            clproffiles(jtype,:) = cn_plchltotfbfiles
+         ENDIF
+         IF (ln_pchltot) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'pchltot'
+            clproffiles(jtype,:) = cn_pchltotfbfiles
+         ENDIF
+         IF (ln_pno3) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'pno3'
+            clproffiles(jtype,:) = cn_pno3fbfiles
+         ENDIF
+         IF (ln_psi4) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'psi4'
+            clproffiles(jtype,:) = cn_psi4fbfiles
+         ENDIF
+         IF (ln_ppo4) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'ppo4'
+            clproffiles(jtype,:) = cn_ppo4fbfiles
+         ENDIF
+         IF (ln_pdic) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'pdic'
+            clproffiles(jtype,:) = cn_pdicfbfiles
+         ENDIF
+         IF (ln_palk) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'palk'
+            clproffiles(jtype,:) = cn_palkfbfiles
+         ENDIF
+         IF (ln_pph) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'pph'
+            clproffiles(jtype,:) = cn_pphfbfiles
+         ENDIF
+         IF (ln_po2) THEN
+            jtype = jtype + 1
+            cobstypesprof(jtype) = 'po2'
+            clproffiles(jtype,:) = cn_po2fbfiles
+         ENDIF
+
+         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles )
+
+      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
+            cobstypessurf(jtype) = 'sla'
+            clsurffiles(jtype,:) = cn_slafbfiles
+         ENDIF
+         IF (ln_sst) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'sst'
+            clsurffiles(jtype,:) = cn_sstfbfiles
+         ENDIF
+         IF (ln_sic) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'sic'
+            clsurffiles(jtype,:) = cn_sicfbfiles
+         ENDIF
+         IF (ln_sss) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'sss'
+            clsurffiles(jtype,:) = cn_sssfbfiles
+         ENDIF
+         IF (ln_slchltot) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slchltot'
+            clsurffiles(jtype,:) = cn_slchltotfbfiles
+         ENDIF
+         IF (ln_slchldia) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slchldia'
+            clsurffiles(jtype,:) = cn_slchldiafbfiles
+         ENDIF
+         IF (ln_slchlnon) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slchlnon'
+            clsurffiles(jtype,:) = cn_slchlnonfbfiles
+         ENDIF
+         IF (ln_slchldin) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slchldin'
+            clsurffiles(jtype,:) = cn_slchldinfbfiles
+         ENDIF
+         IF (ln_slchlmic) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slchlmic'
+            clsurffiles(jtype,:) = cn_slchlmicfbfiles
+         ENDIF
+         IF (ln_slchlnan) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slchlnan'
+            clsurffiles(jtype,:) = cn_slchlnanfbfiles
+         ENDIF
+         IF (ln_slchlpic) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slchlpic'
+            clsurffiles(jtype,:) = cn_slchlpicfbfiles
+         ENDIF
+         IF (ln_schltot) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'schltot'
+            clsurffiles(jtype,:) = cn_schltotfbfiles
+         ENDIF
+         IF (ln_slphytot) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slphytot'
+            clsurffiles(jtype,:) = cn_slphytotfbfiles
+         ENDIF
+         IF (ln_slphydia) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slphydia'
+            clsurffiles(jtype,:) = cn_slphydiafbfiles
+         ENDIF
+         IF (ln_slphynon) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'slphynon'
+            clsurffiles(jtype,:) = cn_slphynonfbfiles
+         ENDIF
+         IF (ln_sspm) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'sspm'
+            clsurffiles(jtype,:) = cn_sspmfbfiles
+         ENDIF
+         IF (ln_sfco2) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'sfco2'
+            clsurffiles(jtype,:) = cn_sfco2fbfiles
+         ENDIF
+         IF (ln_spco2) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'spco2'
+            clsurffiles(jtype,:) = cn_spco2fbfiles
+         ENDIF
+
+         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles )
+
+         DO jtype = 1, nsurftypes
+
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
+               IF ( nn_2dint_sla == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
                ELSE
-                  WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', &
-                     TRIM(profbfiles(ji))
+                  n2dint_type  = nn_2dint_sla
                ENDIF
-               WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji)
-            END DO
-         ENDIF
-         IF (ln_sladt) THEN
-            DO ji = 1, jnumslaact
-               WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', &
-                  TRIM(slafilesact(ji))
-            END DO
-            DO ji = 1, jnumslapas
-               WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', &
-                  TRIM(slafilespas(ji))
-            END DO
-         ENDIF
-         IF (ln_slafb) THEN
-            DO ji = 1, jnumslafb
-               WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', &
-                  TRIM(slafbfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_ghrsst) THEN
-            DO ji = 1, jnumsst
-               WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', &
-                  TRIM(sstfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_sstfb) THEN
-            DO ji = 1, jnumsstfb
-               WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', &
-                  TRIM(sstfbfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_seaice) THEN
-            DO ji = 1, jnumseaice
-               WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', &
-                  TRIM(seaicefiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velavcur) THEN
-            DO ji = 1, jnumvelavcur
-               WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', &
-                  TRIM(velavcurfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velhrcur) THEN
-            DO ji = 1, jnumvelhrcur
-               WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', &
-                  TRIM(velhrcurfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velavadcp) THEN
-            DO ji = 1, jnumvelavadcp
-               WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', &
-                  TRIM(velavadcpfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velhradcp) THEN
-            DO ji = 1, jnumvelhradcp
-               WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', &
-                  TRIM(velhradcpfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velfb) THEN
-            DO ji = 1, jnumvelfb
-               IF (ln_velfb_av(ji)) THEN
-                  WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', &
-                     TRIM(velfbfiles(ji))
+               ztype_avglamscl = rn_sla_avglamscl
+               ztype_avgphiscl = rn_sla_avgphiscl
+               ltype_fp_indegs = ln_sla_fp_indegs
+               ltype_night     = .FALSE.
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
+               IF ( nn_2dint_sst == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
                ELSE
-                  WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', &
-                     TRIM(velfbfiles(ji))
+                  n2dint_type  = nn_2dint_sst
                ENDIF
-            END DO
-         ENDIF
-         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini
-         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend
-         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint
-         WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint
-         WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea
-         WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc
-         WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr
-         WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff
-         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias
-         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis
-         WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes
+               ztype_avglamscl = rn_sst_avglamscl
+               ztype_avgphiscl = rn_sst_avgphiscl
+               ltype_fp_indegs = ln_sst_fp_indegs
+               ltype_night     = ln_sstnight
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
+               IF ( nn_2dint_sic == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
+               ELSE
+                  n2dint_type  = nn_2dint_sic
+               ENDIF
+               ztype_avglamscl = rn_sic_avglamscl
+               ztype_avgphiscl = rn_sic_avgphiscl
+               ltype_fp_indegs = ln_sic_fp_indegs
+               ltype_night     = .FALSE.
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
+               IF ( nn_2dint_sss == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
+               ELSE
+                  n2dint_type  = nn_2dint_sss
+               ENDIF
+               ztype_avglamscl = rn_sss_avglamscl
+               ztype_avgphiscl = rn_sss_avgphiscl
+               ltype_fp_indegs = ln_sss_fp_indegs
+               ltype_night     = .FALSE.
+            ELSE
+               n2dint_type     = nn_2dint_default
+               ztype_avglamscl = rn_default_avglamscl
+               ztype_avgphiscl = rn_default_avgphiscl
+               ltype_fp_indegs = ln_default_fp_indegs
+               ltype_night     = .FALSE.
+            ENDIF
+            
+            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), &
+               &                    nn_2dint_default, n2dint_type,                 &
+               &                    ztype_avglamscl, ztype_avgphiscl,              &
+               &                    ltype_fp_indegs, ltype_night,                  &
+               &                    n2dintsurf, ravglamscl, ravgphiscl,            &
+               &                    lfpindegs, llnightav )
+
+         END DO
 
       ENDIF
-      
+
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+
+
+      !-----------------------------------------------------------------------
+      ! Obs operator parameter checking and initialisations
+      !-----------------------------------------------------------------------
+
       IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
          CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
@@ -464,505 +686,154 @@
       ENDIF
 
-      CALL obs_typ_init
-      
-      CALL mppmap_init
-      
-      ! Parameter control
-#if defined key_diaobs
-      IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
-         & ( .NOT. ln_vel3d ).AND.                                         &
-         & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
-         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN
-         IF(lwp) WRITE(numout,cform_war)
-         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
-            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.'
-         nwarn = nwarn + 1
-      ENDIF
-#endif
-
-      CALL obs_grid_setup( )
-      IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
+      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
          CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
             &                    ' is not available')
       ENDIF
-      IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
-         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
+
+      IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN
+         CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', &
             &                    ' is not available')
       ENDIF
 
+      CALL obs_typ_init
+
+      CALL mppmap_init
+
+      CALL obs_grid_setup( )
+
       !-----------------------------------------------------------------------
       ! Depending on switches read the various observation types
       !-----------------------------------------------------------------------
-      !  - Temperature/salinity profiles
-
-      IF ( ln_t3d .OR. ln_s3d ) THEN
-
-         ! Set the number of variables for profiles to 2 (T and S)
-         nprofvars = 2
-         ! Set the number of extra variables for profiles to 1 (insitu temp).
-         nprofextr = 1
-
-         ! Count how may insitu data sets we have and allocate data.
-         jprofset = 0
-         IF ( ln_ena ) jprofset = jprofset + 1
-         IF ( ln_cor ) jprofset = jprofset + 1
-         IF ( ln_profb ) jprofset = jprofset + jnumprofb
-         nprofsets = jprofset
-         IF ( nprofsets > 0 ) THEN
-            ALLOCATE(ld_enact(nprofsets))
-            ALLOCATE(profdata(nprofsets))
-            ALLOCATE(prodatqc(nprofsets))
-         ENDIF
-
-         jprofset = 0
-          
-         ! ENACT insitu data
-
-         IF ( ln_ena ) THEN
-
-            jprofset = jprofset + 1
+
+      IF ( nproftypes > 0 ) THEN
+
+         ALLOCATE(profdata(nproftypes))
+         ALLOCATE(profdataqc(nproftypes))
+         ALLOCATE(nvarsprof(nproftypes))
+         ALLOCATE(nextrprof(nproftypes))
+
+         DO jtype = 1, nproftypes
+
+            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
+               nvarsprof(jtype) = 2
+               nextrprof(jtype) = 1
+               ALLOCATE(llvar(nvarsprof(jtype)))
+               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
+               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
+               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
+               llvar(1)       = ln_t3d
+               llvar(2)       = ln_s3d
+               zglam(:,:,1)   = glamt(:,:)
+               zglam(:,:,2)   = glamt(:,:)
+               zgphi(:,:,1)   = gphit(:,:)
+               zgphi(:,:,2)   = gphit(:,:)
+               zmask(:,:,:,1) = tmask(:,:,:)
+               zmask(:,:,:,2) = tmask(:,:,:)
+            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
+               nvarsprof(jtype) = 2
+               nextrprof(jtype) = 2
+               ALLOCATE(llvar(nvarsprof(jtype)))
+               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
+               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
+               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
+               llvar(1)       = ln_vel3d
+               llvar(2)       = ln_vel3d
+               zglam(:,:,1)   = glamu(:,:)
+               zglam(:,:,2)   = glamv(:,:)
+               zgphi(:,:,1)   = gphiu(:,:)
+               zgphi(:,:,2)   = gphiv(:,:)
+               zmask(:,:,:,1) = umask(:,:,:)
+               zmask(:,:,:,2) = vmask(:,:,:)
+            ELSE
+               nvarsprof(jtype) = 1
+               nextrprof(jtype) = 0
+               ALLOCATE(llvar(nvarsprof(jtype)))
+               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
+               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
+               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
+               llvar(1)       = .TRUE.
+               zglam(:,:,1)   = glamt(:,:)
+               zgphi(:,:,1)   = gphit(:,:)
+               zmask(:,:,:,1) = tmask(:,:,:)
+            ENDIF
+
+            !Read in profile or profile obs types
+            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       &
+               &               clproffiles(jtype,1:ifilesprof(jtype)), &
+               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
+               &               rn_dobsini, rn_dobsend, llvar, &
+               &               ln_ignmis, ln_s_at_t, .FALSE., &
+               &               kdailyavtypes = nn_profdavtypes )
+
+            DO jvar = 1, nvarsprof(jtype)
+               CALL obs_prof_staend( profdata(jtype), jvar )
+            END DO
+
+            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
+               &               llvar, &
+               &               jpi, jpj, jpk, &
+               &               zmask, zglam, zgphi,  &
+               &               ln_nea, ln_bound_reject, &
+               &               kdailyavtypes = nn_profdavtypes )
             
-            ld_enact(jprofset) = .TRUE.
-
-            CALL obs_rea_pro_dri( 1, profdata(jprofset),          &
-               &                  jnumenact, enactfiles(1:jnumenact), &
-               &                  nprofvars, nprofextr,        &
-               &                  nitend-nit000+2,             &
-               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
-               &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
-               &                  kdailyavtypes = endailyavtypes )
-
-            DO jvar = 1, 2
-
-               CALL obs_prof_staend( profdata(jprofset), jvar )
-
-            END DO
-
-            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-               &              ln_t3d, ln_s3d, ln_nea, &
-               &              kdailyavtypes=endailyavtypes )
-            
-         ENDIF
-
-         ! Coriolis insitu data
-
-         IF ( ln_cor ) THEN
-           
-            jprofset = jprofset + 1
-
-            ld_enact(jprofset) = .FALSE.
-
-            CALL obs_rea_pro_dri( 2, profdata(jprofset),          &
-               &                  jnumcorio, coriofiles(1:jnumcorio), &
-               &                  nprofvars, nprofextr,        &
-               &                  nitend-nit000+2,             &
-               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
-               &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
-
-            DO jvar = 1, 2
-
-               CALL obs_prof_staend( profdata(jprofset), jvar )
-
-            END DO
-
-            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-                 &            ln_t3d, ln_s3d, ln_nea )
-            
-         ENDIF
- 
-         ! Feedback insitu data
-
-         IF ( ln_profb ) THEN
-           
-            DO jset = 1, jnumprofb
-               
-               jprofset = jprofset + 1
-               ld_enact (jprofset) = ln_profb_ena(jset)
-
-               CALL obs_rea_pro_dri( 0, profdata(jprofset),          &
-                  &                  1, profbfiles(jset:jset), &
-                  &                  nprofvars, nprofextr,        &
-                  &                  nitend-nit000+2,             &
-                  &                  dobsini, dobsend, ln_t3d, ln_s3d, &
-                  &                  ln_ignmis, ln_s_at_t, &
-                  &                  ld_enact(jprofset).AND.&
-                  &                  ln_profb_enatim(jset), &
-                  &                  .FALSE., kdailyavtypes = endailyavtypes )
-               
-               DO jvar = 1, 2
-                  
-                  CALL obs_prof_staend( profdata(jprofset), jvar )
-                  
+            DEALLOCATE( llvar )
+            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zglam )
+            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zgphi )
+            CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
+
+         END DO
+
+         DEALLOCATE( ifilesprof, clproffiles )
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         ALLOCATE(surfdata(nsurftypes))
+         ALLOCATE(surfdataqc(nsurftypes))
+         ALLOCATE(nvarssurf(nsurftypes))
+         ALLOCATE(nextrsurf(nsurftypes))
+
+         DO jtype = 1, nsurftypes
+
+            nvarssurf(jtype) = 1
+            nextrsurf(jtype) = 0
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
+
+            !Read in surface obs types
+            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
+               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
+               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
+               &               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), 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
+               jnumsstbias = 0
+               DO jfile = 1, jpmaxnfiles
+                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
+                     &  jnumsstbias = jnumsstbias + 1
                END DO
-               
-               IF ( ld_enact(jprofset) ) THEN
-                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-                     &              ln_t3d, ln_s3d, ln_nea, &
-                     &              kdailyavtypes = endailyavtypes )
-               ELSE
-                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-                     &              ln_t3d, ln_s3d, ln_nea )
+               IF ( jnumsstbias == 0 ) THEN
+                  CALL ctl_stop("ln_sstbias set but no bias files to read in")    
                ENDIF
-               
-            END DO
-
-         ENDIF
+
+               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
+                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
+
+            ENDIF
+
+         END DO
+
+         DEALLOCATE( ifilessurf, clsurffiles )
 
       ENDIF
 
-      !  - Sea level anomalies
-      IF ( ln_sla ) THEN
-        ! Set the number of variables for sla to 1
-         nslavars = 1
-
-         ! Set the number of extra variables for sla to 2
-         nslaextr = 2
-         
-         ! Set the number of sla data sets to 2
-         nslasets = 0
-         IF ( ln_sladt ) THEN
-            nslasets = nslasets + 2
-         ENDIF
-         IF ( ln_slafb ) THEN
-            nslasets = nslasets + jnumslafb
-         ENDIF
-         
-         ALLOCATE(sladata(nslasets))
-         ALLOCATE(sladatqc(nslasets))
-         sladata(:)%nsurf=0
-         sladatqc(:)%nsurf=0
-
-         nslasets = 0
-
-         ! AVISO SLA data
-
-         IF ( ln_sladt ) THEN
-
-            ! Active SLA observations
-            
-            nslasets = nslasets + 1
-            
-            CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
-               &              slafilesact(1:jnumslaact), &
-               &              nslavars, nslaextr, nitend-nit000+2, &
-               &              dobsini, dobsend, ln_ignmis, .FALSE. )
-            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
-               &              ln_sla, ln_nea )
-            
-            ! Passive SLA observations
-            
-            nslasets = nslasets + 1
-            
-            CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
-               &              slafilespas(1:jnumslapas), &
-               &              nslavars, nslaextr, nitend-nit000+2, &
-               &              dobsini, dobsend, ln_ignmis, .FALSE. )
-            
-            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
-               &              ln_sla, ln_nea )
-
-         ENDIF
-         
-         ! Feedback SLA data
-
-         IF ( ln_slafb ) THEN
-
-            DO jset = 1, jnumslafb
-            
-               nslasets = nslasets + 1
-            
-               CALL obs_rea_sla( 0, sladata(nslasets), 1, &
-                  &              slafbfiles(jset:jset), &
-                  &              nslavars, nslaextr, nitend-nit000+2, &
-                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
-               CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
-                  &              ln_sla, ln_nea )
-
-            END DO               
-
-         ENDIF
-         
-         CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
-            
-         ! read in altimeter bias
-         
-         IF ( ln_altbias ) THEN     
-            CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
-         ENDIF
-     
-      ENDIF
-
-      !  - Sea surface height
-      IF ( ln_ssh ) THEN
-         IF(lwp) WRITE(numout,*) ' SSH currently not available'
-      ENDIF
-
-      !  - Sea surface temperature
-      IF ( ln_sst ) THEN
-
-         ! Set the number of variables for sst to 1
-         nsstvars = 1
-
-         ! Set the number of extra variables for sst to 0
-         nsstextr = 0
-
-         nsstsets = 0
-
-         IF (ln_reysst) nsstsets = nsstsets + 1
-         IF (ln_ghrsst) nsstsets = nsstsets + 1
-         IF ( ln_sstfb ) THEN
-            nsstsets = nsstsets + jnumsstfb
-         ENDIF
-
-         ALLOCATE(sstdata(nsstsets))
-         ALLOCATE(sstdatqc(nsstsets))
-         ALLOCATE(ld_sstnight(nsstsets))
-         sstdata(:)%nsurf=0
-         sstdatqc(:)%nsurf=0    
-         ld_sstnight(:)=.false.
-
-         nsstsets = 0
-
-         IF (ln_reysst) THEN
-
-            nsstsets = nsstsets + 1
-
-            ld_sstnight(nsstsets) = ln_sstnight
-
-            CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
-               &                  nsstvars, nsstextr, &
-               &                  nitend-nit000+2, dobsini, dobsend )
-            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
-               &              ln_nea )
-
-        ENDIF
-        
-        IF (ln_ghrsst) THEN
-        
-            nsstsets = nsstsets + 1
-
-            ld_sstnight(nsstsets) = ln_sstnight
-          
-            CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
-               &              sstfiles(1:jnumsst), &
-               &              nsstvars, nsstextr, nitend-nit000+2, &
-               &              dobsini, dobsend, ln_ignmis, .FALSE. )
-            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
-               &              ln_nea )
-
-        ENDIF
-               
-         ! Feedback SST data
-
-         IF ( ln_sstfb ) THEN
-
-            DO jset = 1, jnumsstfb
-            
-               nsstsets = nsstsets + 1
-
-               ld_sstnight(nsstsets) = ln_sstnight
-            
-               CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
-                  &              sstfbfiles(jset:jset), &
-                  &              nsstvars, nsstextr, nitend-nit000+2, &
-                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
-               CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
-                  &              ln_sst, ln_nea )
-
-            END DO               
-
-         ENDIF
-
-      ENDIF
-
-      !  - Sea surface salinity
-      IF ( ln_sss ) THEN
-         IF(lwp) WRITE(numout,*) ' SSS currently not available'
-      ENDIF
-
-      !  - Sea Ice Concentration
-      
-      IF ( ln_seaice ) THEN
-
-         ! Set the number of variables for seaice to 1
-         nseaicevars = 1
-
-         ! Set the number of extra variables for seaice to 0
-         nseaiceextr = 0
-         
-         ! Set the number of data sets to 1
-         nseaicesets = 1
-
-         ALLOCATE(seaicedata(nseaicesets))
-         ALLOCATE(seaicedatqc(nseaicesets))
-         seaicedata(:)%nsurf=0
-         seaicedatqc(:)%nsurf=0
-
-         CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
-            &                 seaicefiles(1:jnumseaice), &
-            &                 nseaicevars, nseaiceextr, nitend-nit000+2, &
-            &                 dobsini, dobsend, ln_ignmis, .FALSE. )
-
-         CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
-            &                 ln_seaice, ln_nea )
- 
-      ENDIF
-
-      IF (ln_vel3d) THEN
-
-         ! Set the number of variables for profiles to 2 (U and V)
-         nvelovars = 2
-
-         ! Set the number of extra variables for profiles to 2 to store 
-         ! rotation parameters
-         nveloextr = 2
-
-         jveloset = 0
-         
-         IF ( ln_velavcur ) jveloset = jveloset + 1
-         IF ( ln_velhrcur ) jveloset = jveloset + 1
-         IF ( ln_velavadcp ) jveloset = jveloset + 1
-         IF ( ln_velhradcp ) jveloset = jveloset + 1
-         IF (ln_velfb) jveloset = jveloset + jnumvelfb
-
-         nvelosets = jveloset
-         IF ( nvelosets > 0 ) THEN
-            ALLOCATE( velodata(nvelosets) )
-            ALLOCATE( veldatqc(nvelosets) )
-            ALLOCATE( ld_velav(nvelosets) )
-         ENDIF
-         
-         jveloset = 0
-         
-         ! Daily averaged data
-
-         IF ( ln_velavcur ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .TRUE.
-            
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
-               &                  velavcurfiles(1:jnumvelavcur), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhrcur ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .FALSE.
-               
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
-               &                  velhrcurfiles(1:jnumvelhrcur), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         ! Daily averaged data
-
-         IF ( ln_velavadcp ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .TRUE.
-            
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
-               &                  velavadcpfiles(1:jnumvelavadcp), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhradcp ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .FALSE.
-               
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
-               &                  velhradcpfiles(1:jnumvelhradcp), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         IF ( ln_velfb ) THEN
-
-            DO jset = 1, jnumvelfb
-            
-               jveloset = jveloset + 1
-
-               ld_velav(jveloset) = ln_velfb_av(jset)
-               
-               CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
-                  &                  velfbfiles(jset:jset), &
-                  &                  nvelovars, nveloextr, &
-                  &                  nitend-nit000+2,              &
-                  &                  dobsini, dobsend, ln_ignmis, &
-                  &                  ld_velav(jveloset), &
-                  &                  .FALSE. )
-               
-               DO jvar = 1, 2
-                  CALL obs_prof_staend( velodata(jveloset), jvar )
-               END DO
-               
-               CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-                  &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-
-
-            END DO
-            
-         ENDIF
-
-      ENDIF
-     
    END SUBROUTINE dia_obs_init
 
@@ -974,13 +845,9 @@
       !!
       !! ** Method  : Call the observation operators on each time step to
-      !!              compute the model equivalent of the following date:
-      !!               - T profiles
-      !!               - S profiles
-      !!               - Sea surface height (referenced to a mean)
-      !!               - Sea surface temperature
-      !!               - Sea surface salinity
-      !!               - Velocity component (U,V) profiles
-      !!
-      !! ** Action  : 
+      !!              compute the model equivalent of the following data:
+      !!               - Profile data, currently T/S or U/V
+      !!               - Surface data, currently SST, SLA or sea-ice concentration.
+      !!
+      !! ** Action  :
       !!
       !! History :
@@ -991,43 +858,85 @@
       !!        !  07-04  (G. Smith) Generalized surface operators
       !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
+      !!        !  15-08  (M. Martin) Combined surface/profile routines.
       !!----------------------------------------------------------------------
       !! * Modules used
-      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
-         & rdt,           &                       
-         & gdept_1d,       &             
-         & tmask, umask, vmask                            
-      USE phycst, ONLY : &              ! Physical constants
-         & rday                         
-      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
-         & tsn,  &             
-         & un, vn,  &
+      USE phycst, ONLY : &         ! Physical constants
+#if defined key_fabm
+         & rt0,          &
+#endif
+         & rday
+      USE oce, ONLY : &            ! Ocean dynamics and tracers variables
+         & tsn,       &
+         & un,        &
+         & vn,        &
          & sshn
 #if defined  key_lim3
-      USE ice, ONLY : &                     ! LIM Ice model variables
+      USE ice, ONLY : &            ! LIM3 Ice model variables
          & frld
 #endif
 #if defined key_lim2
-      USE ice_2, ONLY : &                     ! LIM Ice model variables
+      USE ice_2, ONLY : &          ! LIM2 Ice model variables
          & frld
 #endif
+#if defined key_cice
+      USE sbc_oce, ONLY : fr_i     ! ice fraction
+#endif
+#if defined key_top
+      USE trc, ONLY :  &           ! Biogeochemical state variables
+         & trn
+#endif
+#if defined key_hadocc
+      USE par_hadocc               ! HadOCC parameters
+      USE trc, ONLY :  &
+         & HADOCC_CHL, &
+         & HADOCC_FCO2, &
+         & HADOCC_PCO2, &
+         & HADOCC_FILL_FLT
+      USE had_bgc_const, ONLY: c2n_p
+#elif defined key_medusa
+      USE par_medusa               ! MEDUSA parameters
+      USE sms_medusa, ONLY: &
+         & xthetapn, &
+         & xthetapd
+#if defined key_roam
+      USE sms_medusa, ONLY: &
+         & f2_pco2w, &
+         & f2_fco2w, &
+         & f3_pH
+#endif
+#elif defined key_fabm
+      USE par_fabm                 ! FABM parameters
+      USE fabm, ONLY: &
+         & fabm_get_interior_diagnostic_data
+#endif
+#if defined key_spm
+      USE par_spm, ONLY: &         ! Sediment parameters
+         & jp_spm
+#endif
+
       IMPLICIT NONE
 
       !! * Arguments
-      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
+      INTEGER, INTENT(IN) :: kstp  ! Current timestep
       !! * Local declarations
-      INTEGER :: idaystp                ! Number of timesteps per day
-      INTEGER :: jprofset               ! Profile data set loop variable
-      INTEGER :: jslaset                ! SLA data set loop variable
-      INTEGER :: jsstset                ! SST data set loop variable
-      INTEGER :: jseaiceset             ! sea ice data set loop variable
-      INTEGER :: jveloset               ! velocity profile data loop variable
-      INTEGER :: jvar                   ! Variable number    
-#if ! defined key_lim2 && ! defined key_lim3
-      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
-#endif
-      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
- 
-#if ! defined key_lim2 && ! defined key_lim3
-      CALL wrk_alloc(jpi,jpj,frld) 
+      INTEGER :: idaystp           ! Number of timesteps per day
+      INTEGER :: jtype             ! Data loop variable
+      INTEGER :: jvar              ! Variable number
+      INTEGER :: ji, jj, jk        ! Loop counters
+      REAL(wp) :: tiny             ! small number
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
+         & zprofvar                ! Model values for variables in a prof ob
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
+         & zprofmask               ! Mask associated with zprofvar
+      REAL(wp), POINTER, DIMENSION(:,:) :: &
+         & zsurfvar, &             ! Model values equivalent to surface ob.
+         & zsurfmask               ! Mask associated with surface variable
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zglam,    &             ! Model longitudes for prof variables
+         & zgphi                   ! Model latitudes for prof variables
+      LOGICAL :: llog10            ! Perform log10 transform of variable
+#if defined key_fabm
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & pco2_3d                 ! 3D pCO2 from FABM
 #endif
 
@@ -1036,4 +945,5 @@
          WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
          WRITE(numout,*) '~~~~~~~'
+         CALL FLUSH(numout)
       ENDIF
 
@@ -1041,82 +951,556 @@
 
       !-----------------------------------------------------------------------
-      ! No LIM => frld == 0.0_wp
-      !-----------------------------------------------------------------------
-#if ! defined key_lim2 && ! defined key_lim3
-      frld(:,:) = 0.0_wp
-#endif
-      !-----------------------------------------------------------------------
-      ! Depending on switches call various observation operators
-      !-----------------------------------------------------------------------
-
-      !  - Temperature/salinity profiles
-      IF ( ln_t3d .OR. ln_s3d ) THEN
-         DO jprofset = 1, nprofsets
-            IF ( ld_enact(jprofset) ) THEN
-               CALL obs_pro_opt( prodatqc(jprofset),                     &
-                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   &
-                  &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &
-                  &              gdept_1d, tmask, n1dint, n2dint,        &
-                  &              kdailyavtypes = endailyavtypes )
-            ELSE
-               CALL obs_pro_opt( prodatqc(jprofset),                     &
-                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   &
-                  &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &
-                  &              gdept_1d, tmask, n1dint, n2dint              )
+      ! Call the profile and surface observation operators
+      !-----------------------------------------------------------------------
+
+      IF ( nproftypes > 0 ) THEN
+
+         DO jtype = 1, nproftypes
+
+            ! Allocate local work arrays
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
+            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
+            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
+            
+            ! Defaults which might change
+            DO jvar = 1, profdataqc(jtype)%nvar
+               zprofmask(:,:,:,jvar) = tmask(:,:,:)
+               zglam(:,:,jvar)       = glamt(:,:)
+               zgphi(:,:,jvar)       = gphit(:,:)
+            END DO
+
+            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
+
+            CASE('prof')
+               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
+               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
+
+            CASE('vel')
+               zprofvar(:,:,:,1) = un(:,:,:)
+               zprofvar(:,:,:,2) = vn(:,:,:)
+               zprofmask(:,:,:,1) = umask(:,:,:)
+               zprofmask(:,:,:,2) = vmask(:,:,:)
+               zglam(:,:,1) = glamu(:,:)
+               zglam(:,:,2) = glamv(:,:)
+               zgphi(:,:,1) = gphiu(:,:)
+               zgphi(:,:,2) = gphiv(:,:)
+
+            CASE('plchltot')
+#if defined key_hadocc
+               ! Chlorophyll from HadOCC
+               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
+#elif defined key_medusa
+               ! Add non-diatom and diatom chlorophyll from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
+#elif defined key_fabm
+               ! Add all chlorophyll groups from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
+                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
+#else
+               CALL ctl_stop( ' Trying to run plchltot observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               ! Take the log10 where we can, otherwise exclude
+               tiny = 1.0e-20
+               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt )
+                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:))
+               ELSEWHERE
+                  zprofvar(:,:,:,:)  = obfillflt
+                  zprofmask(:,:,:,:) = 0
+               END WHERE
+               ! Mask out model below any excluded values,
+               ! to avoid interpolation issues
+               DO jvar = 1, profdataqc(jtype)%nvar
+                 DO jj = 1, jpj
+                    DO ji = 1, jpi
+                       depth_loop: DO jk = 1, jpk
+                          IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN
+                             zprofmask(ji,jj,jk:jpk,jvar) = 0
+                             EXIT depth_loop
+                          ENDIF
+                       END DO depth_loop
+                    END DO
+                 END DO
+              END DO
+
+            CASE('pchltot')
+#if defined key_hadocc
+               ! Chlorophyll from HadOCC
+               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
+#elif defined key_medusa
+               ! Add non-diatom and diatom chlorophyll from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
+#elif defined key_fabm
+               ! Add all chlorophyll groups from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
+                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
+#else
+               CALL ctl_stop( ' Trying to run pchltot observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('pno3')
+#if defined key_hadocc
+               ! Dissolved inorganic nitrogen from HadOCC
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut)
+#elif defined key_medusa
+               ! Dissolved inorganic nitrogen from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdin)
+#elif defined key_fabm
+               ! Nitrate from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n)
+#else
+               CALL ctl_stop( ' Trying to run pno3 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('psi4')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run psi4 observation operator', &
+                  &           ' but HadOCC does not simulate silicate' )
+#elif defined key_medusa
+               ! Silicate from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpsil)
+#elif defined key_fabm
+               ! Silicate from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s)
+#else
+               CALL ctl_stop( ' Trying to run psi4 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('ppo4')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
+                  &           ' but HadOCC does not simulate phosphate' )
+#elif defined key_medusa
+               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
+                  &           ' but MEDUSA does not simulate phosphate' )
+#elif defined key_fabm
+               ! Phosphate from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p)
+#else
+               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('pdic')
+#if defined key_hadocc
+               ! Dissolved inorganic carbon from HadOCC
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic)
+#elif defined key_medusa
+               ! Dissolved inorganic carbon from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdic)
+#elif defined key_fabm
+               ! Dissolved inorganic carbon from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
+#else
+               CALL ctl_stop( ' Trying to run pdic observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('palk')
+#if defined key_hadocc
+               ! Alkalinity from HadOCC
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk)
+#elif defined key_medusa
+               ! Alkalinity from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpalk)
+#elif defined key_fabm
+               ! Alkalinity from ERSEM
+               zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3a)
+#else
+               CALL ctl_stop( ' Trying to run palk observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('pph')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run pph observation operator', &
+                  &           ' but HadOCC has no pH diagnostic defined' )
+#elif defined key_medusa && defined key_roam
+               ! pH from MEDUSA
+               zprofvar(:,:,:,1) = f3_pH(:,:,:)
+#elif defined key_fabm
+               ! pH from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph)
+#else
+               CALL ctl_stop( ' Trying to run pph observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('po2')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run po2 observation operator', &
+                  &           ' but HadOCC does not simulate oxygen' )
+#elif defined key_medusa
+               ! Oxygen from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)
+#elif defined key_fabm
+               ! Oxygen from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o)
+#else
+               CALL ctl_stop( ' Trying to run po2 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE DEFAULT
+               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
+
+            END SELECT
+
+            DO jvar = 1, profdataqc(jtype)%nvar
+               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
+                  &               nit000, idaystp, jvar,                   &
+                  &               zprofvar(:,:,:,jvar),                    &
+                  &               fsdept(:,:,:), fsdepw(:,:,:),            & 
+                  &               zprofmask(:,:,:,jvar),                   &
+                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &
+                  &               nn_1dint, nn_2dint_default,              &
+                  &               kdailyavtypes = nn_profdavtypes )
+            END DO
+
+            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
+            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
+            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
+            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
+
+         END DO
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         !Allocate local work arrays
+         CALL wrk_alloc( jpi, jpj, zsurfvar )
+         CALL wrk_alloc( jpi, jpj, zsurfmask )
+#if defined key_fabm
+         CALL wrk_alloc( jpi, jpj, jpk, pco2_3d )
+#endif
+
+         DO jtype = 1, nsurftypes
+
+            !Defaults which might be changed
+            zsurfmask(:,:) = tmask(:,:,1)
+            llog10 = .FALSE.
+
+            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
+            CASE('sst')
+               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
+            CASE('sla')
+               zsurfvar(:,:) = sshn(:,:)
+            CASE('sss')
+               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
+            CASE('sic')
+               IF ( kstp == 0 ) THEN
+                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
+                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
+                        &           'time-step but some obs are valid then.' )
+                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
+                        &           ' sea-ice obs will be missed'
+                  ENDIF
+                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
+                     &                        surfdataqc(jtype)%nsstp(1)
+                  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
+
+            CASE('slchltot')
+#if defined key_hadocc
+               ! Surface chlorophyll from HadOCC
+               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
+#elif defined key_medusa
+               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
+               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
+#elif defined key_fabm
+               ! Add all surface chlorophyll groups from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
+                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
+#else
+               CALL ctl_stop( ' Trying to run slchltot observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slchldia')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slchldia observation operator', &
+                  &           ' but HadOCC does not explicitly simulate diatoms' )
+#elif defined key_medusa
+               ! Diatom surface chlorophyll from MEDUSA
+               zsurfvar(:,:) = trn(:,:,1,jpchd)
+#elif defined key_fabm
+               ! Diatom surface chlorophyll from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1)
+#else
+               CALL ctl_stop( ' Trying to run slchldia observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slchlnon')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
+                  &           ' but HadOCC does not explicitly simulate non-diatoms' )
+#elif defined key_medusa
+               ! Non-diatom surface chlorophyll from MEDUSA
+               zsurfvar(:,:) = trn(:,:,1,jpchn)
+#elif defined key_fabm
+               ! Add all non-diatom surface chlorophyll groups from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
+                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
+#else
+               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slchldin')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slchldin observation operator', &
+                  &           ' but HadOCC does not explicitly simulate dinoflagellates' )
+#elif defined key_medusa
+               CALL ctl_stop( ' Trying to run slchldin observation operator', &
+                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' )
+#elif defined key_fabm
+               ! Dinoflagellate surface chlorophyll from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
+#else
+               CALL ctl_stop( ' Trying to run slchldin observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slchlmic')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
+                  &           ' but HadOCC does not explicitly simulate microphytoplankton' )
+#elif defined key_medusa
+               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
+                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' )
+#elif defined key_fabm
+               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
+#else
+               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slchlnan')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
+                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' )
+#elif defined key_medusa
+               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
+                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' )
+#elif defined key_fabm
+               ! Nanophytoplankton surface chlorophyll from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2)
+#else
+               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slchlpic')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
+                  &           ' but HadOCC does not explicitly simulate picophytoplankton' )
+#elif defined key_medusa
+               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
+                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' )
+#elif defined key_fabm
+               ! Picophytoplankton surface chlorophyll from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3)
+#else
+               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('schltot')
+#if defined key_hadocc
+               ! Surface chlorophyll from HadOCC
+               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
+#elif defined key_medusa
+               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
+               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
+#elif defined key_fabm
+               ! Add all surface chlorophyll groups from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
+                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
+#else
+               CALL ctl_stop( ' Trying to run schltot observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('slphytot')
+#if defined key_hadocc
+               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio
+               zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
+#elif defined key_medusa
+               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA
+               ! multiplied by C:N ratio for each
+               zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
+#elif defined key_fabm
+               ! Add all surface phytoplankton carbon groups from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
+                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
+#else
+               CALL ctl_stop( ' Trying to run slphytot observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slphydia')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slphydia observation operator', &
+                  &           ' but HadOCC does not explicitly simulate diatoms' )
+#elif defined key_medusa
+               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
+               zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd
+#elif defined key_fabm
+               ! Diatom surface phytoplankton carbon from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c)
+#else
+               CALL ctl_stop( ' Trying to run slphydia observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('slphynon')
+#if defined key_hadocc
+               CALL ctl_stop( ' Trying to run slphynon observation operator', &
+                  &           ' but HadOCC does not explicitly simulate non-diatoms' )
+#elif defined key_medusa
+               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
+               zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn
+#elif defined key_fabm
+               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
+                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
+#else
+               CALL ctl_stop( ' Trying to run slphynon observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               llog10 = .TRUE.
+
+            CASE('sspm')
+#if defined key_spm
+               zsurfvar(:,:) = 0.0
+               DO jn = 1, jp_spm
+                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
+               END DO
+#else
+               CALL ctl_stop( ' Trying to run sspm observation operator', &
+                  &           ' but no spm model appears to have been defined' )
+#endif
+
+            CASE('sfco2')
+#if defined key_hadocc
+               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
+               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
+                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
+                  zsurfvar(:,:) = obfillflt
+                  zsurfmask(:,:) = 0
+                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
+                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
+               ENDIF
+#elif defined key_medusa && defined key_roam
+               zsurfvar(:,:) = f2_fco2w(:,:)
+#elif defined key_fabm
+               ! First, get pCO2 from FABM
+               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
+               zsurfvar(:,:) = pco2_3d(:,:,1)
+               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
+               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
+               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
+               ! and
+               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
+               ! Marine Chemistry, 2: 203-215.
+               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
+               ! not explicitly included - atmospheric pressure is not necessarily available so this is
+               ! the best assumption.
+               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
+               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
+               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
+               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
+               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + &
+                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
+                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
+                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
+                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
+                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
+#else
+               CALL ctl_stop( ' Trying to run sfco2 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE('spco2')
+#if defined key_hadocc
+               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
+               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
+                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
+                  zsurfvar(:,:) = obfillflt
+                  zsurfmask(:,:) = 0
+                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
+                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
+               ENDIF
+#elif defined key_medusa && defined key_roam
+               zsurfvar(:,:) = f2_pco2w(:,:)
+#elif defined key_fabm
+               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
+               zsurfvar(:,:) = pco2_3d(:,:,1)
+#else
+               CALL ctl_stop( ' Trying to run spco2 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+            CASE DEFAULT
+
+               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
+
+            END SELECT
+            
+            IF ( llog10 ) THEN
+               ! Take the log10 where we can, otherwise exclude
+               tiny = 1.0e-20
+               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
+                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:))
+               ELSEWHERE
+                  zsurfvar(:,:)  = obfillflt
+                  zsurfmask(:,:) = 0
+               END WHERE
             ENDIF
+
+            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
+               &               nit000, idaystp, zsurfvar, zsurfmask,    &
+               &               n2dintsurf(jtype), llnightav(jtype),     &
+               &               ravglamscl(jtype), ravgphiscl(jtype),     &
+               &               lfpindegs(jtype) )
+
          END DO
+
+         CALL wrk_dealloc( jpi, jpj, zsurfvar )
+         CALL wrk_dealloc( jpi, jpj, zsurfmask )
+#if defined key_fabm
+         CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d )
+#endif
+
       ENDIF
 
-      !  - Sea surface anomaly
-      IF ( ln_sla ) THEN
-         DO jslaset = 1, nslasets
-            CALL obs_sla_opt( sladatqc(jslaset),            &
-               &              kstp, jpi, jpj, nit000, sshn, &
-               &              tmask(:,:,1), n2dint )
-         END DO         
-      ENDIF
-
-      !  - Sea surface temperature
-      IF ( ln_sst ) THEN
-         DO jsstset = 1, nsstsets
-            CALL obs_sst_opt( sstdatqc(jsstset),                &
-               &              kstp, jpi, jpj, nit000, idaystp,  &
-               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  &
-               &              n2dint, ld_sstnight(jsstset) )
-         END DO
-      ENDIF
-
-      !  - Sea surface salinity
-      IF ( ln_sss ) THEN
-         IF(lwp) WRITE(numout,*) ' SSS currently not available'
-      ENDIF
-
-#if defined key_lim2 || defined key_lim3
-      IF ( ln_seaice ) THEN
-         DO jseaiceset = 1, nseaicesets
-            CALL obs_seaice_opt( seaicedatqc(jseaiceset),      &
-               &              kstp, jpi, jpj, nit000, 1.-frld, &
-               &              tmask(:,:,1), n2dint )
-         END DO
-      ENDIF      
-#endif
-
-      !  - Velocity profiles
-      IF ( ln_vel3d ) THEN
-         DO jveloset = 1, nvelosets
-           ! zonal component of velocity
-           CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
-              &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
-                             n1dint, n2dint, ld_velav(jveloset) )
-         END DO
-      ENDIF
-
-#if ! defined key_lim2 && ! defined key_lim3
-      CALL wrk_dealloc(jpi,jpj,frld) 
-#endif
-
    END SUBROUTINE dia_obs
-  
-   SUBROUTINE dia_obs_wri 
+
+   SUBROUTINE dia_obs_wri
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_obs_wri  ***
@@ -1126,5 +1510,5 @@
       !! ** Method  : Call observation diagnostic output routines
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!
       !! History :
@@ -1134,258 +1518,76 @@
       !!        !  07-03  (K. Mogensen) General handling of profiles
       !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
+      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
       !!----------------------------------------------------------------------
+      !! * Modules used
+      USE obs_rot_vel          ! Rotation of velocities
+
       IMPLICIT NONE
 
       !! * Local declarations
-
-      INTEGER :: jprofset                 ! Profile data set loop variable
-      INTEGER :: jveloset                 ! Velocity data set loop variable
-      INTEGER :: jslaset                  ! SLA data set loop variable
-      INTEGER :: jsstset                  ! SST data set loop variable
-      INTEGER :: jseaiceset               ! Sea Ice data set loop variable
-      INTEGER :: jset
-      INTEGER :: jfbini
-      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
-      CHARACTER(LEN=10) :: cdtmp
+      INTEGER :: jtype                    ! Data set loop variable
+      INTEGER :: jo, jvar, jk
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
+         & zu, &
+         & zv
+
       !-----------------------------------------------------------------------
       ! Depending on switches call various observation output routines
       !-----------------------------------------------------------------------
 
-      !  - Temperature/salinity profiles
-
-      IF( ln_t3d .OR. ln_s3d ) THEN
-
-         ! Copy data from prodatqc to profdata structures
-         DO jprofset = 1, nprofsets
-
-            CALL obs_prof_decompress( prodatqc(jprofset), &
-                 &                    profdata(jprofset), .TRUE., numout )
+      IF ( nproftypes > 0 ) THEN
+
+         DO jtype = 1, nproftypes
+
+            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
+
+               ! For velocity data, rotate the model velocities to N/S, E/W
+               ! using the compressed data structure.
+               ALLOCATE( &
+                  & zu(profdataqc(jtype)%nvprot(1)), &
+                  & zv(profdataqc(jtype)%nvprot(2))  &
+                  & )
+
+               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
+
+               DO jo = 1, profdataqc(jtype)%nprof
+                  DO jvar = 1, 2
+                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
+
+                        IF ( jvar == 1 ) THEN
+                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
+                        ELSE
+                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
+                        ENDIF
+
+                     END DO
+                  END DO
+               END DO
+
+               DEALLOCATE( zu )
+               DEALLOCATE( zv )
+
+            END IF
+
+            CALL obs_prof_decompress( profdataqc(jtype), &
+               &                      profdata(jtype), .TRUE., numout )
+
+            CALL obs_wri_prof( profdata(jtype) )
 
          END DO
 
-         ! Write the profiles.
-
-         jprofset = 0
-
-         ! ENACT insitu data
-
-         IF ( ln_ena ) THEN
-           
-            jprofset = jprofset + 1
-
-            CALL obs_wri_p3d( 'enact', profdata(jprofset) )
-
-         ENDIF
-
-         ! Coriolis insitu data
-
-         IF ( ln_cor ) THEN
-            
-            jprofset = jprofset + 1
-
-            CALL obs_wri_p3d( 'corio', profdata(jprofset) )
-            
-         ENDIF
-         
-         ! Feedback insitu data
-
-         IF ( ln_profb ) THEN
-
-            jfbini = jprofset + 1
-
-            DO jprofset = jfbini, nprofsets
-               
-               jset = jprofset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'profb_',jset
-               CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
-
-            END DO
-
-         ENDIF
-
       ENDIF
 
-      !  - Sea surface anomaly
-      IF ( ln_sla ) THEN
-
-         ! Copy data from sladatqc to sladata structures
-         DO jslaset = 1, nslasets
-
-              CALL obs_surf_decompress( sladatqc(jslaset), &
-                 &                    sladata(jslaset), .TRUE., numout )
+      IF ( nsurftypes > 0 ) THEN
+
+         DO jtype = 1, nsurftypes
+
+            CALL obs_surf_decompress( surfdataqc(jtype), &
+               &                      surfdata(jtype), .TRUE., numout )
+
+            CALL obs_wri_surf( surfdata(jtype) )
 
          END DO
 
-         jslaset = 0 
-
-         ! Write the AVISO SLA data
-
-         IF ( ln_sladt ) THEN
-            
-            jslaset = 1
-            CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
-            jslaset = 2
-            CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
-
-         ENDIF
-
-         IF ( ln_slafb ) THEN
-            
-            jfbini = jslaset + 1
-
-            DO jslaset = jfbini, nslasets
-               
-               jset = jslaset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
-               CALL obs_wri_sla( cdtmp, sladata(jslaset) )
-
-            END DO
-
-         ENDIF
-
-      ENDIF
-
-      !  - Sea surface temperature
-      IF ( ln_sst ) THEN
-
-         ! Copy data from sstdatqc to sstdata structures
-         DO jsstset = 1, nsstsets
-     
-              CALL obs_surf_decompress( sstdatqc(jsstset), &
-                 &                    sstdata(jsstset), .TRUE., numout )
-
-         END DO
-
-         jsstset = 0 
-
-         ! Write the AVISO SST data
-
-         IF ( ln_reysst ) THEN
-            
-            jsstset = jsstset + 1
-            CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
-
-         ENDIF
-
-         IF ( ln_ghrsst ) THEN
-            
-            jsstset = jsstset + 1
-            CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
-
-         ENDIF
-
-         IF ( ln_sstfb ) THEN
-            
-            jfbini = jsstset + 1
-
-            DO jsstset = jfbini, nsstsets
-               
-               jset = jsstset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
-               CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
-
-            END DO
-
-         ENDIF
-
-      ENDIF
-
-      !  - Sea surface salinity
-      IF ( ln_sss ) THEN
-         IF(lwp) WRITE(numout,*) ' SSS currently not available'
-      ENDIF
-
-      !  - Sea Ice Concentration
-      IF ( ln_seaice ) THEN
-
-         ! Copy data from seaicedatqc to seaicedata structures
-         DO jseaiceset = 1, nseaicesets
-
-              CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
-                 &                    seaicedata(jseaiceset), .TRUE., numout )
-
-         END DO
-
-         ! Write the Sea Ice data
-         DO jseaiceset = 1, nseaicesets
-      
-            WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
-            CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
-
-         END DO
-
-      ENDIF
-      
-      ! Velocity data
-      IF( ln_vel3d ) THEN
-
-         ! Copy data from veldatqc to velodata structures
-         DO jveloset = 1, nvelosets
-
-            CALL obs_prof_decompress( veldatqc(jveloset), &
-                 &                    velodata(jveloset), .TRUE., numout )
-
-         END DO
-
-         ! Write the profiles.
-
-         jveloset = 0
-
-         ! Daily averaged data
-
-         IF ( ln_velavcur ) THEN
-            
-            jveloset = jveloset + 1
-
-            CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
-
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhrcur ) THEN
-            
-            jveloset = jveloset + 1
-
-            CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
-
-         ENDIF
-
-         ! Daily averaged data
-
-         IF ( ln_velavadcp ) THEN
-            
-            jveloset = jveloset + 1
-
-            CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
-
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhradcp ) THEN
-            
-            jveloset = jveloset + 1
-            
-            CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
-               
-         ENDIF
-
-         ! Feedback velocity data
-
-         IF ( ln_velfb ) THEN
-
-            jfbini = jveloset + 1
-
-            DO jveloset = jfbini, nvelosets
-               
-               jset = jveloset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
-               CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
-
-            END DO
-
-         ENDIF
-         
       ENDIF
 
@@ -1405,30 +1607,15 @@
       !!
       !!----------------------------------------------------------------------
-      !! obs_grid deallocation
+      ! obs_grid deallocation
       CALL obs_grid_deallocate
 
-      !! diaobs deallocation
-      IF ( nprofsets > 0 ) THEN
-          DEALLOCATE(ld_enact, &
-                  &  profdata, &
-                  &  prodatqc)
-      END IF
-      IF ( ln_sla ) THEN
-          DEALLOCATE(sladata, &
-                  &  sladatqc)
-      END IF
-      IF ( ln_seaice ) THEN
-          DEALLOCATE(sladata, &
-                  &  sladatqc)
-      END IF
-      IF ( ln_sst ) THEN
-          DEALLOCATE(sstdata, &
-                  &  sstdatqc)
-      END IF
-      IF ( ln_vel3d ) THEN
-          DEALLOCATE(ld_velav, &
-                  &  velodata, &
-                  &  veldatqc)
-      END IF
+      ! diaobs deallocation
+      IF ( nproftypes > 0 ) &
+         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
+
+      IF ( nsurftypes > 0 ) &
+         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
+         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
+
    END SUBROUTINE dia_obs_dealloc
 
@@ -1436,10 +1623,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ini_date  ***
-      !!          
-      !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
       !!
       !! History :
@@ -1452,6 +1639,4 @@
       USE phycst, ONLY : &            ! Physical constants
          & rday
-!      USE daymod, ONLY : &            ! Time variables
-!         & nmonth_len           
       USE dom_oce, ONLY : &           ! Ocean space and time domain variables
          & rdt
@@ -1460,5 +1645,5 @@
 
       !! * Arguments
-      REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS
 
       !! * Local declarations
@@ -1468,13 +1653,13 @@
       INTEGER :: ihou
       INTEGER :: imin
-      INTEGER :: imday         ! Number of days in month.
-      REAL(KIND=wp) :: zdayfrc ! Fraction of day
-
-      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
-
-      !!----------------------------------------------------------------------
-      !! Initial date initialization (year, month, day, hour, minute)
-      !! (This assumes that the initial date is for 00z))
-      !!----------------------------------------------------------------------
+      INTEGER :: imday       ! Number of days in month.
+      INTEGER, DIMENSION(12) :: &
+         &       imonth_len  ! Length in days of the months of the current year
+      REAL(wp) :: zdayfrc    ! Fraction of day
+
+      !----------------------------------------------------------------------
+      ! Initial date initialization (year, month, day, hour, minute)
+      ! (This assumes that the initial date is for 00z))
+      !----------------------------------------------------------------------
       iyea =   ndate0 / 10000
       imon = ( ndate0 - iyea * 10000 ) / 100
@@ -1483,7 +1668,7 @@
       imin = 0
 
-      !!----------------------------------------------------------------------
-      !! Compute number of days + number of hours + min since initial time
-      !!----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      ! Compute number of days + number of hours + min since initial time
+      !----------------------------------------------------------------------
       iday = iday + ( nit000 -1 ) * rdt / rday
       zdayfrc = ( nit000 -1 ) * rdt / rday
@@ -1492,10 +1677,10 @@
       imin = int( (zdayfrc * 24 - ihou) * 60 )
 
-      !!-----------------------------------------------------------------------
-      !! Convert number of days (iday) into a real date
-      !!----------------------------------------------------------------------
+      !-----------------------------------------------------------------------
+      ! Convert number of days (iday) into a real date
+      !----------------------------------------------------------------------
 
       CALL calc_month_len( iyea, imonth_len )
-      
+
       DO WHILE ( iday > imonth_len(imon) )
          iday = iday - imonth_len(imon)
@@ -1508,7 +1693,7 @@
       END DO
 
-      !!----------------------------------------------------------------------
-      !! Convert it into YYYYMMDD.HHMMSS format.
-      !!----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      ! Convert it into YYYYMMDD.HHMMSS format.
+      !----------------------------------------------------------------------
       ddobsini = iyea * 10000_dp + imon * 100_dp + &
          &       iday + ihou * 0.01_dp + imin * 0.0001_dp
@@ -1520,10 +1705,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE fin_date  ***
-      !!          
-      !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format
       !!
       !! History :
@@ -1535,6 +1720,4 @@
       USE phycst, ONLY : &            ! Physical constants
          & rday
-!      USE daymod, ONLY : &            ! Time variables
-!         & nmonth_len                
       USE dom_oce, ONLY : &           ! Ocean space and time domain variables
          & rdt
@@ -1543,5 +1726,5 @@
 
       !! * Arguments
-      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
 
       !! * Local declarations
@@ -1551,9 +1734,9 @@
       INTEGER :: ihou
       INTEGER :: imin
-      INTEGER :: imday         ! Number of days in month.
-      REAL(KIND=wp) :: zdayfrc       ! Fraction of day
-         
-      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
-            
+      INTEGER :: imday       ! Number of days in month.
+      INTEGER, DIMENSION(12) :: &
+         &       imonth_len  ! Length in days of the months of the current year
+      REAL(wp) :: zdayfrc    ! Fraction of day
+
       !-----------------------------------------------------------------------
       ! Initial date initialization (year, month, day, hour, minute)
@@ -1565,5 +1748,5 @@
       ihou = 0
       imin = 0
-      
+
       !-----------------------------------------------------------------------
       ! Compute number of days + number of hours + min since initial time
@@ -1580,5 +1763,5 @@
 
       CALL calc_month_len( iyea, imonth_len )
-      
+
       DO WHILE ( iday > imonth_len(imon) )
          iday = iday - imonth_len(imon)
@@ -1598,4 +1781,119 @@
 
     END SUBROUTINE fin_date
-    
+
+    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles )
+
+       INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
+       INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
+       INTEGER, DIMENSION(ntypes), INTENT(OUT) :: &
+          &                   ifiles      ! Out number of files for each type
+       CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: &
+          &                   cobstypes   ! List of obs types
+       CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: &
+          &                   cfiles      ! List of files for all types
+
+       !Local variables
+       INTEGER :: jfile
+       INTEGER :: jtype
+
+       DO jtype = 1, ntypes
+
+          ifiles(jtype) = 0
+          DO jfile = 1, jpmaxnfiles
+             IF ( trim(cfiles(jtype,jfile)) /= '' ) &
+                       ifiles(jtype) = ifiles(jtype) + 1
+          END DO
+
+          IF ( ifiles(jtype) == 0 ) THEN
+               CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//   &
+                  &           ' set to true but no files available to read' )
+          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 DO
+
+    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=8), 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 >= 0) .AND. (n2dint_type <= 6) ) THEN
+          n2dint(jtype) = n2dint_type
+       ELSE IF ( n2dint_type == -1 ) THEN
+          n2dint(jtype) = n2dint_default
+       ELSE
+          CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
+            &                    ' is not available')
+       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
