Index: /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9185)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9186)
@@ -1212,4 +1212,13 @@
    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
@@ -1226,34 +1235,54 @@
    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_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
+   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_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_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 9185)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 9186)
@@ -164,4 +164,13 @@
          & 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
 
@@ -188,4 +197,13 @@
       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
@@ -227,4 +245,7 @@
          &            ln_slchlpic, ln_schltot,                        &
          &            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,          &
@@ -242,4 +263,8 @@
          &            cn_schltotfbfiles, 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,            &
@@ -285,4 +310,13 @@
       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
@@ -335,4 +369,13 @@
          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
@@ -375,9 +418,11 @@
       !-----------------------------------------------------------------------
 
-      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
+      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_sspm, ln_sfco2, ln_spco2 /) )
+         &                  ln_sspm,     ln_sfco2,    ln_spco2 /) )
 
       IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
@@ -400,11 +445,56 @@
          IF (ln_t3d .OR. ln_s3d) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', &
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof', &
                &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
          ENDIF
          IF (ln_vel3d) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', &
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel', &
                &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_plchltot) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'plchltot', &
+               &                   cn_plchltotfbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_pchltot) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pchltot', &
+               &                   cn_pchltotfbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_pno3) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pno3', &
+               &                   cn_pno3fbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_psi4) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'psi4', &
+               &                   cn_psi4fbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_ppo4) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'ppo4', &
+               &                   cn_ppo4fbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_pdic) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pdic', &
+               &                   cn_pdicfbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_palk) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'palk', &
+               &                   cn_palkfbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_pph) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pph', &
+               &                   cn_pphfbfiles, ifilesprof, cobstypesprof, clproffiles )
+         ENDIF
+         IF (ln_po2) THEN
+            jtype = jtype + 1
+            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'po2', &
+               &                   cn_po2fbfiles, ifilesprof, cobstypesprof, clproffiles )
          ENDIF
 
@@ -615,6 +705,6 @@
          DO jtype = 1, nproftypes
 
-            nvarsprof(jtype) = 2
             IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
+               nvarsprof(jtype) = 2
                nextrprof(jtype) = 1
                llvar1 = ln_t3d
@@ -626,6 +716,6 @@
                zgphi2 = gphit
                zmask2 = tmask
-            ENDIF
-            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
+            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
+               nvarsprof(jtype) = 2
                nextrprof(jtype) = 2
                llvar1 = ln_vel3d
@@ -637,4 +727,15 @@
                zgphi2 = gphiv
                zmask2 = vmask
+            ELSE
+               nvarsprof(jtype) = 1
+               nextrprof(jtype) = 0
+               llvar1 = .TRUE.
+               llvar2 = .FALSE.
+               zglam1 = glamt
+               zgphi1 = gphit
+               zmask1 = tmask
+               zglam2 = glamt
+               zgphi2 = gphit
+               zmask2 = tmask
             ENDIF
 
@@ -763,19 +864,20 @@
 #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
+      USE par_hadocc
 #elif defined key_medusa && defined key_foam_medusa
-      USE trc, ONLY :  &           ! MEDUSA chlorophyll, fCO2 and pCO2
+      USE trc, ONLY :  &           ! MEDUSA variables
          & trn
-      USE par_medusa, ONLY: &
-         & jpchn, &
-         & jpchd
+      USE par_medusa
 #if defined key_roam
       USE sms_medusa, ONLY: &
          & f2_pco2w, &
-         & f2_fco2w
+         & f2_fco2w, &
+         & f3_pH
 #endif
 #elif defined key_fabm
@@ -864,7 +966,187 @@
                zgphi1(:,:) = gphiu(:,:)
                zgphi2(:,:) = gphiv(:,:)
