Index: /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9305)
+++ /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9306)
@@ -1201,8 +1201,27 @@
    ln_vel3d   = .false.             ! Logical switch for velocity observations
    ln_sss     = .false.             ! Logical swithc for SSS observations
-   ln_logchl  = .false.             ! Logical switch for log(Chl) obs
-   ln_spm     = .false.             ! Logical switch for SPM obs
-   ln_fco2    = .false.             
-   ln_pco2    = .false.
+   ln_slchltot = .false.            ! Logical switch for surface total              log10(chlorophyll) obs
+   ln_slchldia = .false.            ! Logical switch for surface diatom             log10(chlorophyll) obs
+   ln_slchlnon = .false.            ! Logical switch for surface non-diatom         log10(chlorophyll) obs
+   ln_slchldin = .false.            ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs
+   ln_slchlmic = .false.            ! Logical switch for surface microphytoplankton log10(chlorophyll) obs
+   ln_slchlnan = .false.            ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs
+   ln_slchlpic = .false.            ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs
+   ln_schltot  = .false.            ! Logical switch for surface total              chlorophyll        obs
+   ln_slphytot = .false.            ! Logical switch for surface total      log10(phytoplankton carbon) obs
+   ln_slphydia = .false.            ! Logical switch for surface diatom     log10(phytoplankton carbon) obs
+   ln_slphynon = .false.            ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs
+   ln_sspm     = .false.            ! Logical switch for surface suspended particulate matter obs
+   ln_sfco2    = .false.            ! Logical switch for surface fugacity         of carbon dioxide obs
+   ln_spco2    = .false.            ! Logical switch for surface partial pressure of carbon dioxide obs
+   ln_plchltot = .false.            ! Logical switch for profile total log10(chlorophyll) obs
+   ln_pchltot  = .false.            ! Logical switch for profile total chlorophyll obs
+   ln_pno3     = .false.            ! Logical switch for profile nitrate obs
+   ln_psi4     = .false.            ! Logical switch for profile silicate obs
+   ln_ppo4     = .false.            ! Logical switch for profile phosphate obs
+   ln_pdic     = .false.            ! Logical switch for profile dissolved inorganic carbon obs
+   ln_palk     = .false.            ! Logical switch for profile alkalinity obs
+   ln_pph      = .false.            ! Logical switch for profile pH obs
+   ln_po2      = .false.            ! Logical switch for profile dissolved oxygen obs
    ln_altbias = .false.             ! Logical switch for altimeter bias correction
    ln_sstbias = .false.             ! Logical switch for SST bias correction
@@ -1213,37 +1232,63 @@
    ln_s_at_t  = .false.             ! Logical switch for computing model S at T obs if not there
    ln_sstnight = .false.            ! Logical switch for calculating night-time average for SST obs
-   ln_sla_fp_indegs = .true.        ! Logical: T=> averaging footpring is in degrees, F=> in metres
+   ln_default_fp_indegs = .true.    ! Logical: T=> averaging footprint is in degrees, F=> in metres
+   ln_sla_fp_indegs = .true.        ! Logical: T=> averaging footprint is in degrees, F=> in metres
    ln_sst_fp_indegs = .true.
    ln_sss_fp_indegs = .true.
    ln_sic_fp_indegs = .true.
 ! All of the *files* variables below are arrays. Use namelist_cfg to add more files
-   cn_profbfiles = 'profiles_01.nc'    ! Profile feedback input observation file names
-   cn_slafbfiles = 'sla_01.nc'         ! SLA feedback input observation file names
-   cn_sstfbfiles = 'sst_01.nc'         ! SST feedback input observation file names
-   cn_sicfbfiles = 'sic_01.nc'         ! SIC feedback input observation file names
-   cn_velfbfiles = 'vel_01.nc'         ! Velocity feedback input observation file names
-   cn_sssfbfiles = 'sss_01.nc'         ! SSS feedback input observation file names
-   cn_altbiasfile = 'altbias.nc'       ! Altimeter bias input file name
-   cn_sstbiasfiles = 'sstbias.nc'      ! SST bias input file names
-   cn_gridsearchfile = 'gridsearch.nc' ! Grid search file name
-   rn_gridsearchres = 0.5              ! Grid search resolution
-   rn_sla_avglamscl = 0.               ! E/W diameter of SLA observation footprint (metres/degrees)
-   rn_sla_avgphiscl = 0.               ! N/S diameter of SLA observation footprint (metres/degrees)
-   rn_sst_avglamscl = 0.
-   rn_sst_avgphiscl = 0.
-   rn_sss_avglamscl = 0.
-   rn_sss_avgphiscl = 0.
-   rn_sic_avglamscl = 0.
-   rn_sic_avgphiscl = 0.
-   nn_1dint = 0                        ! Type of vertical interpolation method
-   nn_2dint = 0                        ! Default horizontal interpolation method
-   nn_2dint_sla = 0                    ! Horizontal interpolation method for SLA
-   nn_2dint_sst = 0                    ! Horizontal interpolation method for SST
-   nn_2dint_sss = 0                    ! Horizontal interpolation method for SSS
-   nn_2dint_sic = 0                    ! Horizontal interpolation method for SIC
-   nn_msshc = 0                        ! MSSH correction scheme
-   rn_mdtcorr = 1.61                   ! MDT  correction
-   rn_mdtcutoff = 65.0                 ! MDT cutoff for computed correction
-   nn_profdavtypes = -1                ! Profile daily average types - array
+   cn_profbfiles = 'profiles_01.nc'      ! Profile feedback input observation file names
+   cn_slafbfiles = 'sla_01.nc'           ! SLA feedback input observation file names
+   cn_sstfbfiles = 'sst_01.nc'           ! SST feedback input observation file names
+   cn_sicfbfiles = 'sic_01.nc'           ! SIC feedback input observation file names
+   cn_velfbfiles = 'vel_01.nc'           ! Velocity feedback input observation file names
+   cn_sssfbfiles = 'sss_01.nc'           ! SSS feedback input observation file names
+   cn_slchltotfbfiles = 'slchltot_01.nc' ! Surface total              log10(chlorophyll) obs file names
+   cn_slchldiafbfiles = 'slchldia_01.nc' ! Surface diatom             log10(chlorophyll) obs file names
+   cn_slchlnonfbfiles = 'slchlnon_01.nc' ! Surface non-diatom         log10(chlorophyll) obs file names
+   cn_slchldinfbfiles = 'slchldin_01.nc' ! Surface dinoflagellate     log10(chlorophyll) obs file names
+   cn_slchlmicfbfiles = 'slchlmic_01.nc' ! Surface microphytoplankton log10(chlorophyll) obs file names
+   cn_slchlnanfbfiles = 'slchlnan_01.nc' ! Surface nanophytoplankton  log10(chlorophyll) obs file names
+   cn_slchlpicfbfiles = 'slchlpic_01.nc' ! Surface picophytoplankton  log10(chlorophyll) obs file names
+   cn_schltotfbfiles  = 'schltot_01.nc'  ! Surface total              chlorophyll        obs file names
+   cn_slphytotfbfiles = 'slphytot_01.nc' ! Surface total      log10(phytoplankton carbon) obs file names
+   cn_slphydiafbfiles = 'slphydia_01.nc' ! Surface diatom     log10(phytoplankton carbon) obs file names
+   cn_slphynonfbfiles = 'slphynon_01.nc' ! Surface non-diatom log10(phytoplankton carbon) obs file names
+   cn_sspmfbfiles     = 'sspm_01.nc'     ! Surface suspended particulate matter obs file names
+   cn_sfco2fbfiles    = 'sfco2_01.nc'    ! Surface fugacity         of carbon dioxide obs file names
+   cn_spco2fbfiles    = 'spco2_01.nc'    ! Surface partial pressure of carbon dioxide obs file names
+   cn_plchltotfbfiles = 'plchltot_01.nc' ! Profile total log10(chlorophyll) obs file names
+   cn_pchltotfbfiles  = 'pchltot_01.nc'  ! Profile total chlorophyll obs file names
+   cn_pno3fbfiles     = 'pno3_01.nc'     ! Profile nitrate obs file names
+   cn_psi4fbfiles     = 'psi4_01.nc'     ! Profile silicate obs file names
+   cn_ppo4fbfiles     = 'ppo4_01.nc'     ! Profile phosphate obs file names
+   cn_pdicfbfiles     = 'pdic_01.nc'     ! Profile dissolved inorganic carbon obs file names
+   cn_palkfbfiles     = 'palk_01.nc'     ! Profile alkalinity obs file names
+   cn_pphfbfiles      = 'pph_01.nc'      ! Profile pH obs file names
+   cn_po2fbfiles      = 'po2_01.nc'      ! Profile dissolved oxygen obs file names
+   cn_altbiasfile = 'altbias.nc'         ! Altimeter bias input file name
+   cn_sstbiasfiles = 'sstbias.nc'        ! SST bias input file names
+   cn_gridsearchfile = 'gridsearch.nc'   ! Grid search file name
+   rn_gridsearchres = 0.5                ! Grid search resolution
+   rn_default_avglamscl = 0.             ! Default E/W diameter of observation footprint (metres/degrees)
+   rn_default_avgphiscl = 0.             ! Default N/S diameter of observation footprint (metres/degrees)
+   rn_sla_avglamscl = 0.                 ! E/W diameter of SLA observation footprint (metres/degrees)
+   rn_sla_avgphiscl = 0.                 ! N/S diameter of SLA observation footprint (metres/degrees)
+   rn_sst_avglamscl = 0.                 ! E/W diameter of SST observation footprint (metres/degrees)
+   rn_sst_avgphiscl = 0.                 ! N/S diameter of SST observation footprint (metres/degrees)
+   rn_sss_avglamscl = 0.                 ! E/W diameter of SSS observation footprint (metres/degrees)
+   rn_sss_avgphiscl = 0.                 ! N/S diameter of SSS observation footprint (metres/degrees)
+   rn_sic_avglamscl = 0.                 ! E/W diameter of SIC observation footprint (metres/degrees)
+   rn_sic_avgphiscl = 0.                 ! N/S diameter of SIC observation footprint (metres/degrees)
+   nn_1dint = 0                          ! Type of vertical interpolation method
+   nn_2dint_default = 0                  ! Default horizontal interpolation method
+   nn_2dint_sla = -1                     ! Horizontal interpolation method for SLA
+   nn_2dint_sst = -1                     ! Horizontal interpolation method for SST
+   nn_2dint_sss = -1                     ! Horizontal interpolation method for SSS
+   nn_2dint_sic = -1                     ! Horizontal interpolation method for SIC
+   nn_msshc = 0                          ! MSSH correction scheme
+   rn_mdtcorr = 1.61                     ! MDT  correction
+   rn_mdtcutoff = 65.0                   ! MDT cutoff for computed correction
+   nn_profdavtypes = -1                  ! Profile daily average types - array
    
 /
Index: /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 9305)
+++ /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 9306)
@@ -45,27 +45,30 @@
    !! * Module variables
    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_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres
-   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres
-   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres
-   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
-
-   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres)
-   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres)
-   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres)
-   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres)
-   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres)
-   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres)
-   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres)
-   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres)
-
-   INTEGER :: nn_1dint       !: Vertical interpolation method
-   INTEGER :: nn_2dint       !: Default horizontal interpolation method
-   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method 
-   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method 
-   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method 
-   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method 
+      &       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) :: &
@@ -95,5 +98,5 @@
       & profdataqc           !: Profile data after quality control
 
