Index: /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 14165)
@@ -1202,5 +1202,6 @@
    ln_fbd     = .false.             ! Logical switch for Sea Ice Freeboard observations
    ln_vel3d   = .false.             ! Logical switch for velocity observations
-   ln_sss     = .false.             ! Logical swithc for SSS observations
+   ln_sss     = .false.             ! Logical switch for SSS observations
+   ln_ssv     = .false.             ! Logical switch for SSV (surface velocity) observations   
    ln_slchltot = .false.            ! Logical switch for surface total              log10(chlorophyll) obs
    ln_slchldia = .false.            ! Logical switch for surface diatom             log10(chlorophyll) obs
@@ -1241,4 +1242,5 @@
    ln_sst_fp_indegs = .true.
    ln_sss_fp_indegs = .true.
+   ln_ssv_fp_indegs = .true.   
    ln_sic_fp_indegs = .true.
    ln_sit_fp_indegs = .true.
@@ -1253,4 +1255,5 @@
    cn_velfbfiles = 'vel_01.nc'           ! Velocity feedback input observation file names
    cn_sssfbfiles = 'sss_01.nc'           ! SSS feedback input observation file names
+   cn_ssvfbfiles = 'ssv_01.nc'           ! SSV feedback input observation file names   
    cn_slchltotfbfiles = 'slchltot_01.nc' ! Surface total              log10(chlorophyll) obs file names
    cn_slchldiafbfiles = 'slchldia_01.nc' ! Surface diatom             log10(chlorophyll) obs file names
@@ -1289,4 +1292,6 @@
    rn_sss_avglamscl = 0.                 ! E/W diameter of SSS observation footprint (metres/degrees)
    rn_sss_avgphiscl = 0.                 ! N/S diameter of SSS observation footprint (metres/degrees)
+   rn_ssv_avglamscl = 0.                 ! E/W diameter of SSV observation footprint (metres/degrees)
+   rn_ssv_avgphiscl = 0.                 ! N/S diameter of SSV observation footprint (metres/degrees)
    rn_sic_avglamscl = 0.                 ! E/W diameter of SIC observation footprint (metres/degrees)
    rn_sic_avgphiscl = 0.                 ! N/S diameter of SIC observation footprint (metres/degrees)
@@ -1300,4 +1305,5 @@
    nn_2dint_sst = -1                     ! Horizontal interpolation method for SST
    nn_2dint_sss = -1                     ! Horizontal interpolation method for SSS
+   nn_2dint_ssv = -1                     ! Horizontal interpolation method for SSV   
    nn_2dint_sic = -1                     ! Horizontal interpolation method for SIC
    nn_2dint_sit = -1                     ! Horizontal interpolation method for SIT
Index: /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 14165)
@@ -51,7 +51,8 @@
    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_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
@@ -68,4 +69,6 @@
    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
@@ -83,4 +86,5 @@
    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)
@@ -174,4 +178,5 @@
          & 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
@@ -213,4 +218,5 @@
       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
@@ -267,9 +273,11 @@
       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_vel3d,                               &
+         &            ln_sss, ln_ssv, ln_vel3d,                       &
          &            ln_slchltot, ln_slchldia, ln_slchlnon,          &
          &            ln_slchldin, ln_slchlmic, ln_slchlnan,          &
@@ -287,10 +295,11 @@
          &            ln_time_mean_sla_bkg, ln_default_fp_indegs,     &
          &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
-         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
-         &            ln_sit_fp_indegs, ln_fbd_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_velfbfiles, cn_sssfbfiles, cn_ssvfbfiles,    &
          &            cn_slchltotfbfiles, cn_slchldiafbfiles,         &
          &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         &
@@ -312,4 +321,5 @@
          &            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,             &
@@ -317,6 +327,7 @@
          &            nn_1dint, nn_2dint_default,                     &
          &            nn_2dint_sla, nn_2dint_sst,                     &
-         &            nn_2dint_sss, nn_2dint_sic,                     &
-         &            nn_2dint_sit, nn_2dint_fbd,                     &
+         &            nn_2dint_sss, nn_2dint_ssv,                     &
+         &            nn_2dint_sic, nn_2dint_sit,                     &
+         &            nn_2dint_fbd,                                   &
          &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
          &            nn_profdavtypes
@@ -335,4 +346,5 @@
       cn_velfbfiles(:)      = ''
       cn_sssfbfiles(:)      = ''
+      cn_ssvfbfiles(:)      = ''      
       cn_slchltotfbfiles(:) = ''
       cn_slchldiafbfiles(:) = ''
