Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 5997)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 5998)
@@ -6,6 +6,4 @@
    !!======================================================================
 
-   !!----------------------------------------------------------------------
-   !!   'key_diaobs' : Switch on the observation diagnostic computation
    !!----------------------------------------------------------------------
    !!   dia_obs_init : Reading and prepare observations
@@ -16,5 +14,5 @@
    !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
    !!----------------------------------------------------------------------
-   !! * Modules used   
+   !! * Modules used
    USE wrk_nemo                 ! Memory Allocation
    USE par_kind                 ! Precision variables
@@ -22,12 +20,8 @@
    USE par_oce
    USE dom_oce                  ! Ocean space and time domain variables
-   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch
-   USE obs_read_prof            ! Reading and allocation of observations (Coriolis)
-   USE obs_read_sla             ! Reading and allocation of SLA observations  
-   USE obs_read_sst             ! Reading and allocation of SST observations  
+   USE obs_read_prof            ! Reading and allocation of profile obs
+   USE obs_read_surf            ! Reading and allocation of surface obs
    USE obs_sstbias              ! Bias correction routine for SST 
    USE obs_readmdt              ! Reading and allocation of MDT for SLA.
-   USE obs_read_seaice          ! Reading and allocation of Sea Ice observations  
-   USE obs_read_vel             ! Reading and allocation of velocity component observations
    USE obs_prep                 ! Preparation of obs. (grid search etc).
    USE obs_oper                 ! Observation operators
@@ -36,9 +30,5 @@
    USE obs_read_altbias         ! Bias treatment for altimeter
    USE obs_profiles_def         ! Profile data definitions
-   USE obs_profiles             ! Profile data storage
    USE obs_surf_def             ! Surface data definitions
-   USE obs_sla                  ! SLA data storage
-   USE obs_sst                  ! SST data storage
-   USE obs_seaice               ! Sea Ice data storage
    USE obs_types                ! Definitions for observation types
    USE mpp_map                  ! MPP mapping
@@ -55,66 +45,31 @@
       &   calc_date           ! Compute the date of a timestep
 
-   !! * Shared Module variables
-   LOGICAL, PUBLIC, PARAMETER :: &
-#if defined key_diaobs
-      & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics
-#else
-      & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics
-#endif
-
    !! * Module variables
-   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles
-   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles
-   LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set
-   LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set
-   LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles
-   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies 
-   LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files
-   LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files
-   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature
-   LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature
-   LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data
-   LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files
-   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration
-   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations
-   LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data 
-   LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data 
-   LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data 
-   LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data 
-   LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files
-   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
-   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
-   LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations
-   LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land
-   LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias  
-   LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files
-   LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations
-   LOGICAL, PUBLIC :: ln_sstbias     !: Logical switch for bias corection of SST 
+   LOGICAL, PUBLIC :: ln_diaobs   !: Logical switch for the obs operator
+   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs
    
-   REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS
-   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS
-  
-   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method
-   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method 
-
+   INTEGER :: nn_1dint       !: Vertical interpolation method
+   INTEGER :: nn_2dint       !: Horizontal interpolation method
    INTEGER, DIMENSION(imaxavtypes) :: &
-      & endailyavtypes !: ENACT data types which are daily average
-
-   INTEGER, PARAMETER :: MaxNumFiles = 1000
-   LOGICAL, DIMENSION(MaxNumFiles) :: &
-      & ln_profb_ena, & !: Is the feedback files from ENACT data ?
-   !                    !: If so use endailyavtypes
-      & ln_profb_enatim !: Change tim for 820 enact data set.
-
-   INTEGER, DIMENSION(MaxNumFiles), PUBLIC :: sstbias_type !SST bias type    
-      
-   LOGICAL, DIMENSION(MaxNumFiles) :: &
-      & ln_velfb_av   !: Is the velocity feedback files daily average?
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_enact     !: Profile data is ENACT so use endailyavtypes
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_velav     !: Velocity data is daily averaged
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_sstnight  !: SST observation corresponds to night mean
+      & nn_profdavtypes      !: Profile data types representing a daily average
+   INTEGER :: nproftypes     !: Number of profile obs types
+   INTEGER :: nsurftypes     !: Number of surface obs types
+   INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      & nvarsprof, &         !: Number of profile variables
+      & nvarssurf            !: Number of surface variables
+   INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      & nextrprof, &         !: Number of profile extra variables
+      & nextrsurf            !: Number of surface extra variables
+   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:),:: sstbias_type !SST bias type    
+   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
+      & surfdata, &          !: Initial surface data
+      & surfdataqc           !: Surface data after quality control
+   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
+      & profdata, &          !: Initial profile data
+      & profdataqc           !: Profile data after quality control
+
+   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
+      & cobstypesprof, &     !: Profile obs types
+      & cobstypessurf        !: Surface obs types
 
    !!----------------------------------------------------------------------
@@ -144,4 +99,5 @@
       !!        !  07-03  (K. Mogensen) General handling of profiles
       !!        !  14-08  (J.While) Incorporated SST bias correction  
+      !!        !  15-02  (M. Martin) Simplification of namelist and code
       !!----------------------------------------------------------------------
 
@@ -149,70 +105,74 @@
 
       !! * Local declarations
-      CHARACTER(len=128) :: enactfiles(MaxNumFiles)
-      CHARACTER(len=128) :: coriofiles(MaxNumFiles)
-      CHARACTER(len=128) :: profbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: sstfiles(MaxNumFiles)      
-      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: sstbias_files(MaxNumFiles) 
-      CHARACTER(len=128) :: slafilesact(MaxNumFiles)      
-      CHARACTER(len=128) :: slafilespas(MaxNumFiles)      
-      CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: seaicefiles(MaxNumFiles)           
-      CHARACTER(len=128) :: velcurfiles(MaxNumFiles)  
-      CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)    
-      CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
-      CHARACTER(LEN=128) :: reysstname
-      CHARACTER(LEN=12)  :: reysstfmt
-      CHARACTER(LEN=128) :: bias_file
-      CHARACTER(LEN=20)  :: datestr=" ", timestr=" "
-      NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       &
-         &            ln_sla, ln_sladt, ln_slafb,                     &
-         &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       &
-         &            enactfiles, coriofiles, profbfiles,             &
-         &            slafilesact, slafilespas, slafbfiles,           &
-         &            sstfiles, sstfbfiles,                           &
-         &            ln_seaice, seaicefiles,                         &
-         &            dobsini, dobsend, n1dint, n2dint,               &
-         &            nmsshc, mdtcorr, mdtcutoff,                     &
-         &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    &
-         &            ln_sstnight,                                    &
+      INTEGER, PARAMETER :: &
+         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+         & ifilesprof, &         ! Number of profile files
+         & ifilessurf            ! Number of surface files
+      INTEGER :: ios             ! Local integer output status for namelist read
+      INTEGER :: jtype           ! Counter for obs types
+      INTEGER :: jvar            ! Counter for variables
+      INTEGER :: jfile           ! Counter for files
+
+      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
+         & cn_profbfiles, &      ! T/S profile input filenames
+         & cn_sstfbfiles, &      ! Sea surface temperature input filenames
+         & cn_slafbfiles, &      ! Sea level anomaly input filenames
+         & cn_sicfbfiles, &      ! Seaice concentration input filenames
+         & cn_velfbfiles         ! Velocity profile input filenames
+      CHARACTER(LEN=128) :: &
+         & cn_altbiasfile        ! Altimeter bias input filename
+         & cn_sstbias_files      ! Altimeter bias input filenames
+      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
+         & clproffiles, &        ! Profile filenames
+         & clsurffiles           ! Surface filenames
+
+      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
+      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
+      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies 
+      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
+      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
+      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
+      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
+      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
+      LOGICAL :: ln_sstbias     !: Logical switch for bias corection of SST 
+      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
+      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
+      LOGICAL :: llvar1          ! Logical for profile variable 1
+      LOGICAL :: llvar2          ! Logical for profile variable 1
+      LOGICAL :: llnightav       ! Logical for calculating night-time averages
+
+      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
+      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
+      REAL(wp), POINTER, DIMENSION(:,:) :: &
+         & zglam1, &             ! Model longitudes for profile variable 1
+         & zglam2                ! Model longitudes for profile variable 2
+      REAL(wp), POINTER, DIMENSION(:,:) :: &
+         & zgphi1, &             ! Model latitudes for profile variable 1
+         & zgphi2                ! Model latitudes for profile variable 2
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zmask1, &             ! Model land/sea mask associated with variable 1
+         & zmask2                ! Model land/sea mask associated with variable 2
+
+      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
+         &            ln_sst, ln_sic, ln_vel3d,                       &
+         &            ln_altbias, ln_nea, ln_grid_global,             &
          &            ln_grid_search_lookup,                          &