-   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
+   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
       & cobstypesprof, &     !: Profile obs types
       & cobstypessurf        !: Surface obs types
@@ -141,16 +144,36 @@
       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_logchlfbfiles, &   ! Log(Chl) input filenames
-         & cn_spmfbfiles,    &   ! Sediment input filenames
-         & cn_fco2fbfiles,   &   ! fco2 input filenames
-         & cn_pco2fbfiles,   &   ! pco2 input filenames
+         & 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
 
@@ -166,8 +189,27 @@
       LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs
       LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
-      LOGICAL :: ln_logchl       ! Logical switch for log(Chl) obs
-      LOGICAL :: ln_spm          ! Logical switch for sediment obs
-      LOGICAL :: ln_fco2         ! Logical switch for fco2 obs
-      LOGICAL :: ln_pco2         ! Logical switch for pco2 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
@@ -180,29 +222,37 @@
       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 :: llvar1          ! Logical for profile variable 1
-      LOGICAL :: llvar2          ! Logical for profile variable 1
-
-      REAL(wp), POINTER, DIMENSION(:,:) :: &
-         & zglam1, &             ! Model longitudes for profile variable 1
-         & zglam2                ! Model longitudes for profile variable 2
-      REAL(wp), POINTER, DIMENSION(:,:) :: &
-         & zgphi1, &             ! Model latitudes for profile variable 1
-         & zgphi2                ! Model latitudes for profile variable 2
+      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(:,:,:) :: &
-         & zmask1, &             ! Model land/sea mask associated with variable 1
-         & zmask2                ! Model land/sea mask associated with variable 2
+         & 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_logchl, ln_spm, ln_fco2, ln_pco2,            &
+         &            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_sstnight, ln_default_fp_indegs,              &
          &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
          &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
@@ -210,14 +260,24 @@
          &            cn_sstfbfiles, cn_sicfbfiles,                   &
          &            cn_velfbfiles, cn_sssfbfiles,                   &
-         &            cn_logchlfbfiles, cn_spmfbfiles,                &
-         &            cn_fco2fbfiles, cn_pco2fbfiles,                 &
+         &            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,                             &
+         &            nn_1dint, nn_2dint_default,                     &
          &            nn_2dint_sla, nn_2dint_sst,                     &
          &            nn_2dint_sss, nn_2dint_sic,                     &
@@ -225,11 +285,4 @@
          &            nn_profdavtypes
 
-      CALL wrk_alloc( jpi, jpj, zglam1 )
-      CALL wrk_alloc( jpi, jpj, zglam2 )
-      CALL wrk_alloc( jpi, jpj, zgphi1 )
-      CALL wrk_alloc( jpi, jpj, zgphi2 )
-      CALL wrk_alloc( jpi, jpj, jpk, zmask1 )
-      CALL wrk_alloc( jpi, jpj, jpk, zmask2 )
-
       !-----------------------------------------------------------------------
       ! Read namelist parameters
@@ -237,16 +290,35 @@
 
       ! Some namelist arrays need initialising
-      cn_profbfiles(:)    = ''
-      cn_slafbfiles(:)    = ''
-      cn_sstfbfiles(:)    = ''
-      cn_sicfbfiles(:)    = ''
-      cn_velfbfiles(:)    = ''
-      cn_sssfbfiles(:)    = ''
-      cn_logchlfbfiles(:) = ''
-      cn_spmfbfiles(:)    = ''
-      cn_fco2fbfiles(:)   = ''
-      cn_pco2fbfiles(:)   = ''
-      cn_sstbiasfiles(:)  = ''
-      nn_profdavtypes(:)  = -1
+      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 )
@@ -286,8 +358,27 @@
          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 log(Chl) observations              ln_logchl = ', ln_logchl
-         WRITE(numout,*) '             Logical switch for SPM observations                      ln_spm = ', ln_spm
-         WRITE(numout,*) '             Logical switch for FCO2 observations                    ln_fco2 = ', ln_fco2
-         WRITE(numout,*) '             Logical switch for PCO2 observations                    ln_pco2 = ', ln_pco2
+         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
@@ -297,5 +388,21 @@
          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,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint
+         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
@@ -314,7 +421,12 @@
       !-----------------------------------------------------------------------
 
-      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
-      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, &
-         &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) )
+      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
@@ -337,12 +449,59 @@
          IF (ln_t3d .OR. ln_s3d) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', &
-               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
+            cobstypesprof(jtype) = 'prof'
+            clproffiles(jtype,:) = cn_profbfiles
          ENDIF
          IF (ln_vel3d) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', &
-               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
-         ENDIF
+            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
@@ -363,88 +522,153 @@
          IF (ln_sla) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', &
-               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      &
-               &                  nn_2dint, nn_2dint_sla,             &
-               &                  rn_sla_avglamscl, rn_sla_avgphiscl, &
-               &                  ln_sla_fp_indegs, .FALSE.,          &
-               &                  n2dintsurf, ravglamscl, ravgphiscl, &
-               &                  lfpindegs, llnightav )
+            cobstypessurf(jtype) = 'sla'
+            clsurffiles(jtype,:) = cn_slafbfiles
          ENDIF
          IF (ln_sst) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', &
-               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      &
-               &                  nn_2dint, nn_2dint_sst,             &
-               &                  rn_sst_avglamscl, rn_sst_avgphiscl, &
-               &                  ln_sst_fp_indegs, ln_sstnight,      &
-               &                  n2dintsurf, ravglamscl, ravgphiscl, &
-               &                  lfpindegs, llnightav )
-         ENDIF
-#if defined key_lim2 || defined key_lim3 || defined key_cice
+            cobstypessurf(jtype) = 'sst'
+            clsurffiles(jtype,:) = cn_sstfbfiles
+         ENDIF
          IF (ln_sic) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', &
-               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      &
-               &                  nn_2dint, nn_2dint_sic,             &
-               &                  rn_sic_avglamscl, rn_sic_avgphiscl, &
-               &                  ln_sic_fp_indegs, .FALSE.,          &
-               &                  n2dintsurf, ravglamscl, ravgphiscl, &
-               &                  lfpindegs, llnightav )
-         ENDIF
-#endif
+            cobstypessurf(jtype) = 'sic'
+            clsurffiles(jtype,:) = cn_sicfbfiles
+         ENDIF
          IF (ln_sss) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', &
-               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      &
-               &                  nn_2dint, nn_2dint_sss,             &
-               &                  rn_sss_avglamscl, rn_sss_avgphiscl, &
-               &                  ln_sss_fp_indegs, .FALSE.,          &
-               &                  n2dintsurf, ravglamscl, ravgphiscl, &
-               &                  lfpindegs, llnightav )
-         ENDIF
-
-         IF (ln_logchl) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', &
-               &                   cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'logchl',         &
-               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
-               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
+            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
+                  n2dint_type  = nn_2dint_sla
+               ENDIF
+               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
+                  n2dint_type  = nn_2dint_sst
+               ENDIF
+               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 )
-         ENDIF
-
-         IF (ln_spm) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm   ', &
-               &                   cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'spm   ',         &
-               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
-               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
-               &                    lfpindegs, llnightav )
-         ENDIF
-
-         IF (ln_fco2) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2  ', &
-               &                   cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'fco2  ',         &
-               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
-               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
-               &                    lfpindegs, llnightav )
-         ENDIF
-
-         IF (ln_pco2) THEN
-            jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2  ', &
-               &                   cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'pco2  ',         &
-               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
-               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
-               &                    lfpindegs, llnightav )
-         ENDIF
+
+         END DO
 
       ENDIF
@@ -467,6 +691,6 @@
       ENDIF
 
-      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) 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
@@ -491,26 +715,45 @@
          DO jtype = 1, nproftypes
 
-            nvarsprof(jtype) = 2
             IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
+               nvarsprof(jtype) = 2
                nextrprof(jtype) = 1
-               llvar1 = ln_t3d
-               llvar2 = ln_s3d
-               zglam1 = glamt
-               zgphi1 = gphit
-               zmask1 = tmask
-               zglam2 = glamt
-               zgphi2 = gphit
-               zmask2 = tmask
-            ENDIF
-            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
+               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
-               llvar1 = ln_vel3d
-               llvar2 = ln_vel3d
-               zglam1 = glamu
-               zgphi1 = gphiu
-               zmask1 = umask
-               zglam2 = glamv
-               zgphi2 = gphiv
-               zmask2 = vmask
+               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
 
@@ -519,5 +762,5 @@
                &               clproffiles(jtype,1:ifilesprof(jtype)), &
                &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
-               &               rn_dobsini, rn_dobsend, llvar1, llvar2, &
+               &               rn_dobsini, rn_dobsend, llvar, &
                &               ln_ignmis, ln_s_at_t, .FALSE., &
                &               kdailyavtypes = nn_profdavtypes )
@@ -528,9 +771,14 @@
 
             CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
-               &               llvar1, llvar2, &
+               &               llvar, &
                &               jpi, jpj, jpk, &
-               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
+               &               zmask, zglam, zgphi,  &
                &               ln_nea, ln_bound_reject, &
                &               kdailyavtypes = nn_profdavtypes )
+            
+            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
@@ -587,11 +835,4 @@
 
       ENDIF
-
-      CALL wrk_dealloc( jpi, jpj, zglam1 )
-      CALL wrk_dealloc( jpi, jpj, zglam2 )
-      CALL wrk_dealloc( jpi, jpj, zgphi1 )
-      CALL wrk_dealloc( jpi, jpj, zgphi2 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )
 
    END SUBROUTINE dia_obs_init
@@ -639,19 +880,24 @@
 #endif
 #if defined key_hadocc
-      USE trc, ONLY :  &           ! HadOCC chlorophyll, fCO2 and pCO2
+      USE trc, ONLY :  &           ! HadOCC variables
+         & trn, &
          & HADOCC_CHL, &
          & HADOCC_FCO2, &
          & HADOCC_PCO2, &
          & HADOCC_FILL_FLT
-#elif defined key_medusa && defined key_foam_medusa
-      USE trc, ONLY :  &           ! MEDUSA chlorophyll, fCO2 and pCO2
+      USE par_hadocc
+      USE had_bgc_const, ONLY: c2n_p
+#elif defined key_medusa && defined key_foam_medusa
+      USE trc, ONLY :  &           ! MEDUSA variables
          & trn
-      USE par_medusa, ONLY: &
-         & jpchn, &
-         & jpchd
+      USE par_medusa
+      USE sms_medusa, ONLY: &
+         & xthetapn, &
+         & xthetapd
 #if defined key_roam
       USE sms_medusa, ONLY: &
          & f2_pco2w, &
-         & f2_fco2w
+         & f2_fco2w, &
+         & f3_pH
 #endif
 #elif defined key_fabm
@@ -674,33 +920,18 @@
       INTEGER :: jtype             ! Data loop variable
       INTEGER :: jvar              ! Variable number