@@ -400,4 +412,5 @@
          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
@@ -435,4 +448,5 @@
          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
@@ -477,5 +491,6 @@
          &                  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,     &
+      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,  &
@@ -611,4 +626,9 @@
             cobstypessurf(jtype) = 'sss'
             clsurffiles(jtype,:) = cn_sssfbfiles
+         ENDIF
+         IF (ln_ssv) THEN
+            jtype = jtype + 1
+            cobstypessurf(jtype) = 'ssv'
+            clsurffiles(jtype,:) = cn_ssvfbfiles
          ENDIF
          IF (ln_slchltot) THEN
@@ -751,4 +771,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
@@ -934,4 +964,7 @@
                nvarssurf(jtype) = 1
                nextrsurf(jtype) = 2
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN
+               nvarssurf(jtype) = 2
+               nextrsurf(jtype) = 0
             ELSE
                nvarssurf(jtype) = 1
@@ -940,5 +973,23 @@
             
             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'
@@ -955,4 +1006,7 @@
             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'
@@ -994,5 +1048,8 @@
                &               llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars )
 
-            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject, ln_seaicetypes )
+            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
@@ -1024,4 +1081,7 @@
             
             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
@@ -1104,6 +1164,4 @@
 #elif defined key_fabm
       USE par_fabm                 ! FABM parameters
-      USE fabm, ONLY: &
-         & fabm_get_interior_diagnostic_data
 #endif
 #if defined key_spm
@@ -1128,11 +1186,11 @@
       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
          & zprofmask               ! Mask associated with zprofvar
-      REAL(wp), POINTER, DIMENSION(:,:) :: &
+      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 for prof variables
-         & zgphi                   ! Model latitudes for prof variables
+         & zglam,    &             ! Model longitudes
+         & zgphi                   ! Model latitudes
       LOGICAL :: llog10            ! Perform log10 transform of variable
 #if defined key_fabm
@@ -1316,5 +1374,5 @@
 #elif defined key_fabm
                ! Alkalinity from ERSEM
-               zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
+               zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ta)
 #else
                CALL ctl_stop( ' Trying to run palk observation operator', &
@@ -1331,5 +1389,5 @@
 #elif defined key_fabm
                ! pH from ERSEM
-               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph)
+               zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ph)
 #else
                CALL ctl_stop( ' Trying to run pph observation operator', &
@@ -1380,9 +1438,5 @@
 
       IF ( nsurftypes > 0 ) THEN
-
-         !Allocate local work arrays
-         CALL wrk_alloc( jpi, jpj, zsurfvar )
-         CALL wrk_alloc( jpi, jpj, zsurfclim )         
-         CALL wrk_alloc( jpi, jpj, zsurfmask )
+         
 #if defined key_fabm
          CALL wrk_alloc( jpi, jpj, jpk, fabm_3d )
@@ -1391,18 +1445,38 @@
          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
-            zsurfmask(:,:) = tmask(:,:,1)
-            zsurfclim(:,:) = 0._wp          
+            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(:,:) = tsn(:,:,1,jp_tem)
-               IF ( ln_output_clim ) zsurfclim(:,:) = tclim(:,:,1)
+               zsurfvar(:,:,1) = tsn(:,:,1,jp_tem)
+               IF ( ln_output_clim ) zsurfclim(:,:,1) = tclim(:,:,1)
             CASE('sla')
-               zsurfvar(:,:) = sshn(:,:)
+               zsurfvar(:,:,1) = sshn(:,:)
             CASE('sss')
-               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
-               IF ( ln_output_clim ) zsurfclim(:,:) = sclim(:,:,1)              
+               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
@@ -1415,7 +1489,7 @@
                ELSE
 #if defined key_cice
-                  zsurfvar(:,:) = fr_i(:,:)
+                  zsurfvar(:,:,1) = fr_i(:,:)
 #elif defined key_lim2 || defined key_lim3
-                  zsurfvar(:,:) = 1._wp - frld(:,:)
+                  zsurfvar(:,:,1) = 1._wp - frld(:,:)
 #else
                CALL ctl_stop( ' Trying to run sea-ice concentration observation operator', &
@@ -1434,5 +1508,5 @@
                ELSE       
 #if defined key_cice
-                  zsurfvar(:,:) = thick_i(:,:)
+                  zsurfvar(:,:,1) = thick_i(:,:)
 #elif defined key_lim2 || defined key_lim3
                   CALL ctl_stop( ' No sea-ice thickness observation operator defined for LIM model' )
@@ -1465,11 +1539,11 @@
 #if defined key_hadocc
                ! Surface chlorophyll from HadOCC
-               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
+               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1)
 #elif defined key_medusa
                ! Add non-diatom and diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
+               zsurfvar(:,:,1) = 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) + &
+               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
@@ -1485,8 +1559,8 @@
 #elif defined key_medusa
                ! Diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchd)