-         &            grid_search_file, grid_search_res,              &
-         &            ln_grid_global, bias_file, ln_altbias,          &
-         &            endailyavtypes, ln_s_at_t, ln_profb_ena,        &
-         &            ln_vel3d, ln_velavcur, velavcurfiles,           &
-         &            ln_velhrcur, velhrcurfiles,                     &
-         &            ln_velavadcp, velavadcpfiles,                   &
-         &            ln_velhradcp, velhradcpfiles,                   &
-         &            ln_velfb, velfbfiles, ln_velfb_av,              &
-         &            ln_profb_enatim, ln_ignmis, ln_cl4,             &
-         &            ln_sstbias, sstbias_files          
-
-      INTEGER :: jprofset
-      INTEGER :: jveloset
-      INTEGER :: jvar
-      INTEGER :: jnumenact
-      INTEGER :: jnumcorio
-      INTEGER :: jnumprofb
-      INTEGER :: jnumslaact
-      INTEGER :: jnumslapas
-      INTEGER :: jnumslafb
-      INTEGER :: jnumsst
-      INTEGER :: jnumsstfb
-      INTEGER :: jnumsstbias 
-      INTEGER :: jnumseaice
-      INTEGER :: jnumvelavcur
-      INTEGER :: jnumvelhrcur  
-      INTEGER :: jnumvelavadcp
-      INTEGER :: jnumvelhradcp   
-      INTEGER :: jnumvelfb
-      INTEGER :: ji
-      INTEGER :: jset
-      INTEGER :: ios                 ! Local integer output status for namelist read
-      LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
+         &            ln_ignmis, ln_s_at_t, ln_sstnight,              &
+         &            cn_profbfiles, cn_slafbfiles,                   &
+         &            cn_sstfbfiles, cn_sicfbfiles,                   &
+         &            cn_velfbfiles, cn_altbiasfile,                  &
+         &            cn_gridsearchfile, rn_gridsearchres,            &
+         &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     &
+         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
+         &            nn_profdavtypes, ln_sstbias, sstbias_files
+
+      INTEGER :: jnumsstbias  !TG - Is this still needed
+      CALL wrk_alloc( jpi, jpj, zglam1 )
+      CALL wrk_alloc( jpi, jpj, zglam2 )
+      CALL wrk_alloc( jpi, jpj, zgphi1 )
+      CALL wrk_alloc( jpi, jpj, zgphi2 )
+      CALL wrk_alloc( jpi, jpj, jpk, zmask1 )
+      CALL wrk_alloc( jpi, jpj, jpk, zmask2 )
 
       !-----------------------------------------------------------------------
@@ -221,83 +181,41 @@
       
       !Initalise all values in namelist arrays
-      enactfiles(:) = ''
-      coriofiles(:) = ''
-      profbfiles(:) = ''
-      slafilesact(:) = ''
-      slafilespas(:) = ''
-      slafbfiles(:) = ''
-      sstfiles(:)   = ''
-      sstfbfiles(:) = ''
-      seaicefiles(:) = ''
-      velcurfiles(:) = ''
-      veladcpfiles(:) = ''
-      velavcurfiles(:) = ''
-      velhrcurfiles(:) = ''
-      velavadcpfiles(:) = ''
-      velhradcpfiles(:) = ''
-      velfbfiles(:) = ''
-      velcurfiles(:) = ''
-      veladcpfiles(:) = ''
-      endailyavtypes(:) = -1
-      endailyavtypes(1) = 820
-      ln_profb_ena(:) = .FALSE.
-      ln_profb_enatim(:) = .TRUE.
-      ln_velfb_av(:) = .FALSE.
-      ln_ignmis = .FALSE.      
-
-      CALL ini_date( dobsini )
-      CALL fin_date( dobsend )
- 
-      ! Read Namelist namobs : control observation diagnostics
-      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation
+      ALLOCATE(sstbias_type(jpmaxnumfiles))
+      ! Some namelist arrays need initialising
+      cn_profbfiles(:) = ''
+      cn_slafbfiles(:) = ''
+      cn_sstfbfiles(:) = ''
+      cn_sicfbfiles(:) = ''
+      cn_velfbfiles(:) = ''
+      cn_sstbias_files(:) = ''
+      nn_profdavtypes(:) = -1
+
+      CALL ini_date( rn_dobsini )
+      CALL fin_date( rn_dobsend )
+
+      ! Read namelist namobs : control observation diagnostics
+      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
       READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
 
-      REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation
+      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
       READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namobs )
 
-      ! Count number of files for each type
-      IF (ln_ena) THEN
-         lmask(:) = .FALSE.
-         WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
-         jnumenact = COUNT(lmask)
-      ENDIF
-      IF (ln_cor) THEN
-         lmask(:) = .FALSE.
-         WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
-         jnumcorio = COUNT(lmask)
-      ENDIF
-      IF (ln_profb) THEN
-         lmask(:) = .FALSE.
-         WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumprofb = COUNT(lmask)
-      ENDIF
-      IF (ln_sladt) THEN
-         lmask(:) = .FALSE.
-         WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
-         jnumslaact = COUNT(lmask)
-         lmask(:) = .FALSE.
-         WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
-         jnumslapas = COUNT(lmask)
-      ENDIF
-      IF (ln_slafb) THEN
-         lmask(:) = .FALSE.
-         WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumslafb = COUNT(lmask)
-         lmask(:) = .FALSE.
-      ENDIF
-      IF (ln_ghrsst) THEN
-         lmask(:) = .FALSE.
-         WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
-         jnumsst = COUNT(lmask)
-      ENDIF      
-      IF (ln_sstfb) THEN
-         lmask(:) = .FALSE.
-         WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumsstfb = COUNT(lmask)
-         lmask(:) = .FALSE.
-      ENDIF
+      IF ( .NOT. ln_diaobs ) THEN
+         IF(lwp) WRITE(numout,cform_war)
+         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs'
+         RETURN
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Set up list of observation types to be used
+      ! and the files associated with each type
+      !-----------------------------------------------------------------------
+
+      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
+      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) )
+
       IF (ln_sstbias) THEN 
          lmask(:) = .FALSE. 
