Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 15764)
@@ -21,5 +21,6 @@
    !!   dyn_asm_inc    : Apply the dynamic (u and v) increments
    !!   ssh_asm_inc    : Apply the SSH increment
-   !!   seaice_asm_inc : Apply the seaice increment
+   !!   seaice_asm_inc : Apply the sea ice concentration increment
+   !!   sit_asm_inc    : Apply the sea ice thickness increment
    !!----------------------------------------------------------------------
    USE wrk_nemo         ! Memory Allocation
@@ -49,5 +50,7 @@
    PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
    PUBLIC   ssh_asm_inc    !: Apply the SSH increment
-   PUBLIC   seaice_asm_inc !: Apply the seaice increment
+   PUBLIC   seaice_asm_inc !: Apply the seaice concentration increment
+   PUBLIC   sit_asm_inc    !: Apply the seaice thickness increment
+   PUBLIC   bgc_asm_inc    !: Apply the biogeochemistry increments
 
 #if defined key_asminc
@@ -57,4 +60,5 @@
 #endif
    LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
+   LOGICAL, PUBLIC :: ln_avgbkg = .FALSE.      !: No output of the mean background state fields
    LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
    LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
@@ -62,5 +66,7 @@
    LOGICAL, PUBLIC :: ln_dyninc = .FALSE.      !: No dynamics (u and v) assimilation increments
    LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
-   LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment
+   LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE.   !: No sea ice concentration increment
+   LOGICAL, PUBLIC :: ln_sitinc = .FALSE.      !: No sea ice thickness increment
+   LOGICAL, PUBLIC :: lk_bgcinc = .FALSE.      !: No biogeochemistry increments
    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -85,4 +91,5 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   sit_bkginc            ! Increment to the background sea ice thickness
 
    !! * Substitutions
@@ -127,5 +134,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace
       !!
-      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
+      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg, ln_balwri,                &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
          &                 ln_asmdin, ln_asmiau,                           &
@@ -138,4 +145,5 @@
       !-----------------------------------------------------------------------
       ln_seaiceinc = .FALSE.
+      ln_sitinc = .FALSE.
       ln_temnofreeze = .FALSE.
 
@@ -154,11 +162,14 @@
          WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
+         WRITE(numout,*) '   Namelist nam_asminc : set assimilation increment parameters'
          WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
+         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg
+         WRITE(numout,*) '      Logical switch for writing out balancing increments      ln_balwri = ', ln_balwri
          WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
          WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
          WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
+         WRITE(numout,*) '      Logical switch for applying SIC increments               ln_seaiceinc = ', ln_seaiceinc
+         WRITE(numout,*) '      Logical switch for applying SIT increments               ln_sitinc = ', ln_sitinc
          WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
-         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
          WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
          WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
@@ -216,6 +227,8 @@
 
       IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
-           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
-         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
+         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. &
+         &        ( ln_sitinc ).OR.( lk_bgcinc ) )) &
+         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_sitinc and ln_(bgc-variable)inc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -226,6 +239,7 @@
 
       IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
-         &                     )  &
-         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
+         & .AND.( .NOT. ln_sitinc ).AND.( .NOT. lk_bgcinc ) )  &
+         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_sitinc and ln_(bgc-variable)inc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -325,23 +339,34 @@
       !--------------------------------------------------------------------
 
-      ALLOCATE( t_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( s_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( u_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( v_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( ssh_bkginc(jpi,jpj)   )
-      ALLOCATE( seaice_bkginc(jpi,jpj))
+      IF ( ln_trainc ) THEN
+         ALLOCATE( t_bkginc(jpi,jpj,jpk) )
+         ALLOCATE( s_bkginc(jpi,jpj,jpk) )
+         t_bkginc(:,:,:) = 0.0
+         s_bkginc(:,:,:) = 0.0
+      ENDIF
+      IF ( ln_dyninc ) THEN 
+         ALLOCATE( u_bkginc(jpi,jpj,jpk) )
+         ALLOCATE( v_bkginc(jpi,jpj,jpk) )
+         u_bkginc(:,:,:) = 0.0
+         v_bkginc(:,:,:) = 0.0
+      ENDIF
+      IF ( ln_sshinc ) THEN
+         ALLOCATE( ssh_bkginc(jpi,jpj)   )
+         ssh_bkginc(:,:) = 0.0
+      ENDIF
+      IF ( ln_seaiceinc ) THEN 
+         ALLOCATE( seaice_bkginc(jpi,jpj))
+         seaice_bkginc(:,:) = 0.0
+      ENDIF
+      IF ( ln_sitinc ) THEN 
+         ALLOCATE( sit_bkginc(jpi,jpj))
+         sit_bkginc(:,:) = 0.0
+      ENDIF
 #if defined key_asminc
       ALLOCATE( ssh_iau(jpi,jpj)      )
-#endif
-      t_bkginc(:,:,:) = 0.0
-      s_bkginc(:,:,:) = 0.0
-      u_bkginc(:,:,:) = 0.0
-      v_bkginc(:,:,:) = 0.0
-      ssh_bkginc(:,:) = 0.0
-      seaice_bkginc(:,:) = 0.0
-#if defined key_asminc
       ssh_iau(:,:)    = 0.0
 #endif
-      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
+      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
+         &  .OR.( ln_sitinc ).OR.( lk_bgcinc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -406,4 +431,13 @@
             ! to allow for differences in masks
             WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
+         ENDIF 
+
+         IF ( ln_sitinc ) THEN 
+            CALL iom_get( inum, jpdom_autoglo, 'bckinsit', sit_bkginc, 1 ) 
+            ! Apply the masks 
+            sit_bkginc(:,:) = sit_bkginc(:,:) * tmask(:,:,1) 
+            ! Set missing increments to 0.0 rather than 1e+20 
+            ! to allow for differences in masks 
+            WHERE( ABS( sit_bkginc(:,:) ) > 1.0e+10 ) sit_bkginc(:,:) = 0.0 
          ENDIF
 
@@ -766,4 +800,5 @@
       ! Perhaps the following call should be in step
       IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
+      IF   ( ln_sitinc  )      CALL sit_asm_inc ( kt )      ! apply sea ice thickness increment
       !
    END SUBROUTINE tra_asm_inc
@@ -923,4 +958,92 @@
    END SUBROUTINE ssh_asm_inc
 
+   SUBROUTINE sit_asm_inc( kt, kindic )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE sit_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the sea ice thickness assimilation increment.
+      !!
+      !! ** Method  : Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      !
+      INTEGER, INTENT(in)           ::   kt   ! Current time step
+      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
+      !
+      INTEGER  ::   it
+      REAL(wp) ::   zincwgt   ! IAU weight for current time step
+      !!----------------------------------------------------------------------
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
+            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
+            ! EF: Actually CICE is expecting a tendency so is divided by rdt below
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'sit_asm_inc : sea ice thick IAU at time step = ', &
+                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE
+            ndsit_da(:,:) = sit_bkginc(:,:) * zincwgt / rdt
+#endif
+
+            IF ( kt == nitiaufin_r ) THEN
+               DEALLOCATE( sit_bkginc )
+            ENDIF
+
+         ELSE
+
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thickness : CICE case. Zero ice increment tendency into CICE
+            ndsit_da(:,:) = 0.0_wp
+#endif
+
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE
+           ndsit_da(:,:) = sit_bkginc(:,:) / rdt
+#endif
+           IF ( .NOT. PRESENT(kindic) ) THEN
+              DEALLOCATE( sit_bkginc )
+           END IF
+
+         ELSE
+
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thicnkness : CICE case. Zero ice thickness increment tendency into CICE 
+            ndsit_da(:,:) = 0.0_wp
+#endif
+         
+         ENDIF
+
+      ENDIF
+
+   END SUBROUTINE sit_asm_inc
 
    SUBROUTINE seaice_asm_inc( kt, kindic )
@@ -928,5 +1051,5 @@
       !!                    ***  ROUTINE seaice_asm_inc  ***
       !!          
-      !! ** Purpose : Apply the sea ice assimilation increment.
+      !! ** Purpose : Apply the sea ice concentration assimilation increment.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 15764)
@@ -45,7 +45,4 @@
    USE agrif_opa_interp ! agrif
 #endif
-#if defined key_asminc   
-   USE asminc          ! Assimilation increment
-#endif
 
    IMPLICIT NONE
@@ -459,12 +456,6 @@
                 &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
       ENDIF
-#if defined key_asminc
-      !                                         ! Include the IAU weighted SSH increment
-      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
-         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
-      ENDIF
-#endif
-      !                                   !* Fill boundary data arrays for AGRIF
-      !                                   ! ------------------------------------
+      !                                   !* Fill boundary data arrays with AGRIF
+      !                                   ! -------------------------------------
 #if defined key_agrif
          IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 15764)
@@ -22,4 +22,5 @@
    USE obs_read_surf            ! Reading and allocation of surface obs
    USE obs_readmdt              ! Reading and allocation of MDT for SLA.
+   USE obs_readsnowdepth        ! Get model snow depth for conversion of freeboard to ice thickness
    USE obs_prep                 ! Preparation of obs. (grid search etc).
    USE obs_oper                 ! Observation operators
@@ -33,4 +34,5 @@
    USE mpp_map                  ! MPP mapping
    USE lib_mpp                  ! For ctl_warn/stop
+   USE tradmp                   ! For climatological temperature & salinity
 
    IMPLICIT NONE
@@ -44,32 +46,49 @@
 
    !! * 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_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)
- 
+   LOGICAL, PUBLIC :: & 
+      &       lk_diaobs = .TRUE.   !: Include this for backwards compatibility at NEMO 3.6. 
+   LOGICAL :: ln_diaobs            !: Logical switch for the obs operator 
+   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs 
+   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 
+   LOGICAL :: ln_sla_fp_indegs     !: T=> SLA obs footprint size specified in degrees, F=> in metres 
+   LOGICAL :: ln_sst_fp_indegs     !: T=> SST obs footprint size specified in degrees, F=> in metres 
+   LOGICAL :: ln_sss_fp_indegs     !: T=> SSS obs footprint size specified in degrees, F=> in metres 
+   LOGICAL :: ln_ssv_fp_indegs     !: T=> SSV obs footprint size specified in degrees, F=> in metres    
+   LOGICAL :: ln_sic_fp_indegs     !: T=> SIC obs footprint size specified in degrees, F=> in metres 
+   LOGICAL :: ln_sit_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres 
+   LOGICAL :: ln_fbd_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres 
+   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and writing T/S climatology 
+   LOGICAL :: ln_time_mean_sla_bkg !: Logical switch for applying time mean of SLA background to remove tidal signal 
+
+   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_ssv_avglamscl     !: E/W diameter of SSV observation footprint 
+   REAL(wp) :: rn_ssv_avgphiscl     !: N/S diameter of SSV observation footprint 
+   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of SIC observation footprint 
+   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of SIC observation footprint 
+   REAL(wp) :: rn_sit_avglamscl     !: E/W diameter of SIT observation footprint 
+   REAL(wp) :: rn_sit_avgphiscl     !: N/S diameter of SIT observation footprint 
+   REAL(wp) :: rn_fbd_avglamscl     !: E/W diameter of FBD observation footprint 
+   REAL(wp) :: rn_fbd_avgphiscl     !: N/S diameter of FBD observation footprint 
+   REAL(wp), PUBLIC :: & 
+      &        MeanPeriodHours = 24. + (5./6.) !: Meaning period for surface data. 
+
+
+   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_ssv     !: SSV horizontal interpolation method (-1 = default)    
+   INTEGER :: nn_2dint_sic     !: SIC horizontal interpolation method (-1 = default) 
+   INTEGER :: nn_2dint_sit     !: SIT horizontal interpolation method (-1 = default) 
+   INTEGER :: nn_2dint_fbd     !: FBD horizontal interpolation method (-1 = default)
+
    INTEGER, DIMENSION(imaxavtypes) :: &
       & nn_profdavtypes      !: Profile data types representing a daily average
@@ -130,157 +149,186 @@
       !!        !  15-02  (M. Martin) Simplification of namelist and code
       !!----------------------------------------------------------------------
+#if defined key_cice 
+USE sbc_oce, ONLY : &        ! CICE variables 
+   & thick_s                 ! snow depth for freeboard conversion  
+#endif
 
       IMPLICIT NONE
 
       !! * Local declarations
-      INTEGER, PARAMETER :: &
-         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
-      INTEGER, DIMENSION(:), ALLOCATABLE :: &
-         & ifilesprof, &         ! Number of profile files
-         & ifilessurf            ! Number of surface files
-      INTEGER :: ios             ! Local integer output status for namelist read
-      INTEGER :: jtype           ! Counter for obs types
-      INTEGER :: jvar            ! Counter for variables
-      INTEGER :: jfile           ! Counter for files
-      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply
-      INTEGER :: n2dint_type     ! Local version of nn_2dint*
-
-      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
-         & cn_profbfiles,      & ! T/S profile input filenames
-         & cn_sstfbfiles,      & ! Sea surface temperature input filenames
-         & cn_slafbfiles,      & ! Sea level anomaly input filenames
-         & cn_sicfbfiles,      & ! Seaice concentration input filenames
-         & cn_velfbfiles,      & ! Velocity profile input filenames
-         & cn_sssfbfiles,      & ! Sea surface salinity input filenames
-         & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames
-         & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames
-         & cn_slchlnonfbfiles, & ! Surface non-diatom         log10(chlorophyll) input filenames
-         & cn_slchldinfbfiles, & ! Surface dinoflagellate     log10(chlorophyll) input filenames
-         & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames
-         & cn_slchlnanfbfiles, & ! Surface nanophytoplankton  log10(chlorophyll) input filenames
-         & cn_slchlpicfbfiles, & ! Surface picophytoplankton  log10(chlorophyll) input filenames
-         & cn_schltotfbfiles,  & ! Surface total              chlorophyll        input filenames
-         & cn_slphytotfbfiles, & ! Surface total      log10(phytoplankton carbon) input filenames
-         & cn_slphydiafbfiles, & ! Surface diatom     log10(phytoplankton carbon) input filenames
-         & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames
-         & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames
-         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames
-         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames
-         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames
-         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames
-         & cn_pno3fbfiles,     & ! Profile nitrate input filenames
-         & cn_psi4fbfiles,     & ! Profile silicate input filenames
-         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames
-         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames
-         & cn_palkfbfiles,     & ! Profile alkalinity input filenames
-         & cn_pphfbfiles,      & ! Profile pH input filenames
-         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames
-         & cn_sstbiasfiles       ! SST bias input filenames
-
-      CHARACTER(LEN=128) :: &
-         & cn_altbiasfile        ! Altimeter bias input filename
-
-
-      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
-      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
-      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies 
-      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
-      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
-      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs
-      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
-      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs
-      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs
-      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs
-      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs
-      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs
-      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs
-      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs
-      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs
-      LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs
-      LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs
-      LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs
-      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs
-      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs
-      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs
-      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs
-      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs
-      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs
-      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs
-      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs
-      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs
-      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs
-      LOGICAL :: ln_pph          ! Logical switch for profile pH obs
-      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs
-      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
-      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
-      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST
-      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
-      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
-      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary
-
-      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
-      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
-
-      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl
-      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl
-
-      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
-         & clproffiles, &        ! Profile filenames
-         & clsurffiles           ! Surface filenames
-
-      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read
-      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
-      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables)
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & zglam                 ! Model longitudes for profile variables
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & zgphi                 ! Model latitudes for profile variables
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
-         & zmask                 ! Model land/sea mask associated with variables
-
-
-      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
-         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
-         &            ln_slchltot, ln_slchldia, ln_slchlnon,          &
-         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          &
-         &            ln_slchlpic, ln_schltot,                        &
-         &            ln_slphytot, ln_slphydia, ln_slphynon,          &
-         &            ln_sspm,     ln_sfco2,    ln_spco2,             &
-         &            ln_plchltot, ln_pchltot,  ln_pno3,              &
-         &            ln_psi4,     ln_ppo4,     ln_pdic,              &
-         &            ln_palk,     ln_pph,      ln_po2,               &
-         &            ln_altbias, ln_sstbias, ln_nea,                 &
-         &            ln_grid_global, ln_grid_search_lookup,          &
-         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
-         &            ln_sstnight, ln_default_fp_indegs,              &
-         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
-         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
-         &            cn_profbfiles, cn_slafbfiles,                   &
-         &            cn_sstfbfiles, cn_sicfbfiles,                   &
-         &            cn_velfbfiles, cn_sssfbfiles,                   &
-         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         &
-         &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         &
-         &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         &
-         &            cn_slchlpicfbfiles, cn_schltotfbfiles,          &
-         &            cn_slphytotfbfiles, cn_slphydiafbfiles,         &
-         &            cn_slphynonfbfiles, cn_sspmfbfiles,             &
-         &            cn_sfco2fbfiles, cn_spco2fbfiles,               &
-         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          &
-         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, &
-         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  &
-         &            cn_po2fbfiles,                                  &
-         &            cn_sstbiasfiles, cn_altbiasfile,                &
-         &            cn_gridsearchfile, rn_gridsearchres,            &
-         &            rn_dobsini, rn_dobsend,                         &
-         &            rn_default_avglamscl, rn_default_avgphiscl,     &
-         &            rn_sla_avglamscl, rn_sla_avgphiscl,             &
-         &            rn_sst_avglamscl, rn_sst_avgphiscl,             &
-         &            rn_sss_avglamscl, rn_sss_avgphiscl,             &
-         &            rn_sic_avglamscl, rn_sic_avgphiscl,             &
-         &            nn_1dint, nn_2dint_default,                     &
-         &            nn_2dint_sla, nn_2dint_sst,                     &
-         &            nn_2dint_sss, nn_2dint_sic,                     &
-         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
+      INTEGER, PARAMETER :: & 
+         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type 
+      INTEGER, DIMENSION(:), ALLOCATABLE :: & 
+         & ifilesprof, &         ! Number of profile files 
+         & ifilessurf            ! Number of surface files 
+      INTEGER :: ios             ! Local integer output status for namelist read 
+      INTEGER :: jtype           ! Counter for obs types 
+      INTEGER :: jvar            ! Counter for variables 
+      INTEGER :: jfile           ! Counter for files 
+      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
+      INTEGER :: n2dint_type     ! Local version of nn_2dint* 
+
+      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
+         & cn_profbfiles,      & ! T/S profile input filenames 
+         & cn_sstfbfiles,      & ! Sea surface temperature input filenames 
+         & cn_slafbfiles,      & ! Sea level anomaly input filenames 
+         & cn_sicfbfiles,      & ! Seaice concentration input filenames 
+         & cn_sitfbfiles,      & ! Seaice thickness input filenames 
+         & cn_fbdfbfiles,      & ! Seaice freeboard input filenames 
+         & cn_velfbfiles,      & ! Velocity profile input filenames 
+         & cn_sssfbfiles,      & ! Sea surface salinity input filenames 
+         & cn_ssvfbfiles,      & ! Sea surface velocity 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_skd490fbfiles,   & ! Surface Kd490 input filenames 
+         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames 
+         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames 
+         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 
+         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames 
+         & cn_pno3fbfiles,     & ! Profile nitrate input filenames 
+         & cn_psi4fbfiles,     & ! Profile silicate input filenames 
+         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames 
+         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames 
+         & cn_palkfbfiles,     & ! Profile alkalinity input filenames 
+         & cn_pphfbfiles,      & ! Profile pH input filenames 
+         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames 
+         & cn_sstbiasfiles       ! SST bias input filenames 
+
+      CHARACTER(LEN=128) :: & 
+         & cn_altbiasfile        ! Altimeter bias input filename 
+
+
+      LOGICAL :: ln_seaicetypes = .FALSE.          ! Logical switch indicating data type is sea ice 
+      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
+      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles 
+      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies  
+      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
+      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
+      LOGICAL :: ln_sit          ! Logical switch for sea ice thickness 
+      LOGICAL :: ln_fbd          ! Logical switch for sea ice freeboard 
+      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
+      LOGICAL :: ln_ssv          ! Logical switch for sea surface velocity obs       
+      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
+      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs 
+      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs 
+      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs 
+      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs 
+      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 
+      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs 
+      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs 
+      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs 
+      LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs 
+      LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs 
+      LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 
+      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs 
+      LOGICAL :: ln_skd490       ! Logical switch for surface Kd490 
+      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs 
+      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs 
+      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs 
+      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs 
+      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs 
+      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs 
+      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs 
+      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs 
+      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs 
+      LOGICAL :: ln_pph          ! Logical switch for profile pH obs 
+      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs 
+      LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
+      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
+      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST 
+      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
+      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
+      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 
+
+      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
+      REAL(dp) :: rn_dobsend     ! Obs window end date YYYYMMDD.HHMMSS 
+
+      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 
+      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 
+
+      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
+         & clproffiles, &        ! Profile filenames 
+         & clsurffiles           ! Surface filenames 
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars  ! Expected variable names 
+
+      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) 
+      LOGICAL :: ltype_clim      ! Local version of ln_output_clim 
+
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
+         & zglam                 ! Model longitudes for profile variables 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
+         & zgphi                 ! Model latitudes for profile variables 
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
+         & zmask                 ! Model land/sea mask associated with variables 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
+         & zmask_surf            ! Surface model land/sea mask associated with variables 
+
+
+      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
+         &            ln_sst, ln_sic, ln_sit, ln_fbd,                 & 
+         &            ln_sss, ln_ssv, ln_vel3d,                       & 
+         &            ln_slchltot, ln_slchldia, ln_slchlnon,          & 
+         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          & 
+         &            ln_slchlpic, ln_schltot,                        & 
+         &            ln_slphytot, ln_slphydia, ln_slphynon,          & 
+         &            ln_sspm,     ln_sfco2,    ln_spco2,             & 
+         &            ln_skd490,                                      & 
+         &            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_output_clim,                   & 
+         &            ln_time_mean_sla_bkg, ln_default_fp_indegs,     & 
+         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
+         &            ln_sss_fp_indegs, ln_ssv_fp_indegs,             & 
+         &            ln_sic_fp_indegs, ln_sit_fp_indegs,             & 
+         &            ln_fbd_fp_indegs,                               & 
+         &            cn_profbfiles, cn_slafbfiles,                   & 
+         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
+         &            cn_sitfbfiles, cn_fbdfbfiles,                   & 
+         &            cn_velfbfiles, cn_sssfbfiles, cn_ssvfbfiles,    & 
+         &            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_skd490fbfiles,                               & 
+         &            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_ssv_avglamscl, rn_ssv_avgphiscl,             &          
+         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
+         &            rn_sit_avglamscl, rn_sit_avgphiscl,             & 
+         &            rn_fbd_avglamscl, rn_fbd_avgphiscl,             & 
+         &            nn_1dint, nn_2dint_default,                     & 
+         &            nn_2dint_sla, nn_2dint_sst,                     & 
+         &            nn_2dint_sss, nn_2dint_ssv,                     & 
+         &            nn_2dint_sic, nn_2dint_sit,                     & 
+         &            nn_2dint_fbd,                                   & 
+         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
          &            nn_profdavtypes
 