+
+            CASE('plchltot')
+#if defined key_hadocc
+               ! Chlorophyll from HadOCC
+               zprofvar1(:,:,:) = HADOCC_CHL(:,:,:)
+#elif defined key_medusa && defined key_foam_medusa
+               ! Add non-diatom and diatom chlorophyll from MEDUSA
+               zprofvar1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
+#elif defined key_fabm
+               ! Add all chlorophyll groups from ERSEM
+               zprofvar1(:,:,:) = 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
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               ! Take the log10 where we can, otherwise exclude
+               tiny = 1.0e-20
+               WHERE(zprofvar1(:,:,:) > tiny .AND. zprofvar1(:,:,:) /= obfillflt )
+                  zprofvar1(:,:,:)  = LOG10(zprofvar1(:,:,:))
+               ELSEWHERE
+                  zprofvar1(:,:,:)  = obfillflt
+                  zprofmask1(:,:,:) = 0
+               END WHERE
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            CASE('pchltot')
+#if defined key_hadocc
+               ! Chlorophyll from HadOCC
+               zprofvar1(:,:,:) = HADOCC_CHL(:,:,:)
+#elif defined key_medusa && defined key_foam_medusa
+               ! Add non-diatom and diatom chlorophyll from MEDUSA
+               zprofvar1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
+#elif defined key_fabm
+               ! Add all chlorophyll groups from ERSEM
+               zprofvar1(:,:,:) = 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
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            CASE('pno3')
+#if defined key_hadocc
+               ! Dissolved inorganic nitrogen from HadOCC
+               zprofvar1(:,:,:) = trn(:,:,:,jp_had_nut)
+#elif defined key_medusa && defined key_foam_medusa
+               ! Dissolved inorganic nitrogen from MEDUSA
+               zprofvar1(:,:,:) = trn(:,:,:,jpdin)
+#elif defined key_fabm
+               ! Nitrate from ERSEM
+               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n3n)
+#else
+               CALL ctl_stop( ' Trying to run pno3 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            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
+               zprofvar1(:,:,:) = trn(:,:,:,jpsil)
+#elif defined key_fabm
+               ! Silicate from ERSEM
+               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n5s)
+#else
+               CALL ctl_stop( ' Trying to run psi4 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            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
+               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n1p)
+#else
+               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            CASE('pdic')
+#if defined key_hadocc
+               ! Dissolved inorganic carbon from HadOCC
+               zprofvar1(:,:,:) = trn(:,:,:,jp_had_dic)
+#elif defined key_medusa && defined key_foam_medusa
+               ! Dissolved inorganic carbon from MEDUSA
+               zprofvar1(:,:,:) = trn(:,:,:,jpdic)
+#elif defined key_fabm
+               ! Dissolved inorganic carbon from ERSEM
+               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3c)
+#else
+               CALL ctl_stop( ' Trying to run pdic observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            CASE('palk')
+#if defined key_hadocc
+               ! Alkalinity from HadOCC
+               zprofvar1(:,:,:) = trn(:,:,:,jp_had_alk)
+#elif defined key_medusa && defined key_foam_medusa
+               ! Alkalinity from MEDUSA
+               zprofvar1(:,:,:) = trn(:,:,:,jpalk)
+#elif defined key_fabm
+               ! Alkalinity from ERSEM
+               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3a)
+#else
+               CALL ctl_stop( ' Trying to run palk observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            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
+               zprofvar1(:,:,:) = f3_pH(:,:,:)
+#elif defined key_fabm
+               ! pH from ERSEM
+               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3ph)
+#else
+               CALL ctl_stop( ' Trying to run pph observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
+            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
+               zprofvar1(:,:,:) = trn(:,:,:,jpoxy)
+#elif defined key_fabm
+               ! Oxygen from ERSEM
+               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o2o)
+#else
+               CALL ctl_stop( ' Trying to run po2 observation operator', &
+                  &           ' but no biogeochemical model appears to have been defined' )
+#endif
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+
             CASE DEFAULT
                CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
             END SELECT
+            
+            IF ( ( TRIM(cobstypesprof(jtype)) /= 'prof' ) .AND. ( TRIM(cobstypesprof(jtype)) /= 'vel' ) ) THEN
+               zprofvar2(:,:,:)  = zprofvar1(:,:,:)
+               zprofmask2(:,:,:) = zprofmask1(:,:,:)
+               zglam2(:,:)       = zglam1(:,:)
+               zgphi2(:,:)       = zgphi1(:,:)
+            ENDIF
 
             CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
Index: /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 9185)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 9186)
@@ -231,5 +231,7 @@
                   DO ji = 1, jpi
                      prodatqc%vdmean(ji,jj,jk,1) = 0.0