@@ -306,37 +224,88 @@
          lmask(:) = .FALSE. 
       ENDIF      
-      IF (ln_seaice) THEN
-         lmask(:) = .FALSE.
-         WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
-         jnumseaice = COUNT(lmask)
-      ENDIF
-      IF (ln_velavcur) THEN
-         lmask(:) = .FALSE.
-         WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelavcur = COUNT(lmask)
-      ENDIF
-      IF (ln_velhrcur) THEN
-         lmask(:) = .FALSE.
-         WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelhrcur = COUNT(lmask)
-      ENDIF
-      IF (ln_velavadcp) THEN
-         lmask(:) = .FALSE.
-         WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelavadcp = COUNT(lmask)
-      ENDIF
-      IF (ln_velhradcp) THEN
-         lmask(:) = .FALSE.
-         WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelhradcp = COUNT(lmask)
-      ENDIF
-      IF (ln_velfb) THEN
-         lmask(:) = .FALSE.
-         WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
-         jnumvelfb = COUNT(lmask)
-         lmask(:) = .FALSE.
-      ENDIF
-      
-      ! Control print
+
+      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
+         IF(lwp) WRITE(numout,cform_war)
+         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
+            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', &
+            &                    ' are set to .FALSE. so turning off calls to dia_obs'
+         nwarn = nwarn + 1
+         ln_diaobs = .FALSE.
+         RETURN
+      ENDIF
+
+      IF ( nproftypes > 0 ) THEN
+
+         ALLOCATE( cobstypesprof(nproftypes) )
+         ALLOCATE( ifilesprof(nproftypes) )
+         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
+
+         jtype = 0
+         IF (ln_t3d .OR. ln_s3d) THEN
+            jtype = jtype + 1
+            clproffiles(jtype,:) = cn_profbfiles(:)
+            cobstypesprof(jtype) = 'prof  '
+            ifilesprof(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clproffiles(jtype,jfile)) /= '' ) &
+                  ifilesprof(jtype) = ifilesprof(jtype) + 1
+            END DO
+         ENDIF
+         IF (ln_vel3d) THEN
+            jtype = jtype + 1
+            clproffiles(jtype,:) = cn_velfbfiles(:)
+            cobstypesprof(jtype) = 'vel   '
+            ifilesprof(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clproffiles(jtype,jfile)) /= '' ) &
+                  ifilesprof(jtype) = ifilesprof(jtype) + 1
+            END DO
+         ENDIF
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         ALLOCATE( cobstypessurf(nsurftypes) )
+         ALLOCATE( ifilessurf(nsurftypes) )
+         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
+
+         jtype = 0
+         IF (ln_sla) THEN
+            jtype = jtype + 1
+            clsurffiles(jtype,:) = cn_slafbfiles(:)
+            cobstypessurf(jtype) = 'sla   '
+            ifilessurf(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
+                  ifilessurf(jtype) = ifilessurf(jtype) + 1
+            END DO
+         ENDIF
+         IF (ln_sst) THEN
+            jtype = jtype + 1
+            clsurffiles(jtype,:) = cn_sstfbfiles(:)
+            cobstypessurf(jtype) = 'sst   '
+            ifilessurf(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
+                  ifilessurf(jtype) = ifilessurf(jtype) + 1
+            END DO
+         ENDIF
+#if defined key_lim2 || defined key_lim3
+         IF (ln_sic) THEN
+            jtype = jtype + 1
+            clsurffiles(jtype,:) = cn_sicfbfiles(:)
+            cobstypessurf(jtype) = 'sic   '
+            ifilessurf(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
+                  ifilessurf(jtype) = ifilessurf(jtype) + 1
+            END DO
+         ENDIF
+#endif
+
+      ENDIF
+
+      !Write namelist settings to stdout
       IF(lwp) THEN
          WRITE(numout,*)
@@ -344,139 +313,55 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
-         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d
-         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d
-         WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena
-         WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor
-         WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb
-         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla
-         WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt
-         WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb
-         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh
-         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst
-         WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst    
-         WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst
-         WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb
+         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
+         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
+         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
+         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
+         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
+         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
+         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global
          WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias 
-         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight
-         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss
-         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice
-         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d
-         WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur
-         WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur
-         WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp
-         WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp
-         WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb
-         WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
-         WRITE(numout,*) &
-   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup
+         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup
          IF (ln_grid_search_lookup) &
-            WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file
-         IF (ln_ena) THEN
-            DO ji = 1, jnumenact
-               WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', &
-                  TRIM(enactfiles(ji))
+            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
+         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
+         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
+         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
+         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint
+         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
+         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
+         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
+         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
+         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
+         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
+         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
+         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
+         WRITE(numout,*) '          Number of profile obs types: ',nproftypes
+
+         IF ( nproftypes > 0 ) THEN
+            DO jtype = 1, nproftypes
+               DO jfile = 1, ifilesprof(jtype)
+                  WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', &
+                     TRIM(clproffiles(jtype,jfile))
+               END DO
             END DO
          ENDIF
-         IF (ln_cor) THEN
-            DO ji = 1, jnumcorio
-               WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', &
-                  TRIM(coriofiles(ji))
+
+         WRITE(numout,*)'          Number of surface obs types: ',nsurftypes
+         IF ( nsurftypes > 0 ) THEN
+            DO jtype = 1, nsurftypes
+               DO jfile = 1, ifilessurf(jtype)
+                  WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', &
+                     TRIM(clsurffiles(jtype,jfile))
+               END DO
             END DO
          ENDIF
-         IF (ln_profb) THEN
-            DO ji = 1, jnumprofb
-               IF (ln_profb_ena(ji)) THEN
-                  WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', &
-                     TRIM(profbfiles(ji))
-               ELSE
-                  WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', &
-                     TRIM(profbfiles(ji))
-               ENDIF
-               WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji)
-            END DO
-         ENDIF
-         IF (ln_sladt) THEN
-            DO ji = 1, jnumslaact
-               WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', &
-                  TRIM(slafilesact(ji))
-            END DO
-            DO ji = 1, jnumslapas
-               WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', &
-                  TRIM(slafilespas(ji))
-            END DO
-         ENDIF
-         IF (ln_slafb) THEN
-            DO ji = 1, jnumslafb
-               WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', &
-                  TRIM(slafbfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_ghrsst) THEN
-            DO ji = 1, jnumsst
-               WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', &
-                  TRIM(sstfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_sstfb) THEN
-            DO ji = 1, jnumsstfb
-               WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', &
-                  TRIM(sstfbfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_seaice) THEN
-            DO ji = 1, jnumseaice
-               WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', &
-                  TRIM(seaicefiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velavcur) THEN
-            DO ji = 1, jnumvelavcur
-               WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', &
-                  TRIM(velavcurfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velhrcur) THEN
-            DO ji = 1, jnumvelhrcur
-               WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', &
-                  TRIM(velhrcurfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velavadcp) THEN
-            DO ji = 1, jnumvelavadcp
-               WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', &
-                  TRIM(velavadcpfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velhradcp) THEN
-            DO ji = 1, jnumvelhradcp
-               WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', &
-                  TRIM(velhradcpfiles(ji))
-            END DO
-         ENDIF
-         IF (ln_velfb) THEN
-            DO ji = 1, jnumvelfb
-               IF (ln_velfb_av(ji)) THEN
-                  WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', &
-                     TRIM(velfbfiles(ji))
-               ELSE
-                  WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', &
-                     TRIM(velfbfiles(ji))
-               ENDIF
-            END DO
-         ENDIF
-         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini
-         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend
-         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint
-         WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint
-         WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea
-         WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc
-         WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr
-         WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff
-         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias
-         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis
-         WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes
-
-      ENDIF
-      
+         WRITE(numout,*) '~~~~~~~~~~~~'
+
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Obs operator parameter checking and initialisations
+      !-----------------------------------------------------------------------
+
       IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
          CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
@@ -484,313 +369,101 @@
       ENDIF
 
-      CALL obs_typ_init
-      
-      CALL mppmap_init
-      
-      ! Parameter control
-#if defined key_diaobs
-      IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
-         & ( .NOT. ln_vel3d ).AND.                                         &
-         & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
-         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN
-         IF(lwp) WRITE(numout,cform_war)
-         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
-            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.'
-         nwarn = nwarn + 1
-      ENDIF
-#endif
-
-      CALL obs_grid_setup( )
-      IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
+      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
          CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
             &                    ' is not available')
       ENDIF
-      IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
+
+      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN
          CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
             &                    ' is not available')
       ENDIF
 