@@ -289,41 +337,44 @@
       !-----------------------------------------------------------------------
 
-      ! Some namelist arrays need initialising
-      cn_profbfiles(:)      = ''
-      cn_slafbfiles(:)      = ''
-      cn_sstfbfiles(:)      = ''
-      cn_sicfbfiles(:)      = ''
-      cn_velfbfiles(:)      = ''
-      cn_sssfbfiles(:)      = ''
-      cn_slchltotfbfiles(:) = ''
-      cn_slchldiafbfiles(:) = ''
-      cn_slchlnonfbfiles(:) = ''
-      cn_slchldinfbfiles(:) = ''
-      cn_slchlmicfbfiles(:) = ''
-      cn_slchlnanfbfiles(:) = ''
-      cn_slchlpicfbfiles(:) = ''
-      cn_schltotfbfiles(:)  = ''
-      cn_slphytotfbfiles(:) = ''
-      cn_slphydiafbfiles(:) = ''
-      cn_slphynonfbfiles(:) = ''
-      cn_sspmfbfiles(:)     = ''
-      cn_sfco2fbfiles(:)    = ''
-      cn_spco2fbfiles(:)    = ''
-      cn_plchltotfbfiles(:) = ''
-      cn_pchltotfbfiles(:)  = ''
-      cn_pno3fbfiles(:)     = ''
-      cn_psi4fbfiles(:)     = ''
-      cn_ppo4fbfiles(:)     = ''
-      cn_pdicfbfiles(:)     = ''
-      cn_palkfbfiles(:)     = ''
-      cn_pphfbfiles(:)      = ''
-      cn_po2fbfiles(:)      = ''
-      cn_sstbiasfiles(:)    = ''
-      nn_profdavtypes(:)    = -1
-
-      CALL ini_date( rn_dobsini )
-      CALL fin_date( rn_dobsend )
-
-      ! Read namelist namobs : control observation diagnostics
+      cn_profbfiles(:)      = '' 
+      cn_slafbfiles(:)      = '' 
+      cn_sstfbfiles(:)      = '' 
+      cn_sicfbfiles(:)      = '' 
+      cn_sitfbfiles(:)      = '' 
+      cn_fbdfbfiles(:)      = '' 
+      cn_velfbfiles(:)      = '' 
+      cn_sssfbfiles(:)      = '' 
+      cn_ssvfbfiles(:)      = ''       
+      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_skd490fbfiles(:)   = '' 
+      cn_sfco2fbfiles(:)    = '' 
+      cn_spco2fbfiles(:)    = '' 
+      cn_plchltotfbfiles(:) = '' 
+      cn_pchltotfbfiles(:)  = '' 
+      cn_pno3fbfiles(:)     = '' 
+      cn_psi4fbfiles(:)     = '' 
+      cn_ppo4fbfiles(:)     = '' 
+      cn_pdicfbfiles(:)     = '' 
+      cn_palkfbfiles(:)     = '' 
+      cn_pphfbfiles(:)      = '' 
+      cn_po2fbfiles(:)      = '' 
+      cn_sstbiasfiles(:)    = '' 
+      nn_profdavtypes(:)    = -1 
+
+      CALL ini_date( rn_dobsini ) 
+      CALL fin_date( rn_dobsend ) 
+
+      ! Read namelist namobs : control observation diagnostics 
       REWIND( numnam_ref )   ! Namelist namobs in reference namelist
       READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
@@ -351,68 +402,83 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
-         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
-         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
-         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
-         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
-         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
-         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
-         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss
-         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot
-         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia
-         WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon
-         WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin
-         WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic
-         WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan
-         WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic
-         WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot
-         WRITE(numout,*) '             Logical switch for surface total log(phyC) obs      ln_slphytot = ', ln_slphytot
-         WRITE(numout,*) '             Logical switch for surface diatom log(phyC) obs     ln_slphydia = ', ln_slphydia
-         WRITE(numout,*) '             Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon
-         WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm
-         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2
-         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2
-         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot
-         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot
-         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3
-         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4
-         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4
-         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic
-         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk
-         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph
-         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2
-         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global
+         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d 
+         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d 
+         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla 
+         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst 
+         WRITE(numout,*) '             Logical switch for SIC observations                      ln_sic = ', ln_sic 
+         WRITE(numout,*) '             Logical switch for SIT observations                      ln_sit = ', ln_sit 
+         WRITE(numout,*) '             Logical switch for FBD observations                      ln_fbd = ', ln_fbd 
+         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 SSV observations                      ln_ssv = ', ln_ssv          
+         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 Kd490 observations         ln_skd490 = ', ln_skd490 
+         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2 
+         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2 
+         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot 
+         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot 
+         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3 
+         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4 
+         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4 
+         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic 
+         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk 
+         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph 
+         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2 
+         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
          WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
          IF (ln_grid_search_lookup) &
             WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
-         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
-         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
-         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
-         WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default
-         WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla
-         WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst
-         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss
-         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic
-         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl
-         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl
-         WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs
-         WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl
-         WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl
-         WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs
-         WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl
-         WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl
-         WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs
-         WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl
-         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl
-         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs
-         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
-         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject
-         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
-         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
-         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
-         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
-         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias
-         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
-         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
-         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
+            WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini 
+            WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
+            WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
+            WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default 
+            WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla 
+            WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst 
+            WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss 
+            WRITE(numout,*) '             Type of horizontal interpolation method for SSV    nn_2dint_ssv = ', nn_2dint_ssv          
+            WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
+            WRITE(numout,*) '             Type of horizontal interpolation method for SIT    nn_2dint_sit = ', nn_2dint_sit 
+            WRITE(numout,*) '             Type of horizontal interpolation method for FBD    nn_2dint_fbd = ', nn_2dint_fbd 
+            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,*) '             SIT E/W diameter of obs footprint              rn_sit_avglamscl = ', rn_sit_avglamscl 
+            WRITE(numout,*) '             SIT N/S diameter of obs footprint              rn_sit_avgphiscl = ', rn_sit_avgphiscl 
+            WRITE(numout,*) '             SIT obs footprint in deg [T] or m [F]          ln_sit_fp_indegs = ', ln_sit_fp_indegs 
+            WRITE(numout,*) '             FBD E/W diameter of obs footprint              rn_fbd_avglamscl = ', rn_fbd_avglamscl 
+            WRITE(numout,*) '             FBD N/S diameter of obs footprint              rn_fbd_avgphiscl = ', rn_fbd_avgphiscl 
+            WRITE(numout,*) '             FBD obs footprint in deg [T] or m [F]          ln_fbd_fp_indegs = ', ln_fbd_fp_indegs 
+            WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
+            WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
+            WRITE(numout,*) '             MSSH correction scheme nn_msshc = ', nn_msshc 
+            WRITE(numout,*) '             MDT  correction rn_mdtcorr = ', rn_mdtcorr 
+            WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
+            WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
+            WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
+            WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
+            WRITE(numout,*) '             Daily average types nn_profdavtypes = ', nn_profdavtypes 
+            WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
+            WRITE(numout,*) '             Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim 
+            WRITE(numout,*) '             Logical switch for time-mean of SLA        ln_time_mean_sla_bkg = ', ln_time_mean_sla_bkg
       ENDIF
       !-----------------------------------------------------------------------
@@ -421,12 +487,13 @@
       !-----------------------------------------------------------------------
 
-      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 /) )
+      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_sit, ln_fbd,             & 
+         &                  ln_sss, ln_ssv,                                     & 
+         &                  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_skd490,   ln_sfco2,    ln_spco2 /) ) 
 
       IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
@@ -439,4 +506,12 @@
       ENDIF
 
+      IF ( ln_output_clim .AND. ( .NOT. ln_tradmp ) ) THEN 
+         IF(lwp) WRITE(numout,cform_war) 
+         IF(lwp) WRITE(numout,*) ' ln_output_clim is true, but ln_tradmp is false', & 
+            &                    ' so climatological T/S not available and will not be output' 
+         nwarn = nwarn + 1 
+         ln_output_clim = .FALSE. 
+      ENDIF
+
       IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes
       IF ( nproftypes > 0 ) THEN
@@ -535,4 +610,14 @@
             clsurffiles(jtype,:) = cn_sicfbfiles
          ENDIF
+         IF (ln_sit) THEN 
+            jtype = jtype + 1 
+            cobstypessurf(jtype) = 'sit' 
+            clsurffiles(jtype,:) = cn_sitfbfiles 
+         ENDIF 
+         IF (ln_fbd) THEN 
+            jtype = jtype + 1 
+            cobstypessurf(jtype) = 'fbd' 
+            clsurffiles(jtype,:) = cn_fbdfbfiles 
+         ENDIF
          IF (ln_sss) THEN
             jtype = jtype + 1
@@ -540,4 +625,9 @@
             clsurffiles(jtype,:) = cn_sssfbfiles
          ENDIF
+         IF (ln_ssv) THEN 
+            jtype = jtype + 1 
+            cobstypessurf(jtype) = 'ssv' 
+            clsurffiles(jtype,:) = cn_ssvfbfiles 
+         ENDIF
          IF (ln_slchltot) THEN
             jtype = jtype + 1
@@ -599,4 +689,9 @@
             cobstypessurf(jtype) = 'sspm'
             clsurffiles(jtype,:) = cn_sspmfbfiles
+         ENDIF
+         IF (ln_skd490) THEN 
+            jtype = jtype + 1 
+            cobstypessurf(jtype) = 'skd490' 
+            clsurffiles(jtype,:) = cn_skd490fbfiles 
          ENDIF
          IF (ln_sfco2) THEN
@@ -645,4 +740,24 @@
                ltype_fp_indegs = ln_sic_fp_indegs
                ltype_night     = .FALSE.
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 
+               IF ( nn_2dint_sit == -1 ) THEN 
+                  n2dint_type  = nn_2dint_default 
+               ELSE 
+                  n2dint_type  = nn_2dint_sit 
+               ENDIF 
+               ztype_avglamscl = rn_sit_avglamscl 
+               ztype_avgphiscl = rn_sit_avgphiscl 
+               ltype_fp_indegs = ln_sit_fp_indegs 
+               ltype_night     = .FALSE. 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN 
+               IF ( nn_2dint_fbd == -1 ) THEN 
+                  n2dint_type  = nn_2dint_default 
+               ELSE 
+                  n2dint_type  = nn_2dint_fbd 
+               ENDIF 
+               ztype_avglamscl = rn_fbd_avglamscl 
+               ztype_avgphiscl = rn_fbd_avgphiscl 
+               ltype_fp_indegs = ln_fbd_fp_indegs 
+               ltype_night     = .FALSE.
             ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
                IF ( nn_2dint_sss == -1 ) THEN
@@ -654,4 +769,14 @@
                ztype_avgphiscl = rn_sss_avgphiscl
                ltype_fp_indegs = ln_sss_fp_indegs
+               ltype_night     = .FALSE.
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
+               IF ( nn_2dint_ssv == -1 ) THEN 
+                  n2dint_type  = nn_2dint_default 
+               ELSE 
+                  n2dint_type  = nn_2dint_ssv 
+               ENDIF 
+               ztype_avglamscl = rn_ssv_avglamscl 
+               ztype_avgphiscl = rn_ssv_avgphiscl 
+               ltype_fp_indegs = ln_ssv_fp_indegs 
                ltype_night     = .FALSE.
             ELSE
@@ -715,8 +840,12 @@
          DO jtype = 1, nproftypes
 
+            ltype_clim = .FALSE.
+
             IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
                nvarsprof(jtype) = 2
                nextrprof(jtype) = 1
+               IF ( ln_output_clim ) ltype_clim = .TRUE.
                ALLOCATE(llvar(nvarsprof(jtype)))
+               ALLOCATE(clvars(nvarsprof(jtype)))
                CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
                CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
@@ -734,4 +863,5 @@
                nextrprof(jtype) = 2
                ALLOCATE(llvar(nvarsprof(jtype)))
+               ALLOCATE(clvars(nvarsprof(jtype)))
                CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
                CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
@@ -739,4 +869,6 @@
                llvar(1)       = ln_vel3d
                llvar(2)       = ln_vel3d
+               clvars(1)      = 'UVEL' 
+               clvars(2)      = 'VVEL'
                zglam(:,:,1)   = glamu(:,:)
                zglam(:,:,2)   = glamv(:,:)
@@ -749,4 +881,5 @@
                nextrprof(jtype) = 0
                ALLOCATE(llvar(nvarsprof(jtype)))
+               ALLOCATE(clvars(nvarsprof(jtype)))
                CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
                CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
@@ -756,4 +889,23 @@
                zgphi(:,:,1)   = gphit(:,:)
                zmask(:,:,:,1) = tmask(:,:,:)
+               IF ( TRIM(cobstypesprof(jtype)) == 'plchltot' ) THEN 
+                  clvars(1) = 'PLCHLTOT' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pchltot' )  THEN 
+                  clvars(1) = 'PCHLTOT' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pno3' ) THEN 
+                  clvars(1) = 'PNO3' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'psi4' ) THEN 
+                  clvars(1) = 'PSI4' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'ppo4' ) THEN 
+                  clvars(1) = 'PPO4' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pdic' ) THEN 
+                  clvars(1) = 'PDIC' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'palk' ) THEN 
+                  clvars(1) = 'PALK' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pph' ) THEN 
+                  clvars(1) = 'PPH' 
+               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'po2' ) THEN 
+                  clvars(1) = 'PO2' 
+               ENDIF
             ENDIF
 
@@ -763,5 +915,5 @@
                &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
                &               rn_dobsini, rn_dobsend, llvar, &
-               &               ln_ignmis, ln_s_at_t, .FALSE., &
+               &               ln_ignmis, ln_s_at_t, .FALSE., ltype_clim, clvars, &
                &               kdailyavtypes = nn_profdavtypes )
 
@@ -797,39 +949,139 @@
          DO jtype = 1, nsurftypes
 
-            nvarssurf(jtype) = 1
-            nextrsurf(jtype) = 0
-            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
+            ltype_clim = .FALSE. 
+            IF ( ln_output_clim .AND. & 
+               & ( ( TRIM(cobstypessurf(jtype)) == 'sst' ) .OR. & 
+               &   ( TRIM(cobstypessurf(jtype)) == 'sss' ) ) ) & 
+               & ltype_clim = .TRUE. 
+
+            IF ( (TRIM(cobstypessurf(jtype)) == 'sla') .OR. & 
+               & (TRIM(cobstypessurf(jtype)) == 'sit') .OR. & 
+               & (TRIM(cobstypessurf(jtype)) == 'fbd') ) THEN 
+               nvarssurf(jtype) = 1 
+               nextrsurf(jtype) = 2 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
+               nvarssurf(jtype) = 2 
+               nextrsurf(jtype) = 0 
+            ELSE 
+               nvarssurf(jtype) = 1 
+               nextrsurf(jtype) = 0 
+            ENDIF 
+
+            ALLOCATE( clvars( nvarssurf(jtype) ) ) 
+            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zglam ) 
+            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zgphi ) 
+            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zmask_surf ) 
+  
+            IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
+               zglam(:,:,1) = glamu(:,:) 
+               zglam(:,:,2) = glamv(:,:) 
+               zgphi(:,:,1) = gphiu(:,:) 
+               zgphi(:,:,2) = gphiv(:,:) 
+               zmask_surf(:,:,1) = umask(:,:,1) 
+               zmask_surf(:,:,2) = vmask(:,:,1) 
+            ELSE             
+               DO jvar = 1, nvarssurf(jtype) 
+                  zglam(:,:,jvar) = glamt(:,:) 
+                  zgphi(:,:,jvar) = gphit(:,:) 
+                  zmask_surf(:,:,jvar) = tmask(:,:,1)                                        
+               END DO
+            ENDIF
+
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
+               clvars(1) = 'SLA' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 
+               clvars(1) = 'SST' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
+               clvars(1) = 'ICECONC' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN 
+               clvars(1) = 'FBD' 
+               ln_seaicetypes = .TRUE. 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 
+               clvars(1) = 'SIT' 
+               ln_seaicetypes = .TRUE. 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
+               clvars(1) = 'SSS' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
+               clvars(1) = 'UVEL' 
+               clvars(2) = 'VVEL'            
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchltot' ) THEN 
+               clvars(1) = 'SLCHLTOT' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldia' ) THEN 
+               clvars(1) = 'SLCHLDIA' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnon' ) THEN 
+               clvars(1) = 'SLCHLNON' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldin' ) THEN 
+               clvars(1) = 'SLCHLDIN' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlmic' ) THEN 
+               clvars(1) = 'SLCHLMIC' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnan' ) THEN 
+               clvars(1) = 'SLCHLNAN' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlpic' ) THEN 
+               clvars(1) = 'SLCHLPIC' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'schltot' ) THEN 
+               clvars(1) = 'SCHLTOT' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphytot' ) THEN 
+               clvars(1) = 'SLPHYTOT' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphydia' ) THEN 
+               clvars(1) = 'SLPHYDIA' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphynon' ) THEN 
+               clvars(1) = 'SLPHYNON' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sspm' ) THEN 
+               clvars(1) = 'SSPM' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'skd490' ) THEN 
+               clvars(1) = 'SKD490' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sfco2' ) THEN 
+               clvars(1) = 'SFCO2' 
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'spco2' ) THEN 
+               clvars(1) = 'SPCO2' 
+            ENDIF
 
             !Read in surface obs types