-                     prodatqc%vdmean(ji,jj,jk,2) = 0.0
+                     IF ( prodatqc%nvar == 2 ) THEN
+                        prodatqc%vdmean(ji,jj,jk,2) = 0.0
+                     ENDIF
                   END DO
                END DO
@@ -243,7 +245,9 @@
                   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)
+                  IF ( prodatqc%nvar == 2 ) THEN
+                     ! Increment field 2 for computing daily mean
+                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
+                        &                        + pvar2(ji,jj,jk)
+                  ENDIF
                END DO
             END DO
@@ -260,6 +264,8 @@
                      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
+                     IF ( prodatqc%nvar == 2 ) THEN
+                        prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
+                           &                        * zdaystp
+                     ENDIF
                   END DO
                END DO
@@ -297,12 +303,14 @@
          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)
+         IF ( prodatqc%nvar == 2 ) THEN
+            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)
+         ENDIF
       END DO
 
@@ -316,8 +324,10 @@
       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 )
+      IF ( prodatqc%nvar == 2 ) THEN
+         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 )
+      ENDIF
 
       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 
@@ -334,6 +344,8 @@
          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 )
+         IF ( prodatqc%nvar == 2 ) THEN
+            CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
+               &                  prodatqc%vdmean(:,:,:,2), zinm2 )
+         ENDIF
 
       ENDIF
@@ -378,10 +390,12 @@
          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 )
- 
+         IF ( prodatqc%nvar == 2 ) THEN
+            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
          ENDIF
 
@@ -512,13 +526,51 @@
          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
+         IF ( prodatqc%nvar == 2 ) THEN
+            ! 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 
@@ -526,115 +578,79 @@
                   iend = prodatqc%npvend(jobs,2) 
                   inum_obs = iend - ista + 1 
-                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
-
+                  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 
 
-                        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)) 
-       
+                        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 
 
-               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 
+                  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 
-
-            !------------------------------------------------------------- 
-            ! 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 
 
Index: /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 9185)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 9186)
@@ -403,6 +403,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                         )
+      IF ( ld_var2 ) THEN
+         CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), &
+            &              profdata%nqc,     igrdobs                         )
+      ENDIF
 
       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
@@ -441,25 +443,27 @@
       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 )
+      IF ( ld_var2 ) THEN
+         ! 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 )
+      ENDIF
 
       ! -----------------------------------------------------------------------
@@ -535,30 +539,37 @@
          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)
+         IF ( ld_var2 ) THEN
+            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)
+         ENDIF
 
          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)//'     '
+         IF ( ld_var2 ) THEN
+            WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', &
+               &                               '   '//prodatqc%cvars(1)//'     ', &
+               &                               '   '//prodatqc%cvars(2)//'     '
+         ELSE
+            WRITE(numout,'(10X,A,5X,A,5X,A)')'Time step','Profiles', &
+               &                               '   '//prodatqc%cvars(1)//'     '
+         ENDIF
          WRITE(numout,998)
       ENDIF
@@ -588,7 +599,12 @@
          DO jstp = nit000 - 1, nitend
             inrc = jstp - nit000 + 2
-            WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
-               &                    prodatqc%nvstpmpp(inrc,1), &
-               &                    prodatqc%nvstpmpp(inrc,2)
+            IF ( ld_var2 ) THEN
+               WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
+                  &                    prodatqc%nvstpmpp(inrc,1), &
+                  &                    prodatqc%nvstpmpp(inrc,2)
+            ELSE
+               WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
+                  &                    prodatqc%nvstpmpp(inrc,1)
+            ENDIF
          END DO
       ENDIF
Index: /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 9185)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 9186)
@@ -219,7 +219,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
 
@@ -320,12 +320,18 @@
             ALLOCATE( iobsj1(inowin) )
             ALLOCATE( iproc1(inowin) )
-            ALLOCATE( iobsi2(inowin) )
-            ALLOCATE( iobsj2(inowin) )
-            ALLOCATE( iproc2(inowin) )
+            IF ( kvars == 2 ) THEN
+               ALLOCATE( iobsi2(inowin) )
+               ALLOCATE( iobsj2(inowin) )
+               ALLOCATE( iproc2(inowin) )
+            ENDIF
             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
+               IF ( kvars == 2 ) THEN
+                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                     & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               ELSE
+                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE
+               ENDIF
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -347,4 +353,7 @@
                CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, &
                   &                  iproc2, 'V' )