-      INTEGER :: ji, jj            ! Loop counters
-      REAL(wp) :: tiny                  ! small number
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & zprofvar1, &            ! Model values for 1st variable in a prof ob
-         & zprofvar2               ! Model values for 2nd variable in a prof ob
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & zprofmask1, &           ! Mask associated with zprofvar1
-         & zprofmask2              ! Mask associated with zprofvar2
+      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(:,:) :: &
-         & zglam1,    &            ! Model longitudes for prof variable 1
-         & zglam2,    &            ! Model longitudes for prof variable 2
-         & zgphi1,    &            ! Model latitudes for prof variable 1
-         & zgphi2                  ! Model latitudes for prof variable 2
-
-
-      !Allocate local work arrays
-      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
-      CALL wrk_alloc( jpi, jpj, zsurfvar )
-      CALL wrk_alloc( jpi, jpj, zsurfmask )
-      CALL wrk_alloc( jpi, jpj, zglam1 )
-      CALL wrk_alloc( jpi, jpj, zglam2 )
-      CALL wrk_alloc( jpi, jpj, zgphi1 )
-      CALL wrk_alloc( jpi, jpj, zgphi2 )
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zglam,    &             ! Model longitudes for prof variables
+         & zgphi                   ! Model latitudes for prof variables
+      LOGICAL :: llog10            ! Perform log10 transform of variable
+
 
       IF(lwp) THEN
@@ -721,35 +952,212 @@
          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')
-               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
-               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zprofmask2(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zglam2(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
-               zgphi2(:,:) = gphit(:,:)
+               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
+               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
+
             CASE('vel')
-               zprofvar1(:,:,:) = un(:,:,:)
-               zprofvar2(:,:,:) = vn(:,:,:)
-               zprofmask1(:,:,:) = umask(:,:,:)
-               zprofmask2(:,:,:) = vmask(:,:,:)
-               zglam1(:,:) = glamu(:,:)
-               zglam2(:,:) = glamv(:,:)
-               zgphi1(:,:) = gphiu(:,:)
-               zgphi2(:,:) = gphiv(:,:)
+               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 && defined key_foam_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_chl1) + trn(:,:,:,jp_fabm_chl2) + &
+                  &                trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,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 && defined key_foam_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_chl1) + trn(:,:,:,jp_fabm_chl2) + &
+                  &                trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,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 && defined key_foam_medusa
+               ! Dissolved inorganic nitrogen from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdin)
+#elif defined key_fabm
+               ! Nitrate from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,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 && defined key_foam_medusa
+               ! Silicate from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpsil)
+#elif defined key_fabm
+               ! Silicate from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,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 && defined key_foam_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_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 && defined key_foam_medusa
+               ! Dissolved inorganic carbon from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdic)
+#elif defined key_fabm
+               ! Dissolved inorganic carbon from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,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 && defined key_foam_medusa
+               ! Alkalinity from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpalk)
+#elif defined key_fabm
+               ! Alkalinity from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,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_foam_medusa
+               ! pH from MEDUSA
+               zprofvar(:,:,:,1) = f3_pH(:,:,:)
+#elif defined key_fabm
+               ! pH from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,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 && defined key_foam_medusa
+               ! Oxygen from MEDUSA
+               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)
+#elif defined key_fabm
+               ! Oxygen from ERSEM
+               zprofvar(:,:,:,1) = trn(:,:,:,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
 
-            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
-               &               nit000, idaystp,                         &
-               &               zprofvar1, zprofvar2,                    &
-               &               fsdept(:,:,:), fsdepw(:,:,:),            & 
-               &               zprofmask1, zprofmask2,                  &
-               &               zglam1, zglam2, zgphi1, zgphi2,          &
-               &               nn_1dint, nn_2dint,                      &
-               &               kdailyavtypes = nn_profdavtypes )
+            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
@@ -758,4 +1166,8 @@
 
       IF ( nsurftypes > 0 ) THEN
+
+         !Allocate local work arrays
+         CALL wrk_alloc( jpi, jpj, zsurfvar )
+         CALL wrk_alloc( jpi, jpj, zsurfmask )
 
          DO jtype = 1, nsurftypes
@@ -763,4 +1175,5 @@
             !Defaults which might be changed
             zsurfmask(:,:) = tmask(:,:,1)
+            llog10 = .FALSE.
 
             SELECT CASE ( TRIM(cobstypessurf(jtype)) )
@@ -793,27 +1206,186 @@
                ENDIF
 
-            CASE('logchl')
-#if defined key_hadocc
-               zsurfvar(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
+            CASE('slchltot')
+#if defined key_hadocc
+               ! Surface chlorophyll from HadOCC
+               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
 #elif defined key_medusa && defined key_foam_medusa
                ! Add non-diatom and diatom surface chlorophyll from MEDUSA
                zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
 #elif defined key_fabm
-               chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
-               zsurfvar(:,:) = chl_3d(:,:,1)
-#else
-               CALL ctl_stop( ' Trying to run logchl observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               zsurfmask(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
-               ! 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
-            CASE('spm')
+               ! Add all surface chlorophyll groups from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) + trn(:,:,1,jp_fabm_chl2) + &
+                  &            trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,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 && defined key_foam_medusa
+               ! Diatom surface chlorophyll from MEDUSA
+               zsurfvar(:,:) = trn(:,:,1,jpchd)
+#elif defined key_fabm
+               ! Diatom surface chlorophyll from ERSEM
+               zsurfvar(:,:) = trn(:,:,1,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 && defined key_foam_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_chl2) + &
+                  &            trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,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 && defined key_foam_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_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 && defined key_foam_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_chl1) + trn(:,:,1,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 && defined key_foam_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_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 && defined key_foam_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_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 && defined key_foam_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_chl1) + trn(:,:,1,jp_fabm_chl2) + &
+                  &            trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,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 && defined key_foam_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_p1c) + trn(:,:,1,jp_fabm_p2c) + &
+                  &            trn(:,:,1,jp_fabm_p3c) + trn(:,:,1,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 && defined key_foam_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_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 && defined key_foam_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_p2c) + &
+                  &            trn(:,:,1,jp_fabm_p3c) + trn(:,:,1,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
@@ -822,8 +1394,9 @@
                END DO
 #else
-               CALL ctl_stop( ' Trying to run spm observation operator', &
+               CALL ctl_stop( ' Trying to run sspm observation operator', &
                   &           ' but no spm model appears to have been defined' )
 #endif
-            CASE('fco2')
+
+            CASE('sfco2')
 #if defined key_hadocc
                zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
@@ -861,8 +1434,9 @@
                   &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
 #else
-               CALL ctl_stop( ' Trying to run fco2 observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-            CASE('pco2')
+               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
@@ -880,5 +1454,5 @@
                zsurfvar(:,:) = pco2_3d(:,:,1)
 #else
-               CALL ctl_stop( ' Trying to run pCO2 observation operator', &
+               CALL ctl_stop( ' Trying to run spco2 observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
@@ -889,4 +1463,15 @@
 
             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,       &
@@ -898,16 +1483,8 @@
          END DO
 
+         CALL wrk_dealloc( jpi, jpj, zsurfvar )
+         CALL wrk_dealloc( jpi, jpj, zsurfmask )
+
       ENDIF
-
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
-      CALL wrk_dealloc( jpi, jpj, zsurfvar )
-      CALL wrk_dealloc( jpi, jpj, zsurfmask )
-      CALL wrk_dealloc( jpi, jpj, zglam1 )
-      CALL wrk_dealloc( jpi, jpj, zglam2 )
-      CALL wrk_dealloc( jpi, jpj, zgphi1 )
-      CALL wrk_dealloc( jpi, jpj, zgphi2 )
 
    END SUBROUTINE dia_obs
@@ -960,5 +1537,5 @@
                   & )
 
-               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
+               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
 
                DO jo = 1, profdataqc(jtype)%nprof
@@ -1193,43 +1770,40 @@
     END SUBROUTINE fin_date
 
-    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
-       &                         cfilestype, 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, INTENT(IN) :: jtype       ! Index of the current type of obs
-    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                   ifiles      ! Out appended number of files for this type
-
-    CHARACTER(len=6), INTENT(IN) :: ctypein 
-    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
-       &                   cfilestype  ! In list of files for this obs type
-    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                   cobstypes   ! Out appended list of obs types
-    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
-       &                   cfiles      ! Out appended list of files for all types
-
-    !Local variables
-    INTEGER :: jfile
-
-    cfiles(jtype,:) = cfilestype(:)
-    cobstypes(jtype) = ctypein
-    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(ctypein)//   &
-            &           ' 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))
+    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
-    ENDIF
 
     END SUBROUTINE obs_settypefiles
@@ -1242,67 +1816,70 @@
                &                  lfpindegs, lavnight )
 
-    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
-    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
-    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
-    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
-    REAL(wp), INTENT(IN) :: &
-       &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
-       &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
-    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
-    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
-    CHARACTER(len=6), INTENT(IN) :: ctypein 
-
-    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                    n2dint 
-    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                    ravglamscl, ravgphiscl
-    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                    lfpindegs, lavnight
-
-    lavnight(jtype) = lavnight_type
-
-    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
-       n2dint(jtype) = n2dint_type
-    ELSE
-       n2dint(jtype) = n2dint_default
-    ENDIF
-
-    ! For averaging observation footprints set options for size of footprint 
-    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
-       IF ( ravglamscl_type > 0._wp ) THEN
-          ravglamscl(jtype) = ravglamscl_type
+       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( 'Incorrect value set for averaging footprint '// &
-                         'scale (ravglamscl) for observation type '//TRIM(ctypein) )      
+          CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
+            &                    ' is not available')
        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) )      
+       ! 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
 
-       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)'
+       ! 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
-    ENDIF
 
     END SUBROUTINE obs_setinterpopts
Index: /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 9305)
+++ /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 9306)
@@ -60,9 +60,9 @@
 
 
-   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          &
-      &                     kit000, kdaystp,                      &
-      &                     pvar1, pvar2, pgdept, pgdepw,         &
-      &                     pmask1, pmask2,                       &  
-      &                     plam1, plam2, pphi1, pphi2,           &
+   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
+      &                     kit000, kdaystp, kvar,       &
+      &                     pvar, pgdept, pgdepw,        &
+      &                     pmask,                       &  
+      &                     plam, pphi,                  &
       &                     k1dint, k2dint, kdailyavtypes )
 
@@ -134,14 +134,11 @@
       INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header)
       INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
+      INTEGER, INTENT(IN) :: kvar    ! Number of variable in prodatqc
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
-         & pvar1,    &               ! Model field 1
-         & pvar2,    &               ! Model field 2
-         & pmask1,   &               ! Land-sea mask 1
-         & pmask2                    ! Land-sea mask 2
+         & pvar,   &                 ! Model field for variable
+         & pmask                     ! Land-sea mask for variable
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & plam1,    &               ! Model longitudes for variable 1
-         & plam2,    &               ! Model longitudes for variable 2
-         & pphi1,    &               ! Model latitudes for variable 1
-         & pphi2                     ! Model latitudes for variable 2
+         & plam,   &                 ! Model longitudes for variable
+         & pphi                      ! Model latitudes for variable
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
          & pgdept, &                 ! Model array of depth T levels 
@@ -166,8 +163,6 @@
          & idailyavtypes
       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi1, &
-         & igrdi2, &
-         & igrdj1, &
-         & igrdj2
+         & igrdi, &
+         & igrdj
       INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
 
@@ -176,27 +171,19 @@
       REAL(KIND=wp) :: zdaystp
       REAL(KIND=wp), DIMENSION(kpk) :: &
-         & zobsmask1, &
-         & zobsmask2, &
          & zobsk,    &
          & zobs2k
       REAL(KIND=wp), DIMENSION(2,2,1) :: &
          & zweig1, &
-         & zweig2, &
          & zweig
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
-         & zmask1, &
-         & zmask2, &
-         & zint1,  &
-         & zint2,  &
-         & zinm1,  &
-         & zinm2,  &
+         & zmask,  &
+         & zint,   &
+         & zinm,   &
          & zgdept, & 
          & zgdepw
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zglam1, &
-         & zglam2, &
-         & zgphi1, &
-         & zgphi2
-      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2   
+         & zglam,  &
+         & zgphi
+      REAL(KIND=wp), DIMENSION(1) :: zmsk
       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
 