-            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
-               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
-               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
-               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
-
-            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
-
-            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
-               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
-               IF ( ln_altbias ) &
-                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
-            ENDIF
-
-            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
-               jnumsstbias = 0
-               DO jfile = 1, jpmaxnfiles
-                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
-                     &  jnumsstbias = jnumsstbias + 1
-               END DO
-               IF ( jnumsstbias == 0 ) THEN
-                  CALL ctl_stop("ln_sstbias set but no bias files to read in")    
+            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 
+               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
+               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
+               &               rn_dobsini, rn_dobsend, MeanPeriodHours, ln_ignmis, .FALSE., & 
+               &               llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars ) 
+ 
+            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), & 
+               &               jpi, jpj, &             
+               &               zmask_surf, zglam, zgphi, & 
+               &               ln_nea, ln_bound_reject, ln_seaicetypes ) 
+ 
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
+               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
+               IF ( ln_altbias ) & 
+                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
+            ENDIF 
+ 
+#if defined key_cice 
+            IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN 
+               CALL obs_rea_snowdepth( surfdataqc(jtype), n2dintsurf(jtype), thick_s(:,:) ) 
+            ENDIF  
+#endif 
+ 
+            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
+               jnumsstbias = 0 
+               DO jfile = 1, jpmaxnfiles 
+                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
+                     &  jnumsstbias = jnumsstbias + 1 
+               END DO 
+               IF ( jnumsstbias == 0 ) THEN 
+                  CALL ctl_stop("ln_sstbias set but no bias files to read in")     
+ 
                ENDIF
 
-               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
-                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
-
-            ENDIF
-
-         END DO
-
+               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
+                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
+ 
+            ENDIF 
+             
+            DEALLOCATE( clvars ) 
+            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zglam ) 
+            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zgphi ) 
+            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zmask_surf ) 
+ 
+         END DO 
+ 
          DEALLOCATE( ifilessurf, clsurffiles )
 
@@ -879,38 +1131,38 @@
          & frld
 #endif
-#if defined key_cice
-      USE sbc_oce, ONLY : fr_i     ! ice fraction
-#endif
-#if defined key_top
-      USE trc, ONLY :  &           ! Biogeochemical state variables
-         & trn
-#endif
-#if defined key_hadocc
-      USE par_hadocc               ! HadOCC parameters
-      USE trc, ONLY :  &
-         & HADOCC_CHL, &
-         & HADOCC_FCO2, &
-         & HADOCC_PCO2, &
-         & HADOCC_FILL_FLT
-      USE had_bgc_const, ONLY: c2n_p
-#elif defined key_medusa
-      USE par_medusa               ! MEDUSA parameters
-      USE sms_medusa, ONLY: &
-         & xthetapn, &
-         & xthetapd
-#if defined key_roam
-      USE sms_medusa, ONLY: &
-         & f2_pco2w, &
-         & f2_fco2w, &
-         & f3_pH
-#endif
-#elif defined key_fabm
-      USE par_fabm                 ! FABM parameters
-      USE fabm, ONLY: &
-         & fabm_get_interior_diagnostic_data
-#endif
-#if defined key_spm
-      USE par_spm, ONLY: &         ! Sediment parameters
-         & jp_spm
+#if defined key_cice 
+      USE sbc_oce, ONLY : &        ! CICE variables 
+         & fr_i,          &        ! ice fraction 
+         & thick_i                 ! ice thickness 
+#endif 
+#if defined key_top 
+      USE trc, ONLY :  &           ! Biogeochemical state variables 
+         & trn 
+#endif 
+#if defined key_hadocc 
+      USE par_hadocc               ! HadOCC parameters 
+      USE trc, ONLY :  & 
+         & HADOCC_CHL, & 
+         & HADOCC_FCO2, & 
+         & HADOCC_PCO2, & 
+         & HADOCC_FILL_FLT 
+      USE had_bgc_const, ONLY: c2n_p 
+#elif defined key_medusa 
+      USE par_medusa               ! MEDUSA parameters 
+      USE sms_medusa, ONLY: & 
+         & xthetapn, & 
+         & xthetapd 
+#if defined key_roam 
+      USE sms_medusa, ONLY: & 
+         & f2_pco2w, & 
+         & f2_fco2w, & 
+         & f3_pH 
+#endif 
+#elif defined key_fabm 
+      USE par_fabm                 ! FABM parameters 
+#endif 
+#if defined key_spm 
+      USE par_spm, ONLY: &         ! Sediment parameters 
+         & jp_spm 
 #endif
 
@@ -920,23 +1172,26 @@
       INTEGER, INTENT(IN) :: kstp  ! Current timestep
       !! * Local declarations
-      INTEGER :: idaystp           ! Number of timesteps per day
-      INTEGER :: jtype             ! Data loop variable
-      INTEGER :: jvar              ! Variable number
-      INTEGER :: ji, jj, jk        ! Loop counters
-      REAL(wp) :: tiny             ! small number
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
-         & zprofvar                ! Model values for variables in a prof ob
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
-         & zprofmask               ! Mask associated with zprofvar
-      REAL(wp), POINTER, DIMENSION(:,:) :: &
-         & zsurfvar, &             ! Model values equivalent to surface ob.
-         & zsurfmask               ! Mask associated with surface variable
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & zglam,    &             ! Model longitudes for prof variables
-         & zgphi                   ! Model latitudes for prof variables
-      LOGICAL :: llog10            ! Perform log10 transform of variable
-#if defined key_fabm
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & pco2_3d                 ! 3D pCO2 from FABM
+      INTEGER :: idaystp           ! Number of timesteps per day 
+      INTEGER :: imeanstp          ! Number of timesteps for sla averaging 
+      INTEGER :: jtype             ! Data loop variable 
+      INTEGER :: jvar              ! Variable number 
+      INTEGER :: ji, jj, jk        ! Loop counters 
+      REAL(wp) :: tiny             ! small number 
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
+         & zprofvar, &             ! Model values for variables in a prof ob 
+         & zprofclim               ! Climatology 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. 
+         & zsurfclim, &            ! Climatology values for variables in a surface ob. 
+         & zsurfmask               ! Mask associated with surface variable 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
+         & zglam,    &             ! Model longitudes 
+         & zgphi                   ! Model latitudes 
+      LOGICAL :: llog10            ! Perform log10 transform of variable 
+#if defined key_fabm 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
+         & fabm_3d                 ! 3D variable from FABM 
 #endif
 
@@ -954,549 +1209,663 @@
       !-----------------------------------------------------------------------
 