+            ELSE
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, &
+                  &                  iproc1, 'T' )
             ENDIF
 
@@ -352,6 +361,10 @@
             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
+               IF ( kvars == 2 ) THEN
+                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                     & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               ELSE
+                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE
+               ENDIF
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -360,20 +373,29 @@
                   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')
+                  IF ( kvars == 2 ) THEN
+                     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')
+                     ENDIF
                   ENDIF
                ENDIF
             END DO
-            DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 )
+            DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1 )
+            IF ( kvars == 2 ) THEN
+               DEALLOCATE( iobsi2, iobsj2, iproc2 )
+            ENDIF
 
             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
+               IF ( kvars == 2 ) THEN
+                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                     & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               ELSE
+                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE
+               ENDIF
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -407,12 +429,21 @@
                      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.
+                     IF ( kvars == 2 ) 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
+                           llvalprof = .TRUE.
+                        ENDIF
+                     ELSE
+                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                           &    ldvar1 ) ) THEN
+                           ip3dt = ip3dt + 1
+                           llvalprof = .TRUE.
+                        ENDIF
                      ENDIF
                   END DO loop_p_count
@@ -438,6 +469,10 @@
          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
+            IF ( kvars == 2 ) THEN
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            ELSE
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE
+            ENDIF
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -453,6 +488,10 @@
          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
+            IF ( kvars == 2 ) THEN
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            ELSE
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE
+            ENDIF
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -501,7 +540,11 @@
          ji = iprofidx(iindx(jk))
 
-            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+         IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+         IF ( kvars == 2 ) THEN
             IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
                & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+         ELSE
+            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE
+         ENDIF
 
          IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
@@ -519,6 +562,10 @@
 
             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 ( kvars == 2 ) THEN
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
+                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            ELSE
+               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE
+            ENDIF
 
             loop_prof : DO ij = 1, inpfiles(jj)%nlev
@@ -535,10 +582,12 @@
                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
-
+               IF ( kvars == 2 ) THEN
+                  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
                ENDIF
 
@@ -575,6 +624,8 @@
                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)
+               IF ( kvars == 2 ) THEN
+                  profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2)
+                  profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2)
+               ENDIF
 
                ! Profile WMO number
@@ -616,13 +667,23 @@
                   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
+                     IF ( kvars == 2 ) 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
+                        ENDIF
                      ELSE
-                        CYCLE
+                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
+                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                           &    ldvar1 ) ) THEN
+                           ip3dt = ip3dt + 1
+                        ELSE
+                           CYCLE
+                        ENDIF
                      ENDIF
 
@@ -693,60 +754,62 @@
                   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
+                  IF ( kvars == 2 ) THEN
+                     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
-
-                     ! 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
 
@@ -764,9 +827,13 @@
 
       CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp )
-      CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp )
+      IF ( kvars == 2 ) THEN
+         CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp )
+      ENDIF
       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 )
+      IF ( kvars == 2 ) THEN
+         CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 )
+      ENDIF
 
       !-----------------------------------------------------------------------
@@ -795,33 +862,39 @@
             & '---------------------------------------------------------------'
          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
-         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,*) 
+         IF ( kvars == 2 ) THEN
+            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
+            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
       ENDIF
 
       IF (ldsatt) THEN
          profdata%nvprot(1)    = ip3dt
-         profdata%nvprot(2)    = ip3dt
          profdata%nvprotmpp(1) = ip3dtmpp
-         profdata%nvprotmpp(2) = ip3dtmpp
+         IF ( kvars == 2 ) THEN
+            profdata%nvprot(2)    = ip3dt
+            profdata%nvprotmpp(2) = ip3dtmpp
+         ENDIF
       ELSE
          profdata%nvprot(1)    = ivar1t
-         profdata%nvprot(2)    = ivar2t
          profdata%nvprotmpp(1) = ivar1tmpp
-         profdata%nvprotmpp(2) = ivar2tmpp
+         IF ( kvars == 2 ) THEN
+            profdata%nvprot(2)    = ivar2t
+            profdata%nvprotmpp(2) = ivar2tmpp
+         ENDIF
       ENDIF
       profdata%nprof        = iprof
Index: /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 9185)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 9186)
@@ -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)