@@ -230,6 +217,5 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     prodatqc%vdmean(ji,jj,jk,1) = 0.0
-                     prodatqc%vdmean(ji,jj,jk,2) = 0.0
+                     prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
                   END DO
                END DO
@@ -240,10 +226,7 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! Increment field 1 for computing daily mean
-                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                     &                        + pvar1(ji,jj,jk)
-                  ! Increment field 2 for computing daily mean
-                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                     &                        + pvar2(ji,jj,jk)
+                  ! Increment field for computing daily mean
+                  prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
+                     &                           + pvar(ji,jj,jk)
                END DO
             END DO
@@ -258,8 +241,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                        &                        * zdaystp
-                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                        &                        * zdaystp
+                     prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
+                        &                           * zdaystp
                   END DO
                END DO
@@ -271,16 +252,10 @@
       ! Get the data for interpolation
       ALLOCATE( &
-         & igrdi1(2,2,ipro),      &
-         & igrdi2(2,2,ipro),      &
-         & igrdj1(2,2,ipro),      &
-         & igrdj2(2,2,ipro),      &
-         & zglam1(2,2,ipro),      &
-         & zglam2(2,2,ipro),      &
-         & zgphi1(2,2,ipro),      &
-         & zgphi2(2,2,ipro),      &
-         & zmask1(2,2,kpk,ipro),  &
-         & zmask2(2,2,kpk,ipro),  &
-         & zint1(2,2,kpk,ipro),   &
-         & zint2(2,2,kpk,ipro),   &
+         & igrdi(2,2,ipro),       &
+         & igrdj(2,2,ipro),       &
+         & zglam(2,2,ipro),       &
+         & zgphi(2,2,ipro),       &
+         & zmask(2,2,kpk,ipro),   &
+         & zint(2,2,kpk,ipro),    &
          & zgdept(2,2,kpk,ipro),  & 
          & zgdepw(2,2,kpk,ipro)   & 
@@ -289,20 +264,12 @@
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
          iobs = jobs - prodatqc%nprofup
-         igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1
-         igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1
-         igrdj1(1,2,iobs) = prodatqc%mj(jobs,1)
-         igrdi1(2,1,iobs) = prodatqc%mi(jobs,1)
-         igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1)
-         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1)
-         igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1
-         igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1
-         igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1
-         igrdj2(1,2,iobs) = prodatqc%mj(jobs,2)
-         igrdi2(2,1,iobs) = prodatqc%mi(jobs,2)
-         igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1
-         igrdi2(2,2,iobs) = prodatqc%mi(jobs,2)
-         igrdj2(2,2,iobs) = prodatqc%mj(jobs,2)
+         igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1
+         igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1
+         igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1
+         igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar)
+         igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar)
+         igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1
+         igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar)
+         igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar)
       END DO
 
@@ -311,29 +278,19 @@
       zgdepw(:,:,:,:) = 0.0
 
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 )
-      
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 )
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 )
-
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam )
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar,   zint )
+
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 
 
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
 
-         ALLOCATE( &
-            & zinm1(2,2,kpk,ipro),  &
-            & zinm2(2,2,kpk,ipro)   &
-            & )
-
-         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, &
-            &                  prodatqc%vdmean(:,:,:,1), zinm1 )
-         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
-            &                  prodatqc%vdmean(:,:,:,2), zinm2 )
+         ALLOCATE( zinm(2,2,kpk,ipro) )
+
+         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &
+            &                  prodatqc%vdmean(:,:,:,kvar), zinm )
 
       ENDIF
@@ -370,21 +327,13 @@
          ! Horizontal weights 
          ! Masked values are calculated later.  
-         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
+         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
 
             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
-               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), &
-               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 )
-
-         ENDIF
-
-         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
-
-            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
-               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), &
-               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 )
- 
-         ENDIF
-
-         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
+               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
+               &                   zmask(:,:,1,iobs), zweig1, zmsk )
+
+         ENDIF
+
+         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
 
             zobsk(:) = obfillflt
@@ -396,6 +345,6 @@
 
                   ! vertically interpolate all 4 corners 
-                  ista = prodatqc%npvsta(jobs,1) 
-                  iend = prodatqc%npvend(jobs,1) 
+                  ista = prodatqc%npvsta(jobs,kvar) 
+                  iend = prodatqc%npvend(jobs,kvar) 
                   inum_obs = iend - ista + 1 
                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
@@ -406,20 +355,20 @@
                         IF ( k1dint == 1 ) THEN 
                            CALL obs_int_z1d_spl( kpk, & 
-                              &     zinm1(iin,ijn,:,iobs), & 
+                              &     zinm(iin,ijn,:,iobs), & 
                               &     zobs2k, zgdept(iin,ijn,:,iobs), & 
-                              &     zmask1(iin,ijn,:,iobs)) 
+                              &     zmask(iin,ijn,:,iobs)) 
                         ENDIF 
        
                         CALL obs_level_search(kpk, & 
                            &    zgdept(iin,ijn,:,iobs), & 
-                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                           &    inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 
                            &    iv_indic) 
 
                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
-                           &    prodatqc%var(1)%vdep(ista:iend), & 
-                           &    zinm1(iin,ijn,:,iobs), & 
+                           &    prodatqc%var(kvar)%vdep(ista:iend), & 
+                           &    zinm(iin,ijn,:,iobs), & 
                            &    zobs2k, interp_corner(iin,ijn,:), & 
                            &    zgdept(iin,ijn,:,iobs), & 
-                           &    zmask1(iin,ijn,:,iobs)) 
+                           &    zmask(iin,ijn,:,iobs)) 
        
                      ENDDO 
@@ -433,6 +382,6 @@
      
                ! vertically interpolate all 4 corners 
-               ista = prodatqc%npvsta(jobs,1) 
-               iend = prodatqc%npvend(jobs,1) 
+               ista = prodatqc%npvsta(jobs,kvar) 
+               iend = prodatqc%npvend(jobs,kvar) 
                inum_obs = iend - ista + 1 
                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
@@ -442,7 +391,7 @@
                      IF ( k1dint == 1 ) THEN 
                         CALL obs_int_z1d_spl( kpk, & 
-                           &    zint1(iin,ijn,:,iobs),& 
+                           &    zint(iin,ijn,:,iobs),& 
                            &    zobs2k, zgdept(iin,ijn,:,iobs), & 
-                           &    zmask1(iin,ijn,:,iobs)) 
+                           &    zmask(iin,ijn,:,iobs)) 
   
                      ENDIF 
@@ -450,13 +399,13 @@
                      CALL obs_level_search(kpk, & 
                          &        zgdept(iin,ijn,:,iobs),& 
-                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                         &        inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 
                          &        iv_indic) 
 
                      CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
-                         &          prodatqc%var(1)%vdep(ista:iend),     & 
-                         &          zint1(iin,ijn,:,iobs),            & 
+                         &          prodatqc%var(kvar)%vdep(ista:iend),     & 
+                         &          zint(iin,ijn,:,iobs),            & 
                          &          zobs2k,interp_corner(iin,ijn,:), & 
                          &          zgdept(iin,ijn,:,iobs),         & 
-                         &          zmask1(iin,ijn,:,iobs) )      
+                         &          zmask(iin,ijn,:,iobs) )      
          
                   ENDDO 
@@ -482,17 +431,17 @@
                   DO ijn=1,2 
      
-                     depth_loop1: DO ik=kpk,2,-1 
-                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
+                     depth_loop: DO ik=kpk,2,-1 
+                        IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
                             
                            zweig(iin,ijn,1) = &  
                               & zweig1(iin,ijn,1) * & 
                               & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
-                              &  - prodatqc%var(1)%vdep(iend)),0._wp) 
+                              &  - prodatqc%var(kvar)%vdep(iend)),0._wp) 
                             
-                           EXIT depth_loop1 
+                           EXIT depth_loop
 
                         ENDIF 
 
-                     ENDDO depth_loop1 
+                     ENDDO depth_loop
      
                   ENDDO 
@@ -500,9 +449,9 @@
    
                CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
-                  &              prodatqc%var(1)%vmod(iend:iend) ) 
+                  &              prodatqc%var(kvar)%vmod(iend:iend) ) 
 
                   ! Set QC flag for any observations found below the bottom
                   ! needed as the check here is more strict than that in obs_prep
-               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
+               IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
  
             ENDDO 
@@ -510,149 +459,16 @@
             DEALLOCATE(interp_corner,iv_indic) 
           
-         ENDIF 
-
-         ! For the second variable
-         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
-
-            zobsk(:) = obfillflt
-
-            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
-
-               IF ( idayend == 0 )  THEN
-                  ! Daily averaged data
-
-                  ! vertically interpolate all 4 corners 
-                  ista = prodatqc%npvsta(jobs,2) 
-                  iend = prodatqc%npvend(jobs,2) 
-                  inum_obs = iend - ista + 1 
-                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
-
-                  DO iin=1,2 
-                     DO ijn=1,2 
-
-                        IF ( k1dint == 1 ) THEN 
-                           CALL obs_int_z1d_spl( kpk, & 
-                              &     zinm2(iin,ijn,:,iobs), & 
-                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
-                              &     zmask2(iin,ijn,:,iobs)) 
-                        ENDIF 
-       
-                        CALL obs_level_search(kpk, & 
-                           &    zgdept(iin,ijn,:,iobs), & 
-                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
-                           &    iv_indic) 
-
-                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
-                           &    prodatqc%var(2)%vdep(ista:iend), & 
-                           &    zinm2(iin,ijn,:,iobs), & 
-                           &    zobs2k, interp_corner(iin,ijn,:), & 
-                           &    zgdept(iin,ijn,:,iobs), & 
-                           &    zmask2(iin,ijn,:,iobs)) 
-       
-                     ENDDO 
-                  ENDDO 
-
-               ENDIF !idayend
-
-            ELSE   
-
-               ! Point data 
-     
-               ! vertically interpolate all 4 corners 
-               ista = prodatqc%npvsta(jobs,2) 
-               iend = prodatqc%npvend(jobs,2) 
-               inum_obs = iend - ista + 1 
-               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
-               DO iin=1,2  
-                  DO ijn=1,2 
-                    
-                     IF ( k1dint == 1 ) THEN 
-                        CALL obs_int_z1d_spl( kpk, & 
-                           &    zint2(iin,ijn,:,iobs),& 
-                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
-                           &    zmask2(iin,ijn,:,iobs)) 
-  
-                     ENDIF 
-       
-                     CALL obs_level_search(kpk, & 
-                         &        zgdept(iin,ijn,:,iobs),& 
-                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
-                         &        iv_indic) 
-
-                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
-                         &          prodatqc%var(2)%vdep(ista:iend),     & 
-                         &          zint2(iin,ijn,:,iobs),            & 
-                         &          zobs2k,interp_corner(iin,ijn,:), & 
-                         &          zgdept(iin,ijn,:,iobs),         & 
-                         &          zmask2(iin,ijn,:,iobs) )      
-         
-                  ENDDO 
-               ENDDO 
-             
-            ENDIF 
-
-            !------------------------------------------------------------- 
-            ! Compute the horizontal interpolation for every profile level 
-            !------------------------------------------------------------- 
-             
-            DO ikn=1,inum_obs 
-               iend=ista+ikn-1
-                  
-               zweig(:,:,1) = 0._wp 
-   
-               ! This code forces the horizontal weights to be  
-               ! zero IF the observation is below the bottom of the  
-               ! corners of the interpolation nodes, Or if it is in  
-               ! the mask. This is important for observations near  
-               ! steep bathymetry 
-               DO iin=1,2 
-                  DO ijn=1,2 
-     
-                     depth_loop2: DO ik=kpk,2,-1 
-                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
-                            
-                           zweig(iin,ijn,1) = &  
-                              & zweig2(iin,ijn,1) * & 
-                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
-                              &  - prodatqc%var(2)%vdep(iend)),0._wp) 
-                            
-                           EXIT depth_loop2 
-
-                        ENDIF 
-
-                     ENDDO depth_loop2 
-     
-                  ENDDO 
-               ENDDO 
-   
-               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
-                  &              prodatqc%var(2)%vmod(iend:iend) ) 
-
-                  ! Set QC flag for any observations found below the bottom
-                  ! needed as the check here is more strict than that in obs_prep
-               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
- 
-            ENDDO 
- 
-            DEALLOCATE(interp_corner,iv_indic) 
-          
-         ENDIF 
+         ENDIF
 
       ENDDO
 
       ! Deallocate the data for interpolation