-      IF ( nproftypes > 0 ) THEN
-
-         DO jtype = 1, nproftypes
-
-            ! Allocate local work arrays
-            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
-            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
-            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
-            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
-            
-            ! Defaults which might change
-            DO jvar = 1, profdataqc(jtype)%nvar
-               zprofmask(:,:,:,jvar) = tmask(:,:,:)
-               zglam(:,:,jvar)       = glamt(:,:)
-               zgphi(:,:,jvar)       = gphit(:,:)
-            END DO
-
-            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
-
-            CASE('prof')
-               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
-               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
-
-            CASE('vel')
-               zprofvar(:,:,:,1) = un(:,:,:)
-               zprofvar(:,:,:,2) = vn(:,:,:)
-               zprofmask(:,:,:,1) = umask(:,:,:)
-               zprofmask(:,:,:,2) = vmask(:,:,:)
-               zglam(:,:,1) = glamu(:,:)
-               zglam(:,:,2) = glamv(:,:)
-               zgphi(:,:,1) = gphiu(:,:)
-               zgphi(:,:,2) = gphiv(:,:)
-
-            CASE('plchltot')
-#if defined key_hadocc
-               ! Chlorophyll from HadOCC
-               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
-#elif defined key_medusa
-               ! Add non-diatom and diatom chlorophyll from MEDUSA
-               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
-#elif defined key_fabm
-               ! Add all chlorophyll groups from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
-                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
-#else
-               CALL ctl_stop( ' Trying to run plchltot observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               ! Take the log10 where we can, otherwise exclude
-               tiny = 1.0e-20
-               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt )
-                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:))
-               ELSEWHERE
-                  zprofvar(:,:,:,:)  = obfillflt
-                  zprofmask(:,:,:,:) = 0
-               END WHERE
-               ! Mask out model below any excluded values,
-               ! to avoid interpolation issues
-               DO jvar = 1, profdataqc(jtype)%nvar
-                 DO jj = 1, jpj
-                    DO ji = 1, jpi
-                       depth_loop: DO jk = 1, jpk
-                          IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN
-                             zprofmask(ji,jj,jk:jpk,jvar) = 0
-                             EXIT depth_loop
-                          ENDIF
-                       END DO depth_loop
-                    END DO
-                 END DO
-              END DO
-
-            CASE('pchltot')
-#if defined key_hadocc
-               ! Chlorophyll from HadOCC
-               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
-#elif defined key_medusa
-               ! Add non-diatom and diatom chlorophyll from MEDUSA
-               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
-#elif defined key_fabm
-               ! Add all chlorophyll groups from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
-                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
-#else
-               CALL ctl_stop( ' Trying to run pchltot observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('pno3')
-#if defined key_hadocc
-               ! Dissolved inorganic nitrogen from HadOCC
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut)
-#elif defined key_medusa
-               ! Dissolved inorganic nitrogen from MEDUSA
-               zprofvar(:,:,:,1) = trn(:,:,:,jpdin)
-#elif defined key_fabm
-               ! Nitrate from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n)
-#else
-               CALL ctl_stop( ' Trying to run pno3 observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('psi4')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run psi4 observation operator', &
-                  &           ' but HadOCC does not simulate silicate' )
-#elif defined key_medusa
-               ! Silicate from MEDUSA
-               zprofvar(:,:,:,1) = trn(:,:,:,jpsil)
-#elif defined key_fabm
-               ! Silicate from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s)
-#else
-               CALL ctl_stop( ' Trying to run psi4 observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('ppo4')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
-                  &           ' but HadOCC does not simulate phosphate' )
-#elif defined key_medusa
-               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
-                  &           ' but MEDUSA does not simulate phosphate' )
-#elif defined key_fabm
-               ! Phosphate from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p)
-#else
-               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('pdic')
-#if defined key_hadocc
-               ! Dissolved inorganic carbon from HadOCC
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic)
-#elif defined key_medusa
-               ! Dissolved inorganic carbon from MEDUSA
-               zprofvar(:,:,:,1) = trn(:,:,:,jpdic)
-#elif defined key_fabm
-               ! Dissolved inorganic carbon from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
-#else
-               CALL ctl_stop( ' Trying to run pdic observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('palk')
-#if defined key_hadocc
-               ! Alkalinity from HadOCC
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk)
-#elif defined key_medusa
-               ! Alkalinity from MEDUSA
-               zprofvar(:,:,:,1) = trn(:,:,:,jpalk)
-#elif defined key_fabm
-               ! Alkalinity from ERSEM
-               zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
-#else
-               CALL ctl_stop( ' Trying to run palk observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('pph')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run pph observation operator', &
-                  &           ' but HadOCC has no pH diagnostic defined' )
-#elif defined key_medusa && defined key_roam
-               ! pH from MEDUSA
-               zprofvar(:,:,:,1) = f3_pH(:,:,:)
-#elif defined key_fabm
-               ! pH from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph)
-#else
-               CALL ctl_stop( ' Trying to run pph observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('po2')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run po2 observation operator', &
-                  &           ' but HadOCC does not simulate oxygen' )
-#elif defined key_medusa
-               ! Oxygen from MEDUSA
-               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)
-#elif defined key_fabm
-               ! Oxygen from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o)
-#else
-               CALL ctl_stop( ' Trying to run po2 observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE DEFAULT
-               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
-
-            END SELECT
-
-            DO jvar = 1, profdataqc(jtype)%nvar
-               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
-                  &               nit000, idaystp, jvar,                   &
-                  &               zprofvar(:,:,:,jvar),                    &
-                  &               fsdept(:,:,:), fsdepw(:,:,:),            & 
-                  &               zprofmask(:,:,:,jvar),                   &
-                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &
-                  &               nn_1dint, nn_2dint_default,              &
-                  &               kdailyavtypes = nn_profdavtypes )
-            END DO
-
-            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
-            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
-            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
-            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
-
-         END DO
-
-      ENDIF
-
-      IF ( nsurftypes > 0 ) THEN
-
-         !Allocate local work arrays
-         CALL wrk_alloc( jpi, jpj, zsurfvar )
-         CALL wrk_alloc( jpi, jpj, zsurfmask )
-#if defined key_fabm
-         CALL wrk_alloc( jpi, jpj, jpk, pco2_3d )
-#endif
-
-         DO jtype = 1, nsurftypes
-
-            !Defaults which might be changed
-            zsurfmask(:,:) = tmask(:,:,1)
-            llog10 = .FALSE.
-
-            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
-            CASE('sst')
-               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
-            CASE('sla')
-               zsurfvar(:,:) = sshn(:,:)
-            CASE('sss')
-               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
-            CASE('sic')
-               IF ( kstp == 0 ) THEN
-                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
-                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
-                        &           'time-step but some obs are valid then.' )
-                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
-                        &           ' sea-ice obs will be missed'
-                  ENDIF
-                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
-                     &                        surfdataqc(jtype)%nsstp(1)
-                  CYCLE
-               ELSE
-#if defined key_cice
-                  zsurfvar(:,:) = fr_i(:,:)
-#elif defined key_lim2 || defined key_lim3
-                  zsurfvar(:,:) = 1._wp - frld(:,:)
-#else
-               CALL ctl_stop( ' Trying to run sea-ice observation operator', &
-                  &           ' but no sea-ice model appears to have been defined' )
-#endif
-               ENDIF
-
-            CASE('slchltot')
-#if defined key_hadocc
-               ! Surface chlorophyll from HadOCC
-               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
-#elif defined key_medusa
-               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
-#elif defined key_fabm
-               ! Add all surface chlorophyll groups from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
-                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
-#else
-               CALL ctl_stop( ' Trying to run slchltot observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slchldia')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slchldia observation operator', &
-                  &           ' but HadOCC does not explicitly simulate diatoms' )
-#elif defined key_medusa
-               ! Diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchd)
-#elif defined key_fabm
-               ! Diatom surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1)
-#else
-               CALL ctl_stop( ' Trying to run slchldia observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slchlnon')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
-                  &           ' but HadOCC does not explicitly simulate non-diatoms' )
-#elif defined key_medusa
-               ! Non-diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchn)
-#elif defined key_fabm
-               ! Add all non-diatom surface chlorophyll groups from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
-                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
-#else
-               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slchldin')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slchldin observation operator', &
-                  &           ' but HadOCC does not explicitly simulate dinoflagellates' )
-#elif defined key_medusa
-               CALL ctl_stop( ' Trying to run slchldin observation operator', &
-                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' )
-#elif defined key_fabm
-               ! Dinoflagellate surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
-#else
-               CALL ctl_stop( ' Trying to run slchldin observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slchlmic')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
-                  &           ' but HadOCC does not explicitly simulate microphytoplankton' )
-#elif defined key_medusa
-               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
-                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' )
-#elif defined key_fabm
-               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
-#else
-               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slchlnan')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
-                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' )
-#elif defined key_medusa
-               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
-                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' )
-#elif defined key_fabm
-               ! Nanophytoplankton surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2)
-#else
-               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slchlpic')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
-                  &           ' but HadOCC does not explicitly simulate picophytoplankton' )
-#elif defined key_medusa
-               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
-                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' )
-#elif defined key_fabm
-               ! Picophytoplankton surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3)
-#else
-               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('schltot')
-#if defined key_hadocc
-               ! Surface chlorophyll from HadOCC
-               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
-#elif defined key_medusa
-               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
-#elif defined key_fabm
-               ! Add all surface chlorophyll groups from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
-                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
-#else
-               CALL ctl_stop( ' Trying to run schltot observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('slphytot')
-#if defined key_hadocc
-               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio
-               zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
-#elif defined key_medusa
-               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA
-               ! multiplied by C:N ratio for each
-               zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
-#elif defined key_fabm
-               ! Add all surface phytoplankton carbon groups from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
-                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
-#else
-               CALL ctl_stop( ' Trying to run slphytot observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slphydia')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slphydia observation operator', &
-                  &           ' but HadOCC does not explicitly simulate diatoms' )
-#elif defined key_medusa
-               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
-               zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd
-#elif defined key_fabm
-               ! Diatom surface phytoplankton carbon from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c)
-#else
-               CALL ctl_stop( ' Trying to run slphydia observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('slphynon')
-#if defined key_hadocc
-               CALL ctl_stop( ' Trying to run slphynon observation operator', &
-                  &           ' but HadOCC does not explicitly simulate non-diatoms' )
-#elif defined key_medusa
-               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
-               zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn
-#elif defined key_fabm
-               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
-                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
-#else
-               CALL ctl_stop( ' Trying to run slphynon observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-               llog10 = .TRUE.
-
-            CASE('sspm')
-#if defined key_spm
-               zsurfvar(:,:) = 0.0
-               DO jn = 1, jp_spm
-                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
-               END DO
-#else
-               CALL ctl_stop( ' Trying to run sspm observation operator', &
-                  &           ' but no spm model appears to have been defined' )
-#endif
-
-            CASE('sfco2')
-#if defined key_hadocc
-               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
-               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
-                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
-                  zsurfvar(:,:) = obfillflt
-                  zsurfmask(:,:) = 0
-                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
-                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
-               ENDIF
-#elif defined key_medusa && defined key_roam
-               zsurfvar(:,:) = f2_fco2w(:,:)
-#elif defined key_fabm
-               ! First, get pCO2 from FABM
-               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
-               zsurfvar(:,:) = pco2_3d(:,:,1)
-               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
-               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
-               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
-               ! and
-               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
-               ! Marine Chemistry, 2: 203-215.
-               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
-               ! not explicitly included - atmospheric pressure is not necessarily available so this is
-               ! the best assumption.
-               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
-               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
-               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
-               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
-               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + &
-                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
-                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
-                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
-                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
-                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
-#else
-               CALL ctl_stop( ' Trying to run sfco2 observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE('spco2')
-#if defined key_hadocc
-               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
-               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
-                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
-                  zsurfvar(:,:) = obfillflt
-                  zsurfmask(:,:) = 0
-                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
-                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
-               ENDIF
-#elif defined key_medusa && defined key_roam
-               zsurfvar(:,:) = f2_pco2w(:,:)
-#elif defined key_fabm
-               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
-               zsurfvar(:,:) = pco2_3d(:,:,1)
-#else
-               CALL ctl_stop( ' Trying to run spco2 observation operator', &
-                  &           ' but no biogeochemical model appears to have been defined' )
-#endif
-
-            CASE DEFAULT
-
-               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
-
-            END SELECT
-            
-            IF ( llog10 ) THEN
-               ! Take the log10 where we can, otherwise exclude
-               tiny = 1.0e-20
-               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
-                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:))
-               ELSEWHERE
-                  zsurfvar(:,:)  = obfillflt
-                  zsurfmask(:,:) = 0
-               END WHERE
-            ENDIF
-
-            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
-               &               nit000, idaystp, zsurfvar, zsurfmask,    &
-               &               n2dintsurf(jtype), llnightav(jtype),     &
-               &               ravglamscl(jtype), ravgphiscl(jtype),     &
-               &               lfpindegs(jtype) )
-
-         END DO
-
-         CALL wrk_dealloc( jpi, jpj, zsurfvar )
-         CALL wrk_dealloc( jpi, jpj, zsurfmask )
-#if defined key_fabm
-         CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d )
-#endif
-
-      ENDIF
+      IF ( nproftypes > 0 ) THEN 
+ 
+         DO jtype = 1, nproftypes 
+ 
+            ! Allocate local work arrays 
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  ) 
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 
+            CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam     ) 
+            CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi     ) 
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim )     
+                               
+            ! Defaults which might change 
+            DO jvar = 1, profdataqc(jtype)%nvar 
+               zprofmask(:,:,:,jvar) = tmask(:,:,:) 
+               zglam(:,:,jvar)       = glamt(:,:) 
+               zgphi(:,:,jvar)       = gphit(:,:) 
+               zprofclim(:,:,:,jvar) = 0._wp 
+            END DO 
+ 
+            SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
+ 
+            CASE('prof') 
+               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 
+               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 
+               IF ( ln_output_clim ) THEN           
+                  zprofclim(:,:,:,1) = tclim(:,:,:) 
+                  zprofclim(:,:,:,2) = sclim(:,:,:) 
+               ENDIF 
+                
+            CASE('vel') 
+               zprofvar(:,:,:,1) = un(:,:,:) 
+               zprofvar(:,:,:,2) = vn(:,:,:) 
+               zprofmask(:,:,:,1) = umask(:,:,:) 
+               zprofmask(:,:,:,2) = vmask(:,:,:) 
+               zglam(:,:,1) = glamu(:,:) 
+               zglam(:,:,2) = glamv(:,:) 
+               zgphi(:,:,1) = gphiu(:,:) 
+               zgphi(:,:,2) = gphiv(:,:) 
+ 
+            CASE('plchltot') 
+#if defined key_hadocc 
+               ! Chlorophyll from HadOCC 
+               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
+#elif defined key_medusa 
+               ! Add non-diatom and diatom chlorophyll from MEDUSA 
+               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
+#elif defined key_fabm 
+               ! Add all chlorophyll groups from ERSEM 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
+                  & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
+#else 
+               CALL ctl_stop( ' Trying to run plchltot observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               ! Take the log10 where we can, otherwise exclude 
+               tiny = 1.0e-20 
+               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 
+                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:)) 
+               ELSEWHERE 
+                  zprofvar(:,:,:,:)  = obfillflt 
+                  zprofmask(:,:,:,:) = 0 
+               END WHERE 
+               ! Mask out model below any excluded values, 
+               ! to avoid interpolation issues 
+               DO jvar = 1, profdataqc(jtype)%nvar 
+                 DO jj = 1, jpj 
+                    DO ji = 1, jpi 
+                       depth_loop: DO jk = 1, jpk 
+                          IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 
+                             zprofmask(ji,jj,jk:jpk,jvar) = 0 
+                             EXIT depth_loop 
+                          ENDIF 
+                       END DO depth_loop 
+                    END DO 
+                 END DO 
+              END DO 
+ 
+            CASE('pchltot') 
+#if defined key_hadocc 
+               ! Chlorophyll from HadOCC 
+               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
+#elif defined key_medusa 
+               ! Add non-diatom and diatom chlorophyll from MEDUSA 
+               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
+#elif defined key_fabm 
+               ! Add all chlorophyll groups from ERSEM 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
+                  & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
+#else 
+               CALL ctl_stop( ' Trying to run pchltot observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('pno3') 
+#if defined key_hadocc 
+               ! Dissolved inorganic nitrogen from HadOCC 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 
+#elif defined key_medusa 
+               ! Dissolved inorganic nitrogen from MEDUSA 
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 
+#elif defined key_fabm 
+               ! Nitrate from ERSEM 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 
+#else 
+               CALL ctl_stop( ' Trying to run pno3 observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('psi4') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
+                  &           ' but HadOCC does not simulate silicate' ) 
+#elif defined key_medusa 
+               ! Silicate from MEDUSA 
+               zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 
+#elif defined key_fabm 
+               ! Silicate from ERSEM 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 
+#else 
+               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('ppo4') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
+                  &           ' but HadOCC does not simulate phosphate' ) 
+#elif defined key_medusa 
+               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
+                  &           ' but MEDUSA does not simulate phosphate' ) 
+#elif defined key_fabm 
+               ! Phosphate from ERSEM 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 
+#else 
+               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('pdic') 
+#if defined key_hadocc 
+               ! Dissolved inorganic carbon from HadOCC 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 
+#elif defined key_medusa 
+               ! Dissolved inorganic carbon from MEDUSA 
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 
+#elif defined key_fabm 
+               ! Dissolved inorganic carbon from ERSEM 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 
+#else 
+               CALL ctl_stop( ' Trying to run pdic observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('palk') 
+#if defined key_hadocc 
+               ! Alkalinity from HadOCC 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 
+#elif defined key_medusa 
+               ! Alkalinity from MEDUSA 
+               zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 
+#elif defined key_fabm 
+               ! Alkalinity from ERSEM 
+               zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ta) 
+#else 
+               CALL ctl_stop( ' Trying to run palk observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('pph') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run pph observation operator', & 
+                  &           ' but HadOCC has no pH diagnostic defined' ) 
+#elif defined key_medusa && defined key_roam 
+               ! pH from MEDUSA 
+               zprofvar(:,:,:,1) = f3_pH(:,:,:) 
+#elif defined key_fabm 
+               ! pH from ERSEM 
+               zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ph) 
+#else 
+               CALL ctl_stop( ' Trying to run pph observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('po2') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run po2 observation operator', & 
+                  &           ' but HadOCC does not simulate oxygen' ) 
+#elif defined key_medusa 
+               ! Oxygen from MEDUSA 
+               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 
+#elif defined key_fabm 
+               ! Oxygen from ERSEM 
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 
+#else 
+               CALL ctl_stop( ' Trying to run po2 observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE DEFAULT 
+               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
+ 
+            END SELECT 
+ 
+            DO jvar = 1, profdataqc(jtype)%nvar 
+               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
+                  &               nit000, idaystp, jvar,                   & 
+                  &               zprofvar(:,:,:,jvar),                    & 
+                  &               zprofclim(:,:,:,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     ) 
+            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim  )             
+ 
+         END DO 
+ 
+      ENDIF 
+ 
+      IF ( nsurftypes > 0 ) THEN 
+          
+#if defined key_fabm 
+         CALL wrk_alloc( jpi, jpj, jpk, fabm_3d ) 
+#endif 
+ 
+         DO jtype = 1, nsurftypes 
+ 
+            !Allocate local work arrays 
+            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar  ) 
+            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim )          
+            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask ) 
+            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam     ) 
+            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi     ) 
+ 
+            !Defaults which might be changed 
+            DO jvar = 1, surfdataqc(jtype)%nvar             
+               zsurfmask(:,:,jvar) = tmask(:,:,1) 
+               zsurfclim(:,:,jvar) = 0._wp 
+               zglam(:,:,jvar) = glamt(:,:) 
+               zgphi(:,:,jvar) = gphit(:,:)    
+            END DO              
+            llog10 = .FALSE. 
+ 
+            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
+            CASE('sst') 
+               zsurfvar(:,:,1) = tsn(:,:,1,jp_tem) 
+               IF ( ln_output_clim ) zsurfclim(:,:,1) = tclim(:,:,1) 
+            CASE('sla') 
+               zsurfvar(:,:,1) = sshn(:,:) 
+            CASE('sss') 
+               zsurfvar(:,:,1) = tsn(:,:,1,jp_sal) 
+               IF ( ln_output_clim ) zsurfclim(:,:,1) = sclim(:,:,1)               
+            CASE('ssv') 
+               zsurfvar(:,:,1) = un(:,:,1) 
+               zsurfvar(:,:,2) = vn(:,:,1) 
+               zsurfmask(:,:,1) = umask(:,:,1) 
+               zsurfmask(:,:,2) = vmask(:,:,1)    
+               zglam(:,:,1) = glamu(:,:) 
+               zglam(:,:,2) = glamv(:,:)                
+               zgphi(:,:,1) = gphiu(:,:) 
+               zgphi(:,:,2) = gphiv(:,:)                   
+            CASE('sic') 
+               IF ( kstp == 0 ) THEN 
+                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
+                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
+                        &           'time-step but some obs are valid then.' ) 
+                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
+                        &           ' sea-ice concentration obs will be missed' 
+                  ENDIF 
+               ELSE 
+#if defined key_cice 
+                  zsurfvar(:,:,1) = fr_i(:,:) 
+#elif defined key_lim2 || defined key_lim3 
+                  zsurfvar(:,:,1) = 1._wp - frld(:,:) 
+#else 
+               CALL ctl_stop( ' Trying to run sea-ice concentration observation operator', & 
+                  &           ' but no sea-ice model appears to have been defined' ) 
+#endif 
+               ENDIF 
+ 
+            CASE('sit') 
+               IF ( kstp == 0 ) THEN 
+                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
+                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
+                        &           'time-step but some obs are valid then.' ) 
+                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
+                        &           ' sea-ice thickness obs will be missed and QC flag set to 4' 
+                  ENDIF                                  
+               ELSE        
+#if defined key_cice 
+                  zsurfvar(:,:,1) = thick_i(:,:) 
+#elif defined key_lim2 || defined key_lim3 
+                  CALL ctl_stop( ' No sea-ice thickness observation operator defined for LIM model' ) 
+#else 
+                  CALL ctl_stop( ' Trying to run sea-ice thickness observation operator', & 
+                     &           ' but no sea-ice model appears to have been defined' ) 
+#endif 
+               ENDIF 
+ 
+            CASE('fbd') 
+               IF ( kstp == 0 ) THEN 
+                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
+                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
+                        &           'time-step but some obs are valid then.' ) 
+                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
+                        &           ' sea-ice freeboard obs will be missed and QC flag set to 4' 
+                  ENDIF                                  
+               ELSE        
+#if defined key_cice 
+                  zsurfvar(:,:) = thick_i(:,:) 
+#elif defined key_lim2 || defined key_lim3 
+                  CALL ctl_stop( ' No sea-ice freeboard observation operator defined for LIM model' ) 
+#else 
+                  CALL ctl_stop( ' Trying to run sea-ice freeboard observation operator', & 
+                     &           ' but no sea-ice model appears to have been defined' ) 
+#endif 
+               ENDIF 
+ 
+            CASE('slchltot') 
+#if defined key_hadocc 
+               ! Surface chlorophyll from HadOCC 
+               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1) 
+#elif defined key_medusa 
+               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
+               zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
+#elif defined key_fabm 
+               ! Add all surface chlorophyll groups from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
+                  & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
+#else 
+               CALL ctl_stop( ' Trying to run slchltot observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slchldia') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slchldia observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate diatoms' ) 
+#elif defined key_medusa 
+               ! Diatom surface chlorophyll from MEDUSA 
+               zsurfvar(:,:,1) = trn(:,:,1,jpchd) 
+#elif defined key_fabm 
+               ! Diatom surface chlorophyll from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 
+#else 
+               CALL ctl_stop( ' Trying to run slchldia observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slchlnon') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
+#elif defined key_medusa 
+               ! Non-diatom surface chlorophyll from MEDUSA 
+               zsurfvar(:,:,1) = trn(:,:,1,jpchn) 
+#elif defined key_fabm 
+               ! Add all non-diatom surface chlorophyll groups from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
+                  & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
+#else 
+               CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slchldin') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate dinoflagellates' ) 
+#elif defined key_medusa 
+               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
+                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' ) 
+#elif defined key_fabm 
+               ! Dinoflagellate surface chlorophyll from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
+#else 
+               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slchlmic') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate microphytoplankton' ) 
+#elif defined key_medusa 
+               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
+                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' ) 
+#elif defined key_fabm 
+               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
+#else 
+               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slchlnan') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' ) 
+#elif defined key_medusa 
+               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
+                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 
+#elif defined key_fabm 
+               ! Nanophytoplankton surface chlorophyll from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 
+#else 
+               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slchlpic') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate picophytoplankton' ) 
+#elif defined key_medusa 
+               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
+                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' ) 
+#elif defined key_fabm 
+               ! Picophytoplankton surface chlorophyll from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 
+#else 
+               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('schltot') 
+#if defined key_hadocc 
+               ! Surface chlorophyll from HadOCC 
+               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1) 
+#elif defined key_medusa 
+               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
+               zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
+#elif defined key_fabm 
+               ! Add all surface chlorophyll groups from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
+                  & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
+#else 
+               CALL ctl_stop( ' Trying to run schltot observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('slphytot') 
+#if defined key_hadocc 
+               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_had_phy) * c2n_p 
+#elif defined key_medusa 
+               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 
+               ! multiplied by C:N ratio for each 
+               zsurfvar(:,:,1) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 
+#elif defined key_fabm 
+               ! Add all surface phytoplankton carbon groups from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
+                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
+#else 
+               CALL ctl_stop( ' Trying to run slphytot observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slphydia') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slphydia observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate diatoms' ) 
+#elif defined key_medusa 
+               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
+               zsurfvar(:,:,1) = trn(:,:,1,jpphd) * xthetapd 
+#elif defined key_fabm 
+               ! Diatom surface phytoplankton carbon from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 
+#else 
+               CALL ctl_stop( ' Trying to run slphydia observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('slphynon') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run slphynon observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
+#elif defined key_medusa 
+               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
+               zsurfvar(:,:,1) = trn(:,:,1,jpphn) * xthetapn 
+#elif defined key_fabm 
+               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
+                  & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
+#else 
+               CALL ctl_stop( ' Trying to run slphynon observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+               llog10 = .TRUE. 
+ 
+            CASE('sspm') 
+#if defined key_spm 
+               zsurfvar(:,:,1) = 0.0 
+               DO jn = 1, jp_spm 
+                  zsurfvar(:,:,1) = zsurfvar(:,:,1) + trn(:,:,1,jn)   ! sum SPM sizes 
+               END DO 
+#else 
+               CALL ctl_stop( ' Trying to run sspm observation operator', & 
+                  &           ' but no spm model appears to have been defined' ) 
+#endif 
+ 
+            CASE('skd490') 
+#if defined key_hadocc 
+               CALL ctl_stop( ' Trying to run skd490 observation operator', & 
+                  &           ' but HadOCC does not explicitly simulate Kd490' ) 
+#elif defined key_medusa 
+               CALL ctl_stop( ' Trying to run skd490 observation operator', & 
+                  &           ' but MEDUSA does not explicitly simulate Kd490' ) 
+#elif defined key_fabm 
+               ! light_Kd_band3 diagnostic variable if using spectral optical model 
+               ! light_xEPS diagnostic variable if using standard ERSEM light model 
+               IF ( jp_fabm_kd490 /= -1 ) THEN 
+                  fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_kd490) 
+               ELSEIF ( jp_fabm_xeps /= -1 ) THEN 
+                  fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_xeps) 
+               ELSE 
+                  CALL ctl_stop( ' Trying to run skd490 observation operator', & 
+                     &           ' but cannot access Kd490 from ERSEM' ) 
+               ENDIF 
+               zsurfvar(:,:,1) = fabm_3d(:,:,1) 
+#else 
+               CALL ctl_stop( ' Trying to run skd490 observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('sfco2') 
+#if defined key_hadocc 
+               zsurfvar(:,:,1) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
+               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 
+                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
+                  zsurfvar(:,:) = obfillflt 
+                  zsurfmask(:,:,1) = 0 
+                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
+                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
+               ENDIF 
+#elif defined key_medusa && defined key_roam 
+               zsurfvar(:,:,1) = f2_fco2w(:,:) 
+#elif defined key_fabm 
+               ! First, get pCO2 from FABM 
+               fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc) 
+               zsurfvar(:,:,1) = fabm_3d(:,:,1) 
+               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
+               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
+               ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 
+               ! and 
+               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 
+               ! Marine Chemistry, 2: 203-215. 
+               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 
+               ! not explicitly included - atmospheric pressure is not necessarily available so this is 
+               ! the best assumption. 
+               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 
+               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 
+               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
+               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
+               zsurfvar(:,:,1) = zsurfvar(:,:,1) * EXP((-1636.75 + & 
+                  &              12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
+                  &              0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
+                  &              0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
+                  &              2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
+                  &              (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
+#else 
+               CALL ctl_stop( ' Trying to run sfco2 observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE('spco2') 
+#if defined key_hadocc 
+               zsurfvar(:,:,1) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
+               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 
+                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
+                  zsurfvar(:,:,1) = obfillflt 
+                  zsurfmask(:,:,1) = 0 
+                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
+                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
+               ENDIF 
+#elif defined key_medusa && defined key_roam 
+               zsurfvar(:,:,1) = f2_pco2w(:,:) 
+#elif defined key_fabm 
+               fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc) 
+               zsurfvar(:,:,1) = fabm_3d(:,:,1) 
+#else 
+               CALL ctl_stop( ' Trying to run spco2 observation operator', & 
+                  &           ' but no biogeochemical model appears to have been defined' ) 
+#endif 
+ 
+            CASE DEFAULT 
+ 
+               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 
+ 
+            END SELECT 
+             
+            IF ( llog10 ) THEN 
+               ! Take the log10 where we can, otherwise exclude 
+               tiny = 1.0e-20 
+               WHERE(zsurfvar(:,:,1) > tiny .AND. zsurfvar(:,:,1) /= obfillflt ) 
+                  zsurfvar(:,:,1)  = LOG10(zsurfvar(:,:,1)) 
+               ELSEWHERE 
+                  zsurfvar(:,:,1)  = obfillflt 
+                  zsurfmask(:,:,1) = 0 
+               END WHERE 
+            ENDIF 
+ 
+            DO jvar = 1, surfdataqc(jtype)%nvar 
+ 
+               IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND. & 
+                     &  ln_time_mean_sla_bkg ) THEN 
+                  !Number of time-steps in meaning period 
+                  imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt ) 
+                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
+                     &               nit000, idaystp, jvar,                   & 
+                     &               zsurfvar(:,:,jvar),                      & 
+                     &               zsurfclim(:,:,jvar),                     & 
+                     &               zsurfmask(:,:,jvar),                     & 
+                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                      
+                     &               n2dintsurf(jtype), llnightav(jtype),     & 
+                     &               ravglamscl(jtype), ravgphiscl(jtype),    & 
+                     &               lfpindegs(jtype), kmeanstp = imeanstp ) 
+ 
+               ELSE 
+                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
+                     &               nit000, idaystp, jvar,                   & 
+                     &               zsurfvar(:,:,jvar),                      & 
+                     &               zsurfclim(:,:,jvar),                     & 
+                     &               zsurfmask(:,:,jvar),                     & 
+                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                                           
+                     &               n2dintsurf(jtype), llnightav(jtype),     & 
+                     &               ravglamscl(jtype), ravgphiscl(jtype),    & 
+                     &               lfpindegs(jtype) ) 
+               ENDIF 
+ 
+            END DO 
+ 
+            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar  ) 
+            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim )                   
+            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask ) 
+            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam     ) 
+            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi     )          
+ 
+         END DO 
+#if defined key_fabm 
+         CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d ) 
+#endif 
+ 
+      ENDIF 
 
    END SUBROUTINE dia_obs
@@ -1549,5 +1918,5 @@
                   & )
 
-               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
+               CALL obs_rotvel_pro( profdataqc(jtype), nn_2dint_default, zu, zv )
 
                DO jo = 1, profdataqc(jtype)%nprof
@@ -1582,4 +1951,23 @@
 
          DO jtype = 1, nsurftypes
+
+            IF ( TRIM(cobstypessurf(jtype)) == 'vel' ) THEN 
+               ! For velocity data, rotate the model velocities to N/S, E/W 
+               ! using the compressed data structure. 
+               ALLOCATE( & 
+                  & zu(surfdataqc(jtype)%nsurf), & 
+                  & zv(surfdataqc(jtype)%nsurf)  & 
+                  & ) 
+
+               CALL obs_rotvel_surf( surfdataqc(jtype), nn_2dint_default, zu, zv ) 
+
+               DO jo = 1, surfdataqc(jtype)%nsurf 
+                  surfdataqc(jtype)%rmod(jo,1) = zu(jo) 
+                  surfdataqc(jtype)%rmod(jo,2) = zv(jo) 
+               END DO 
+
+               DEALLOCATE( zu ) 
+               DEALLOCATE( zv ) 
+            END IF 
 
             CALL obs_surf_decompress( surfdataqc(jtype), &
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90	(revision 15764)
@@ -13,5 +13,8 @@
       !! ** Method  : Straightforward search
       !!
-      !! ** Action  : 
+      !! ** Action  : Will return level associated with T-point below the obs 
+      !!              depth, except when observation is in the top box will  
+      !!              return level 2. Also, if obs depth greater than depth  
+      !!              of last wet T-point (kpk-1) will return level kpk.
       !!
       !! History :
@@ -43,6 +46,6 @@
       DO ji = 1, kobs 
          kobsk(ji) = 1
-         depk: DO jk = 2, kgrd
-            IF ( pgrddep(jk) >= pobsdep(ji) ) EXIT depk
+         depk: DO jk = 2, kgrd-1 
+            IF ( pgrddep(jk) > pobsdep(ji) ) EXIT depk
          END DO depk
          kobsk(ji) = jk
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 15764)
@@ -37,4 +37,8 @@
    USE obs_grid,      ONLY : & 
       & obs_level_search     
+#if defined key_cice 
+   USE ice_constants, ONLY : &    ! For conversion from sea ice freeboard to thickness 
+      & rhos, rhoi, rhow 
+#endif
 
    IMPLICIT NONE