+      CALL obs_typ_init
+
+      CALL mppmap_init
+
+      CALL obs_grid_setup( )
+
       !-----------------------------------------------------------------------
       ! Depending on switches read the various observation types
       !-----------------------------------------------------------------------
-      !  - Temperature/salinity profiles
-
-      IF ( ln_t3d .OR. ln_s3d ) THEN
-
-         ! Set the number of variables for profiles to 2 (T and S)
-         nprofvars = 2
-         ! Set the number of extra variables for profiles to 1 (insitu temp).
-         nprofextr = 1
-
-         ! Count how may insitu data sets we have and allocate data.
-         jprofset = 0
-         IF ( ln_ena ) jprofset = jprofset + 1
-         IF ( ln_cor ) jprofset = jprofset + 1
-         IF ( ln_profb ) jprofset = jprofset + jnumprofb
-         nprofsets = jprofset
-         IF ( nprofsets > 0 ) THEN
-            ALLOCATE(ld_enact(nprofsets))
-            ALLOCATE(profdata(nprofsets))
-            ALLOCATE(prodatqc(nprofsets))
-         ENDIF
-
-         jprofset = 0
-          
-         ! ENACT insitu data
-
-         IF ( ln_ena ) THEN
-
-            jprofset = jprofset + 1
-            
-            ld_enact(jprofset) = .TRUE.
-
-            CALL obs_rea_pro_dri( 1, profdata(jprofset),          &
-               &                  jnumenact, enactfiles(1:jnumenact), &
-               &                  nprofvars, nprofextr,        &
-               &                  nitend-nit000+2,             &
-               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
-               &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
-               &                  kdailyavtypes = endailyavtypes )
-
-            DO jvar = 1, 2
-
-               CALL obs_prof_staend( profdata(jprofset), jvar )
-
+
+      IF ( nproftypes > 0 ) THEN
+
+         ALLOCATE(profdata(nproftypes))
+         ALLOCATE(profdataqc(nproftypes))
+         ALLOCATE(nvarsprof(nproftypes))
+         ALLOCATE(nextrprof(nproftypes))
+
+         DO jtype = 1, nproftypes
+
+            nvarsprof(jtype) = 2
+            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
+               nextrprof(jtype) = 1
+               llvar1 = ln_t3d
+               llvar2 = ln_s3d
+               zglam1 = glamt
+               zgphi1 = gphit
+               zmask1 = tmask
+               zglam2 = glamt
+               zgphi2 = gphit
+               zmask2 = tmask
+            ENDIF
+            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
+               nextrprof(jtype) = 2
+               llvar1 = ln_vel3d
+               llvar2 = ln_vel3d
+               zglam1 = glamu
+               zgphi1 = gphiu
+               zmask1 = umask
+               zglam2 = glamv
+               zgphi2 = gphiv
+               zmask2 = vmask
+            ENDIF
+
+            !Read in profile or profile obs types
+            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       &
+               &               clproffiles(jtype,1:ifilesprof(jtype)), &
+               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
+               &               rn_dobsini, rn_dobsend, llvar1, llvar2, &
+               &               ln_ignmis, ln_s_at_t, .FALSE., &
+               &               kdailyavtypes = nn_profdavtypes )
+
+            DO jvar = 1, nvarsprof(jtype)
+               CALL obs_prof_staend( profdata(jtype), jvar )
             END DO
 