-      DEALLOCATE( &
-         & igrdi1, &
-         & igrdi2, &
-         & igrdj1, &
-         & igrdj2, &
-         & zglam1, &
-         & zglam2, &
-         & zgphi1, &
-         & zgphi2, &
-         & zmask1, &
-         & zmask2, &
-         & zint1,  &
-         & zint2,  &
+      DEALLOCATE(  &
+         & igrdi,  &
+         & igrdj,  &
+         & zglam,  &
+         & zgphi,  &
+         & zmask,  &
+         & zint,   &
          & zgdept, &
          & zgdepw  &
@@ -661,11 +477,10 @@
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
-         DEALLOCATE( &
-            & zinm1,  &
-            & zinm2   &
-            & )
+         DEALLOCATE( zinm )
       ENDIF
 
-      prodatqc%nprofup = prodatqc%nprofup + ipro 
+      IF ( kvar == prodatqc%nvar ) THEN
+         prodatqc%nprofup = prodatqc%nprofup + ipro 
+      ENDIF
 
    END SUBROUTINE obs_prof_opt
Index: /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 9305)
+++ /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 9306)
@@ -255,7 +255,7 @@
 
 
-   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var1, ld_var2, &
+   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, &
       &                     kpi, kpj, kpk, &
-      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  &
+      &                     zmask, pglam, pgphi,  &
       &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff )
 
@@ -284,6 +284,6 @@
       TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data
       TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening
-      LOGICAL, INTENT(IN) :: ld_var1              ! Observed variables switches
-      LOGICAL, INTENT(IN) :: ld_var2
+      LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: &
+         & ld_var                                 ! Observed variables switches
       LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land
       LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary
@@ -291,12 +291,9 @@
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes                          ! Types for daily averages
-      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
-         & zmask1, &
-         & zmask2
-      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & pglam1, &
-         & pglam2, &
-         & pgphi1, &
-         & pgphi2
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: &
+         & zmask
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: &
+         & pglam, &
+         & pgphi
       INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
 
@@ -309,37 +306,31 @@
       INTEGER :: imin0
       INTEGER :: icycle       ! Current assimilation cycle
-                              ! Counters for observations that are
-      INTEGER :: iotdobs      !  - outside time domain
-      INTEGER :: iosdv1obs    !  - outside space domain (variable 1)
-      INTEGER :: iosdv2obs    !  - outside space domain (variable 2)
-      INTEGER :: ilanv1obs    !  - within a model land cell (variable 1)
-      INTEGER :: ilanv2obs    !  - within a model land cell (variable 2)
-      INTEGER :: inlav1obs    !  - close to land (variable 1)
-      INTEGER :: inlav2obs    !  - close to land (variable 2)
-      INTEGER :: ibdyv1obs    !  - boundary (variable 1) 
-      INTEGER :: ibdyv2obs    !  - boundary (variable 2)      
-      INTEGER :: igrdobs      !  - fail the grid search
-      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa
-      INTEGER :: iuvchkv      !
-                              ! Global counters for observations that are
-      INTEGER :: iotdobsmpp   !  - outside time domain
-      INTEGER :: iosdv1obsmpp !  - outside space domain (variable 1)
-      INTEGER :: iosdv2obsmpp !  - outside space domain (variable 2)
-      INTEGER :: ilanv1obsmpp !  - within a model land cell (variable 1)
-      INTEGER :: ilanv2obsmpp !  - within a model land cell (variable 2)
-      INTEGER :: inlav1obsmpp !  - close to land (variable 1)
-      INTEGER :: inlav2obsmpp !  - close to land (variable 2)
-      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1) 
-      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)      
-      INTEGER :: igrdobsmpp   !  - fail the grid search
-      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa
-      INTEGER :: iuvchkvmpp   !
+                                                       ! Counters for observations that are
+      INTEGER                           :: iotdobs     !  - outside time domain
+      INTEGER, DIMENSION(profdata%nvar) :: iosdvobs    !  - outside space domain
+      INTEGER, DIMENSION(profdata%nvar) :: ilanvobs    !  - within a model land cell
+      INTEGER, DIMENSION(profdata%nvar) :: inlavobs    !  - close to land
+      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs    !  - boundary   
+      INTEGER                           :: igrdobs     !  - fail the grid search
+      INTEGER                           :: iuvchku     !  - reject UVEL if VVEL rejected
+      INTEGER                           :: iuvchkv     !  - reject VVEL if UVEL rejected
+                                                       ! Global counters for observations that are
+      INTEGER                           :: iotdobsmpp  !  - outside time domain
+      INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp !  - outside space domain
+      INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp !  - within a model land cell
+      INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp !  - close to land
+      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp !  - boundary
+      INTEGER :: igrdobsmpp                            !  - fail the grid search
+      INTEGER :: iuvchkumpp                            !  - reject UVEL if VVEL rejected
+      INTEGER :: iuvchkvmpp                            !  - reject VVEL if UVEL rejected
       TYPE(obs_prof_valid) ::  llvalid      ! Profile selection 
       TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
-         & llvvalid           ! var1,var2 selection 
+         & llvvalid           ! vars selection 
       INTEGER :: jvar         ! Variable loop variable
       INTEGER :: jobs         ! Obs. loop variable
       INTEGER :: jstp         ! Time loop variable
       INTEGER :: inrc         ! Time index variable
+      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line
+      CHARACTER(LEN=256) :: cout2  ! Diagnostic output line
 
       IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...'
@@ -355,18 +346,14 @@
       icycle = no     ! Assimilation cycle
 
-      ! Diagnotics counters for various failures.
-
-      iotdobs   = 0
-      igrdobs   = 0
-      iosdv1obs = 0
-      iosdv2obs = 0
-      ilanv1obs = 0
-      ilanv2obs = 0
-      inlav1obs = 0
-      inlav2obs = 0
-      ibdyv1obs = 0
-      ibdyv2obs = 0
-      iuvchku   = 0
-      iuvchkv   = 0
+      ! Diagnostics counters for various failures.
+
+      iotdobs     = 0
+      igrdobs     = 0
+      iosdvobs(:) = 0
+      ilanvobs(:) = 0
+      inlavobs(:) = 0
+      ibdyvobs(:) = 0
+      iuvchku     = 0
+      iuvchkv     = 0
 
 
@@ -401,8 +388,8 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), &
-         &              profdata%nqc,     igrdobs                         )
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), &
-         &              profdata%nqc,     igrdobs                         )
+      DO jvar = 1, profdata%nvar
+         CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,jvar), profdata%mj(:,jvar), &
+            &              profdata%nqc,     igrdobs )
+      END DO
 
       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
@@ -419,47 +406,26 @@
       ! -----------------------------------------------------------------------
 
-      ! Variable 1
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   &
-         &                 profdata%npvsta(:,1),  profdata%npvend(:,1), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          &
-         &                 profdata%var(1)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(1)%vdep,                        &
-         &                 pglam1,                pgphi1,               &
-         &                 gdept_1d,              zmask1,               &
-         &                 profdata%nqc,          profdata%var(1)%nvqc, &
-         &                 iosdv1obs,             ilanv1obs,            &
-         &                 inlav1obs,             ld_nea,               &
-         &                 ibdyv1obs,             ld_bound_reject,      &
-         &                 iqc_cutoff       )
-
-      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp )
-      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp )
-      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp )
-      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp )
-
-      ! Variable 2
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   &
-         &                 profdata%npvsta(:,2),  profdata%npvend(:,2), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          & 
-         &                 profdata%var(2)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(2)%vdep,                        &
-         &                 pglam2,                pgphi2,               &
-         &                 gdept_1d,              zmask2,               &
-         &                 profdata%nqc,          profdata%var(2)%nvqc, &
-         &                 iosdv2obs,             ilanv2obs,            &
-         &                 inlav2obs,             ld_nea,               &
-         &                 ibdyv2obs,             ld_bound_reject,      &
-         &                 iqc_cutoff       )
-
-      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp )
-      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp )
-      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp )
-      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp )
+      DO jvar = 1, profdata%nvar
+         CALL obs_coo_spc_3d( profdata%nprof,          profdata%nvprot(jvar),   &
+            &                 profdata%npvsta(:,jvar), profdata%npvend(:,jvar), &
+            &                 jpi,                     jpj,                     &
+            &                 jpk,                                              &
+            &                 profdata%mi,             profdata%mj,             &
+            &                 profdata%var(jvar)%mvk,                           &
+            &                 profdata%rlam,           profdata%rphi,           &
+            &                 profdata%var(jvar)%vdep,                          &
+            &                 pglam(:,:,jvar),         pgphi(:,:,jvar),         &
+            &                 gdept_1d,                zmask(:,:,:,jvar),       &
+            &                 profdata%nqc,            profdata%var(jvar)%nvqc, &
+            &                 iosdvobs(jvar),          ilanvobs(jvar),          &
+            &                 inlavobs(jvar),          ld_nea,                  &
+            &                 ibdyvobs(jvar),          ld_bound_reject,         &
+            &                 iqc_cutoff       )
+
+         CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) )
+      END DO
 
       ! -----------------------------------------------------------------------
@@ -512,54 +478,44 @@
       
          WRITE(numout,*)
-         WRITE(numout,*) ' Profiles outside time domain                     = ', &
+         WRITE(numout,*) ' Profiles outside time domain                       = ', &
             &            iotdobsmpp