+               zsurfvar(:,:,1) = trn(:,:,1,jpchd)
 #elif defined key_fabm
                ! Diatom surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1)
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1)
 #else
                CALL ctl_stop( ' Trying to run slchldia observation operator', &
@@ -1501,8 +1575,8 @@
 #elif defined key_medusa
                ! Non-diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchn)
+               zsurfvar(:,:,1) = 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) + &
+               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
@@ -1521,5 +1595,5 @@
 #elif defined key_fabm
                ! Dinoflagellate surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
 #else
                CALL ctl_stop( ' Trying to run slchldin observation operator', &
@@ -1537,5 +1611,5 @@
 #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)
+               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', &
@@ -1553,5 +1627,5 @@
 #elif defined key_fabm
                ! Nanophytoplankton surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2)
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2)
 #else
                CALL ctl_stop( ' Trying to run slchlnan observation operator', &
@@ -1569,5 +1643,5 @@
 #elif defined key_fabm
                ! Picophytoplankton surface chlorophyll from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3)
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3)
 #else
                CALL ctl_stop( ' Trying to run slchlpic observation operator', &
@@ -1579,12 +1653,12 @@
 #if defined key_hadocc
                ! Surface chlorophyll from HadOCC
-               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
+               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1)
 #elif defined key_medusa
                ! Add non-diatom and diatom surface chlorophyll from MEDUSA
-               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
+               zsurfvar(:,:,1) = 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)
+               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', &
@@ -1595,12 +1669,12 @@
 #if defined key_hadocc
                ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio
-               zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
+               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(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
+               zsurfvar(:,:,1) = (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) + &
+               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
@@ -1616,8 +1690,8 @@
 #elif defined key_medusa
                ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
-               zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd
+               zsurfvar(:,:,1) = trn(:,:,1,jpphd) * xthetapd
 #elif defined key_fabm
                ! Diatom surface phytoplankton carbon from ERSEM
-               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c)
+               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c)
 #else
                CALL ctl_stop( ' Trying to run slphydia observation operator', &
@@ -1632,9 +1706,9 @@
 #elif defined key_medusa
                ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
-               zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn
+               zsurfvar(:,:,1) = 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)
+               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', &
@@ -1645,7 +1719,7 @@
             CASE('sspm')
 #if defined key_spm
-               zsurfvar(:,:) = 0.0
+               zsurfvar(:,:,1) = 0.0
                DO jn = 1, jp_spm
-                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
+                  zsurfvar(:,:,1) = zsurfvar(:,:,1) + trn(:,:,1,jn)   ! sum SPM sizes
                END DO
 #else
@@ -1662,7 +1736,15 @@
                   &           ' but MEDUSA does not explicitly simulate Kd490' )
 #elif defined key_fabm
-               ! light_xEPS diagnostic variable
-               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_xeps)
-               zsurfvar(:,:) = fabm_3d(:,:,1)
+               ! 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', &
@@ -1672,18 +1754,18 @@
             CASE('sfco2')
 #if defined key_hadocc
-               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from 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(:,:) = 0
+                  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(:,:) = f2_fco2w(:,:)
+               zsurfvar(:,:,1) = f2_fco2w(:,:)
 #elif defined key_fabm
                ! First, get pCO2 from FABM
-               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
-               zsurfvar(:,:) = fabm_3d(:,:,1)
+               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
@@ -1699,10 +1781,10 @@
                ! 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)))
+               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', &
@@ -1712,17 +1794,17 @@
             CASE('spco2')
 #if defined key_hadocc
-               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
+               zsurfvar(:,:,1) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
                IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
                   & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
-                  zsurfvar(:,:) = obfillflt
-                  zsurfmask(:,:) = 0
+                  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(:,:) = f2_pco2w(:,:)
-#elif defined key_fabm
-               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
-               zsurfvar(:,:) = fabm_3d(:,:,1)
+               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', &
@@ -1739,37 +1821,49 @@
                ! Take the log10 where we can, otherwise exclude
                tiny = 1.0e-20
-               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
-                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:))
+               WHERE(zsurfvar(:,:,1) > tiny .AND. zsurfvar(:,:,1) /= obfillflt )
+                  zsurfvar(:,:,1)  = LOG10(zsurfvar(:,:,1))
                ELSEWHERE