@@ -60,9 +64,9 @@
 
 
-   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
-      &                     kit000, kdaystp, kvar,       &
-      &                     pvar, pgdept, pgdepw,        &
-      &                     pmask,                       &  
-      &                     plam, pphi,                  &
+   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 
+      &                     kit000, kdaystp, kvar,       & 
+      &                     pvar, pclim,                 & 
+      &                     pgdept, pgdepw, pmask,       &   
+      &                     plam, pphi,                  & 
       &                     k1dint, k2dint, kdailyavtypes )
 
@@ -137,4 +141,5 @@
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
          & pvar,   &                 ! Model field for variable
+         & pclim,  &                 ! Climatology field for variable
          & pmask                     ! Land-sea mask for variable
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
@@ -172,5 +177,6 @@
       REAL(KIND=wp), DIMENSION(kpk) :: &
          & zobsk,    &
-         & zobs2k
+         & zobs2k,   & 
+         & zclm2k
       REAL(KIND=wp), DIMENSION(2,2,1) :: &
          & zweig1, &
@@ -178,4 +184,5 @@
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
          & zmask,  &
+         & zclim,  &
          & zint,   &
          & zinm,   &
@@ -187,4 +194,5 @@
       REAL(KIND=wp), DIMENSION(1) :: zmsk
       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
+      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim
 
       LOGICAL :: ld_dailyav
@@ -262,4 +270,6 @@
          & )
 
+      IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) )
+
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
          iobs = jobs - prodatqc%nprofup
@@ -285,4 +295,8 @@
       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 ) 
+
+      IF ( prodatqc%lclim ) THEN 
+         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim )             
+      ENDIF
 
       ! At the end of the day also get interpolated means
@@ -349,4 +363,5 @@
                   inum_obs = iend - ista + 1 
                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
+                  IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
 
                   DO iin=1,2 
@@ -358,4 +373,11 @@
                               &     zobs2k, zgdept(iin,ijn,:,iobs), & 
                               &     zmask(iin,ijn,:,iobs)) 
+
+                           IF ( prodatqc%lclim ) THEN 
+                              CALL obs_int_z1d_spl( kpk, &  
+                                 &     zclim(iin,ijn,:,iobs), &  
+                                 &     zclm2k, zgdept(iin,ijn,:,iobs), &  
+                                 &     zmask(iin,ijn,:,iobs))  
+                           ENDIF
                         ENDIF 
        
@@ -372,4 +394,12 @@
                            &    zmask(iin,ijn,:,iobs)) 
        
+                        IF ( prodatqc%lclim ) THEN 
+                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
+                              &    prodatqc%var(kvar)%vdep(ista:iend), &  
+                              &    zclim(iin,ijn,:,iobs), &  
+                              &    zclm2k, interp_corner_clim(iin,ijn,:), &  
+                              &    zgdept(iin,ijn,:,iobs), &  
+                              &    zmask(iin,ijn,:,iobs))  
+                        ENDIF
                      ENDDO 
                   ENDDO 
@@ -386,4 +416,5 @@
                inum_obs = iend - ista + 1 
                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
+               IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
                DO iin=1,2  
                   DO ijn=1,2 
@@ -395,4 +426,10 @@
                            &    zmask(iin,ijn,:,iobs)) 
   
+                        IF ( prodatqc%lclim ) THEN 
+                           CALL obs_int_z1d_spl( kpk, &  
+                              &    zclim(iin,ijn,:,iobs),&  
+                              &    zclm2k, zgdept(iin,ijn,:,iobs), &  
+                              &    zmask(iin,ijn,:,iobs))  
+                        ENDIF
                      ENDIF 
        
@@ -408,5 +445,13 @@
                          &          zgdept(iin,ijn,:,iobs),         & 
                          &          zmask(iin,ijn,:,iobs) )      
-         
+
+                     IF ( prodatqc%lclim ) THEN 
+                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
+                            &          prodatqc%var(kvar)%vdep(ista:iend),     &  
+                            &          zclim(iin,ijn,:,iobs),  &  
+                            &          zclm2k,interp_corner_clim(iin,ijn,:), &  
+                            &          zgdept(iin,ijn,:,iobs), &  
+                            &          zmask(iin,ijn,:,iobs) )    
+                     ENDIF
                   ENDDO 
                ENDDO 
@@ -451,4 +496,9 @@
                   &              prodatqc%var(kvar)%vmod(iend:iend) ) 
 
+               IF ( prodatqc%lclim ) THEN 
+                  CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), &  
+                     &              prodatqc%var(kvar)%vclm(iend:iend) ) 
+               ENDIF
+
                   ! Set QC flag for any observations found below the bottom
                   ! needed as the check here is more strict than that in obs_prep
@@ -458,4 +508,5 @@
  
             DEALLOCATE(interp_corner,iv_indic) 
+            IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim )
           
          ENDIF
@@ -475,4 +526,6 @@
          & )
 
+      IF ( prodatqc%lclim ) DEALLOCATE( zclim )
+
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
@@ -486,8 +539,10 @@
    END SUBROUTINE obs_prof_opt
 
-   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
-      &                     kit000, kdaystp, psurf, psurfmask,   &
-      &                     k2dint, ldnightav, plamscl, pphiscl, &
-      &                     lindegrees )
+   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 
+      &                     kit000, kdaystp, kvar, & 
+      &                     psurf, pclim, psurfmask, & 
+      &                     plam, pphi, & 
+      &                     k2dint, ldnightav, plamscl, pphiscl, & 
+      &                     lindegrees, kmeanstp )
 
       !!-----------------------------------------------------------------------
@@ -522,4 +577,5 @@
       !!      ! 15-02 (M. Martin) Combined routine for surface types
       !!      ! 17-03 (M. Martin) Added horizontal averaging options
+      !!      ! 20-08 (M. Martin) Added surface velocity options
       !!-----------------------------------------------------------------------
 
@@ -538,8 +594,16 @@
                                        !   (kit000-1 = restart time)
       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
+      INTEGER, INTENT(IN) :: kvar      ! Number of variables in surfdataqc
       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
+      INTEGER, INTENT(IN), OPTIONAL :: & 
+                             kmeanstp  ! Number of time steps for the time meaning 
+                                       ! Averaging is triggered if present and greater than one
       REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
          & psurf,  &                   ! Model surface field
+         & pclim,  &                   ! Climatological surface field
          & psurfmask                   ! Land-sea mask
+      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
+         & plam,   &                   ! Model longitudes for variable 
+         & pphi                        ! Model latitudes for variable
       LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
       REAL(KIND=wp), INTENT(IN) :: &
@@ -559,4 +623,5 @@
       INTEGER :: imodi, imodj
       INTEGER :: idayend
+      INTEGER :: imeanend
       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
          & igrdi,   &
@@ -569,6 +634,7 @@
       REAL(wp) :: zlam
       REAL(wp) :: zphi
-      REAL(wp), DIMENSION(1) :: zext, zobsmask
+      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm
       REAL(wp) :: zdaystp
+      REAL(wp) :: zmeanstp
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
          & zweig,  &
@@ -577,4 +643,5 @@
          & zsurfm, &
          & zsurftmp, &
+         & zclim,  &
          & zglam,  &
          & zgphi,  &
@@ -586,4 +653,6 @@
          & zouttmp, &
          & zmeanday    ! to compute model sst in region of 24h daylight (pole)
+
+      LOGICAL :: l_timemean
 
       !------------------------------------------------------------------------
@@ -594,4 +663,9 @@
       isurf = surfdataqc%nsstp(inrc)
 
+      l_timemean = .FALSE. 
+      IF ( PRESENT( kmeanstp ) ) THEN 
+         IF ( kmeanstp > 1 ) l_timemean = .TRUE. 
+      ENDIF
+
       ! Work out the maximum footprint size for the 
       ! interpolation/averaging in model grid-points - has to be even.
@@ -599,4 +673,40 @@
       CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
 
+      IF ( l_timemean ) THEN 
+         ! Initialize time mean for first timestep 
+         imeanend = MOD( kt - kit000 + 1, kmeanstp ) 
+         IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend 
+
+         ! Added kt == 0 test to catch restart case  
+         IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN 
+            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 
+            DO jj = 1, jpj 
+               DO ji = 1, jpi 
+                  surfdataqc%vdmean(ji,jj,kvar) = 0.0 
+               END DO 
+            END DO 
+         ENDIF 
+
+         ! On each time-step, increment the field for computing time mean 
+         IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt 
+         DO jj = 1, jpj 
+            DO ji = 1, jpi 
+               surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) & 
+                  &                            + psurf(ji,jj) 
+            END DO 
+         END DO 
+
+         ! Compute the time mean at the end of time period 
+         IF ( imeanend == 0 ) THEN 
+            zmeanstp = 1.0 / REAL( kmeanstp ) 
+            IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp 
+            DO jj = 1, jpj 
+               DO ji = 1, jpi 
+                  surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) & 
+                     &                            * zmeanstp 
+               END DO 
+            END DO 
+         ENDIF 
+      ENDIF !l_timemean
 
       IF ( ldnightav ) THEN
@@ -665,17 +775,24 @@
 
       ALLOCATE( &
-         & zweig(imaxifp,imaxjfp,1),      &
-         & igrdi(imaxifp,imaxjfp,isurf), &
-         & igrdj(imaxifp,imaxjfp,isurf), &
-         & zglam(imaxifp,imaxjfp,isurf), &
-         & zgphi(imaxifp,imaxjfp,isurf), &
-         & zmask(imaxifp,imaxjfp,isurf), &
-         & zsurf(imaxifp,imaxjfp,isurf), &
-         & zsurftmp(imaxifp,imaxjfp,isurf),  &
-         & zglamf(imaxifp+1,imaxjfp+1,isurf), &
-         & zgphif(imaxifp+1,imaxjfp+1,isurf), &
-         & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
-         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
-         & )
+         & zweig(imaxifp,imaxjfp,1),       & 
+         & igrdi(imaxifp,imaxjfp,isurf),   & 
+         & igrdj(imaxifp,imaxjfp,isurf),   & 
+         & zglam(imaxifp,imaxjfp,isurf),   & 
+         & zgphi(imaxifp,imaxjfp,isurf),   & 
+         & zmask(imaxifp,imaxjfp,isurf),   & 
+         & zsurf(imaxifp,imaxjfp,isurf),   & 
+         & zsurftmp(imaxifp,imaxjfp,isurf) & 
+         & ) 
+
+      IF ( k2dint > 4 ) THEN 
+         ALLOCATE( & 
+            & zglamf(imaxifp+1,imaxjfp+1,isurf),  & 
+            & zgphif(imaxifp+1,imaxjfp+1,isurf),  & 
+            & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
+            & igrdjp1(imaxifp+1,imaxjfp+1,isurf)  & 
+            & ) 
+      ENDIF 
+      
+      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) ) 
 
       DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
@@ -695,6 +812,8 @@
                IF ( imodj > jpjglo ) imodj = jpjglo
 
-               igrdip1(ji+1,jj+1,iobs) = imodi
-               igrdjp1(ji+1,jj+1,iobs) = imodj
+               IF ( k2dint > 4 ) THEN 
+                  igrdip1(ji+1,jj+1,iobs) = imodi 
+                  igrdjp1(ji+1,jj+1,iobs) = imodj 
+               ENDIF
                
                IF ( ji >= 1 .AND. jj >= 1 ) THEN
@@ -713,10 +832,29 @@
       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, psurfmask, zmask )
-      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
-         &                  igrdi, igrdj, psurf, zsurf )
-      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
-         &                  igrdip1, igrdjp1, glamf, zglamf )
-      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
-         &                  igrdip1, igrdjp1, gphif, zgphif )
+
+      ! At the end of the averaging period get interpolated means 
+      IF ( l_timemean ) THEN 
+         IF ( imeanend == 0 ) THEN 
+            ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) ) 
+            IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt 
+            CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
+               &                  igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm ) 
+         ENDIF 
+      ELSE 
+         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
+            &                  igrdi, igrdj, psurf, zsurf ) 
+      ENDIF 
+
+      IF ( k2dint > 4 ) THEN          
+         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
+            &                  igrdip1, igrdjp1, glamf, zglamf ) 
+         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
+            &                  igrdip1, igrdjp1, gphif, zgphif ) 
+      ENDIF
+
+      IF ( surfdataqc%lclim ) THEN  
+         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
+            &                  igrdi, igrdj, pclim, zclim ) 
+      ENDIF
 
       ! At the end of the day get interpolated means
@@ -758,48 +896,85 @@
          zphi = surfdataqc%rphi(jobs)
 
-         IF ( ldnightav .AND. idayend == 0 ) THEN
-            ! Night-time averaged data
-            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
-         ELSE
-            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
+         IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN 
+            ! Night-time or N=kmeanstp timestep averaged data 
+            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
+         ELSE 
+            zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
+         ENDIF 
+
+         IF (   ( .NOT. l_timemean ) .OR. &  
+             &  ( l_timemean .AND. imeanend == 0) ) THEN 
+            IF ( k2dint <= 4 ) THEN 
+
+               ! Get weights to interpolate the model value to the observation point 
+               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
+                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
+                  &                   zmask(:,:,iobs), zweig, zobsmask ) 
+
+               ! Interpolate the model value to the observation point  
+               CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
+
+               IF ( surfdataqc%lclim ) THEN   
+                  CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 
+               ENDIF 
+
+
+            ELSE 
+
+               ! Get weights to average the model SLA to the observation footprint 
+               CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
+                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
+                  &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
+                  &                   zmask(:,:,iobs), plamscl, pphiscl, & 
+                  &                   lindegrees, zweig, zobsmask ) 
+
+               ! Average the model SST to the observation footprint 
+               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
+                  &              zweig, zsurftmp(:,:,iobs), zext ) 
+
+               IF ( surfdataqc%lclim ) THEN   
+                  CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
+                     &              zweig, zclim(:,:,iobs), zclm ) 
+               ENDIF 
+
+            ENDIF 
+
+            IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
+               ! ... Remove the MDT from the SSH at the observation point to get the SLA 
+               surfdataqc%rext(jobs,1) = zext(1) 
+               surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
+            ELSE 
+               surfdataqc%rmod(jobs,kvar) = zext(1) 
+            ENDIF 
+
+            IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 
+
+            IF ( zext(1) == obfillflt ) THEN 
+               ! If the observation value is a fill value, set QC flag to bad 
+               surfdataqc%nqc(jobs) = 4 
+            ENDIF          
+         ENDIF 
+
+         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 
+ 
+#if defined key_cice 
+         IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN 
+            ! Convert radar freeboard to true freeboard (add 1/4 snow depth; 1/4 based on ratio of speed of light in vacuum compared to snow (3.0e8 vs 2.4e8 m/s)) 
+            surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1)  
+            surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2) 
+            ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad 
+            IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN 
+               surfdataqc%nqc(jobs) = 4 
+            ENDIF            
+            ! Convert corrected freeboard to ice thickness following Tilling et al. (2016) 
+            surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)*rhow + surfdataqc%rext(jobs,2)*rhos)/ & 
+                                      (rhow - rhoi) 
          ENDIF
 
-         IF ( k2dint <= 4 ) THEN
-
-            ! Get weights to interpolate the model value to the observation point
-            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
-               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
-               &                   zmask(:,:,iobs), zweig, zobsmask )
-
-            ! Interpolate the model value to the observation point 
-            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
-
-         ELSE
-
-            ! Get weights to average the model SLA to the observation footprint
-            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
-               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
-               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
-               &                   zmask(:,:,iobs), plamscl, pphiscl, &
-               &                   lindegrees, zweig, zobsmask )
-
-            ! Average the model SST to the observation footprint
-            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
-               &              zweig, zsurftmp(:,:,iobs),  zext )
-
-         ENDIF
-
-         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
-            ! ... Remove the MDT from the SSH at the observation point to get the SLA
-            surfdataqc%rext(jobs,1) = zext(1)
-            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
-         ELSE
-            surfdataqc%rmod(jobs,1) = zext(1)
-         ENDIF
-         
-         IF ( zext(1) == obfillflt ) THEN
-            ! If the observation value is a fill value, set QC flag to bad
-            surfdataqc%nqc(jobs) = 4
-         ENDIF
+         IF ( zext(1) == obfillflt ) THEN 
+            ! If the observation value is a fill value, set QC flag to bad 
+            surfdataqc%nqc(jobs) = 4 
+         ENDIF 
+#endif defined key_cice
 
       END DO
@@ -814,19 +989,28 @@
          & zmask, &
          & zsurf, &
-         & zsurftmp, &
-         & zglamf, &
-         & zgphif, &
-         & igrdip1,&
-         & igrdjp1 &
+         & zsurftmp &
          & )
 
-      ! At the end of the day also deallocate night-time mean array
-      IF ( idayend == 0 .AND. ldnightav ) THEN
-         DEALLOCATE( &
-            & zsurfm  &
+      IF ( k2dint > 4 ) THEN 
+         DEALLOCATE( &      
+            & zglamf, & 
+            & zgphif, & 
+            & igrdip1,& 
+            & igrdjp1 & 
             & )
       ENDIF
 
-      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
+      IF ( surfdataqc%lclim ) DEALLOCATE( zclim ) 
+
+      ! At the end of the day also deallocate night-time mean array 
+      IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN
+         DEALLOCATE( & 
+            & zsurfm  & 
+            & ) 
+      ENDIF
+
+      IF ( kvar == surfdataqc%nvar ) THEN 
+         surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
+      ENDIF
 
    END SUBROUTINE obs_surf_opt
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 15764)
@@ -16,4 +16,5 @@
    USE par_kind, ONLY : & ! Precision variables
       & wp   
+   USE dom_oce            ! ocean space and time domain
    USE in_out_manager     ! I/O manager
    USE obs_profiles_def   ! Definitions for storage arrays for profiles
@@ -52,6 +53,9 @@
 CONTAINS
 
-   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, &
-                            kqc_cutoff )
+   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, & 
+      &                     kpi, kpj, &    
+      &                     zmask, pglam, pgphi, & 
+      &                     ld_nea, ld_bound_reject, & 
+      &                     ld_seaicetypes, kqc_cutoff )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_pre_sla  ***
@@ -81,6 +85,13 @@
       TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data
       TYPE(obs_surf), INTENT(INOUT) :: surfdataqc  ! Subset of surface data not failing screening
+      INTEGER, INTENT(IN) :: kpi, kpj              ! Local domain sizes       
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: & 
+         & zmask       
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: & 
+         & pglam, & 
+         & pgphi
       LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land
       LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary
+      LOGICAL, INTENT(IN) :: ld_seaicetypes        ! Switch to indicate sea ice data
       INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
       !! * Local declarations
@@ -93,22 +104,25 @@
       INTEGER :: icycle       ! Current assimilation cycle
                               ! Counters for observations that
-      INTEGER :: iotdobs      !  - outside time domain
-      INTEGER :: iosdsobs     !  - outside space domain
-      INTEGER :: ilansobs     !  - within a model land cell
-      INTEGER :: inlasobs     !  - close to land
-      INTEGER :: igrdobs      !  - fail the grid search
-      INTEGER :: ibdysobs     !  - close to open boundary
-                              ! Global counters for observations that
-      INTEGER :: iotdobsmpp     !  - outside time domain
-      INTEGER :: iosdsobsmpp    !  - outside space domain
-      INTEGER :: ilansobsmpp    !  - within a model land cell
-      INTEGER :: inlasobsmpp    !  - close to land
-      INTEGER :: igrdobsmpp     !  - fail the grid search
-      INTEGER :: ibdysobsmpp  !  - close to open boundary
-      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
-         & llvalid            ! SLA data selection
-      INTEGER :: jobs         ! Obs. loop variable
-      INTEGER :: jstp         ! Time loop variable
-      INTEGER :: inrc         ! Time index variable
+      INTEGER                           :: iotdobs      !  - outside time domain 
+      INTEGER, DIMENSION(surfdata%nvar) :: iosdsobs     !  - outside space domain 
+      INTEGER, DIMENSION(surfdata%nvar) :: ilansobs     !  - within a model land cell 
+      INTEGER, DIMENSION(surfdata%nvar) :: inlasobs     !  - close to land 
+      INTEGER, DIMENSION(surfdata%nvar) :: ibdysobs     !  - close to open boundary 
+      INTEGER                           :: igrdobs      !  - fail the grid search       
+                                                        ! Global counters for observations that 
+      INTEGER                           :: iotdobsmpp   !  - outside time domain 
+      INTEGER, DIMENSION(surfdata%nvar) :: iosdsobsmpp  !  - outside space domain 
+      INTEGER, DIMENSION(surfdata%nvar) :: ilansobsmpp  !  - within a model land cell 
+      INTEGER, DIMENSION(surfdata%nvar) :: inlasobsmpp  !  - close to land 
+      INTEGER, DIMENSION(surfdata%nvar) :: ibdysobsmpp  !  - close to open boundary 
+      INTEGER                           :: igrdobsmpp   !  - fail the grid search 
+
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
+         & llvalid            ! SLA data selection 
+      INTEGER :: jobs         ! Obs. loop counter 
+      INTEGER :: jvar         ! Variable loop counter 
+      INTEGER :: jstp         ! Time loop variable 
+      INTEGER :: inrc         ! Time index variable 
+      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line 
 
       IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...'