-         WRITE(numout,*) ' Remaining profiles that failed grid search       = ', &
+         WRITE(numout,*) ' Remaining profiles that failed grid search         = ', &
             &            igrdobsmpp
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain       = ', &
-            &            iosdv1obsmpp
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points             = ', &
-            &            ilanv1obsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',&
-               &            inlav1obsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept)    = ',&
-               &            inlav1obsmpp
-         ENDIF
-         IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
-            WRITE(numout,*) ' U observation rejected since V rejected     = ', &
-               &            iuvchku
-         ENDIF
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',&
-               &            ibdyv1obsmpp
-         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', &
-            &            prodatqc%nvprotmpp(1)
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain       = ', &
-            &            iosdv2obsmpp
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points             = ', &
-            &            ilanv2obsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',&
-               &            inlav2obsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept)    = ',&
-               &            inlav2obsmpp
-         ENDIF
-         IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
-            WRITE(numout,*) ' V observation rejected since U rejected     = ', &
-               &            iuvchkv
-         ENDIF
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',&
-               &            ibdyv2obsmpp
-         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', &
-            &            prodatqc%nvprotmpp(2)
+         DO jvar = 1, profdata%nvar
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain       = ', &
+               &            iosdvobsmpp(jvar)
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points             = ', &
+               &            ilanvobsmpp(jvar)
+            IF (ld_nea) THEN
+               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ',&
+                  &            inlavobsmpp(jvar)
+            ELSE
+               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept)    = ',&
+                  &            inlavobsmpp(jvar)
+            ENDIF
+            IF ( TRIM(profdata%cvars(jvar)) == 'UVEL' ) THEN
+               WRITE(numout,*) ' U observation rejected since V rejected     = ', &
+                  &            iuvchku
+            ELSE IF ( TRIM(profdata%cvars(jvar)) == 'VVEL' ) THEN
+               WRITE(numout,*) ' V observation rejected since U rejected     = ', &
+                  &            iuvchkv
+            ENDIF
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ',&
+                  &            ibdyvobsmpp(jvar)
+            WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted                             = ', &
+               &            prodatqc%nvprotmpp(jvar)
+         END DO
 
          WRITE(numout,*)
          WRITE(numout,*) ' Number of observations per time step :'
          WRITE(numout,*)
-         WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', &
-            &                               '     '//prodatqc%cvars(1)//'     ', &
-            &                               '     '//prodatqc%cvars(2)//'     '
-         WRITE(numout,998)
+         WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles'
+         WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------'
+         DO jvar = 1, prodatqc%nvar
+            WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar))
+            WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------'
+         END DO
+         WRITE(numout,*) cout1
+         WRITE(numout,*) cout2
       ENDIF
       
@@ -588,12 +544,11 @@
          DO jstp = nit000 - 1, nitend
             inrc = jstp - nit000 + 2
-            WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
-               &                    prodatqc%nvstpmpp(inrc,1), &
-               &                    prodatqc%nvstpmpp(inrc,2)
+            WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc)
+            DO jvar = 1, prodatqc%nvar
+               WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar)
+            END DO
+            WRITE(numout,*) cout1
          END DO
       ENDIF
-
-998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
-999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
 
    END SUBROUTINE obs_pre_prof
Index: /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 9305)
+++ /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 9306)
@@ -45,5 +45,5 @@
    SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, &
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                     ldvar1, ldvar2, ldignmis, ldsatt, &
+      &                     ldvar, ldignmis, ldsatt, &
       &                     ldmod, kdailyavtypes )
       !!---------------------------------------------------------------------
@@ -74,6 +74,5 @@
       INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var
       INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
-      LOGICAL, INTENT(IN) :: ldvar1     ! Observed variables switches
-      LOGICAL, INTENT(IN) :: ldvar2
+      LOGICAL, DIMENSION(kvars), INTENT(IN) :: ldvar     ! Observed variables switches
       LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
       LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points
@@ -105,29 +104,22 @@
       INTEGER :: iprof
       INTEGER :: iproftot
-      INTEGER :: ivar1t0
-      INTEGER :: ivar2t0
-      INTEGER :: ivar1t
-      INTEGER :: ivar2t
+      INTEGER, DIMENSION(kvars) :: ivart0
+      INTEGER, DIMENSION(kvars) :: ivart
       INTEGER :: ip3dt
       INTEGER :: ios
       INTEGER :: ioserrcount
-      INTEGER :: ivar1tmpp
-      INTEGER :: ivar2tmpp
+      INTEGER, DIMENSION(kvars) :: ivartmpp
       INTEGER :: ip3dtmpp
       INTEGER :: itype
       INTEGER, DIMENSION(knumfiles) :: &
          & irefdate
-      INTEGER, DIMENSION(ntyp1770+1) :: &
-         & itypvar1,    &
-         & itypvar1mpp, &
-         & itypvar2,    &
-         & itypvar2mpp 
+      INTEGER, DIMENSION(ntyp1770+1,kvars) :: &
+         & itypvar,    &
+         & itypvarmpp
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
+         & iobsi,    &
+         & iobsj,    &
+         & iproc
       INTEGER, DIMENSION(:), ALLOCATABLE :: &
-         & iobsi1,    &
-         & iobsj1,    &
-         & iproc1,    &
-         & iobsi2,    &
-         & iobsj2,    &
-         & iproc2,    &
          & iindx,    &
          & ifileidx, &
@@ -147,4 +139,5 @@
       LOGICAL :: llvalprof
       LOGICAL :: lldavtimset
+      LOGICAL :: llcycle
       TYPE(obfbdata), POINTER, DIMENSION(:) :: &
          & inpfiles
@@ -152,6 +145,5 @@
       ! Local initialization
       iprof = 0
-      ivar1t0 = 0
-      ivar2t0 = 0
+      ivart0(:) = 0
       ip3dt = 0
 
@@ -219,7 +211,7 @@
                &                ldgrid = .TRUE. )
 
-            IF ( inpfiles(jj)%nvar < 2 ) THEN
+            IF ( inpfiles(jj)%nvar /= kvars ) THEN
                CALL ctl_stop( 'Feedback format error: ', &
-                  &           ' less than 2 vars in profile file' )
+                  &           ' unexpected number of vars in profile file' )
             ENDIF
 
@@ -308,6 +300,12 @@
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -317,15 +315,18 @@
             ALLOCATE( zlam(inowin)  )
             ALLOCATE( zphi(inowin)  )
-            ALLOCATE( iobsi1(inowin) )
-            ALLOCATE( iobsj1(inowin) )
-            ALLOCATE( iproc1(inowin) )
-            ALLOCATE( iobsi2(inowin) )
-            ALLOCATE( iobsj2(inowin) )
-            ALLOCATE( iproc2(inowin) )
+            ALLOCATE( iobsi(inowin,kvars) )
+            ALLOCATE( iobsj(inowin,kvars) )
+            ALLOCATE( iproc(inowin,kvars) )
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -336,15 +337,20 @@
             END DO
 
-            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN
-               CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, &
-                  &                  iproc1, 'T' )
-               iobsi2(:) = iobsi1(:)
-               iobsj2(:) = iobsj1(:)
-               iproc2(:) = iproc1(:)
-            ELSEIF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN
-               CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, &
-                  &                  iproc1, 'U' )
-               CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, &
-                  &                  iproc2, 'V' )
+            ! Assume anything other than velocity is on T grid
+            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
+                  &                  iproc(:,1), 'U' )
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), &
+                  &                  iproc(:,2), 'V' )
+            ELSE
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
+                  &                  iproc(:,1), 'T' )
+               IF ( kvars > 1 ) THEN
+                  DO jvar = 2, kvars
+                     iobsi(:,jvar) = iobsi(:,1)
+                     iobsj(:,jvar) = iobsj(:,1)
+                     iproc(:,jvar) = iproc(:,1)
+                  END DO
+               ENDIF
             ENDIF
 
@@ -352,28 +358,43 @@
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
                   inowin = inowin + 1
-                  inpfiles(jj)%iproc(ji,1) = iproc1(inowin)
-                  inpfiles(jj)%iobsi(ji,1) = iobsi1(inowin)
-                  inpfiles(jj)%iobsj(ji,1) = iobsj1(inowin)
-                  inpfiles(jj)%iproc(ji,2) = iproc2(inowin)
-                  inpfiles(jj)%iobsi(ji,2) = iobsi2(inowin)
-                  inpfiles(jj)%iobsj(ji,2) = iobsj2(inowin)
-                  IF ( inpfiles(jj)%iproc(ji,1) /= &
-                     & inpfiles(jj)%iproc(ji,2) ) THEN
-                     CALL ctl_stop( 'Error in obs_read_prof:', &
-                        & 'var1 and var2 observation on different processors')
+                  DO jvar = 1, kvars
+                     inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar)
+                     inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar)
+                     inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar)
+                  END DO
+                  IF ( kvars > 1 ) THEN
+                     DO jvar = 2, kvars
+                        IF ( inpfiles(jj)%iproc(ji,jvar) /= &
+                           & inpfiles(jj)%iproc(ji,1) ) THEN
+                           CALL ctl_stop( 'Error in obs_read_prof:', &
+                              & 'observation on different processors for different vars')
+                        ENDIF
+                     END DO
                   ENDIF
                ENDIF
             END DO
-            DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 )
+            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
 
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -384,37 +405,29 @@
                   ENDIF
                   llvalprof = .FALSE.
-                  IF ( ldvar1 ) THEN
-                     loop_t_count : DO ij = 1,inpfiles(jj)%nlev
-                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
-                           & CYCLE
-                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-                           ivar1t0 = ivar1t0 + 1
-                        ENDIF
-                     END DO loop_t_count
-                  ENDIF
-                  IF ( ldvar2 ) THEN
-                     loop_s_count : DO ij = 1,inpfiles(jj)%nlev
-                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
-                           & CYCLE
-                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-                           ivar2t0 = ivar2t0 + 1
-                        ENDIF
-                     END DO loop_s_count
-                  ENDIF
-                  loop_p_count : DO ij = 1,inpfiles(jj)%nlev
+                  DO jvar = 1, kvars
+                     IF ( ldvar(jvar) ) THEN
+                        DO ij = 1,inpfiles(jj)%nlev
+                           IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
+                              & CYCLE
+                           IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                              & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
+                              ivart0(jvar) = ivart0(jvar) + 1
+                           ENDIF
+                        END DO
+                     ENDIF
+                  END DO
+                  DO ij = 1,inpfiles(jj)%nlev
                      IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                         & CYCLE
-                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &    ldvar1 ) .OR. &
-                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &     ldvar2 ) ) THEN
-                        ip3dt = ip3dt + 1
-                        llvalprof = .TRUE.
-                     ENDIF
-                  END DO loop_p_count
+                     DO jvar = 1, kvars
+                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                           &    ldvar(jvar) ) ) THEN
+                           ip3dt = ip3dt + 1
+                           llvalprof = .TRUE.
+                           EXIT
+                        ENDIF
+                     END DO
+                  END DO
 
                   IF ( llvalprof ) iprof = iprof + 1
@@ -438,6 +451,12 @@
          DO ji = 1, inpfiles(jj)%nobs
             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            llcycle = .TRUE.
+            DO jvar = 1, kvars
+               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                  llcycle = .FALSE.
+                  EXIT
+               ENDIF
+            END DO
+            IF ( llcycle ) CYCLE
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -453,6 +472,12 @@
          DO ji = 1, inpfiles(jj)%nobs
             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            llcycle = .TRUE.
+            DO jvar = 1, kvars
+               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                  llcycle = .FALSE.
+                  EXIT
+               ENDIF
+            END DO
+            IF ( llcycle ) CYCLE
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -470,9 +495,7 @@
       iv3dt(:) = -1
       IF (ldsatt) THEN
-         iv3dt(1) = ip3dt
-         iv3dt(2) = ip3dt
+         iv3dt(:) = ip3dt
       ELSE
-         iv3dt(1) = ivar1t0
-         iv3dt(2) = ivar2t0
+         iv3dt(:) = ivart0(:)
       ENDIF
       CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
@@ -487,11 +510,7 @@
 
       ip3dt = 0
-      ivar1t = 0
-      ivar2t = 0
-      itypvar1   (:) = 0
-      itypvar1mpp(:) = 0
-
-      itypvar2   (:) = 0
-      itypvar2mpp(:) = 0
+      ivart(:) = 0
+      itypvar   (:,:) = 0
+      itypvarmpp(:,:) = 0
 
       ioserrcount = 0