-                  zsurfvar(:,:)  = obfillflt
-                  zsurfmask(:,:) = 0
+                  zsurfvar(:,:,1)  = obfillflt
+                  zsurfmask(:,:,1) = 0
                END WHERE
             ENDIF
 
-            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, zsurfvar,               &
-                  &               zsurfclim, zsurfmask,                    &
-                  &               n2dintsurf(jtype), llnightav(jtype),     &
-                  &               ravglamscl(jtype), ravgphiscl(jtype),    &
-                  &               lfpindegs(jtype), kmeanstp = imeanstp )
-
-
-            ELSE
-               CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
-                  &               nit000, idaystp, zsurfvar,               &
-                  &               zsurfclim, zsurfmask,                    &
-                  &               n2dintsurf(jtype), llnightav(jtype),     &
-                  &               ravglamscl(jtype), ravgphiscl(jtype),    &
-                  &               lfpindegs(jtype) )
-            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
-
-         CALL wrk_dealloc( jpi, jpj, zsurfvar )
-         CALL wrk_dealloc( jpi, jpj, zsurfmask )
 #if defined key_fabm
          CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d )
@@ -1827,5 +1921,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
@@ -1860,4 +1954,26 @@
 
          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/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 14165)
@@ -544,5 +544,7 @@
 
    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
-      &                     kit000, kdaystp, psurf, pclim, psurfmask,   &
+      &                     kit000, kdaystp, kvar, &
+      &                     psurf, pclim, psurfmask, &
+      &                     plam, pphi, &
       &                     k2dint, ldnightav, plamscl, pphiscl, &
       &                     lindegrees, kmeanstp )
@@ -579,4 +581,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      
       !!-----------------------------------------------------------------------
 
@@ -595,4 +598,5 @@
                                        !   (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 :: &
@@ -603,4 +607,7 @@
          & 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) :: &
@@ -670,5 +677,4 @@
       CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
 
-
       IF ( l_timemean ) THEN
          ! Initialize time mean for first timestep
@@ -681,5 +687,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  surfdataqc%vdmean(ji,jj) = 0.0
+                  surfdataqc%vdmean(ji,jj,kvar) = 0.0
                END DO
             END DO
@@ -690,6 +696,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
-                  &                        + psurf(ji,jj)
+               surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
+                  &                            + psurf(ji,jj)
             END DO
          END DO
@@ -701,6 +707,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
-                     &                       * zmeanstp
+                  surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
+                     &                            * zmeanstp
                END DO
             END DO
@@ -729,5 +735,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  surfdataqc%vdmean(ji,jj) = 0.0
+                  surfdataqc%vdmean(ji,jj,kvar) = 0.0
                   zmeanday(ji,jj) = 0.0
                   icount_night(ji,jj) = 0
@@ -743,6 +749,6 @@
             DO ji = 1, jpi
                ! Increment the temperature field for computing night mean and counter
-               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
-                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
+               surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar)  &
+                      &                        + psurf(ji,jj) * REAL( imask_night(ji,jj) )
                zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
                icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
@@ -758,10 +764,10 @@
                   ! Test if "no night" point
                   IF ( icount_night(ji,jj) > 0 ) THEN
-                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
-                       &                        / REAL( icount_night(ji,jj) )
+                     surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
+                       &                            / REAL( icount_night(ji,jj) )
                   ELSE
                      !At locations where there is no night (e.g. poles),
                      ! calculate daily mean instead of night-time mean.
-                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
+                     surfdataqc%vdmean(ji,jj,kvar) = zmeanday(ji,jj) * zdaystp
                   ENDIF
                END DO
@@ -772,19 +778,23 @@
 
       ! Get the data for interpolation
-
       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) )
@@ -793,5 +803,5 @@
          iobs = jobs - surfdataqc%nsurfup
          DO ji = 0, imaxifp
-            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
+            imodi = surfdataqc%mi(jobs,kvar) - int(imaxifp/2) + ji - 1
             
             !Deal with wrap around in longitude
@@ -800,5 +810,5 @@
             
             DO jj = 0, imaxjfp
-               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
+               imodj = surfdataqc%mj(jobs,kvar) - int(imaxjfp/2) + jj - 1
                !If model values are out of the domain to the north/south then
                !set them to be the edge of the domain
@@ -806,6 +816,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
@@ -819,7 +831,7 @@
 
       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
-         &                  igrdi, igrdj, glamt, zglam )
+         &                  igrdi, igrdj, plam, zglam )
       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