-            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-               &              ln_t3d, ln_s3d, ln_nea, &
-               &              kdailyavtypes=endailyavtypes )
-            
-         ENDIF
-
-         ! Coriolis insitu data
-
-         IF ( ln_cor ) THEN
-           
-            jprofset = jprofset + 1
-
-            ld_enact(jprofset) = .FALSE.
-
-            CALL obs_rea_pro_dri( 2, profdata(jprofset),          &
-               &                  jnumcorio, coriofiles(1:jnumcorio), &
-               &                  nprofvars, nprofextr,        &
-               &                  nitend-nit000+2,             &
-               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
-               &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
-
-            DO jvar = 1, 2
-
-               CALL obs_prof_staend( profdata(jprofset), jvar )
-
-            END DO
-
-            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-                 &            ln_t3d, ln_s3d, ln_nea )
-            
-         ENDIF
- 
-         ! Feedback insitu data
-
-         IF ( ln_profb ) THEN
-           
-            DO jset = 1, jnumprofb
-               
-               jprofset = jprofset + 1
-               ld_enact (jprofset) = ln_profb_ena(jset)
-
-               CALL obs_rea_pro_dri( 0, profdata(jprofset),          &
-                  &                  1, profbfiles(jset:jset), &
-                  &                  nprofvars, nprofextr,        &
-                  &                  nitend-nit000+2,             &
-                  &                  dobsini, dobsend, ln_t3d, ln_s3d, &
-                  &                  ln_ignmis, ln_s_at_t, &
-                  &                  ld_enact(jprofset).AND.&
-                  &                  ln_profb_enatim(jset), &
-                  &                  .FALSE., kdailyavtypes = endailyavtypes )
-               
-               DO jvar = 1, 2
-                  
-                  CALL obs_prof_staend( profdata(jprofset), jvar )
-                  
-               END DO
-               
-               IF ( ld_enact(jprofset) ) THEN
-                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-                     &              ln_t3d, ln_s3d, ln_nea, &
-                     &              kdailyavtypes = endailyavtypes )
-               ELSE
-                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
-                     &              ln_t3d, ln_s3d, ln_nea )
-               ENDIF
-               
-            END DO
-
-         ENDIF
-
-      ENDIF
-
-      !  - Sea level anomalies
-      IF ( ln_sla ) THEN
-        ! Set the number of variables for sla to 1
-         nslavars = 1
-
-         ! Set the number of extra variables for sla to 2
-         nslaextr = 2
-         
-         ! Set the number of sla data sets to 2
-         nslasets = 0
-         IF ( ln_sladt ) THEN
-            nslasets = nslasets + 2
-         ENDIF
-         IF ( ln_slafb ) THEN
-            nslasets = nslasets + jnumslafb
-         ENDIF
-         
-         ALLOCATE(sladata(nslasets))
-         ALLOCATE(sladatqc(nslasets))
-         sladata(:)%nsurf=0
-         sladatqc(:)%nsurf=0
-
-         nslasets = 0
-
-         ! AVISO SLA data
-
-         IF ( ln_sladt ) THEN
-
-            ! Active SLA observations
-            
-            nslasets = nslasets + 1
-            
-            CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
-               &              slafilesact(1:jnumslaact), &
-               &              nslavars, nslaextr, nitend-nit000+2, &
-               &              dobsini, dobsend, ln_ignmis, .FALSE. )
-            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
-               &              ln_sla, ln_nea )
-            
-            ! Passive SLA observations
-            
-            nslasets = nslasets + 1
-            
-            CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
-               &              slafilespas(1:jnumslapas), &
-               &              nslavars, nslaextr, nitend-nit000+2, &
-               &              dobsini, dobsend, ln_ignmis, .FALSE. )
-            
-            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
-               &              ln_sla, ln_nea )
-
-         ENDIF
-         
-         ! Feedback SLA data
-
-         IF ( ln_slafb ) THEN
-
-            DO jset = 1, jnumslafb
-            
-               nslasets = nslasets + 1
-            
-               CALL obs_rea_sla( 0, sladata(nslasets), 1, &
-                  &              slafbfiles(jset:jset), &
-                  &              nslavars, nslaextr, nitend-nit000+2, &
-                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
-               CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
-                  &              ln_sla, ln_nea )
-
-            END DO               
-
-         ENDIF
-         
-         CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
-            
-         ! read in altimeter bias
-         
-         IF ( ln_altbias ) THEN     
-            CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
-         ENDIF
-     
-      ENDIF
-
-      !  - Sea surface height
-      IF ( ln_ssh ) THEN
-         IF(lwp) WRITE(numout,*) ' SSH currently not available'
-      ENDIF
-
-      !  - Sea surface temperature
-      IF ( ln_sst ) THEN
-
-         ! Set the number of variables for sst to 1
-         nsstvars = 1
-
-         ! Set the number of extra variables for sst to 0
-         nsstextr = 0
-
-         nsstsets = 0
-
-         IF (ln_reysst) nsstsets = nsstsets + 1
-         IF (ln_ghrsst) nsstsets = nsstsets + 1
-         IF ( ln_sstfb ) THEN
-            nsstsets = nsstsets + jnumsstfb
-         ENDIF
-
-         ALLOCATE(sstdata(nsstsets))
-         ALLOCATE(sstdatqc(nsstsets))
-         ALLOCATE(ld_sstnight(nsstsets))
-         sstdata(:)%nsurf=0
-         sstdatqc(:)%nsurf=0    
-         ld_sstnight(:)=.false.
-
-         nsstsets = 0
-
-         IF (ln_reysst) THEN
-
-            nsstsets = nsstsets + 1
-
-            ld_sstnight(nsstsets) = ln_sstnight
-
-            CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
-               &                  nsstvars, nsstextr, &
-               &                  nitend-nit000+2, dobsini, dobsend )
-            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
-               &              ln_nea )
-
-        ENDIF
-        
-        IF (ln_ghrsst) THEN
-        
-            nsstsets = nsstsets + 1
-
-            ld_sstnight(nsstsets) = ln_sstnight
-          
-            CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
-               &              sstfiles(1:jnumsst), &
-               &              nsstvars, nsstextr, nitend-nit000+2, &
-               &              dobsini, dobsend, ln_ignmis, .FALSE. )
-            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
-               &              ln_nea )
-
-        ENDIF
-               
-         ! Feedback SST data
-
-         IF ( ln_sstfb ) THEN
-
-            DO jset = 1, jnumsstfb
-            
-               nsstsets = nsstsets + 1
-
-               ld_sstnight(nsstsets) = ln_sstnight
-            
-               CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
-                  &              sstfbfiles(jset:jset), &
-                  &              nsstvars, nsstextr, nitend-nit000+2, &
-                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
-               CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
-                  &              ln_sst, ln_nea )
-
-            END DO               
-
-         ENDIF
+            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
+               &               llvar1, llvar2, &
+               &               jpi, jpj, jpk, &
+               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
+               &               ln_nea, kdailyavtypes = nn_profdavtypes )
+
+         END DO
+
+         DEALLOCATE( ifilesprof, clproffiles )
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         ALLOCATE(surfdata(nsurftypes))
+         ALLOCATE(surfdataqc(nsurftypes))
+         ALLOCATE(nvarssurf(nsurftypes))
+         ALLOCATE(nextrsurf(nsurftypes))
+
+         DO jtype = 1, nsurftypes
+
+            nvarssurf(jtype) = 1
+            nextrsurf(jtype) = 0
+            llnightav = .FALSE.
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
+            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight
+
+            !Read in surface obs types
+            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
+               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
+               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
+               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav )
          
          !Read in bias field and correct SST.
@@ -804,195 +477,24 @@
          ENDIF
          
-      ENDIF
-
-      !  - Sea surface salinity
-      IF ( ln_sss ) THEN
-         IF(lwp) WRITE(numout,*) ' SSS currently not available'
-      ENDIF
-
-      !  - Sea Ice Concentration
-      
-      IF ( ln_seaice ) THEN
-
-         ! Set the number of variables for seaice to 1
-         nseaicevars = 1
-
-         ! Set the number of extra variables for seaice to 0
-         nseaiceextr = 0
-         
-         ! Set the number of data sets to 1
-         nseaicesets = 1
-
-         ALLOCATE(seaicedata(nseaicesets))
-         ALLOCATE(seaicedatqc(nseaicesets))
-         seaicedata(:)%nsurf=0
-         seaicedatqc(:)%nsurf=0
-
-         CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
-            &                 seaicefiles(1:jnumseaice), &
-            &                 nseaicevars, nseaiceextr, nitend-nit000+2, &
-            &                 dobsini, dobsend, ln_ignmis, .FALSE. )
-
-         CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
-            &                 ln_seaice, ln_nea )
- 
-      ENDIF
-
-      IF (ln_vel3d) THEN
-
-         ! Set the number of variables for profiles to 2 (U and V)
-         nvelovars = 2
-
-         ! Set the number of extra variables for profiles to 2 to store 
-         ! rotation parameters
-         nveloextr = 2
-
-         jveloset = 0
-         
-         IF ( ln_velavcur ) jveloset = jveloset + 1
-         IF ( ln_velhrcur ) jveloset = jveloset + 1
-         IF ( ln_velavadcp ) jveloset = jveloset + 1
-         IF ( ln_velhradcp ) jveloset = jveloset + 1
-         IF (ln_velfb) jveloset = jveloset + jnumvelfb
-
-         nvelosets = jveloset
-         IF ( nvelosets > 0 ) THEN
-            ALLOCATE( velodata(nvelosets) )
-            ALLOCATE( veldatqc(nvelosets) )
-            ALLOCATE( ld_velav(nvelosets) )
-         ENDIF
-         
-         jveloset = 0
-         
-         ! Daily averaged data
-
-         IF ( ln_velavcur ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .TRUE.
-            
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
-               &                  velavcurfiles(1:jnumvelavcur), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhrcur ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .FALSE.
-               
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
-               &                  velhrcurfiles(1:jnumvelhrcur), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         ! Daily averaged data
-
-         IF ( ln_velavadcp ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .TRUE.
-            
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
-               &                  velavadcpfiles(1:jnumvelavadcp), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhradcp ) THEN
-            
-            jveloset = jveloset + 1
-            
-            ld_velav(jveloset) = .FALSE.
-               
-            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
-               &                  velhradcpfiles(1:jnumvelhradcp), &
-               &                  nvelovars, nveloextr, &
-               &                  nitend-nit000+2,              &
-               &                  dobsini, dobsend, ln_ignmis, &
-               &                  ld_velav(jveloset), &
-               &                  .FALSE. )
-            
-            DO jvar = 1, 2
-               CALL obs_prof_staend( velodata(jveloset), jvar )
-            END DO
-            
-            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-            
-         ENDIF
-
-         IF ( ln_velfb ) THEN
-
-            DO jset = 1, jnumvelfb
-            
-               jveloset = jveloset + 1
-
-               ld_velav(jveloset) = ln_velfb_av(jset)
-               
-               CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
-                  &                  velfbfiles(jset:jset), &
-                  &                  nvelovars, nveloextr, &
-                  &                  nitend-nit000+2,              &
-                  &                  dobsini, dobsend, ln_ignmis, &
-                  &                  ld_velav(jveloset), &
-                  &                  .FALSE. )
-               
-               DO jvar = 1, 2
-                  CALL obs_prof_staend( velodata(jveloset), jvar )
-               END DO
-               
-               CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
-                  &              ln_vel3d, ln_nea, ld_velav(jveloset) )
-
-
-            END DO
-            
-         ENDIF
-
-      ENDIF
-     
+            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea )
+
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
+               CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint )
+               IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile )
+            ENDIF
+
+         END DO
+
+         DEALLOCATE( ifilessurf, clsurffiles )
+
+      ENDIF
+
+      CALL wrk_dealloc( jpi, jpj, zglam1 )
+      CALL wrk_dealloc( jpi, jpj, zglam2 )
+      CALL wrk_dealloc( jpi, jpj, zgphi1 )
+      CALL wrk_dealloc( jpi, jpj, zgphi2 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )
+
    END SUBROUTINE dia_obs_init
 