@@ -140,18 +154,28 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_tim( icycle, &
-         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-         &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, &
-         &              surfdata%nday,    surfdata%nhou, surfdata%nmin, &
-         &              surfdata%nqc,     surfdata%mstp, iotdobs        )
-
-      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
-      
-      ! -----------------------------------------------------------------------
-      ! Check for surface data failing the grid search
-      ! -----------------------------------------------------------------------
-
-      CALL obs_coo_grd( surfdata%nsurf,   surfdata%mi, surfdata%mj, &
-         &              surfdata%nqc,     igrdobs                         )
+      IF ( ld_seaicetypes ) THEN 
+         CALL obs_coo_tim( icycle, & 
+            &              iyea0,   imon0,   iday0,   ihou0, imin0,      & 
+            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
+            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
+            &              surfdata%nqc,     surfdata%mstp, iotdobs,       & 
+            &              ld_seaicetypes = ld_seaicetypes ) 
+      ELSE 
+         CALL obs_coo_tim( icycle, & 
+            &              iyea0,   imon0,   iday0,   ihou0, imin0,      & 
+            &              surfdata%nsurf,   surfdata%nyea, surfdata%nmon, & 
+            &              surfdata%nday,    surfdata%nhou, surfdata%nmin, & 
+            &              surfdata%nqc,     surfdata%mstp, iotdobs        ) 
+      ENDIF 
+
+      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
+
+      ! ----------------------------------------------------------------------- 
+      ! Check for surface data failing the grid search 
+      ! ----------------------------------------------------------------------- 
+      DO jvar = 1, surfdata%nvar 
+         CALL obs_coo_grd( surfdata%nsurf, surfdata%mi(:,jvar), surfdata%mj(:,jvar), & 
+            &              surfdata%nqc,     igrdobs ) 
+      END DO 
 
       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
@@ -161,19 +185,20 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_spc_2d( surfdata%nsurf,              &
-         &                 jpi,          jpj,          &
-         &                 surfdata%mi,   surfdata%mj,   & 
-         &                 surfdata%rlam, surfdata%rphi, &
-         &                 glamt,        gphit,        &
-         &                 tmask(:,:,1), surfdata%nqc,  &
-         &                 iosdsobs,     ilansobs,     &
-         &                 inlasobs,     ld_nea,       &
-         &                 ibdysobs,     ld_bound_reject, &
-         &                 iqc_cutoff                     )
-
-      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
-      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
-      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
-      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
+      DO jvar = 1, surfdata%nvar 
+         CALL obs_coo_spc_2d( surfdata%nsurf,                     & 
+            &                 jpi,          jpj,                  & 
+            &                 surfdata%mi(:,jvar), surfdata%mj(:,jvar), &  
+            &                 surfdata%rlam, surfdata%rphi,       & 
+            &                 pglam(:,:,jvar), pgphi(:,:,jvar),   & 
+            &                 zmask(:,:,jvar), surfdata%nqc(:),   & 
+            &                 iosdsobs(jvar), ilansobs(jvar),     & 
+            &                 inlasobs(jvar),     ld_nea,         & 
+            &                 ibdysobs(jvar), ld_bound_reject,    & 
+            &                 iqc_cutoff                     ) 
+         CALL obs_mpp_sum_integer( iosdsobs(jvar), iosdsobsmpp(jvar) ) 
+         CALL obs_mpp_sum_integer( ilansobs(jvar), ilansobsmpp(jvar) ) 
+         CALL obs_mpp_sum_integer( inlasobs(jvar), inlasobsmpp(jvar) ) 
+         CALL obs_mpp_sum_integer( ibdysobs(jvar), ibdysobsmpp(jvar) ) 
+      END DO 
 
       ! -----------------------------------------------------------------------
@@ -204,32 +229,44 @@
       ! Update the total observation counter array
       
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain                  = ', &
-            &            iotdobsmpp
-         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search    = ', &
-            &            igrdobsmpp
-         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain       = ', &
-            &            iosdsobsmpp
-         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points             = ', &
-            &            ilansobsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', &
-               &            inlasobsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept)    = ', &
-               &            inlasobsmpp
-         ENDIF
-         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', &
-            &            ibdysobsmpp  
-         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', &
-            &            surfdataqc%nsurfmpp
-
-         WRITE(numout,*)
-         WRITE(numout,*) ' Number of observations per time step :'
-         WRITE(numout,*)
-         WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1)
-         WRITE(numout,'(10X,A,5X,A)')'---------','-----------------'
-         CALL FLUSH(numout)
+      IF(lwp) THEN 
+
+         DO jvar = 1, surfdataqc%nvar        
+            IF ( jvar == 1 ) THEN 
+               cout1=TRIM(surfdataqc%cvars(1))                   
+            ELSE 
+               WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdataqc%cvars(jvar))             
+            ENDIF 
+         END DO 
+
+         WRITE(numout,*) 
+         WRITE(numout,*) ' '//TRIM(cout1)//' data outside time domain                  = ', & 
+            &            iotdobsmpp 
+         WRITE(numout,*) ' Remaining '//TRIM(cout1)//' data that failed grid search    = ', & 
+            &            igrdobsmpp 
+
+         DO jvar = 1, surfdataqc%nvar             
+            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data outside space domain       = ', & 
+                &            iosdsobsmpp(jvar) 
+             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data at land points             = ', & 
+                &            ilansobsmpp(jvar) 
+             IF (ld_nea) THEN 
+                WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (removed) = ', & 
+                   &            inlasobsmpp(jvar) 
+             ELSE 
+                WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (kept)    = ', & 
+                   &            inlasobsmpp(jvar) 
+             ENDIF      
+             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near open boundary (removed) = ', & 
+                &            ibdysobsmpp(jvar) 
+          END DO 
+          WRITE(numout,*) ' '//TRIM(cout1)//' data accepted = ', & 
+             &            surfdataqc%nsurfmpp 
+
+         WRITE(numout,*) 
+         WRITE(numout,*) ' Number of observations per time step:' 
+         WRITE(numout,*) 
+         WRITE(numout,'(10X,A,10X,A)')'Time step',TRIM(cout1) 
+         WRITE(numout,'(10X,A,5X,A)')'---------','-----------------' 
+         CALL FLUSH(numout) 
       ENDIF
       
@@ -434,5 +471,5 @@
 
       IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
-         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff )
+         CALL obs_uv_rej_pro( profdata, iuvchku, iuvchkv, iqc_cutoff )
          CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
          CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
@@ -558,5 +595,5 @@
       &                    kobsno,                                        &
       &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
-      &                    kobsqc,  kobsstp, kotdobs                      )
+      &                    kobsqc,  kobsstp, kotdobs, ld_seaicetypes      )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_tim ***
@@ -606,4 +643,5 @@
          & kobsstp          ! Number of time steps up to the 
                             ! observation time
+      LOGICAL, OPTIONAL, INTENT(IN) :: ld_seaicetypes
 
       !! * Local declarations
@@ -620,4 +658,5 @@
       INTEGER :: iskip
       INTEGER :: idaystp
+      INTEGER :: icecount
       REAL(KIND=wp) :: zminstp
       REAL(KIND=wp) :: zhoustp
@@ -713,4 +752,17 @@
             kotdobs = kotdobs + 1
             CYCLE
+         ENDIF
+
+       ! Flag sea ice observations falling on initial timestep 
+         IF ( PRESENT(ld_seaicetypes) ) THEN 
+
+              IF ( ( kobsstp(jobs) == (nit000 - 1) ) ) THEN 
+                 IF (lwp) WRITE(numout,*)( 'Sea-ice not initialised on zeroth '// & 
+                           &    'time-step but SIT observation valid then, flagging '// & 
+                                'in time check subroutine obs_coo_tim.' ) 
+                 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
+                 kotdobs      = kotdobs + 1 
+                 CYCLE 
+              ENDIF 
          ENDIF
 
@@ -1159,8 +1211,10 @@
 #endif 
       REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
+         & zgdept, &
          & zgdepw         
       REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
          & zglam, &           ! Model longitude at grid points
-         & zgphi              ! Model latitude at grid points
+         & zgphi, &           ! Model latitude at grid points 
+         & zbathy             ! Index of deepest wet level at grid points
       INTEGER, DIMENSION(2,2,kprofno) :: &
          & igrdi, &           ! Grid i,j
@@ -1170,4 +1224,5 @@
       INTEGER :: iig, ijg           ! i,j of observation on model grid point.
       INTEGER :: jobs, jobsp, jk, ji, jj
+      REAL(KIND=wp) :: maxdepw
 
       ! Get grid point indices
@@ -1220,7 +1275,10 @@
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
       CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
+      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, REAL(mbathy), zbathy )
       ! Need to know the bathy depth for each observation for sco
       CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), &
          &                  zgdepw )
+      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdept(:,:,:), & 
+         &                  zgdept )
 
       DO jobs = 1, kprofno
@@ -1258,4 +1316,14 @@
          DO jobsp = kpstart(jobs), kpend(jobs)
 
+            ! Calculate max T and W depths of 2x2 grid 
+            maxdepw=zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs) 
+            DO jj = 1, 2 
+               DO ji = 1, 2 
+                  IF ( zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) > maxdepw ) THEN 
+                     maxdepw = zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) 
+                  END IF 
+               END DO 
+            END DO
+
             ! Flag if the observation falls outside the model spatial domain
             IF (       ( pobslam(jobs) < -180.         )       &
@@ -1264,5 +1332,5 @@
                &  .OR. ( pobsphi(jobs) >   90.         )       &
                &  .OR. ( pobsdep(jobsp) < 0.0          )       &
-               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
+               &  .OR. ( pobsdep(jobsp) >= maxdepw ) ) THEN
                kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
                kosdobs = kosdobs + 1
@@ -1327,9 +1395,4 @@
             ENDIF
             
-            ! Set observation depth equal to that of the first model depth
-            IF ( pobsdep(jobsp) <= pdep(1) ) THEN
-               pobsdep(jobsp) = pdep(1)  
-            ENDIF
-            
 #if defined key_bdy
             ! Flag if the observation falls close to the boundary rim
@@ -1405,7 +1468,7 @@
    END SUBROUTINE obs_pro_rej
 
-   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE obs_uv_rej ***
+   SUBROUTINE obs_uv_rej_pro( profdata, knumu, knumv, kqc_cutoff )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE obs_uv_rej_pro ***
       !!
       !! ** Purpose : Reject u if v is rejected and vice versa
@@ -1439,5 +1502,5 @@
             & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
 
-            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
+            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej_pro')
             RETURN
 
@@ -1461,5 +1524,5 @@
       END DO
 
-   END SUBROUTINE obs_uv_rej
+   END SUBROUTINE obs_uv_rej_pro
 
 END MODULE obs_prep
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90	(revision 15764)
@@ -72,5 +72,6 @@
          & vdep,  &       !: Depth coordinate of profile data
          & vobs,  &       !: Profile data
-         & vmod           !: Model counterpart of the profile data vector
+         & vmod,  &       !: Model counterpart of the profile data vector 
+         & vclm           !: Climatological counterpart of the profile data vector
 
       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
@@ -101,4 +102,6 @@
       INTEGER :: npk
       INTEGER :: nprofup  !: Observation counter used in obs_oper
+
+      LOGICAL :: lclim    !: Climatology will be calculated for this structure
 
       ! Bookkeeping arrays with sizes equal to number of variables
@@ -198,5 +201,5 @@
    
    SUBROUTINE obs_prof_alloc( prof,  kvar, kext, kprof,  &
-      &                       ko3dt, kstp, kpi, kpj, kpk )
+      &                       ko3dt, kstp, kpi, kpj, kpk, ldclim )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE obs_prof_alloc  ***
@@ -221,4 +224,5 @@
       INTEGER, INTENT(IN) :: kpj
       INTEGER, INTENT(IN) :: kpk
+      LOGICAL, INTENT(IN) :: ldclim
 
       !!* Local variables
@@ -236,4 +240,6 @@
       prof%npj       = kpj
       prof%npk       = kpk
+
+      prof%lclim     = ldclim
 
       ! Allocate arrays of size number of variables
@@ -503,4 +509,9 @@
             & )
       ENDIF
+      IF (prof%lclim) THEN 
+         ALLOCATE( &  
+            & prof%var(kvar)%vclm(kobs) & 
+            & ) 
+      ENDIF
 
    END SUBROUTINE obs_prof_alloc_var
@@ -537,4 +548,9 @@
          DEALLOCATE( & 
             & prof%var(kvar)%vext  &
+            & ) 
+      ENDIF 
+      IF (prof%lclim) THEN 
+         DEALLOCATE( &  
+            & prof%var(kvar)%vclm  &
             & )
       ENDIF
@@ -630,5 +646,6 @@
             &                 inprof,    invpro,    &
             &                 prof%nstp, prof%npi,  &
-            &                 prof%npj,  prof%npk )
+            &                 prof%npj,  prof%npk,  & 
+            &                 prof%lclim )
       ENDIF
 
@@ -745,4 +762,8 @@
                            &                      prof%var(jvar)%vext(jj,jext)
                      END DO
+                     IF (newprof%lclim) THEN 
+                        newprof%var(jvar)%vclm(invpro(jvar))   = & 
+                           &                           prof%var(jvar)%vclm(jj) 
+                     ENDIF
                   
                      ! nvind is the index of the original variable data
@@ -870,4 +891,7 @@
                      &                        prof%var(jvar)%vext(jj,jext)
                END DO
+               IF (prof%lclim) THEN 
+                  oldprof%var(jvar)%vclm(jl)   = prof%var(jvar)%vclm(jj) 
+               ENDIF
                
             END DO
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 15764)
@@ -46,5 +46,5 @@
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
       &                     ldvar, ldignmis, ldsatt, &
-      &                     ldmod, kdailyavtypes )
+      &                     ldmod, ldclim, cdvars, kdailyavtypes )
       !!---------------------------------------------------------------------
       !!
@@ -78,6 +78,8 @@
       LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points
       LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
+      LOGICAL, INTENT(IN) :: ldclim     ! Set flag to show climatology will be output
       REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
       REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
+      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes                ! Types of daily average observations
@@ -86,5 +88,5 @@
       CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
       CHARACTER(len=8) :: clrefdate
-      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
       INTEGER :: jvar
       INTEGER :: ji
@@ -220,16 +222,20 @@
             ENDIF
 
-            IF ( jj == 1 ) THEN
-               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
-               DO ji = 1, inpfiles(jj)%nvar
-                 clvars(ji) = inpfiles(jj)%cname(ji)
-               END DO
+            IF ( jj == 1 ) THEN 
+               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) ) 
+               DO ji = 1, inpfiles(jj)%nvar 
+                 clvarsin(ji) = inpfiles(jj)%cname(ji) 
+                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN 
+                    CALL ctl_stop( 'Feedback file variables do not match', & 
+                        &           ' expected variable names for this type' ) 
+                 ENDIF 
+               END DO 
             ELSE
-               DO ji = 1, inpfiles(jj)%nvar
-                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
-                     CALL ctl_stop( 'Feedback file variables not consistent', &
-                        &           ' with previous files for this type' )
-                  ENDIF
-               END DO
+               DO ji = 1, inpfiles(jj)%nvar 
+                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN 
+                     CALL ctl_stop( 'Feedback file variables not consistent', & 
+                        &           ' with previous files for this type' ) 
+                  ENDIF 
+               END DO 
             ENDIF
 
@@ -500,5 +506,5 @@
       ENDIF
       CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
-         &                 kstp, jpi, jpj, jpk )
+         &                 kstp, jpi, jpj, jpk, ldclim )
 
       ! * Read obs/positions, QC, all variable and assign to profdata
@@ -506,5 +512,5 @@
       profdata%nprof     = 0
       profdata%nvprot(:) = 0
-      profdata%cvars(:)  = clvars(:)
+      profdata%cvars(:)  = clvarsin(:)
       iprof = 0
 
@@ -696,4 +702,7 @@
                               profdata%var(jvar)%vmod(ivart(jvar)) = &
                                  &                inpfiles(jj)%padd(ij,ji,1,jvar)
+                           ENDIF
+                           IF ( profdata%lclim ) THEN 
+                              profdata%var(jvar)%vclm(ivart(jvar)) = fbrmdi 
                            ENDIF
                            ! Count number of profile var1 data as function of type
@@ -805,5 +814,5 @@
       ! Deallocate temporary data
       !-----------------------------------------------------------------------
-      DEALLOCATE( ifileidx, iprofidx, zdat, clvars )
+      DEALLOCATE( ifileidx, iprofidx, zdat, clvarsin )
 
       !-----------------------------------------------------------------------
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 15764)
@@ -39,6 +39,6 @@
 
    SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
-      &                     kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                     ldignmis, ldmod, ldnightav )
+      &                     kvars, kextr, kstp, ddobsini, ddobsend, MeanPeriodHours, &
+      &                     ldignmis, ldmod, ldnightav, ldclim, ln_time_mean_sla_bkg, cdvars )
       !!---------------------------------------------------------------------
       !!
@@ -71,14 +71,19 @@
       LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
       LOGICAL, INTENT(IN) :: ldnightav  ! Observations represent a night-time average
+      LOGICAL, INTENT(IN) :: ldclim     ! Will include climatology at obs points.
+      LOGICAL, INTENT(IN) :: ln_time_mean_sla_bkg     ! Will reset times to end of averaging period.
       REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
       REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
+      REAL(wp), INTENT(IN) :: MeanPeriodHours ! Averaging period in hours
+      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
 
       !! * Local declarations
       CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
       CHARACTER(len=8) :: clrefdate
-      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
       INTEGER :: ji
       INTEGER :: jj
       INTEGER :: jk
+      INTEGER :: jvar
       INTEGER :: iflag
       INTEGER :: inobf
@@ -103,8 +108,9 @@
          & ityp, &
          & itypmpp
-      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
          & iobsi,    &
          & iobsj,    &
-         & iproc,    &
+         & iproc
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &         
          & iindx,    &
          & ifileidx, &
@@ -122,5 +128,6 @@
       TYPE(obfbdata), POINTER, DIMENSION(:) :: &
          & inpfiles
-
+      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line
+      
       ! Local initialization
       iobs  = 0
@@ -181,4 +188,9 @@
                &                ldgrid = .TRUE. )
 
+            IF ( inpfiles(jj)%nvar /= kvars ) THEN
+               CALL ctl_stop( 'Feedback format error: ', &
+                  &           ' unexpected number of vars in feedback file' )
+            ENDIF
+
             IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
                CALL ctl_stop( 'Model not in input data' )
@@ -187,11 +199,15 @@
 
             IF ( jj == 1 ) THEN
-               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
+               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
                DO ji = 1, inpfiles(jj)%nvar
-                 clvars(ji) = inpfiles(jj)%cname(ji)
+                 clvarsin(ji) = inpfiles(jj)%cname(ji)
+                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN
+                    CALL ctl_stop( 'Feedback file variables do not match', &
+                        &           ' expected variable names for this type' )
+                 ENDIF
                END DO
             ELSE
                DO ji = 1, inpfiles(jj)%nvar
-                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
+                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
                      CALL ctl_stop( 'Feedback file variables not consistent', &
                         &           ' with previous files for this type' )
@@ -265,8 +281,22 @@
 
             IF ( inpfiles(jj)%nobs > 0 ) THEN