-         &                  igrdi, igrdj, gphit, zgphi )
+         &                  igrdi, igrdj, pphi, zgphi )
       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, psurfmask, zmask )
@@ -831,5 +843,5 @@
             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(:,:), zsurfm )
+               &                  igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm )
          ENDIF
       ELSE
@@ -858,5 +870,5 @@
 
          CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
-            &               surfdataqc%vdmean(:,:), zsurfm )
+            &                  surfdataqc%vdmean(:,:,kvar), zsurfm )
 
       ENDIF
@@ -937,5 +949,5 @@
                surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
             ELSE
-               surfdataqc%rmod(jobs,1) = zext(1)
+               surfdataqc%rmod(jobs,kvar) = zext(1)
             ENDIF
 
@@ -985,11 +997,16 @@
          & zmask, &
          & zsurf, &
-         & zsurftmp, &
-         & zglamf, &
-         & zgphif, &
-         & igrdip1,&
-         & igrdjp1 &
+         & zsurftmp &
          & )
 
+      IF ( k2dint > 4 ) THEN
+         DEALLOCATE( &     
+            & zglamf, &
+            & zgphif, &
+            & igrdip1,&
+            & igrdjp1 &
+            & )
+      ENDIF
+            
       IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
 
@@ -1001,5 +1018,7 @@
       ENDIF
 
-      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
+      IF ( kvar == surfdataqc%nvar ) THEN
+         surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
+      ENDIF
 
    END SUBROUTINE obs_surf_opt
Index: /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 14165)
@@ -53,8 +53,11 @@
 CONTAINS
 
-   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, &
-                            ld_seaicetypes, kqc_cutoff )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE obs_pre_sla  ***
+   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, &
+      &                     kpi, kpj, &   
+      &                     zmask, pglam, pgphi, &
+      &                     ld_nea, ld_bound_reject, &
+      &                     ld_seaicetypes, kqc_cutoff )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE obs_pre_surf  ***
       !!
       !! ** Purpose : First level check and screening of surface observations
@@ -82,4 +85,10 @@
       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
@@ -94,24 +103,27 @@
       INTEGER :: imin0
       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
+                                                        ! Counters for observations that are
+      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 variable
+      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...'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
@@ -130,8 +142,8 @@
       iotdobs  = 0
       igrdobs  = 0
-      iosdsobs = 0
-      ilansobs = 0
-      inlasobs = 0
-      ibdysobs = 0 
+      iosdsobs(:) = 0
+      ilansobs(:) = 0
+      inlasobs(:) = 0
+      ibdysobs(:) = 0 
 
       ! Set QC cutoff to optional value if provided
@@ -162,8 +174,9 @@
       ! Check for surface data failing the grid search
       ! -----------------------------------------------------------------------
-
-      CALL obs_coo_grd( surfdata%nsurf,   surfdata%mi, surfdata%mj, &
-         &              surfdata%nqc,     igrdobs                         )
-
+      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 )
 
@@ -172,20 +185,21 @@
       ! -----------------------------------------------------------------------
 
-      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
+           
       ! -----------------------------------------------------------------------
       ! Copy useful data from the surfdata data structure to
@@ -216,29 +230,41 @@
       
       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,*) ' '//surfdataqc%cvars(1)//' data outside time domain                  = ', &
+         WRITE(numout,*) ' '//TRIM(cout1)//' data outside time domain                  = ', &
             &            iotdobsmpp
-         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search    = ', &
+         WRITE(numout,*) ' Remaining '//TRIM(cout1)//' 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
+
+         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',surfdataqc%cvars(1)
+         WRITE(numout,'(10X,A,10X,A)')'Time step',TRIM(cout1)
          WRITE(numout,'(10X,A,5X,A)')'---------','-----------------'
          CALL FLUSH(numout)
@@ -445,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 )
@@ -1457,7 +1483,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
@@ -1513,5 +1539,5 @@
       END DO
 
-   END SUBROUTINE obs_uv_rej
+   END SUBROUTINE obs_uv_rej_pro
 
 END MODULE obs_prep
Index: /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 14165)
@@ -85,4 +85,5 @@
       INTEGER :: jj
       INTEGER :: jk
+      INTEGER :: jvar
       INTEGER :: iflag
       INTEGER :: inobf
@@ -107,8 +108,9 @@
          & ityp, &
          & itypmpp
-      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
          & iobsi,    &
          & iobsj,    &
-         & iproc,    &
+         & iproc
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &         
          & iindx,    &
          & ifileidx, &
@@ -126,5 +128,6 @@
       TYPE(obfbdata), POINTER, DIMENSION(:) :: &
          & inpfiles