@@ -1004,13 +506,9 @@
       !!
       !! ** Method  : Call the observation operators on each time step to
-      !!              compute the model equivalent of the following date:
-      !!               - T profiles
-      !!               - S profiles
-      !!               - Sea surface height (referenced to a mean)
-      !!               - Sea surface temperature
-      !!               - Sea surface salinity
-      !!               - Velocity component (U,V) profiles
-      !!
-      !! ** Action  : 
+      !!              compute the model equivalent of the following data:
+      !!               - Profile data, currently T/S or U/V
+      !!               - Surface data, currently SST, SLA or sea-ice concentration.
+      !!
+      !! ** Action  :
       !!
       !! History :
@@ -1023,9 +521,8 @@
       !!        !  14-08  (J. While) observation operator for profiles in 
       !!                             generalised vertical coordinates
+      !!        !  15-08  (M. Martin) Combined surface/profile routines.
       !!----------------------------------------------------------------------
       !! * Modules used
       USE dom_oce, ONLY : &             ! Ocean space and time domain variables
-         & rdt,           &                       
-         & gdept_1d,       &             
 #if defined key_vvl 
          & gdept_n,       &
@@ -1033,5 +530,4 @@
          & gdept_1d,      &
 #endif                                        
-         & tmask, umask, vmask                            
       USE phycst, ONLY : &              ! Physical constants
          & rday                         
@@ -1039,11 +535,13 @@
          & tsn,  &             
          & un, vn,  &
+      USE phycst, ONLY : &         ! Physical constants
+         & rday
          & sshn
 #if defined  key_lim3
-      USE ice, ONLY : &                     ! LIM Ice model variables
+      USE ice, ONLY : &            ! LIM3 Ice model variables
          & frld
 #endif
 #if defined key_lim2
-      USE ice_2, ONLY : &                     ! LIM Ice model variables
+      USE ice_2, ONLY : &          ! LIM2 Ice model variables
          & frld
 #endif
@@ -1051,18 +549,38 @@
 
       !! * Arguments
-      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
+      INTEGER, INTENT(IN) :: kstp  ! Current timestep
       !! * Local declarations
-      INTEGER :: idaystp                ! Number of timesteps per day
-      INTEGER :: jprofset               ! Profile data set loop variable
-      INTEGER :: jslaset                ! SLA data set loop variable
-      INTEGER :: jsstset                ! SST data set loop variable
-      INTEGER :: jseaiceset             ! sea ice data set loop variable
-      INTEGER :: jveloset               ! velocity profile data loop variable
-      INTEGER :: jvar                   ! Variable number    
+      INTEGER :: idaystp           ! Number of timesteps per day
+      INTEGER :: jtype             ! Data loop variable
+      INTEGER :: jvar              ! Variable number
+      INTEGER :: ji, jj            ! Loop counters
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zprofvar1, &            ! Model values for 1st variable in a prof ob
+         & zprofvar2               ! Model values for 2nd variable in a prof ob
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zprofmask1, &           ! Mask associated with zprofvar1
+         & zprofmask2              ! Mask associated with zprofvar2
+      REAL(wp), POINTER, DIMENSION(:,:) :: &
+         & zsurfvar                ! Model values equivalent to surface ob.
+      REAL(wp), POINTER, DIMENSION(:,:) :: &
+         & zglam1,    &            ! Model longitudes for prof variable 1
+         & zglam2,    &            ! Model longitudes for prof variable 2
+         & zgphi1,    &            ! Model latitudes for prof variable 1
+         & zgphi2                  ! Model latitudes for prof variable 2
 #if ! defined key_lim2 && ! defined key_lim3
-      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
+      REAL(wp), POINTER, DIMENSION(:,:) :: frld
 #endif
-      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
- 
+      LOGICAL :: llnightav        ! Logical for calculating night-time average
+
+      !Allocate local work arrays
+      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
+      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
+      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
+      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
+      CALL wrk_alloc( jpi, jpj, zsurfvar )
+      CALL wrk_alloc( jpi, jpj, zglam1 )
+      CALL wrk_alloc( jpi, jpj, zglam2 )
+      CALL wrk_alloc( jpi, jpj, zgphi1 )
+      CALL wrk_alloc( jpi, jpj, zgphi2 )
 #if ! defined key_lim2 && ! defined key_lim3
       CALL wrk_alloc(jpi,jpj,frld) 
@@ -1084,7 +602,8 @@
 #endif
       !-----------------------------------------------------------------------
-      ! Depending on switches call various observation operators
-      !-----------------------------------------------------------------------
-
+      ! Call the profile and surface observation operators
+      !-----------------------------------------------------------------------
+
+<<<<<<< .working
       !  - Temperature/salinity profiles
       IF ( ln_t3d .OR. ln_s3d ) THEN
@@ -1121,56 +640,111 @@
          END DO
       ENDIF