@@ -501,7 +520,13 @@
          ji = iprofidx(iindx(jk))
 
-            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+         IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+         llcycle = .TRUE.
+         DO jvar = 1, kvars
+            IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+               llcycle = .FALSE.
+               EXIT
+            ENDIF
+         END DO
+         IF ( llcycle ) CYCLE
 
          IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
@@ -519,6 +544,12 @@
 
             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            llcycle = .TRUE.
+            DO jvar = 1, kvars
+               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                  llcycle = .FALSE.
+                  EXIT
+               ENDIF
+            END DO
+            IF ( llcycle ) CYCLE
 
             loop_prof : DO ij = 1, inpfiles(jj)%nlev
@@ -527,19 +558,13 @@
                   & CYCLE
 
-               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-
-                  llvalprof = .TRUE. 
-                  EXIT loop_prof
-
-               ENDIF
-
-               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-
-                  llvalprof = .TRUE. 
-                  EXIT loop_prof
-
-               ENDIF
+               DO jvar = 1, kvars
+                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                     & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
+
+                     llvalprof = .TRUE. 
+                     EXIT loop_prof
+
+                  ENDIF
+               END DO
 
             END DO loop_prof
@@ -573,8 +598,8 @@
 
                ! Coordinate search parameters
-               profdata%mi  (iprof,1) = inpfiles(jj)%iobsi(ji,1)
-               profdata%mj  (iprof,1) = inpfiles(jj)%iobsj(ji,1)
-               profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2)
-               profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2)
+               DO jvar = 1, kvars
+                  profdata%mi  (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar)
+                  profdata%mj  (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar)
+               END DO
 
                ! Profile WMO number
@@ -616,138 +641,84 @@
                   IF (ldsatt) THEN
 
-                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &    ldvar1 ) .OR. &
-                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &   ldvar2 ) ) THEN
-                        ip3dt = ip3dt + 1
-                     ELSE
-                        CYCLE
+                     DO jvar = 1, kvars
+                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                           &    ldvar(jvar) ) ) THEN
+                           ip3dt = ip3dt + 1
+                           EXIT
+                        ELSE IF ( jvar == kvars ) THEN
+                           CYCLE loop_p
+                        ENDIF
+                     END DO
+
+                  ENDIF
+
+                  DO jvar = 1, kvars
+                  
+                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                       &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                       &    ldvar(jvar) ) .OR. ldsatt ) THEN
+
+                        IF (ldsatt) THEN
+
+                           ivart(jvar) = ip3dt
+
+                        ELSE
+
+                           ivart(jvar) = ivart(jvar) + 1
+
+                        ENDIF
+
+                        ! Depth of jvar observation
+                        profdata%var(jvar)%vdep(ivart(jvar)) = &
+                           &                inpfiles(jj)%pdep(ij,ji)
+
+                        ! Depth of jvar observation QC
+                        profdata%var(jvar)%idqc(ivart(jvar)) = &
+                           &                inpfiles(jj)%idqc(ij,ji)
+
+                        ! Depth of jvar observation QC flags
+                        profdata%var(jvar)%idqcf(:,ivart(jvar)) = &
+                           &                inpfiles(jj)%idqcf(:,ij,ji)
+
+                        ! Profile index
+                        profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof
+
+                        ! Vertical index in original profile
+                        profdata%var(jvar)%nvlidx(ivart(jvar)) = ij
+
+                        ! Profile jvar value
+                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
+                           profdata%var(jvar)%vobs(ivart(jvar)) = &
+                              &                inpfiles(jj)%pob(ij,ji,jvar)
+                           IF ( ldmod ) THEN
+                              profdata%var(jvar)%vmod(ivart(jvar)) = &
+                                 &                inpfiles(jj)%padd(ij,ji,1,jvar)
+                           ENDIF
+                           ! Count number of profile var1 data as function of type
+                           itypvar( profdata%ntyp(iprof) + 1, jvar ) = &
+                              & itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1
+                        ELSE
+                           profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi
+                        ENDIF
+
+                        ! Profile jvar qc
+                        profdata%var(jvar)%nvqc(ivart(jvar)) = &
+                           & inpfiles(jj)%ivlqc(ij,ji,jvar)
+
+                        ! Profile jvar qc flags
+                        profdata%var(jvar)%nvqcf(:,ivart(jvar)) = &
+                           & inpfiles(jj)%ivlqcf(:,ij,ji,jvar)
+
+                        ! Profile insitu T value
+                        IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN
+                           profdata%var(jvar)%vext(ivart(jvar),1) = &
+                              &                inpfiles(jj)%pext(ij,ji,1)
+                        ENDIF
+
                      ENDIF
-
-                  ENDIF
-
-                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                    &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                    &    ldvar1 ) .OR. ldsatt ) THEN
-
-                     IF (ldsatt) THEN
-
-                        ivar1t = ip3dt
-
-                     ELSE
-
-                        ivar1t = ivar1t + 1
-
-                     ENDIF
-
-                     ! Depth of var1 observation
-                     profdata%var(1)%vdep(ivar1t) = &
-                        &                inpfiles(jj)%pdep(ij,ji)
-
-                     ! Depth of var1 observation QC
-                     profdata%var(1)%idqc(ivar1t) = &
-                        &                inpfiles(jj)%idqc(ij,ji)
-
-                     ! Depth of var1 observation QC flags
-                     profdata%var(1)%idqcf(:,ivar1t) = &
-                        &                inpfiles(jj)%idqcf(:,ij,ji)
-
-                     ! Profile index
-                     profdata%var(1)%nvpidx(ivar1t) = iprof
-
-                     ! Vertical index in original profile
-                     profdata%var(1)%nvlidx(ivar1t) = ij
-
-                     ! Profile var1 value
-                     IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                        & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-                        profdata%var(1)%vobs(ivar1t) = &
-                           &                inpfiles(jj)%pob(ij,ji,1)
-                        IF ( ldmod ) THEN
-                           profdata%var(1)%vmod(ivar1t) = &
-                              &                inpfiles(jj)%padd(ij,ji,1,1)
-                        ENDIF
-                        ! Count number of profile var1 data as function of type
-                        itypvar1( profdata%ntyp(iprof) + 1 ) = &
-                           & itypvar1( profdata%ntyp(iprof) + 1 ) + 1
-                     ELSE
-                        profdata%var(1)%vobs(ivar1t) = fbrmdi
-                     ENDIF
-
-                     ! Profile var1 qc
-                     profdata%var(1)%nvqc(ivar1t) = &
-                        & inpfiles(jj)%ivlqc(ij,ji,1)
-
-                     ! Profile var1 qc flags
-                     profdata%var(1)%nvqcf(:,ivar1t) = &
-                        & inpfiles(jj)%ivlqcf(:,ij,ji,1)
-
-                     ! Profile insitu T value
-                     IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN
-                        profdata%var(1)%vext(ivar1t,1) = &
-                           &                inpfiles(jj)%pext(ij,ji,1)
-                     ENDIF
-
-                  ENDIF
-
-                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                     &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    .AND. &
-                     &   ldvar2 ) .OR. ldsatt ) THEN
-
-                     IF (ldsatt) THEN
-
-                        ivar2t = ip3dt
-
-                     ELSE
-
-                        ivar2t = ivar2t + 1
-
-                     ENDIF
-
-                     ! Depth of var2 observation
-                     profdata%var(2)%vdep(ivar2t) = &
-                        &                inpfiles(jj)%pdep(ij,ji)
-
-                     ! Depth of var2 observation QC
-                     profdata%var(2)%idqc(ivar2t) = &
-                        &                inpfiles(jj)%idqc(ij,ji)
-
-                     ! Depth of var2 observation QC flags
-                     profdata%var(2)%idqcf(:,ivar2t) = &
-                        &                inpfiles(jj)%idqcf(:,ij,ji)
-
-                     ! Profile index
-                     profdata%var(2)%nvpidx(ivar2t) = iprof
-
-                     ! Vertical index in original profile
-                     profdata%var(2)%nvlidx(ivar2t) = ij
-
-                     ! Profile var2 value
-                  IF (  ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. &
-                    &   ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    )  ) THEN
-                        profdata%var(2)%vobs(ivar2t) = &
-                           &                inpfiles(jj)%pob(ij,ji,2)
-                        IF ( ldmod ) THEN
-                           profdata%var(2)%vmod(ivar2t) = &
-                              &                inpfiles(jj)%padd(ij,ji,1,2)
-                        ENDIF
-                        ! Count number of profile var2 data as function of type
-                        itypvar2( profdata%ntyp(iprof) + 1 ) = &
-                           & itypvar2( profdata%ntyp(iprof) + 1 ) + 1
-                     ELSE
-                        profdata%var(2)%vobs(ivar2t) = fbrmdi
-                     ENDIF
-
-                     ! Profile var2 qc
-                     profdata%var(2)%nvqc(ivar2t) = &
-                        & inpfiles(jj)%ivlqc(ij,ji,2)
-
-                     ! Profile var2 qc flags
-                     profdata%var(2)%nvqcf(:,ivar2t) = &
-                        & inpfiles(jj)%ivlqcf(:,ij,ji,2)
-
-                  ENDIF
+                  
+                  END DO
 
                END DO loop_p
@@ -763,10 +734,12 @@
       !-----------------------------------------------------------------------
 
-      CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp )
-      CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp )
+      DO jvar = 1, kvars
+         CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) )
+      END DO
       CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  )
 
-      CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 )
-      CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 )
+      DO jvar = 1, kvars
+         CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 )
+      END DO
 
       !-----------------------------------------------------------------------
@@ -778,50 +751,33 @@
          WRITE(numout,'(1X,A)') '------------'
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) )
-         WRITE(numout,'(1X,A)') '------------------------'
-         DO ji = 0, ntyp1770
-            IF ( itypvar1mpp(ji+1) > 0 ) THEN
-               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
-                  & cwmonam1770(ji)(1:52),' = ', &
-                  & itypvar1mpp(ji+1)
-            ENDIF
+         DO jvar = 1, kvars
+            WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) )
+            WRITE(numout,'(1X,A)') '------------------------'
+            DO ji = 0, ntyp1770
+               IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN
+                  WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
+                     & cwmonam1770(ji)(1:52),' = ', &
+                     & itypvarmpp(ji+1,jvar)
+               ENDIF
+            END DO
+            WRITE(numout,'(1X,A)') &
+               & '---------------------------------------------------------------'
+            WRITE(numout,'(1X,A55,I8)') &
+               & 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// &
+               & '             = ', ivartmpp(jvar)
+            WRITE(numout,'(1X,A)') &
+               & '---------------------------------------------------------------'
+            WRITE(numout,*) 
          END DO
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// &
-            & '             = ', ivar1tmpp
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) )
-         WRITE(numout,'(1X,A)') '------------------------'
-         DO ji = 0, ntyp1770
-            IF ( itypvar2mpp(ji+1) > 0 ) THEN
-               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
-                  & cwmonam1770(ji)(1:52),' = ', &
-                  & itypvar2mpp(ji+1)
-            ENDIF
+      ENDIF
+
+      IF (ldsatt) THEN
+         profdata%nvprot(:)    = ip3dt
+         profdata%nvprotmpp(:) = ip3dtmpp
+      ELSE
+         DO jvar = 1, kvars
+            profdata%nvprot(jvar)    = ivart(jvar)
+            profdata%nvprotmpp(jvar) = ivartmpp(jvar)
          END DO
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// &
-            & '             = ', ivar2tmpp
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,*) 
-      ENDIF
-
-      IF (ldsatt) THEN
-         profdata%nvprot(1)    = ip3dt
-         profdata%nvprot(2)    = ip3dt
-         profdata%nvprotmpp(1) = ip3dtmpp
-         profdata%nvprotmpp(2) = ip3dtmpp
-      ELSE
-         profdata%nvprot(1)    = ivar1t
-         profdata%nvprot(2)    = ivar2t
-         profdata%nvprotmpp(1) = ivar1tmpp
-         profdata%nvprotmpp(2) = ivar2tmpp
       ENDIF
       profdata%nprof        = iprof