-
+      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line
+      
       ! Local initialization
       iobs  = 0
@@ -278,7 +281,7 @@
 
             IF ( inpfiles(jj)%nobs > 0 ) THEN
-               inpfiles(jj)%iproc = -1
-               inpfiles(jj)%iobsi = -1
-               inpfiles(jj)%iobsj = -1
+               inpfiles(jj)%iproc(:,:) = -1
+               inpfiles(jj)%iobsi(:,:) = -1
+               inpfiles(jj)%iobsj(:,:) = -1
             ENDIF
 
@@ -295,5 +298,5 @@
 	            END DO
             ENDIF   
-
+            
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
@@ -305,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
@@ -318,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
@@ -325,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
@@ -448,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)
@@ -480,5 +503,7 @@
 
                ! 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)
@@ -488,5 +513,7 @@
                ! 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)
@@ -494,10 +521,16 @@
                   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 ) surfdata%rclm(iobs,1) = fbrmdi
+               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
@@ -521,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
@@ -533,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)') &
Index: /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90	(revision 14165)
@@ -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,6 +29,7 @@
    PRIVATE
 
-   PUBLIC obs_rotvel            ! Rotate the observations
-
+   PUBLIC obs_rotvel_pro, &     ! Rotate the profile velocity observations
+      &   obs_rotvel_surf       ! Rotate the surface velocity observations
+      
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -37,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE obs_rotvel( profdata, k2dint, pu, pv )
+   SUBROUTINE obs_rotvel_pro( profdata, k2dint, pu, pv )
       !!---------------------------------------------------------------------
       !!
@@ -228,5 +231,190 @@
       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/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 14165)
@@ -57,6 +57,4 @@
 
       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
@@ -71,4 +69,9 @@
          & 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
@@ -92,5 +95,5 @@
          & rext           !: Extra fields interpolated to observation points
 
-      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
+      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: &
          & vdmean         !: Time averaged of model field
 
@@ -176,6 +179,4 @@
 
       ALLOCATE( &
-         & surf%mi(ksurf),      &
-         & surf%mj(ksurf),      &
          & surf%mt(ksurf),      &
          & surf%nsidx(ksurf),   &
@@ -201,6 +202,8 @@
 
       ALLOCATE( & 
+         & surf%mi(ksurf,kvar),   &
+         & surf%mj(ksurf,kvar),   &
          & surf%robs(ksurf,kvar), &
-         & surf%rmod(ksurf,kvar) &
+         & surf%rmod(ksurf,kvar)  &
          & )   
 
@@ -230,5 +233,5 @@
 
       ALLOCATE( &
-         & surf%vdmean(kpi,kpj) &
+         & surf%vdmean(kpi,kpj,kvar) &
          & )
 
@@ -405,6 +408,6 @@
             insurf = insurf + 1
 
-            newsurf%mi(insurf)    = surf%mi(ji)
-            newsurf%mj(insurf)    = surf%mj(ji)
+            newsurf%mi(insurf,:)  = surf%mi(ji,:)
+            newsurf%mj(insurf,:)  = surf%mj(ji,:)
             newsurf%mt(insurf)    = surf%mt(ji)
             newsurf%nsidx(insurf) = surf%nsidx(ji)
@@ -496,6 +499,6 @@
          jj=surf%nsind(ji)
 
-         oldsurf%mi(jj)    = surf%mi(ji)
-         oldsurf%mj(jj)    = surf%mj(ji)
+         oldsurf%mi(jj,:)  = surf%mi(ji,:)
+         oldsurf%mj(jj,:)  = surf%mj(ji,:)
          oldsurf%mt(jj)    = surf%mt(ji)
          oldsurf%nsidx(jj) = surf%nsidx(ji)
Index: /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 14164)
+++ /branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 14165)
@@ -349,5 +349,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)
@@ -437,10 +437,11 @@
       CHARACTER(LEN=10) :: clfiletype
       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
-      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
-      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
-      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
+      CHARACTER(LEN=ilenlong), DIMENSION(surfdata%nvar) :: cllongname  ! Long name of variable
+      CHARACTER(LEN=ilenunit), DIMENSION(surfdata%nvar) :: clunits     ! Units of variable
+      CHARACTER(LEN=ilengrid), DIMENSION(surfdata%nvar) :: clgrid      ! Grid of variable
       INTEGER :: jo
       INTEGER :: ja
       INTEGER :: je
+      INTEGER :: jv
       INTEGER :: iadd
       INTEGER :: iext
@@ -522,22 +523,22 @@
       CASE('SST')
 