-
-      !  - Sea surface anomaly
-      IF ( ln_sla ) THEN
-         DO jslaset = 1, nslasets
-            CALL obs_sla_opt( sladatqc(jslaset),            &
-               &              kstp, jpi, jpj, nit000, sshn, &
-               &              tmask(:,:,1), n2dint )
-         END DO         
-      ENDIF
-
-      !  - Sea surface temperature
-      IF ( ln_sst ) THEN
-         DO jsstset = 1, nsstsets
-            CALL obs_sst_opt( sstdatqc(jsstset),                &
-               &              kstp, jpi, jpj, nit000, idaystp,  &
-               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  &
-               &              n2dint, ld_sstnight(jsstset) )
+=======
+      IF ( nproftypes > 0 ) THEN
+>>>>>>> .merge-right.r5997
+
+         DO jtype = 1, nproftypes
+
+            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
+            CASE('prof')
+               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
+               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zprofmask2(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zglam2(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+               zgphi2(:,:) = gphit(:,:)
+            CASE('vel')
+               zprofvar1(:,:,:) = un(:,:,:)
+               zprofvar2(:,:,:) = vn(:,:,:)
+               zprofmask1(:,:,:) = umask(:,:,:)
+               zprofmask2(:,:,:) = vmask(:,:,:)
+               zglam1(:,:) = glamu(:,:)
+               zglam2(:,:) = glamv(:,:)
+               zgphi1(:,:) = gphiu(:,:)
+               zgphi2(:,:) = gphiv(:,:)
+            END SELECT
+
+            IF( ln_zco .OR. ln_zps ) THEN 
+               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
+                  &               nit000, idaystp,                         &
+                  &               zprofvar1, zprofvar2,                    &
+                  &               gdept_1d, zprofmask1, zprofmask2,        &
+                  &               zglam1, zglam2, zgphi1, zgphi2,          &
+                  &               nn_1dint, nn_2dint,                      &
+                  &               kdailyavtypes = nn_profdavtypes )
+            ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof')
+               !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION
+               CALL obs_pro_sco_opt( prodatqc(jtype),                    & 
+                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
+                  &              zprofvar1, zprofvar2,                   & 
+                  &              fsdept(:,:,:), fsdepw(:,:,:),           &
+                  &              tmask, nn_1dint, nn_2dint,              & 
+                  &              kdailyavtypes = nn_profdavtypes ) 
+            ELSE
+               ctl_stop('DIA_OBS: Generalised vertical interpolation not'// &
+                         'yet working for velocity date (turn off velocity observations')
+            ENDIF
+
          END DO
-      ENDIF
-
-      !  - Sea surface salinity
-      IF ( ln_sss ) THEN
-         IF(lwp) WRITE(numout,*) ' SSS currently not available'
-      ENDIF
-
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         DO jtype = 1, nsurftypes
+
+            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
+            CASE('sst')
+               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
+               llnightav = ln_sstnight
+            CASE('sla')
+               zsurfvar(:,:) = sshn(:,:)
+               llnightav = .FALSE.
 #if defined key_lim2 || defined key_lim3
-      IF ( ln_seaice ) THEN
-         DO jseaiceset = 1, nseaicesets
-            CALL obs_seaice_opt( seaicedatqc(jseaiceset),      &
-               &              kstp, jpi, jpj, nit000, 1.-frld, &
-               &              tmask(:,:,1), n2dint )
+            CASE('sic')
+               IF ( kstp == 0 ) THEN
+                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
+                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
+                        &           'time-step but some obs are valid then.' )
+                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
+                        &           ' sea-ice obs will be missed'
+                  ENDIF
+                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
+                     &                        surfdataqc(jtype)%nsstp(1)
+                  CYCLE
+               ELSE
+                  zsurfvar(:,:) = 1._wp - frld(:,:)
+               ENDIF
+
+               llnightav = .FALSE.
+#endif
+            END SELECT
+
+            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
+               &               nit000, idaystp, zsurfvar, tmask(:,:,1), &
+               &               nn_2dint, llnightav )
+
          END DO
-      ENDIF      
+
+      ENDIF
+
+      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
+      CALL wrk_dealloc( jpi, jpj, zsurfvar )
+      CALL wrk_dealloc( jpi, jpj, zglam1 )
+      CALL wrk_dealloc( jpi, jpj, zglam2 )
+      CALL wrk_dealloc( jpi, jpj, zgphi1 )
+      CALL wrk_dealloc( jpi, jpj, zgphi2 )
+#if ! defined key_lim2 && ! defined key_lim3
+      CALL wrk_dealloc(jpi,jpj,frld)
 #endif
 
-      !  - Velocity profiles
-      IF ( ln_vel3d ) THEN
-         DO jveloset = 1, nvelosets
-           ! zonal component of velocity
-           CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
-              &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
-                             n1dint, n2dint, ld_velav(jveloset) )
-         END DO
-      ENDIF
-
-#if ! defined key_lim2 && ! defined key_lim3
-      CALL wrk_dealloc(jpi,jpj,frld) 
-#endif
-
    END SUBROUTINE dia_obs
-  
-   SUBROUTINE dia_obs_wri 
+
+   SUBROUTINE dia_obs_wri
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_obs_wri  ***
@@ -1180,5 +754,5 @@
       !! ** Method  : Call observation diagnostic output routines
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!
       !! History :
@@ -1188,258 +762,76 @@
       !!        !  07-03  (K. Mogensen) General handling of profiles
       !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
-      !!----------------------------------------------------------------------
+      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      USE obs_rot_vel          ! Rotation of velocities
+
       IMPLICIT NONE
 
       !! * Local declarations
-
-      INTEGER :: jprofset                 ! Profile data set loop variable
-      INTEGER :: jveloset                 ! Velocity data set loop variable
-      INTEGER :: jslaset                  ! SLA data set loop variable
-      INTEGER :: jsstset                  ! SST data set loop variable
-      INTEGER :: jseaiceset               ! Sea Ice data set loop variable
-      INTEGER :: jset
-      INTEGER :: jfbini
-      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
-      CHARACTER(LEN=10) :: cdtmp
+      INTEGER :: jtype                    ! Data set loop variable
+      INTEGER :: jo, jvar, jk
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
+         & zu, &
+         & zv
+
       !-----------------------------------------------------------------------
       ! Depending on switches call various observation output routines
       !-----------------------------------------------------------------------
 
-      !  - Temperature/salinity profiles
-
-      IF( ln_t3d .OR. ln_s3d ) THEN
-
-         ! Copy data from prodatqc to profdata structures
-         DO jprofset = 1, nprofsets
-
-            CALL obs_prof_decompress( prodatqc(jprofset), &
-                 &                    profdata(jprofset), .TRUE., numout )
+      IF ( nproftypes > 0 ) THEN
+
+         DO jtype = 1, nproftypes
+
+            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
+
+               ! For velocity data, rotate the model velocities to N/S, E/W
+               ! using the compressed data structure.
+               ALLOCATE( &
+                  & zu(profdataqc(jtype)%nvprot(1)), &
+                  & zv(profdataqc(jtype)%nvprot(2))  &
+                  & )
+
+               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
+
+               DO jo = 1, profdataqc(jtype)%nprof
+                  DO jvar = 1, 2
+                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
+
+                        IF ( jvar == 1 ) THEN
+                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
+                        ELSE
+                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
+                        ENDIF
+
+                     END DO
+                  END DO
+               END DO
+
+               DEALLOCATE( zu )
+               DEALLOCATE( zv )
+
+            END IF
+
+            CALL obs_prof_decompress( profdataqc(jtype), &
+               &                      profdata(jtype), .TRUE., numout )
+
+            CALL obs_wri_prof( profdata(jtype) )
 
          END DO
 