@@ -830,14 +786,11 @@
       ! Model level search
       !-----------------------------------------------------------------------
-      IF ( ldvar1 ) THEN
-         CALL obs_level_search( jpk, gdept_1d, &
-            & profdata%nvprot(1), profdata%var(1)%vdep, &
-            & profdata%var(1)%mvk )
-      ENDIF
-      IF ( ldvar2 ) THEN
-         CALL obs_level_search( jpk, gdept_1d, &
-            & profdata%nvprot(2), profdata%var(2)%vdep, &
-            & profdata%var(2)%mvk )
-      ENDIF
+      DO jvar = 1, kvars
+         IF ( ldvar(jvar) ) THEN
+            CALL obs_level_search( jpk, gdept_1d, &
+               & profdata%nvprot(jvar), profdata%var(jvar)%vdep, &
+               & profdata%var(jvar)%mvk )
+         ENDIF
+      END DO
 
       !-----------------------------------------------------------------------
Index: /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 9305)
+++ /branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 9306)
@@ -84,4 +84,7 @@
       CHARACTER(LEN=40) :: clfname
       CHARACTER(LEN=10) :: clfiletype
+      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
+      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
+      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
       INTEGER :: ilevel
       INTEGER :: jvar
@@ -111,5 +114,5 @@
       ! Find maximum level
       ilevel = 0
-      DO jvar = 1, 2
+      DO jvar = 1, profdata%nvar
          ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
       END DO
@@ -175,5 +178,90 @@
          END DO
 
+      CASE('PLCHLTOT')
+
+         clfiletype = 'plchltotfb'
+         cllongname = 'log10(chlorophyll concentration)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('PCHLTOT')
+
+         clfiletype = 'pchltotfb'
+         cllongname = 'chlorophyll concentration'
+         clunits    = 'mg/m3'
+         clgrid     = 'T'
+
+      CASE('PNO3')
+
+         clfiletype = 'pno3fb'
+         cllongname = 'nitrate'
+         clunits    = 'mmol/m3'
+         clgrid     = 'T'
+
+      CASE('PSI4')
+
+         clfiletype = 'psi4fb'
+         cllongname = 'silicate'
+         clunits    = 'mmol/m3'
+         clgrid     = 'T'
+
+      CASE('PPO4')
+
+         clfiletype = 'ppo4fb'
+         cllongname = 'phosphate'
+         clunits    = 'mmol/m3'
+         clgrid     = 'T'
+
+      CASE('PDIC')
+
+         clfiletype = 'pdicfb'
+         cllongname = 'dissolved inorganic carbon'
+         clunits    = 'mmol/m3'
+         clgrid     = 'T'
+
+      CASE('PALK')
+
+         clfiletype = 'palkfb'
+         cllongname = 'alkalinity'
+         clunits    = 'meq/m3'
+         clgrid     = 'T'
+
+      CASE('PPH')
+
+         clfiletype = 'pphfb'
+         cllongname = 'pH'
+         clunits    = '-'
+         clgrid     = 'T'
+
+      CASE('PO2')
+
+         clfiletype = 'po2fb'
+         cllongname = 'dissolved oxygen'
+         clunits    = 'mmol/m3'
+         clgrid     = 'T'
+
       END SELECT
+      
+      IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. &
+         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
+         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
+            &                 1 + iadd, iext, .TRUE. )
+         fbdata%cname(1)      = profdata%cvars(1)
+         fbdata%coblong(1)    = cllongname
+         fbdata%cobunit(1)    = clunits
+         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
+         fbdata%caddunit(1,1) = clunits
+         fbdata%cgrid(:)      = clgrid
+         DO je = 1, iext
+            fbdata%cextname(je) = pext%cdname(je)
+            fbdata%cextlong(je) = pext%cdlong(je,1)
+            fbdata%cextunit(je) = pext%cdunit(je,1)
+         END DO
+         DO ja = 1, iadd
+            fbdata%caddname(1+ja) = padd%cdname(ja)
+            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
+            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
+         END DO
+      ENDIF
 
       fbdata%caddname(1)   = 'Hx'
@@ -228,5 +316,5 @@
             &           krefdate = 19500101 )
          ! Reform the profiles arrays for output
-         DO jvar = 1, 2
+         DO jvar = 1, profdata%nvar
             DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
                ik = profdata%var(jvar)%nvlidx(jk)
@@ -322,4 +410,7 @@
       CHARACTER(LEN=10) :: clfiletype
       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
+      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
+      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
+      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
       INTEGER :: jo
       INTEGER :: ja
@@ -344,4 +435,7 @@
       SELECT CASE ( TRIM(surfdata%cvars(1)) )
       CASE('SLA')
+
+         ! SLA needs special treatment because of MDT, so is all done here
+         ! Other variables are done more generically
 
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
@@ -374,11 +468,138 @@
       CASE('SST')
 
+         clfiletype = 'sstfb'
+         cllongname = 'Sea surface temperature'
+         clunits    = 'Degree centigrade'
+         clgrid     = 'T'
+
+      CASE('ICECONC')
+
+         clfiletype = 'sicfb'
+         cllongname = 'Sea ice'
+         clunits    = 'Fraction'
+         clgrid     = 'T'
+
+      CASE('SSS')
+
+         clfiletype = 'sssfb'
+         cllongname = 'Sea surface salinity'
+         clunits    = 'psu'
+         clgrid     = 'T'
+
+      CASE('SLCHLTOT','LOGCHL','LogChl','logchl')
+
+         clfiletype = 'slchltotfb'
+         cllongname = 'Surface total log10(chlorophyll)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('SLCHLDIA')
+
+         clfiletype = 'slchldiafb'
+         cllongname = 'Surface diatom log10(chlorophyll)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('SLCHLNON')
+
+         clfiletype = 'slchlnonfb'
+         cllongname = 'Surface non-diatom log10(chlorophyll)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('SLCHLDIN')
+
+         clfiletype = 'slchldinfb'
+         cllongname = 'Surface dinoflagellate log10(chlorophyll)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('SLCHLMIC')
+
+         clfiletype = 'slchlmicfb'
+         cllongname = 'Surface microphytoplankton log10(chlorophyll)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('SLCHLNAN')
+
+         clfiletype = 'slchlnanfb'
+         cllongname = 'Surface nanophytoplankton log10(chlorophyll)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('SLCHLPIC')
+
+         clfiletype = 'slchlpicfb'
+         cllongname = 'Surface picophytoplankton log10(chlorophyll)'
+         clunits    = 'log10(mg/m3)'
+         clgrid     = 'T'
+
+      CASE('SCHLTOT')
+
+         clfiletype = 'schltotfb'
+         cllongname = 'Surface total chlorophyll'
+         clunits    = 'mg/m3'
+         clgrid     = 'T'
+
+      CASE('SLPHYTOT')
+
+         clfiletype = 'slphytotfb'
+         cllongname = 'Surface total log10(phytoplankton carbon)'
+         clunits    = 'log10(mmolC/m3)'
+         clgrid     = 'T'
+
+      CASE('SLPHYDIA')
+
+         clfiletype = 'slphydiafb'
+         cllongname = 'Surface diatom log10(phytoplankton carbon)'
+         clunits    = 'log10(mmolC/m3)'
+         clgrid     = 'T'
+
+      CASE('SLPHYNON')
+
+         clfiletype = 'slphynonfb'
+         cllongname = 'Surface non-diatom log10(phytoplankton carbon)'
+         clunits    = 'log10(mmolC/m3)'
+         clgrid     = 'T'
+
+      CASE('SSPM')
+
+         clfiletype = 'sspmfb'
+         cllongname = 'Surface suspended particulate matter'
+         clunits    = 'g/m3'
+         clgrid     = 'T'
+
+      CASE('SFCO2','FCO2','fCO2','fco2')
+
+         clfiletype = 'sfco2fb'
+         cllongname = 'Surface fugacity of carbon dioxide'
+         clunits    = 'uatm'
+         clgrid     = 'T'
+
+      CASE('SPCO2')
+
+         clfiletype = 'spco2fb'
+         cllongname = 'Surface partial pressure of carbon dioxide'
+         clunits    = 'uatm'
+         clgrid     = 'T'
+
+      CASE DEFAULT
+
+         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
+
+      END SELECT
+
+      ! SLA needs special treatment because of MDT, so is done above
+      ! Remaining variables treated more generically
+
+      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
+      
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
             &                 1 + iadd, iext, .TRUE. )
 
-         clfiletype = 'sstfb'
          fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'Sea surface temperature'
-         fbdata%cobunit(1)    = 'Degree centigrade'
+         fbdata%coblong(1)    = cllongname
+         fbdata%cobunit(1)    = clunits
          DO je = 1, iext
             fbdata%cextname(je) = pext%cdname(je)
@@ -386,7 +607,11 @@
             fbdata%cextunit(je) = pext%cdunit(je,1)
          END DO
-         fbdata%caddlong(1,1) = 'Model interpolated SST'
-         fbdata%caddunit(1,1) = 'Degree centigrade'
-         fbdata%cgrid(1)      = 'T'
+         IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
+            fbdata%caddlong(1,1) = 'Model interpolated ICE'
+         ELSE
+            fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1))
+         ENDIF
+         fbdata%caddunit(1,1) = clunits
+         fbdata%cgrid(1)      = clgrid
          DO ja = 1, iadd
             fbdata%caddname(1+ja) = padd%cdname(ja)
@@ -395,148 +620,6 @@
          END DO
 
-      CASE('ICECONC')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'sicfb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'Sea ice'
-         fbdata%cobunit(1)    = 'Fraction'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated ICE'
-         fbdata%caddunit(1,1) = 'Fraction'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE('SSS')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'sssfb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'Sea surface salinity'
-         fbdata%cobunit(1)    = 'psu'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated SSS'
-         fbdata%caddunit(1,1) = 'psu'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE('LOGCHL','LogChl','logchl')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'logchlfb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'logchl concentration'
-         fbdata%cobunit(1)    = 'mg/m3'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated LOGCHL'
-         fbdata%caddunit(1,1) = 'mg/m3'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE('SPM','Spm','spm')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'spmfb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'spm'
-         fbdata%cobunit(1)    = 'g/m3'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated spm'
-         fbdata%caddunit(1,1) = 'g/m3'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE('FCO2','fCO2','fco2')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'fco2fb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'fco2'
-         fbdata%cobunit(1)    = 'uatm'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated fco2'
-         fbdata%caddunit(1,1) = 'uatm'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE('PCO2','pCO2','pco2')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'pco2fb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'pco2'
-         fbdata%cobunit(1)    = 'uatm'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated pco2'
-         fbdata%caddunit(1,1) = 'uatm'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE DEFAULT
-
-         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
-
-      END SELECT
-
+      ENDIF
+      
       fbdata%caddname(1)   = 'Hx'
 