-         clfiletype = 'sstfb'
-         cllongname = 'Sea surface temperature'
-         clunits    = 'Degree centigrade'
-         clgrid     = 'T'
+         clfiletype    = 'sstfb'
+         cllongname(1) = 'Sea surface temperature'
+         clunits(1)    = 'Degree centigrade'
+         clgrid(1)     = 'T'
          
       CASE('ICECONC')
 
-         clfiletype = 'sicfb'
-         cllongname = 'Sea ice concentration'
-         clunits    = 'Fraction'
-         clgrid     = 'T'
+         clfiletype    = 'sicfb'
+         cllongname(1) = 'Sea ice concentration'
+         clunits(1)    = 'Fraction'
+         clgrid(1)     = 'T'
 
       CASE('SIT')
 
-         clfiletype = 'sitfb'
-         cllongname = 'Sea ice thickness'
-         clunits    = 'm'
-         clgrid     = 'T'
+         clfiletype    = 'sitfb'
+         cllongname(1) = 'Sea ice thickness'
+         clunits(1)    = 'm'
+         clgrid(1)     = 'T'
 
       CASE('FBD')
@@ -552,57 +553,67 @@
       CASE('SSS')
 
-         clfiletype = 'sssfb'
-         cllongname = 'Sea surface salinity'
-         clunits    = 'psu'
-         clgrid     = 'T'
+         clfiletype    = 'sssfb'
+         cllongname(1) = 'Sea surface salinity'
+         clunits(1)    = 'psu'
+         clgrid(1)     = 'T'
+
+      CASE('UVEL')
+
+         clfiletype    = 'ssvfb'
+         cllongname(1) = 'Eastward sea surface velocity'
+         clunits(1)    = 'm s-1'
+         clgrid(1)     = 'U'
+         cllongname(2) = 'Northward sea surface velocity'
+         clunits(2)    = 'm s-1'
+         clgrid(2)     = 'V'
          
       CASE('SLCHLTOT')
 
-         clfiletype = 'slchltotfb'
-         cllongname = 'Surface total log10(chlorophyll)'
-         clunits    = 'log10(mg/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slchltotfb'
+         cllongname(1) = 'Surface total log10(chlorophyll)'
+         clunits(1)    = 'log10(mg/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLCHLDIA')
 
-         clfiletype = 'slchldiafb'
-         cllongname = 'Surface diatom log10(chlorophyll)'
-         clunits    = 'log10(mg/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slchldiafb'
+         cllongname(1) = 'Surface diatom log10(chlorophyll)'
+         clunits(1)    = 'log10(mg/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLCHLNON')
 
-         clfiletype = 'slchlnonfb'
-         cllongname = 'Surface non-diatom log10(chlorophyll)'
-         clunits    = 'log10(mg/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slchlnonfb'
+         cllongname(1) = 'Surface non-diatom log10(chlorophyll)'
+         clunits(1)    = 'log10(mg/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLCHLDIN')
 
-         clfiletype = 'slchldinfb'
-         cllongname = 'Surface dinoflagellate log10(chlorophyll)'
-         clunits    = 'log10(mg/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slchldinfb'
+         cllongname(1) = 'Surface dinoflagellate log10(chlorophyll)'
+         clunits(1)    = 'log10(mg/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLCHLMIC')
 
-         clfiletype = 'slchlmicfb'
-         cllongname = 'Surface microphytoplankton log10(chlorophyll)'
-         clunits    = 'log10(mg/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slchlmicfb'
+         cllongname(1) = 'Surface microphytoplankton log10(chlorophyll)'
+         clunits(1)    = 'log10(mg/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLCHLNAN')
 
-         clfiletype = 'slchlnanfb'
-         cllongname = 'Surface nanophytoplankton log10(chlorophyll)'
-         clunits    = 'log10(mg/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slchlnanfb'
+         cllongname(1) = 'Surface nanophytoplankton log10(chlorophyll)'
+         clunits(1)    = 'log10(mg/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLCHLPIC')
 
-         clfiletype = 'slchlpicfb'
-         cllongname = 'Surface picophytoplankton log10(chlorophyll)'
-         clunits    = 'log10(mg/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slchlpicfb'
+         cllongname(1) = 'Surface picophytoplankton log10(chlorophyll)'
+         clunits(1)    = 'log10(mg/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SCHLTOT')
@@ -615,50 +626,50 @@
       CASE('SLPHYTOT')
 