-         ! Write the profiles.
-
-         jprofset = 0
-
-         ! ENACT insitu data
-
-         IF ( ln_ena ) THEN
-           
-            jprofset = jprofset + 1
-
-            CALL obs_wri_p3d( 'enact', profdata(jprofset) )
-
-         ENDIF
-
-         ! Coriolis insitu data
-
-         IF ( ln_cor ) THEN
-            
-            jprofset = jprofset + 1
-
-            CALL obs_wri_p3d( 'corio', profdata(jprofset) )
-            
-         ENDIF
-         
-         ! Feedback insitu data
-
-         IF ( ln_profb ) THEN
-
-            jfbini = jprofset + 1
-
-            DO jprofset = jfbini, nprofsets
-               
-               jset = jprofset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'profb_',jset
-               CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
-
-            END DO
-
-         ENDIF
-
-      ENDIF
-
-      !  - Sea surface anomaly
-      IF ( ln_sla ) THEN
-
-         ! Copy data from sladatqc to sladata structures
-         DO jslaset = 1, nslasets
-
-              CALL obs_surf_decompress( sladatqc(jslaset), &
-                 &                    sladata(jslaset), .TRUE., numout )
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         DO jtype = 1, nsurftypes
+
+            CALL obs_surf_decompress( surfdataqc(jtype), &
+               &                      surfdata(jtype), .TRUE., numout )
+
+            CALL obs_wri_surf( surfdata(jtype) )
 
          END DO
 
-         jslaset = 0 
-
-         ! Write the AVISO SLA data
-
-         IF ( ln_sladt ) THEN
-            
-            jslaset = 1
-            CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
-            jslaset = 2
-            CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
-
-         ENDIF
-
-         IF ( ln_slafb ) THEN
-            
-            jfbini = jslaset + 1
-
-            DO jslaset = jfbini, nslasets
-               
-               jset = jslaset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
-               CALL obs_wri_sla( cdtmp, sladata(jslaset) )
-
-            END DO
-
-         ENDIF
-
-      ENDIF
-
-      !  - Sea surface temperature
-      IF ( ln_sst ) THEN
-
-         ! Copy data from sstdatqc to sstdata structures
-         DO jsstset = 1, nsstsets
-     
-              CALL obs_surf_decompress( sstdatqc(jsstset), &
-                 &                    sstdata(jsstset), .TRUE., numout )
-
-         END DO
-
-         jsstset = 0 
-
-         ! Write the AVISO SST data
-
-         IF ( ln_reysst ) THEN
-            
-            jsstset = jsstset + 1
-            CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
-
-         ENDIF
-
-         IF ( ln_ghrsst ) THEN
-            
-            jsstset = jsstset + 1
-            CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
-
-         ENDIF
-
-         IF ( ln_sstfb ) THEN
-            
-            jfbini = jsstset + 1
-
-            DO jsstset = jfbini, nsstsets
-               
-               jset = jsstset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
-               CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
-
-            END DO
-
-         ENDIF
-
-      ENDIF
-
-      !  - Sea surface salinity
-      IF ( ln_sss ) THEN
-         IF(lwp) WRITE(numout,*) ' SSS currently not available'
-      ENDIF
-
-      !  - Sea Ice Concentration
-      IF ( ln_seaice ) THEN
-
-         ! Copy data from seaicedatqc to seaicedata structures
-         DO jseaiceset = 1, nseaicesets
-
-              CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
-                 &                    seaicedata(jseaiceset), .TRUE., numout )
-
-         END DO
-
-         ! Write the Sea Ice data
-         DO jseaiceset = 1, nseaicesets
-      
-            WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
-            CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
-
-         END DO
-
-      ENDIF
-      
-      ! Velocity data
-      IF( ln_vel3d ) THEN
-
-         ! Copy data from veldatqc to velodata structures
-         DO jveloset = 1, nvelosets
-
-            CALL obs_prof_decompress( veldatqc(jveloset), &
-                 &                    velodata(jveloset), .TRUE., numout )
-
-         END DO
-
-         ! Write the profiles.
-
-         jveloset = 0
-
-         ! Daily averaged data
-
-         IF ( ln_velavcur ) THEN
-            
-            jveloset = jveloset + 1
-
-            CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
-
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhrcur ) THEN
-            
-            jveloset = jveloset + 1
-
-            CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
-
-         ENDIF
-
-         ! Daily averaged data
-
-         IF ( ln_velavadcp ) THEN
-            
-            jveloset = jveloset + 1
-
-            CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
-
-         ENDIF
-
-         ! High frequency data
-
-         IF ( ln_velhradcp ) THEN
-            
-            jveloset = jveloset + 1
-            
-            CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
-               
-         ENDIF
-
-         ! Feedback velocity data
-
-         IF ( ln_velfb ) THEN
-
-            jfbini = jveloset + 1
-
-            DO jveloset = jfbini, nvelosets
-               
-               jset = jveloset - jfbini + 1
-               WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
-               CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
-
-            END DO
-
-         ENDIF
-         
       ENDIF
 
@@ -1459,30 +851,14 @@
       !!
       !!----------------------------------------------------------------------
-      !! obs_grid deallocation
+      ! obs_grid deallocation
       CALL obs_grid_deallocate
 
-      !! diaobs deallocation
-      IF ( nprofsets > 0 ) THEN
-          DEALLOCATE(ld_enact, &
-                  &  profdata, &
-                  &  prodatqc)
-      END IF
-      IF ( ln_sla ) THEN
-          DEALLOCATE(sladata, &
-                  &  sladatqc)
-      END IF
-      IF ( ln_seaice ) THEN
-          DEALLOCATE(sladata, &
-                  &  sladatqc)
-      END IF
-      IF ( ln_sst ) THEN
-          DEALLOCATE(sstdata, &
-                  &  sstdatqc)
-      END IF
-      IF ( ln_vel3d ) THEN
-          DEALLOCATE(ld_velav, &
-                  &  velodata, &
-                  &  veldatqc)
-      END IF
+      ! diaobs deallocation
+      IF ( nproftypes > 0 ) &
+         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
+
+      IF ( nsurftypes > 0 ) &
+         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf )
+
    END SUBROUTINE dia_obs_dealloc
 
@@ -1496,4 +872,6 @@
       !!
       !! ** Action  : Get date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
       !!
       !! History :
@@ -1507,6 +885,4 @@
       USE phycst, ONLY : &            ! Physical constants
          & rday
-!      USE daymod, ONLY : &            ! Time variables
-!         & nmonth_len           
       USE dom_oce, ONLY : &           ! Ocean space and time domain variables
          & rdt
@@ -1524,6 +900,8 @@
       INTEGER :: ihou
       INTEGER :: imin
-      INTEGER :: imday         ! Number of days in month.
-      REAL(KIND=wp) :: zdayfrc ! Fraction of day
+      INTEGER :: imday       ! Number of days in month.
+      INTEGER, DIMENSION(12) :: &
+         &       imonth_len  ! Length in days of the months of the current year
+      REAL(wp) :: zdayfrc    ! Fraction of day
 
       INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
@@ -1554,10 +932,10 @@
       iday = iday + kstp * rdt / rday 
 
-      !!-----------------------------------------------------------------------
-      !! Convert number of days (iday) into a real date
-      !!----------------------------------------------------------------------
+      !-----------------------------------------------------------------------
+      ! Convert number of days (iday) into a real date
+      !----------------------------------------------------------------------
 
       CALL calc_month_len( iyea, imonth_len )
-      
+
       DO WHILE ( iday > imonth_len(imon) )
          iday = iday - imonth_len(imon)
@@ -1570,7 +948,7 @@
       END DO
 
-      !!----------------------------------------------------------------------
-      !! Convert it into YYYYMMDD.HHMMSS format.
-      !!----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      ! Convert it into YYYYMMDD.HHMMSS format.
+      !----------------------------------------------------------------------
       ddobs = iyea * 10000_dp + imon * 100_dp + &
          &    iday + ihou * 0.01_dp + imin * 0.0001_dp
@@ -1626,5 +1004,5 @@
 
       !! * Arguments
-      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
 
       CALL calc_date( nitend, ddobsfin )