-               inpfiles(jj)%iproc = -1
-               inpfiles(jj)%iobsi = -1
-               inpfiles(jj)%iobsj = -1
-            ENDIF
+               inpfiles(jj)%iproc(:,:) = -1
+               inpfiles(jj)%iobsi(:,:) = -1
+               inpfiles(jj)%iobsj(:,:) = -1
+            ENDIF
+
+            !If SLA observations are representing a time mean then set the time
+            !of the obs to the end of that meaning period relative to the start of the run
+            IF ( ln_time_mean_sla_bkg .AND. ( TRIM( clvarsin(1) ) == 'SLA' ) ) THEN
+               DO ji = 1, inpfiles(jj)%nobs
+                  ! Only do this for obs within time window
+                  IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
+                     & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
+	                  inpfiles(jj)%ptim(ji) = &
+	                        & djulini(jj) + (MeanPeriodHours/24.)
+                  ENDIF      
+	            END DO
+            ENDIF   
+            
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
@@ -278,7 +308,7 @@
             ALLOCATE( zlam(inowin)  )
             ALLOCATE( zphi(inowin)  )
-            ALLOCATE( iobsi(inowin) )
-            ALLOCATE( iobsj(inowin) )
-            ALLOCATE( iproc(inowin) )
+            ALLOCATE( iobsi(inowin,kvars) )
+            ALLOCATE( iobsj(inowin,kvars) )
+            ALLOCATE( iproc(inowin,kvars) )
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
@@ -291,5 +321,21 @@
             END DO
 
-            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
+            ! 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
 
             inowin = 0
@@ -298,7 +344,9 @@
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
                   inowin = inowin + 1
-                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
-                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
-                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
+                  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
                ENDIF
             END DO
@@ -359,13 +407,17 @@
          &               iindx   )
 
-      CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj )
+      CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj, ldclim )
 
       ! Read obs/positions, QC, all variable and assign to surfdata
 
       iobs = 0
-      surfdata%cvars(:)  = clvars(:)
+      surfdata%cvars(:)  = clvarsin(:)
       IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'SLA' ) ) THEN
          surfdata%cext(1) = 'SSH'
          surfdata%cext(2) = 'MDT'
+      ENDIF
+      IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'FBD' ) ) THEN
+           surfdata%cext(1) = 'freeboard'
+           surfdata%cext(2) = 'thick_s'
       ENDIF
       IF ( iextr > kextr ) surfdata%cext(iextr) = 'STD'
@@ -417,7 +469,9 @@
 
                ! Coordinate search parameters
-               surfdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
-               surfdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
-
+               DO jvar = 1, kvars
+                  surfdata%mi(iobs,jvar) = inpfiles(jj)%iobsi(ji,jvar)
+                  surfdata%mj(iobs,jvar) = inpfiles(jj)%iobsj(ji,jvar)
+               END DO
+               
                ! WMO number
                surfdata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji)
@@ -449,9 +503,17 @@
 
                ! Observed value
-               surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+               DO jvar = 1, kvars               
+                  surfdata%robs(iobs,jvar) = inpfiles(jj)%pob(1,ji,jvar)
+               END DO
+               IF ( TRIM(surfdata%cvars(1)) == 'FBD' ) THEN
+                   surfdata%rext(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+                   surfdata%rext(iobs,2) = fbrmdi
+               ENDIF
 
                ! Model and MDT is set to fbrmdi unless read from file
                IF ( ldmod ) THEN
-                  surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
+                  DO jvar = 1, kvars                              
+                     surfdata%rmod(iobs,jvar) = inpfiles(jj)%padd(1,ji,1,jvar)
+                  END DO
                   IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
                      surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
@@ -459,8 +521,17 @@
                   ENDIF
                 ELSE
-                  surfdata%rmod(iobs,1) = fbrmdi
+                  DO jvar = 1, kvars                
+                     surfdata%rmod(iobs,jvar) = fbrmdi
+                  END DO
                   IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi
                ENDIF
-
+               
+               ! Initialise climatology if set
+               IF ( surfdata%lclim ) THEN
+                  DO jvar = 1, kvars
+                     surfdata%rclm(iobs,jvar) = fbrmdi
+                  END DO
+               ENDIF
+               
                ! STD (obs error standard deviation) read from file and passed through obs operator
                IF ( iadd_std(jj) /= -1 ) THEN
@@ -483,7 +554,14 @@
       !-----------------------------------------------------------------------
       IF (lwp) THEN
-
+         DO jvar = 1, surfdata%nvar       
+            IF ( jvar == 1 ) THEN
+               cout1=TRIM(surfdata%cvars(1))                  
+            ELSE
+               WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdata%cvars(jvar))            
+            ENDIF
+         END DO
+ 
          WRITE(numout,*)
-         WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data'
+         WRITE(numout,'(1X,A)')TRIM( cout1 )//' data'
          WRITE(numout,'(1X,A)')'--------------'
          DO jj = 1,8
@@ -495,5 +573,5 @@
             & '---------------------------------------------------------------'
          WRITE(numout,'(1X,A,I8)') &
-            & 'Total data for variable '//TRIM( surfdata%cvars(1) )// &
+            & 'Total data for variable '//TRIM( cout1 )// &
             & '           = ', iobsmpp
          WRITE(numout,'(1X,A)') &
@@ -506,5 +584,5 @@
       ! Deallocate temporary data
       !-----------------------------------------------------------------------
-      DEALLOCATE( ifileidx, isurfidx, zdat, clvars )
+      DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin )
 
       !-----------------------------------------------------------------------
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readsnowdepth.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readsnowdepth.F90	(revision 15764)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readsnowdepth.F90	(revision 15764)
@@ -0,0 +1,152 @@
+MODULE obs_readsnowdepth
+   !!======================================================================
+   !!                       ***  MODULE obs_readsnowdepth  ***
+   !! Observation diagnostics: Get the snow depth for freeboard conversion to thickness
+   !!======================================================================
+   !! History :      ! 2018-10 (E. Fiedler) Initial version
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   obs_rea_snowdepth    : Driver for reading MDT
+   !!----------------------------------------------------------------------
+   USE wrk_nemo         ! Memory Allocation
+   USE par_kind         ! Precision variables
+   USE par_oce          ! Domain parameters
+   USE in_out_manager   ! I/O manager
+   USE obs_surf_def     ! Surface observation definitions
+   USE obs_inter_sup    ! Interpolation support routines
+   USE obs_inter_h2d    ! 2D interpolation
+   USE obs_utils        ! Various observation tools
+   USE iom_nf90         ! IOM NetCDF
+   USE netcdf           ! NetCDF library
+   USE lib_mpp          ! MPP library
+   USE dom_oce, ONLY : &                  ! Domain variables
+      &                    tmask, tmask_i, e1t, e2t, gphit, glamt
+   USE obs_const, ONLY :   obfillflt      ! Fillvalue
+   USE oce      , ONLY :   sshn           ! Model variables
+
+
+
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   obs_rea_snowdepth     ! called by dia_obs
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE obs_rea_snowdepth( fbddata, k2dint, thick_s )
+      !!---------------------------------------------------------------------
+      !!
+      !!                   *** ROUTINE obs_rea_snowdepth ***
+      !!
+      !! ** Purpose : Get snowdepth from CICE
+      !!
+      !! ** Method  : 
+      !!
+      !! ** Action  : 
+      !!----------------------------------------------------------------------
+      USE iom
+      !
+      TYPE(obs_surf), INTENT(inout)  ::   fbddata   ! Sea ice freeboard data
+      INTEGER       , INTENT(in)     ::   k2dint    ! ?
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: thick_s ! Model snow depth
+     
+      CHARACTER(LEN=12), PARAMETER ::   cpname  = 'obs_rea_snowdepth'
+
+      INTEGER ::   jobs                ! Obs loop variable
+      INTEGER ::   jpi_thick_s, jpj_thick_s      ! Number of grid point in lat/lon for the snow depth
+      INTEGER ::   iico, ijco          ! Grid point indices
+
+      !
+      REAL(wp), DIMENSION(1)     ::   zext, zobsmask
+      REAL(wp), DIMENSION(2,2,1) ::   zweig
+      !
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zmask, z_thick_s_l, zglam, zgphi
+      INTEGER , DIMENSION(:,:,:), ALLOCATABLE ::   igrdi, igrdj
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) ::  z_thick_s, thick_s_mask
+         
+      REAL(wp) :: zlam, zphi, zfill, zinfill    ! local scalar
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc(jpi,jpj,z_thick_s,thick_s_mask) 
+
+      IF(lwp)WRITE(numout,*) 
+      IF(lwp)WRITE(numout,*) ' obs_rea_snowdepth : Get model snow depth for freeboard conversion to sea ice thickness'
+      IF(lwp)WRITE(numout,*) ' ------------- '
+      CALL FLUSH(numout)
+
+      ! Get ice thickness information
+      z_thick_s = thick_s
+
+      ! Setup mask based on tmask
+      thick_s_mask(:,:) = tmask(:,:,1)
+
+      ! Interpolate the snow depth already on the model grid at the observation point
+  
+      ALLOCATE( &
+         & igrdi(2,2,fbddata%nsurf), &
+         & igrdj(2,2,fbddata%nsurf), &
+         & zglam(2,2,fbddata%nsurf), &
+         & zgphi(2,2,fbddata%nsurf), &
+         & zmask(2,2,fbddata%nsurf), &
+         & z_thick_s_l(2,2,fbddata%nsurf)  &
+         & )
+         
+      DO jobs = 1, fbddata%nsurf
+
+         igrdi(1,1,jobs) = fbddata%mi(jobs)-1
+         igrdj(1,1,jobs) = fbddata%mj(jobs)-1
+         igrdi(1,2,jobs) = fbddata%mi(jobs)-1
+         igrdj(1,2,jobs) = fbddata%mj(jobs)
+         igrdi(2,1,jobs) = fbddata%mi(jobs)
+         igrdj(2,1,jobs) = fbddata%mj(jobs)-1
+         igrdi(2,2,jobs) = fbddata%mi(jobs)
+         igrdj(2,2,jobs) = fbddata%mj(jobs)
+
+      END DO
+
+      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, glamt  , zglam )
+      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, gphit  , zgphi )
+      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, thick_s_mask, zmask )
+      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, z_thick_s  , z_thick_s_l )
+
+      DO jobs = 1, fbddata%nsurf
+            
+         zlam = fbddata%rlam(jobs)
+         zphi = fbddata%rphi(jobs)
+
+         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
+            &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
+            &                   zmask(:,:,jobs), zweig, zobsmask )
+            
+         CALL obs_int_h2d( 1, 1, zweig, z_thick_s_l(:,:,jobs),  zext )
+ 
+         fbddata%rext(jobs,2) = zext(1)
+
+! mark any masked data with a QC flag
+         IF( zobsmask(1) == 0 )   fbddata%nqc(jobs) = IBSET(fbddata%nqc(jobs),15)
+
+         END DO
+         
+      DEALLOCATE( &
+         & igrdi, &
+         & igrdj, &
+         & zglam, &
+         & zgphi, &
+         & zmask, &
+         & z_thick_s_l  &
+         & )
+
+      CALL wrk_dealloc(jpi,jpj,z_thick_s,thick_s_mask) 
+      IF(lwp)WRITE(numout,*) ' ------------- '
+      !
+   END SUBROUTINE obs_rea_snowdepth
+   
+   !!======================================================================
+END MODULE obs_readsnowdepth
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90	(revision 15764)
@@ -6,5 +6,6 @@
 
    !!----------------------------------------------------------------------
-   !!   obs_rotvel : Rotate velocity data into N-S,E-W directorions
+   !!   obs_rotvel_pro :  Rotate profile velocity data into N-S,E-W directions 
+   !!   obs_rotvel_surf : Rotate surface velocity data into N-S,E-W directions
    !!----------------------------------------------------------------------
    !! * Modules used   
@@ -17,4 +18,5 @@
    USE obs_utils                ! For error handling
    USE obs_profiles_def         ! Profile definitions
+   USE obs_surf_def             ! Surface definitions
    USE obs_inter_h2d            ! Horizontal interpolation
    USE obs_inter_sup            ! MPP support routines for interpolation
@@ -27,5 +29,6 @@
    PRIVATE
 
-   PUBLIC obs_rotvel            ! Rotate the observations
+   PUBLIC obs_rotvel_pro, &     ! Rotate the profile velocity observations 
+      &   obs_rotvel_surf       ! Rotate the surface velocity observations
 
    !!----------------------------------------------------------------------
@@ -37,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE obs_rotvel( profdata, k2dint, pu, pv )
+   SUBROUTINE obs_rotvel_pro( profdata, k2dint, pu, pv )
       !!---------------------------------------------------------------------
       !!
@@ -228,5 +231,189 @@
       CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 
 
-   END SUBROUTINE obs_rotvel
+   END SUBROUTINE obs_rotvel_pro
+
+   SUBROUTINE obs_rotvel_surf( surfdata, k2dint, pu, pv ) 
+      !!--------------------------------------------------------------------- 
+      !! 
+      !!                   *** ROUTINE obs_rotvel_surf *** 
+      !! 
+      !! ** Purpose : Rotate surface velocity data into N-S,E-W directorions 
+      !! 
+      !! ** Method  : Interpolation of geo2ocean coefficients on U,V grid 
+      !!              to observation point followed by a similar computations 
+      !!              as in geo2ocean. 
+      !! 
+      !! ** Action  : Review if there is a better way to do this. 
+      !! 
+      !! References :  
+      !! 
+      !! History :   
+      !!      ! :  2009-02 (K. Mogensen) : New routine 
+      !!---------------------------------------------------------------------- 
+      !! * Modules used 
+      !! * Arguments 
+      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Surface data to be read 
+      INTEGER, INTENT(IN) :: k2dint     ! Horizontal interpolation methed 
+      REAL(wp), DIMENSION(*) :: & 
+         & pu, & 
+         & pv 
+      !! * Local declarations 
+      REAL(wp), DIMENSION(2,2,1) :: zweig 
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
+         & zmasku, & 
+         & zmaskv, & 
+         & zcoslu, & 
+         & zsinlu, & 
+         & zcoslv, & 
+         & zsinlv, & 
+         & zglamu, & 
+         & zgphiu, & 
+         & zglamv, & 
+         & zgphiv 
+      REAL(wp), DIMENSION(1) :: & 
+         & zsinu, & 
+         & zcosu, & 
+         & zsinv, & 
+         & zcosv 
+      REAL(wp) :: zsin 
+      REAL(wp) :: zcos 
+      REAL(wp), DIMENSION(1) :: zobsmask 
+      REAL(wp), POINTER, DIMENSION(:,:) :: zsingu,zcosgu,zsingv,zcosgv 
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
+         & igrdiu, & 
+         & igrdju, & 
+         & igrdiv, & 
+         & igrdjv 
+      INTEGER :: ji 
+      INTEGER :: jk 
+
+      CALL wrk_alloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)  
+
+      !----------------------------------------------------------------------- 
+      ! Allocate data for message parsing and interpolation 
+      !----------------------------------------------------------------------- 
+
+      ALLOCATE( & 
+         & igrdiu(2,2,surfdata%nsurf), & 
+         & igrdju(2,2,surfdata%nsurf), & 
+         & zglamu(2,2,surfdata%nsurf), & 
+         & zgphiu(2,2,surfdata%nsurf), & 
+         & zmasku(2,2,surfdata%nsurf), & 
+         & zcoslu(2,2,surfdata%nsurf), & 
+         & zsinlu(2,2,surfdata%nsurf), & 
+         & igrdiv(2,2,surfdata%nsurf), & 
+         & igrdjv(2,2,surfdata%nsurf), & 
+         & zglamv(2,2,surfdata%nsurf), & 
+         & zgphiv(2,2,surfdata%nsurf), & 
+         & zmaskv(2,2,surfdata%nsurf), & 
+         & zcoslv(2,2,surfdata%nsurf), & 
+         & zsinlv(2,2,surfdata%nsurf)  & 
+         & ) 
+
+      !----------------------------------------------------------------------- 
+      ! Receive the angles on the U and V grids. 
+      !----------------------------------------------------------------------- 
+
+      CALL obs_rot( zsingu, zcosgu, zsingv, zcosgv ) 
+
+      DO ji = 1, surfdata%nsurf 
+         igrdiu(1,1,ji) = surfdata%mi(ji,1)-1 
+         igrdju(1,1,ji) = surfdata%mj(ji,1)-1 
+         igrdiu(1,2,ji) = surfdata%mi(ji,1)-1 
+         igrdju(1,2,ji) = surfdata%mj(ji,1) 
+         igrdiu(2,1,ji) = surfdata%mi(ji,1) 
+         igrdju(2,1,ji) = surfdata%mj(ji,1)-1 
+         igrdiu(2,2,ji) = surfdata%mi(ji,1) 
+         igrdju(2,2,ji) = surfdata%mj(ji,1) 
+         igrdiv(1,1,ji) = surfdata%mi(ji,2)-1 
+         igrdjv(1,1,ji) = surfdata%mj(ji,2)-1 
+         igrdiv(1,2,ji) = surfdata%mi(ji,2)-1 
+         igrdjv(1,2,ji) = surfdata%mj(ji,2) 
+         igrdiv(2,1,ji) = surfdata%mi(ji,2) 
+         igrdjv(2,1,ji) = surfdata%mj(ji,2)-1 
+         igrdiv(2,2,ji) = surfdata%mi(ji,2) 
+         igrdjv(2,2,ji) = surfdata%mj(ji,2) 
+      END DO 
+
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
+         &                  glamu, zglamu ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
+         &                  gphiu, zgphiu ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
+         &                  umask(:,:,1), zmasku ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
+         &                  zsingu, zsinlu ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
+         &                  zcosgu, zcoslu ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
+         &                  glamv, zglamv ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
+         &                  gphiv, zgphiv ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
+         &                  vmask(:,:,1), zmaskv ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
+         &                  zsingv, zsinlv ) 
+      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
+         &                  zcosgv, zcoslv ) 
+
+      DO ji = 1, surfdata%nsurf 
+
+         CALL obs_int_h2d_init( 1, 1, k2dint, & 
+            &                   surfdata%rlam(ji), surfdata%rphi(ji), & 
+            &                   zglamu(:,:,ji), zgphiu(:,:,ji),       & 
+            &                   zmasku(:,:,ji), zweig, zobsmask ) 
+
+         CALL obs_int_h2d( 1, 1, zweig, zsinlu(:,:,ji),  zsinu ) 
+
+         CALL obs_int_h2d( 1, 1, zweig, zcoslu(:,:,ji),  zcosu ) 
+
+         CALL obs_int_h2d_init( 1, 1, k2dint, & 
+            &                   surfdata%rlam(ji), surfdata%rphi(ji), & 
+            &                   zglamv(:,:,ji), zgphiv(:,:,ji),       & 
+            &                   zmaskv(:,:,ji), zweig, zobsmask ) 
+
+         CALL obs_int_h2d( 1, 1, zweig, zsinlv(:,:,ji),  zsinv ) 
+
+         CALL obs_int_h2d( 1, 1, zweig, zcoslv(:,:,ji),  zcosv ) 
+
+         ! Assume that the angle at observation point is the  
+         ! mean of u and v cosines/sines 
+
+         zcos = 0.5_wp * ( zcosu(1) + zcosv(1) ) 
+         zsin = 0.5_wp * ( zsinu(1) + zsinv(1) ) 
+
+         IF ( ( surfdata%rmod(ji,1) /= fbrmdi ) .AND. & 
+            & ( surfdata%rmod(ji,2) /= fbrmdi ) ) THEN 
+            pu(ji) = surfdata%rmod(ji,1) * zcos - & 
+               &     surfdata%rmod(ji,2) * zsin 
+            pv(ji) = surfdata%rmod(ji,2) * zcos + & 
+               &     surfdata%rmod(ji,1) * zsin 
+         ELSE 
+            pu(ji) = fbrmdi 
+            pv(ji) = fbrmdi 
+         ENDIF 
+
+      END DO 
+
+      DEALLOCATE( & 
+         & igrdiu, & 
+         & igrdju, & 
+         & zglamu, & 
+         & zgphiu, & 
+         & zmasku, & 
+         & zcoslu, & 
+         & zsinlu, & 
+         & igrdiv, & 
+         & igrdjv, & 
+         & zglamv, & 
+         & zgphiv, & 
+         & zmaskv, & 
+         & zcoslv, & 
+         & zsinlv  & 
+         & ) 
+
+      CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)  
+
+   END SUBROUTINE obs_rotvel_surf 
 
 END MODULE obs_rot_vel
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 15764)
@@ -52,9 +52,9 @@
       INTEGER :: nrec       !: Number of surface observation records in window
 
