Index: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 12610)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 13393)
@@ -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
@@ -67,4 +68,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
@@ -80,4 +83,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)
@@ -169,4 +173,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
@@ -207,4 +212,5 @@
       LOGICAL :: ln_sit          ! Logical switch for sea ice thickness
       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
@@ -261,8 +267,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_sss, ln_vel3d,               &
+         &            ln_sst, ln_sic, ln_sit, ln_sss,                 &
+         &            ln_ssv, ln_vel3d,                               &
          &            ln_slchltot, ln_slchldia, ln_slchlnon,          &
          &            ln_slchldin, ln_slchlmic, ln_slchlnan,          &
@@ -280,9 +289,9 @@
          &            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_sss_fp_indegs, ln_ssv_fp_indegs,             &
+         &            ln_sic_fp_indegs, ln_sit_fp_indegs,             &
          &            cn_profbfiles, cn_slafbfiles,                   &
          &            cn_sstfbfiles, cn_sicfbfiles,                   &
-         &            cn_sitfbfiles,                                  &
+         &            cn_sitfbfiles, cn_ssvfbfiles,                   &
          &            cn_velfbfiles, cn_sssfbfiles,                   &
          &            cn_slchltotfbfiles, cn_slchldiafbfiles,         &
@@ -305,9 +314,11 @@
          &            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,             &
          &            nn_1dint, nn_2dint_default,                     &
          &            nn_2dint_sla, nn_2dint_sst,                     &
-         &            nn_2dint_sss, nn_2dint_sic, nn_2dint_sit,       &
+         &            nn_2dint_sss, nn_2dint_ssv,                     &
+         &            nn_2dint_sic, nn_2dint_sit,                     &
          &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
          &            nn_profdavtypes
@@ -325,4 +336,5 @@
       cn_velfbfiles(:)      = ''
       cn_sssfbfiles(:)      = ''
+      cn_ssvfbfiles(:)      = ''      
       cn_slchltotfbfiles(:) = ''
       cn_slchldiafbfiles(:) = ''
@@ -389,4 +401,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
@@ -424,4 +437,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
@@ -462,5 +476,5 @@
          &                  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_sss,             &
+      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sit, ln_sss, ln_ssv,     &
          &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, &
          &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  &
@@ -591,4 +605,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
@@ -721,4 +740,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
@@ -903,4 +932,7 @@
                nvarssurf(jtype) = 1
                nextrsurf(jtype) = 2
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN
+               nvarssurf(jtype) = 2
+               nextrsurf(jtype) = 0
             ELSE
                nvarssurf(jtype) = 1
@@ -909,5 +941,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'
@@ -921,4 +971,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'
@@ -960,5 +1013,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
@@ -990,4 +1046,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
@@ -1094,11 +1153,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
@@ -1346,9 +1405,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 )
@@ -1357,18 +1412,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
@@ -1381,7 +1456,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', &
@@ -1400,5 +1475,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' )
@@ -1412,11 +1487,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
@@ -1432,8 +1507,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', &
@@ -1448,8 +1523,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
@@ -1468,5 +1543,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', &
@@ -1484,5 +1559,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', &
@@ -1500,5 +1575,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', &
@@ -1516,5 +1591,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', &
@@ -1526,12 +1601,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', &
@@ -1542,12 +1617,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
@@ -1563,8 +1638,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', &
@@ -1579,9 +1654,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', &
@@ -1592,7 +1667,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
@@ -1611,5 +1686,5 @@
                ! light_xEPS diagnostic variable
                fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_xeps)
-               zsurfvar(:,:) = fabm_3d(:,:,1)
+               zsurfvar(:,:,1) = fabm_3d(:,:,1)
 #else
                CALL ctl_stop( ' Trying to run skd490 observation operator', &
@@ -1619,18 +1694,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)
+               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
@@ -1646,10 +1721,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', &
@@ -1659,17 +1734,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(:,:)
+               zsurfvar(:,:,1) = f2_pco2w(:,:)
 #elif defined key_fabm
                fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
-               zsurfvar(:,:) = fabm_3d(:,:,1)
+               zsurfvar(:,:,1) = fabm_3d(:,:,1)
 #else
                CALL ctl_stop( ' Trying to run spco2 observation operator', &
@@ -1686,43 +1761,55 @@
                ! 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
+            
             ! Change label of data from FBD ("freeboard") to SIT ("Sea Ice
             ! Thickness")
             IF ( TRIM(surfdataqc(jtype)%cvars(1)) == 'FBD' ) THEN
-                 surfdata(jtype)%cvars(1) = 'SIT'
+                 surfdataqc(jtype)%cvars(1) = 'SIT'
             ENDIF
 
+            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 )
@@ -1780,5 +1867,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
@@ -1813,4 +1900,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), &