-         clfiletype = 'slphytotfb'
-         cllongname = 'Surface total log10(phytoplankton carbon)'
-         clunits    = 'log10(mmolC/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slphytotfb'
+         cllongname(1) = 'Surface total log10(phytoplankton carbon)'
+         clunits(1)    = 'log10(mmolC/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLPHYDIA')
 
-         clfiletype = 'slphydiafb'
-         cllongname = 'Surface diatom log10(phytoplankton carbon)'
-         clunits    = 'log10(mmolC/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slphydiafb'
+         cllongname(1) = 'Surface diatom log10(phytoplankton carbon)'
+         clunits(1)    = 'log10(mmolC/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SLPHYNON')
 
-         clfiletype = 'slphynonfb'
-         cllongname = 'Surface non-diatom log10(phytoplankton carbon)'
-         clunits    = 'log10(mmolC/m3)'
-         clgrid     = 'T'
+         clfiletype    = 'slphynonfb'
+         cllongname(1) = 'Surface non-diatom log10(phytoplankton carbon)'
+         clunits(1)    = 'log10(mmolC/m3)'
+         clgrid(1)     = 'T'
 
       CASE('SSPM')
 
-         clfiletype = 'sspmfb'
-         cllongname = 'Surface suspended particulate matter'
-         clunits    = 'g/m3'
-         clgrid     = 'T'
+         clfiletype    = 'sspmfb'
+         cllongname(1) = 'Surface suspended particulate matter'
+         clunits(1)    = 'g/m3'
+         clgrid(1)     = 'T'
 
       CASE('SKD490')
 
-         clfiletype = 'skd490fb'
-         cllongname = 'Surface attenuation coefficient of downwelling radiation at 490 nm'
-         clunits    = 'm-1'
-         clgrid     = 'T'
+         clfiletype    = 'skd490fb'
+         cllongname(1) = 'Surface attenuation coefficient of downwelling radiation at 490 nm'
+         clunits(1)    = 'm-1'
+         clgrid(1)     = 'T'
 
       CASE('SFCO2')
 
-         clfiletype = 'sfco2fb'
-         cllongname = 'Surface fugacity of carbon dioxide'
-         clunits    = 'uatm'
-         clgrid     = 'T'
+         clfiletype    = 'sfco2fb'
+         cllongname(1) = 'Surface fugacity of carbon dioxide'
+         clunits(1)    = 'uatm'
+         clgrid(1)     = 'T'
 
       CASE('SPCO2')
 
-         clfiletype = 'spco2fb'
-         cllongname = 'Surface partial pressure of carbon dioxide'
-         clunits    = 'uatm'
-         clgrid     = 'T'
+         clfiletype    = 'spco2fb'
+         cllongname(1) = 'Surface partial pressure of carbon dioxide'
+         clunits(1)    = 'uatm'
+         clgrid(1)     = 'T'
 
       CASE DEFAULT
@@ -673,10 +684,12 @@
       IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
       
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
+         CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, &
             &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )
 
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = cllongname
-         fbdata%cobunit(1)    = clunits
+         DO jv = 1, surfdata%nvar
+            fbdata%cname(jv)      = surfdata%cvars(jv)
+            fbdata%coblong(jv)    = cllongname(jv)
+            fbdata%cobunit(jv)    = clunits(jv)
+         END DO
          DO je = 1, iext
             fbdata%cextname(je) = pext%cdname(je)
@@ -684,33 +697,39 @@
             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 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(jv)
+            fbdata%cgrid(jv)      = clgrid(jv)
+         END DO            
          DO ja = 1, iadd
             fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      ENDIF
-      
+            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
          fbdata%caddname(1+iadd_mdt+iadd_std)   = surfdata%cext(indx_std)
-         fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation'
-         fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1)
+         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'
-         fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,1) = 'Climatology'
-         fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,1) = fbdata%cobunit(1)
-      ENDIF
-      
-      
+         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
+
       WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
 
@@ -719,5 +738,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
 
@@ -733,5 +752,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)
@@ -745,9 +764,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, &
@@ -759,42 +782,45 @@
             &           fbdata%ptim(jo),   &
             &           krefdate = 19500101 )
-                    
-         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
- 
-         ! Additional variables.
-         ! Hx is always the first additional variable
-         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
-         ! MDT is output as an additional variable if SLA obs type
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
-            fbdata%padd(1,jo,1+iadd_mdt,1) = 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,1) = 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,1) = surfdata%rclm(jo,1)
-         ENDIF
-         ! Then other additional variables are output
-         DO ja = 1, iadd
-            fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,1) = &
-               & surfdata%rext(jo,padd%ipoint(ja))
-         END DO
+         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)         