+      LOGICAL :: lclim      !: Climatology will be calculated for this structure
+
       ! Arrays with size equal to the number of surface observations
 
       INTEGER, POINTER, DIMENSION(:) :: &
-         & mi,   &        !: i-th grid coord. for interpolating to surface observation
-         & mj,   &        !: j-th grid coord. for interpolating to surface observation
          & mt,   &        !: time record number for gridded data
          & nsidx,&        !: Surface observation number
@@ -69,4 +69,8 @@
          & ntyp           !: Type of surface observation product
 
+      INTEGER, POINTER, DIMENSION(:,:) :: & 
+         & mi,   &        !: i-th grid coord. for interpolating to surface observation 
+         & mj             !: j-th grid coord. for interpolating to surface observation
+
       CHARACTER(len=8), POINTER, DIMENSION(:) :: &
          & cvars          !: Variable names
@@ -84,10 +88,11 @@
       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
          & robs, &        !: Surface observation 
-         & rmod           !: Model counterpart of the surface observation vector
+         & rmod, &        !: Model counterpart of the surface observation vector 
+         & rclm           !: Climatological counterpart of the surface observation vector
 
       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
          & rext           !: Extra fields interpolated to observation points
 
-      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
+      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: &
          & vdmean         !: Time averaged of model field
 
@@ -111,5 +116,4 @@
 
       ! Is this a gridded product?
-     
       LOGICAL :: lgrid
 
@@ -124,5 +128,5 @@
 CONTAINS
    
-   SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj )
+   SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj, ldclim )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE obs_surf_alloc  ***
@@ -143,4 +147,5 @@
       INTEGER, INTENT(IN) :: kpi     ! Number of 3D grid points
       INTEGER, INTENT(IN) :: kpj
+      LOGICAL, INTENT(IN) :: ldclim
 
       !!* Local variables
@@ -157,4 +162,5 @@
       surf%npi      = kpi
       surf%npj      = kpj
+      surf%lclim    = ldclim
 
       ! Allocate arrays of size number of variables
@@ -171,6 +177,4 @@
 
       ALLOCATE( &
-         & surf%mi(ksurf),      &
-         & surf%mj(ksurf),      &
          & surf%mt(ksurf),      &
          & surf%nsidx(ksurf),   &
@@ -196,8 +200,12 @@
 
       ALLOCATE( & 
+         & surf%mi(ksurf,kvar),   & 
+         & surf%mj(ksurf,kvar),   &
          & surf%robs(ksurf,kvar), &
          & surf%rmod(ksurf,kvar)  &
          & )   
 
+      IF (surf%lclim) ALLOCATE( surf%rclm(ksurf,kvar) )
+
       ! Allocate arrays of number of extra fields at observation points
 
@@ -223,5 +231,5 @@
 
       ALLOCATE( &
-         & surf%vdmean(kpi,kpj) &
+         & surf%vdmean(kpi,kpj,kvar) &
          & )
 
@@ -293,4 +301,5 @@
          & )
 
+      IF (surf%lclim) DEALLOCATE( surf%rclm )
       ! Deallocate arrays of number of extra fields at observation points
 
@@ -371,5 +380,5 @@
       IF ( lallocate ) THEN
          CALL obs_surf_alloc( newsurf,  insurf, surf%nvar, &
-            & surf%nextra, surf%nstp, surf%npi, surf%npj )
+            & surf%nextra, surf%nstp, surf%npi, surf%npj, surf%lclim )
       ENDIF
 
@@ -418,4 +427,5 @@
                newsurf%robs(insurf,jk)  = surf%robs(ji,jk)
                newsurf%rmod(insurf,jk)  = surf%rmod(ji,jk)
+               IF (newsurf%lclim) newsurf%rclm(insurf,jk) = surf%rclm(ji,jk)
                
             END DO
@@ -514,4 +524,5 @@
             oldsurf%robs(jj,jk)  = surf%robs(ji,jk)
             oldsurf%rmod(jj,jk)  = surf%rmod(ji,jk)
+            IF (surf%lclim) oldsurf%rclm(jj,jk)  = surf%rclm(ji,jk)
 
          END DO
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 15764)
@@ -28,5 +28,5 @@
    USE obs_const
    USE obs_mpp              ! MPP support routines for observation diagnostics
-   USE lib_mpp		    ! MPP routines
+   USE lib_mpp              ! MPP routines
 
    IMPLICIT NONE
@@ -95,6 +95,10 @@
       INTEGER :: je
       INTEGER :: iadd
+      INTEGER :: iadd_clm ! 1 if climatology present
       INTEGER :: iext
       REAL(wp) :: zpres
+
+      iadd_clm = 0  
+      IF ( profdata%lclim ) iadd_clm = 1
 
       IF ( PRESENT( padd ) ) THEN
@@ -137,4 +141,10 @@
          fbdata%caddunit(1,1) = 'Degrees centigrade'
          fbdata%caddunit(1,2) = 'PSU'
+         IF ( profdata%lclim ) THEN 
+            fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature' 
+            fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity' 
+            fbdata%caddunit(2,1) = 'Degrees centigrade' 
+            fbdata%caddunit(2,2) = 'PSU' 
+         ENDIF
          fbdata%cgrid(:)      = 'T'
          DO je = 1, iext
@@ -170,4 +180,10 @@
          fbdata%caddunit(1,1) = 'm/s'
          fbdata%caddunit(1,2) = 'm/s'
+         IF ( profdata%lclim ) THEN 
+            fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity' 
+            fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity' 
+            fbdata%caddunit(2,1) = 'm/s' 
+            fbdata%caddunit(2,2) = 'm/s' 
+         ENDIF
          fbdata%cgrid(1)      = 'U' 
          fbdata%cgrid(2)      = 'V'
@@ -252,4 +268,8 @@
          fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
          fbdata%caddunit(1,1) = clunits
+         IF ( profdata%lclim ) THEN 
+            fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname) 
+            fbdata%caddunit(2,1) = clunits 
+         ENDIF
          fbdata%cgrid(:)      = clgrid
          DO je = 1, iext
@@ -266,4 +286,6 @@
 
       fbdata%caddname(1)   = 'Hx'
+
+      IF ( profdata%lclim ) fbdata%caddname(1+iadd_clm)   = 'CLM'
 
       WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
@@ -319,5 +341,4 @@
             DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
                ik = profdata%var(jvar)%nvlidx(jk)
-               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
                fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
                fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
@@ -327,5 +348,5 @@
                   fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
                   fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
-                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111')
+                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111')
                ELSE
                   fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
@@ -333,6 +354,10 @@
                ENDIF
                fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
-               DO ja = 1, iadd
-                  fbdata%padd(ik,jo,1+ja,jvar) = &
+               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 
+               IF ( profdata%lclim ) THEN            
+                  fbdata%padd(ik,jo,1+iadd_clm,jvar) = profdata%var(jvar)%vclm(jk)      
+               ENDIF               
+               DO ja = 1, iadd 
+                  fbdata%padd(ik,jo,1+iadd_clm+ja,jvar) = &
                      & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
                END DO
@@ -416,19 +441,24 @@
       INTEGER :: ja
       INTEGER :: je
+      INTEGER :: jv
       INTEGER :: iadd
       INTEGER :: iext
       INTEGER :: indx_std
       INTEGER :: iadd_std
-
+      INTEGER :: iadd_clm      
+      INTEGER :: iadd_mdt
+
+      IF ( PRESENT( pext ) ) THEN 
+         iext = pext%inum 
+      ELSE 
+         iext = 0 
+      ENDIF 
+
+      ! Set up number of additional variables to be ouput: 
+      ! Hx, CLM, STD, MDT...
       IF ( PRESENT( padd ) ) THEN
          iadd = padd%inum
       ELSE
          iadd = 0
-      ENDIF
-
-      IF ( PRESENT( pext ) ) THEN
-         iext = pext%inum
-      ELSE
-         iext = 0
       ENDIF
 
@@ -444,4 +474,10 @@
       ENDIF
       
+      iadd_clm = 0 
+      IF ( surfdata%lclim ) iadd_clm = 1 
+
+      iadd_mdt = 0 
+      IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_mdt = 1 
+
       CALL init_obfbdata( fbdata )
 
@@ -451,7 +487,9 @@
          ! 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, &
-            &                 2 + iadd + iadd_std, 1 + iext, .TRUE. )
+         ! No climatology for SLA, MDT is our best estimate of that and is already output.
+
+         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
+            &                 1 + iadd_mdt + iadd_std + iadd, & 
+            &                 1 + iext, .TRUE. )
 
          clfiletype = 'slafb'
@@ -493,4 +531,20 @@
          clgrid     = 'T'
 
+      CASE('SIT') 
+
+         clfiletype    = 'sitfb' 
+         cllongname    = 'Sea ice thickness' 
+         clunits       = 'm' 
+         clgrid        = 'T' 
+
+      CASE('FBD') 
+
+         clfiletype = 'fbdfb' 
+         ! Change label from FBD ("freeboard") to SIT ("Sea Ice Thickness") 
+         surfdata%cvars(1) = 'SIT' 
+         cllongname = 'Sea ice thickness' 
+         clunits    = 'm' 
+         clgrid     = 'T' 
+
       CASE('SSS')
 
@@ -500,101 +554,118 @@
          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('UVEL') 
+
+         clfiletype    = 'ssvfb' 
+         cllongname    = 'Eastward sea surface velocity' 
+         clunits       = 'm s-1' 
+         clgrid        = 'U' 
+         cllongname    = 'Northward sea surface velocity' 
+         clunits       = 'm s-1' 
+         clgrid        = 'V' 
+
+      CASE('SLCHLTOT') 
+
+         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('SKD490') 
+
+         clfiletype    = 'skd490fb' 
+         cllongname    = 'Surface attenuation coefficient of downwelling radiation at 490 nm' 
+         clunits       = 'm-1' 
+         clgrid        = 'T' 
+
+      CASE('SFCO2') 
+
+         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
@@ -607,38 +678,52 @@
       ! Remaining variables treated more generically
 
-      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
-      
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd + iadd_std, iext, .TRUE. )
-
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = cllongname
-         fbdata%cobunit(1)    = clunits
-         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
-         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+iadd_std+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
-         END DO
-
+      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 
+
+         CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, & 
+            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. ) 
+
+         DO jv = 1, surfdata%nvar 
+            fbdata%cname(jv)      = surfdata%cvars(jv) 
+            fbdata%coblong(jv)    = cllongname 
+            fbdata%cobunit(jv)    = clunits
+         END DO 
+         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 jv = 1, surfdata%nvar          
+            IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 
+               fbdata%caddlong(1,jv) = 'Model interpolated ICE' 
+            ELSE 
+               fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv)) 
+            ENDIF 
+            fbdata%caddunit(1,jv) = clunits
+            fbdata%cgrid(jv)      = clgrid
+         END DO             
+         DO ja = 1, iadd 
+            fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja) 
+            DO jv = 1, surfdata%nvar                      
+               fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv) 
+               fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv) 
+            END DO 
+         END DO 
       ENDIF
       
       fbdata%caddname(1)   = 'Hx'
-      IF ( indx_std /= -1 ) THEN
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_std = iadd_std + 1
-         fbdata%caddname(1+iadd_std)   = surfdata%cext(indx_std)
-         fbdata%caddlong(1+iadd_std,1) = 'Obs error standard deviation'
-         fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1)
+      IF ( indx_std /= -1 ) THEN 
+         fbdata%caddname(1+iadd_mdt+iadd_std)   = surfdata%cext(indx_std) 
+         DO jv = 1, surfdata%nvar                         
+            fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation' 
+            fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1) 
+         END DO 
+      ENDIF 
+
+      IF ( surfdata%lclim ) THEN 
+         fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm)   = 'CLM' 
+         DO jv = 1, surfdata%nvar                                  
+            fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,jv) = 'Climatology' 
+            fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,jv) = fbdata%cobunit(1) 
+         END DO 
       ENDIF
       
@@ -649,5 +734,5 @@
          WRITE(numout,*)'obs_wri_surf :'
          WRITE(numout,*)'~~~~~~~~~~~~~'
-         WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
+         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
       ENDIF
 
@@ -663,5 +748,5 @@
             fbdata%ioqc(jo)    = 4
             fbdata%ioqcf(1,jo) = 0
-            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
+            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111')
          ELSE
             fbdata%ioqc(jo)    = surfdata%nqc(jo)
@@ -675,9 +760,13 @@
          fbdata%kindex(jo)    = surfdata%nsfil(jo)
          IF (ln_grid_global) THEN
-            fbdata%iobsi(jo,1) = surfdata%mi(jo)
-            fbdata%iobsj(jo,1) = surfdata%mj(jo)
+            DO jv = 1, surfdata%nvar 
+               fbdata%iobsi(jo,jv) = surfdata%mi(jo,jv) 
+               fbdata%iobsj(jo,jv) = surfdata%mj(jo,jv) 
+            END DO
          ELSE
-            fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
-            fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
+            DO jv = 1, surfdata%nvar          
+               fbdata%iobsi(jo,jv) = mig(surfdata%mi(jo,jv)) 
+               fbdata%iobsj(jo,jv) = mjg(surfdata%mj(jo,jv)) 
+            END DO
          ENDIF
          CALL greg2jul( 0, &
@@ -689,33 +778,48 @@
             &           fbdata%ptim(jo),   &
             &           krefdate = 19500101 )
-         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
-         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
          fbdata%pdep(1,jo)     = 0.0
          fbdata%idqc(1,jo)     = 0
          fbdata%idqcf(:,1,jo)  = 0
-         IF ( surfdata%nqc(jo) > 255 ) THEN
-            fbdata%ivqc(jo,1)       = 4
-            fbdata%ivlqc(1,jo,1)    = 4
-            fbdata%ivlqcf(1,1,jo,1) = 0
-            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
-         ELSE
-            fbdata%ivqc(jo,1)       = surfdata%nqc(jo)
-            fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo)
-            fbdata%ivlqcf(:,1,jo,1) = 0
-         ENDIF
-         fbdata%iobsk(1,jo,1)  = 0
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
-         IF ( indx_std /= -1 ) THEN
-            fbdata%padd(1,jo,1+iadd_std,1) = surfdata%rext(jo,indx_std)
-         ENDIF
-         
-         DO ja = 1, iadd
-            fbdata%padd(1,jo,2+iadd_std+ja,1) = &
-               & surfdata%rext(jo,padd%ipoint(ja))
-         END DO
-         DO je = 1, iext
-            fbdata%pext(1,jo,1+je) = &
-               & surfdata%rext(jo,pext%ipoint(je))
+         DO jv = 1, surfdata%nvar  
+            fbdata%pob(1,jo,jv)    = surfdata%robs(jo,jv) 
+
+            IF ( surfdata%nqc(jo) > 255 ) THEN 
+               fbdata%ivqc(jo,jv)       = 4 
+               fbdata%ivlqc(1,jo,jv)    = 4 
+               fbdata%ivlqcf(1,1,jo,jv) = 0 
+               fbdata%ivlqcf(2,1,jo,jv) = IAND(surfdata%nqc(jo),b'0000000011111111') 
+            ELSE 
+               fbdata%ivqc(jo,jv)       = surfdata%nqc(jo) 
+               fbdata%ivlqc(1,jo,jv)    = surfdata%nqc(jo) 
+               fbdata%ivlqcf(:,1,jo,jv) = 0 
+            ENDIF 
+            fbdata%iobsk(1,jo,jv)  = 0 
+
+            ! Additional variables. 
+            ! Hx is always the first additional variable 
+            fbdata%padd(1,jo,1,jv) = surfdata%rmod(jo,jv) 
+            ! MDT is output as an additional variable if SLA obs type 
+            IF ( TRIM(surfdata%cvars(jv)) == 'SLA' ) THEN 
+               fbdata%padd(1,jo,1+iadd_mdt,jv) = surfdata%rext(jo,1) 
+            ENDIF     
+            ! STD is output as an additional variable if available 
+            IF ( indx_std /= -1 ) THEN 
+               fbdata%padd(1,jo,1+iadd_mdt+iadd_std,jv) = surfdata%rext(jo,indx_std) 
+            ENDIF 
+            ! CLM is output as an additional variable if available 
+            IF ( surfdata%lclim ) THEN 
+               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,jv) = surfdata%rclm(jo,jv) 
+            ENDIF 
+            ! Then other additional variables are output 
+            DO ja = 1, iadd 
+               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = & 
+                  & surfdata%rext(jo,padd%ipoint(ja)) 
+            END DO 
+         END DO          
+         ! Extra variables 
+         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)          
+         DO je = 1, iext 
+            fbdata%pext(1,jo,1+je) = & 
+               & surfdata%rext(jo,pext%ipoint(je)) 
          END DO
       END DO
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90	(revision 15764)
@@ -62,29 +62,52 @@
          z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
          z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
-         IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp
-
-         zsum = z1dm + z1dp
          
-         IF ( k1dint == 0 ) THEN
-
-            !-----------------------------------------------------------------
-            !  Linear interpolation
-            !-----------------------------------------------------------------
-            pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &
-               &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum
-
-         ELSEIF ( k1dint == 1 ) THEN
-
-            !-----------------------------------------------------------------
-            ! Cubic spline interpolation
-            !-----------------------------------------------------------------
-            zsum2 = zsum * zsum
-            pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) &
-               &           + z1dp                             * pobsk (kkco(jdep)  ) &
-               &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) &
-               &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) &
-               &             ) / 6.0_wp                                              &
-               &          ) / zsum
-
+         ! Where both levels are masked, return a fill value 
+         IF ( ( pobsmask(kkco(jdep)-1) == 0.0_wp ) .AND. (pobsmask(kkco(jdep)) == 0.0_wp) ) THEN 
+            pobs(jdep)  = 99999. 
+         ELSE 
+          
+            ! Where upper level is masked (e.g., under ice cavity), only use deeper level 
+            ! otherwise where ob is at or above upper level model T-point,  
+            ! use upper model level rather than extrapolate 
+            IF ( pobsmask(kkco(jdep)-1) == 0.0_wp ) THEN 
+               z1dm = 0.0_wp 
+            ELSE IF ( pobsdep(jdep) <= pdep(kkco(jdep)-1) ) THEN 
+               z1dp = 0.0_wp    
+            END IF    
+ 
+            ! Where deeper level is masked (e.g., near sea bed), only use upper level 
+            ! otherwise where ob is at or below deeper level model T-point,  
+            ! use deeper model level rather than extrapolate 
+            IF ( pobsmask(kkco(jdep)) == 0.0_wp ) THEN 
+               z1dp = 0.0_wp 
+            ELSE IF ( pobsdep(jdep) >= pdep(kkco(jdep)) ) THEN 
+               z1dm = 0.0_wp    
+            END IF    
+ 
+            zsum = z1dm + z1dp 
+          
+            IF ( k1dint == 0 ) THEN 
+ 
+               !----------------------------------------------------------------- 
+               !  Linear interpolation 
+               !----------------------------------------------------------------- 
+               pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) & 
+                  &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum 
+ 
+            ELSEIF ( k1dint == 1 ) THEN 
+ 
+               !----------------------------------------------------------------- 
+               ! Cubic spline interpolation 
+               !----------------------------------------------------------------- 
+               zsum2 = zsum * zsum 
+               pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) & 
+                  &           + z1dp                             * pobsk (kkco(jdep)  ) & 
+                  &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) & 
+                  &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) & 
+                  &             ) / 6.0_wp                                              & 
+                  &          ) / zsum 
+ 
+            ENDIF 
          ENDIF
       END DO
Index: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 15754)
+++ branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 15764)
@@ -56,4 +56,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1)
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tclim    !: temperature climatology on each time step(Celcius) 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sclim    !: salinity climatology on each time step (psu)
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -70,5 +73,6 @@
       !!                ***  FUNCTION tra_dmp_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), STAT= tra_dmp_alloc )
+      ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), & 
+         &      tclim(jpi,jpj,jpk) , sclim(jpi,jpj,jpk), STAT=tra_dmp_alloc )
       !
       IF( lk_mpp            )   CALL mpp_sum ( tra_dmp_alloc )
@@ -110,4 +114,7 @@
       !                           !==   input T-S data at kt   ==!
       CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt
+
+      tclim(:,:,:) = zts_dta(:,:,:,jp_tem) 
+      sclim(:,:,:) = zts_dta(:,:,:,jp_sal)
       !
       SELECT CASE ( nn_zdmp )     !==    type of damping   ==!
