Index: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 4990)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 6140)
@@ -25,6 +25,6 @@
    USE netcdf                   ! NetCDF library
    USE obs_oper                 ! Observation operators
-   USE obs_prof_io              ! Profile files I/O (non-FB files)
    USE lib_mpp                  ! For ctl_warn/stop
+   USE obs_fbm                  ! Feedback routines
 
    IMPLICIT NONE
@@ -33,5 +33,5 @@
    PRIVATE
 
-   PUBLIC obs_rea_pro_dri  ! Read the profile observations 
+   PUBLIC obs_rea_prof  ! Read the profile observations 
 
    !!----------------------------------------------------------------------
@@ -42,18 +42,17 @@
 
 CONTAINS
- 
-   SUBROUTINE obs_rea_pro_dri( kformat, &
-      &                        profdata, knumfiles, cfilenames, &
-      &                        kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                        ldt3d, lds3d, ldignmis, ldsatt, ldavtimset, &
-      &                        ldmod, kdailyavtypes )
+
+   SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, &
+      &                     kvars, kextr, kstp, ddobsini, ddobsend, &
+      &                     ldvar1, ldvar2, ldignmis, ldsatt, &
+      &                     ldmod, kdailyavtypes )
       !!---------------------------------------------------------------------
       !!
-      !!                   *** ROUTINE obs_rea_pro_dri ***
+      !!                   *** ROUTINE obs_rea_prof ***
       !!
       !! ** Purpose : Read from file the profile observations
       !!
-      !! ** Method  : Depending on kformat either ENACT, CORIOLIS or
-      !!              feedback data files are read
+      !! ** Method  : Read feedback data in and transform to NEMO internal 
+      !!              profile data structure
       !!
       !! ** Action  : 
@@ -63,31 +62,30 @@
       !! History :  
       !!      ! :  2009-09 (K. Mogensen) : New merged version of old routines
+      !!      ! :  2015-08 (M. Martin) : Merged profile and velocity routines
       !!----------------------------------------------------------------------
-      !! * Modules used
-   
+
       !! * Arguments
-      INTEGER ::  kformat    ! Format of input data
-      !                      ! 1: ENACT
-      !                      ! 2: Coriolis
-      TYPE(obs_prof), INTENT(OUT) ::  profdata     ! Profile data to be read
-      INTEGER, INTENT(IN) :: knumfiles      ! Number of files to read in
+      TYPE(obs_prof), INTENT(OUT) :: &
+         & profdata                     ! Profile data to be read
+      INTEGER, INTENT(IN) :: knumfiles  ! Number of files to read
       CHARACTER(LEN=128), INTENT(IN) ::  &
-         & cfilenames(knumfiles)  ! File names to read in
+         & cdfilenames(knumfiles)        ! File names to read in
       INTEGER, INTENT(IN) :: kvars      ! Number of variables in profdata
-      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var in profdata
-      INTEGER, INTENT(IN) :: kstp        ! Ocean time-step index
-      LOGICAL, INTENT(IN) :: ldt3d       ! Observed variables switches
-      LOGICAL, INTENT(IN) :: lds3d
-      LOGICAL, INTENT(IN) :: ldignmis    ! Ignore missing files
-      LOGICAL, INTENT(IN) :: ldsatt      ! Compute salinity at all temperature points
-      LOGICAL, INTENT(IN) :: ldavtimset  ! Correct time for daily averaged data
-      LOGICAL, INTENT(IN) :: ldmod       ! Initialize model from input data
-      REAL(KIND=dp), INTENT(IN) :: ddobsini    ! Obs. ini time in YYYYMMDD.HHMMSS
-      REAL(KIND=dp), INTENT(IN) :: ddobsend    ! Obs. end time in YYYYMMDD.HHMMSS
+      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var
+      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
+      LOGICAL, INTENT(IN) :: ldvar1     ! Observed variables switches
+      LOGICAL, INTENT(IN) :: ldvar2
+      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
+      LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points
+      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
+      REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
-         & kdailyavtypes
+         & kdailyavtypes                ! Types of daily average observations
 
       !! * Local declarations
-      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_pro_dri'
+      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
+      CHARACTER(len=8) :: clrefdate
+      CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars
       INTEGER :: jvar
       INTEGER :: ji
@@ -105,57 +103,62 @@
       INTEGER :: imin
       INTEGER :: isec
+      INTEGER :: iprof
+      INTEGER :: iproftot
+      INTEGER :: ivar1t0
+      INTEGER :: ivar2t0
+      INTEGER :: ivar1t
+      INTEGER :: ivar2t
+      INTEGER :: ip3dt
+      INTEGER :: ios
+      INTEGER :: ioserrcount
+      INTEGER :: ivar1tmpp
+      INTEGER :: ivar2tmpp
+      INTEGER :: ip3dtmpp
+      INTEGER :: itype
       INTEGER, DIMENSION(knumfiles) :: &
          & irefdate
       INTEGER, DIMENSION(ntyp1770+1) :: &
-         & itypt,    &
-         & ityptmpp, &
-         & ityps,    &
-         & itypsmpp 
-      INTEGER :: it3dtmpp
-      INTEGER :: is3dtmpp
-      INTEGER :: ip3dtmpp
+         & itypvar1,    &
+         & itypvar1mpp, &
+         & itypvar2,    &
+         & itypvar2mpp 
       INTEGER, DIMENSION(:), ALLOCATABLE :: &
-         & iobsi,    &
-         & iobsj,    &
-         & iproc,    &
+         & iobsi1,    &
+         & iobsj1,    &
+         & iproc1,    &
+         & iobsi2,    &
+         & iobsj2,    &
+         & iproc2,    &
          & iindx,    &
          & ifileidx, &
          & iprofidx
-      INTEGER :: itype
       INTEGER, DIMENSION(imaxavtypes) :: &
          & idailyavtypes
+      INTEGER, DIMENSION(kvars) :: &
+         & iv3dt
       REAL(wp), DIMENSION(:), ALLOCATABLE :: &
          & zphi, &
          & zlam
-      real(wp), DIMENSION(:), ALLOCATABLE :: &
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
          & zdat
+      REAL(wp), DIMENSION(knumfiles) :: &
+         & djulini, &
+         & djulend
       LOGICAL :: llvalprof
+      LOGICAL :: lldavtimset
       TYPE(obfbdata), POINTER, DIMENSION(:) :: &
          & inpfiles
-      real(wp), DIMENSION(knumfiles) :: &
-         & djulini, &
-         & djulend
-      INTEGER :: iprof
-      INTEGER :: iproftot
-      INTEGER :: it3dt0
-      INTEGER :: is3dt0
-      INTEGER :: it3dt
-      INTEGER :: is3dt
-      INTEGER :: ip3dt
-      INTEGER :: ios
-      INTEGER :: ioserrcount
-      INTEGER, DIMENSION(kvars) :: &
-         & iv3dt
-      CHARACTER(len=8) :: cl_refdate
-   
+
       ! Local initialization
       iprof = 0
-      it3dt0 = 0
-      is3dt0 = 0
+      ivar1t0 = 0
+      ivar2t0 = 0
       ip3dt = 0
 
       ! Daily average types
+      lldavtimset = .FALSE.
       IF ( PRESENT(kdailyavtypes) ) THEN
          idailyavtypes(:) = kdailyavtypes(:)
+         IF ( ANY (idailyavtypes(:) /= -1) ) lldavtimset = .TRUE.
       ELSE
          idailyavtypes(:) = -1
@@ -163,21 +166,13 @@
 
       !-----------------------------------------------------------------------
-      ! Check data the model part is just with feedback data files
-      !-----------------------------------------------------------------------
-      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
-         CALL ctl_stop( 'Model can only be read from feedback data' )
-         RETURN
-      ENDIF
-
-      !-----------------------------------------------------------------------
       ! Count the number of files needed and allocate the obfbdata type
       !-----------------------------------------------------------------------
-      
+
       inobf = knumfiles
-      
+
       ALLOCATE( inpfiles(inobf) )
 
       prof_files : DO jj = 1, inobf
-          
+
          !---------------------------------------------------------------------
          ! Prints
@@ -186,5 +181,5 @@
             WRITE(numout,*)
             WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', &
-               & TRIM( TRIM( cfilenames(jj) ) )
+               & TRIM( TRIM( cdfilenames(jj) ) )
             WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
             WRITE(numout,*)
@@ -194,25 +189,25 @@
          !  Initialization: Open file and get dimensions only
          !---------------------------------------------------------------------
-         
-         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
+
+         iflag = nf90_open( TRIM( cdfilenames(jj) ), nf90_nowrite, &
             &                      i_file_id )
-         
+
          IF ( iflag /= nf90_noerr ) THEN
 
             IF ( ldignmis ) THEN
                inpfiles(jj)%nobs = 0
-               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
+               CALL ctl_warn( 'File ' // TRIM( cdfilenames(jj) ) // &
                   &           ' not found' )
             ELSE 
-               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
+               CALL ctl_stop( 'File ' // TRIM( cdfilenames(jj) ) // &
                   &           ' not found' )
             ENDIF
 
          ELSE 
-            
+
             !------------------------------------------------------------------
-            !  Close the file since it is opened in read_proffile
+            !  Close the file since it is opened in read_obfbdata
             !------------------------------------------------------------------
-            
+
             iflag = nf90_close( i_file_id )
 
@@ -220,39 +215,31 @@
             !  Read the profile file into inpfiles
             !------------------------------------------------------------------
-            IF ( kformat == 0 ) THEN
-               CALL init_obfbdata( inpfiles(jj) )
-               IF(lwp) THEN
-                  WRITE(numout,*)
-                  WRITE(numout,*)'Reading from feedback file :', &
-                     &           TRIM( cfilenames(jj) )
-               ENDIF
-               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
-                  &                ldgrid = .TRUE. )
-               IF ( inpfiles(jj)%nvar < 2 ) THEN
-                  CALL ctl_stop( 'Feedback format error' )
-                  RETURN
-               ENDIF
-               IF ( TRIM(inpfiles(jj)%cname(1)) /= 'POTM' ) THEN
-                  CALL ctl_stop( 'Feedback format error' )
-                  RETURN
-               ENDIF
-               IF ( TRIM(inpfiles(jj)%cname(2)) /= 'PSAL' ) THEN
-                  CALL ctl_stop( 'Feedback format error' )
-                  RETURN
-               ENDIF
-               IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
-                  CALL ctl_stop( 'Model not in input data' )
-                  RETURN
-               ENDIF
-            ELSEIF ( kformat == 1 ) THEN
-               CALL read_enactfile( TRIM( cfilenames(jj) ), inpfiles(jj), &
-                  &                 numout, lwp, .TRUE. )
-            ELSEIF ( kformat == 2 ) THEN
-               CALL read_coriofile( TRIM( cfilenames(jj) ), inpfiles(jj), &
-                  &                 numout, lwp, .TRUE. )
+            CALL init_obfbdata( inpfiles(jj) )
+            CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
+               &                ldgrid = .TRUE. )
+
+            IF ( inpfiles(jj)%nvar < 2 ) THEN
+               CALL ctl_stop( 'Feedback format error: ', &
+                  &           ' less than 2 vars in profile file' )
+            ENDIF
+
+            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
+               CALL ctl_stop( 'Model not in input data' )
+            ENDIF
+
+            IF ( jj == 1 ) THEN
+               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
+               DO ji = 1, inpfiles(jj)%nvar
+                 clvars(ji) = inpfiles(jj)%cname(ji)
+               END DO
             ELSE
-               CALL ctl_stop( 'File format unknown' )
-            ENDIF
-            
+               DO ji = 1, inpfiles(jj)%nvar
+                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
+                     CALL ctl_stop( 'Feedback file variables not consistent', &
+                        &           ' with previous files for this type' )
+                  ENDIF
+               END DO
+            ENDIF
+
             !------------------------------------------------------------------
             !  Change longitude (-180,180)
@@ -272,7 +259,7 @@
             !  Calculate the date  (change eventually)
             !------------------------------------------------------------------
-            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
-            READ(cl_refdate,'(I8)') irefdate(jj)
-            
+            clrefdate=inpfiles(jj)%cdjuldref(1:8)
+            READ(clrefdate,'(I8)') irefdate(jj)
+
             CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
             CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
@@ -283,26 +270,38 @@
 
             ioserrcount=0
-            IF ( ldavtimset ) THEN
+            IF ( lldavtimset ) THEN
+
+               IF ( ANY ( idailyavtypes(:) /= -1 ) .AND. lwp) THEN
+                  WRITE(numout,*)' Resetting time of daily averaged', &
+                     &           ' observations to the end of the day'
+               ENDIF
+
                DO ji = 1, inpfiles(jj)%nobs
-                  ! 
-                  !  for daily averaged data for example
-                  !  MRB data (itype==820) force the time
-                  !  to be the  end of the day
-                  !
                   READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 900 ) itype
 900               IF ( ios /= 0 ) THEN
-                     itype = 0         ! Set type to zero if there is a problem in the string conversion
-                  ENDIF
-                  IF ( ANY (idailyavtypes == itype ) ) THEN
-                     inpfiles(jj)%ptim(ji) = &
-                     & INT(inpfiles(jj)%ptim(ji)) + 1
-                  ENDIF
+                     ! Set type to zero if there is a problem in the string conversion
+                     itype = 0
+                  ENDIF
+
+                  IF ( ANY ( idailyavtypes(:) == itype ) ) THEN
+                  !  for daily averaged data force the time
+                  !  to be the last time-step of the day, but still within the day.
+                     IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
+                        inpfiles(jj)%ptim(ji) = &
+                           & INT(inpfiles(jj)%ptim(ji)) + 0.9999
+                     ELSE
+                        inpfiles(jj)%ptim(ji) = &
+                           & INT(inpfiles(jj)%ptim(ji)) - 0.0001
+                     ENDIF
+                  ENDIF
+
                END DO
-            ENDIF
-            
+
+            ENDIF
+
             IF ( inpfiles(jj)%nobs > 0 ) THEN
-               inpfiles(jj)%iproc = -1
-               inpfiles(jj)%iobsi = -1
-               inpfiles(jj)%iobsj = -1
+               inpfiles(jj)%iproc(:,:) = -1
+               inpfiles(jj)%iobsi(:,:) = -1
+               inpfiles(jj)%iobsj(:,:) = -1
             ENDIF
             inowin = 0
@@ -318,7 +317,10 @@
             ALLOCATE( zlam(inowin)  )
             ALLOCATE( zphi(inowin)  )
-            ALLOCATE( iobsi(inowin) )
-            ALLOCATE( iobsj(inowin) )
-            ALLOCATE( iproc(inowin) )
+            ALLOCATE( iobsi1(inowin) )
+            ALLOCATE( iobsj1(inowin) )
+            ALLOCATE( iproc1(inowin) )
+            ALLOCATE( iobsi2(inowin) )
+            ALLOCATE( iobsj2(inowin) )
+            ALLOCATE( iproc2(inowin) )
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
@@ -334,5 +336,16 @@
             END DO
 
-            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
+            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, &
+                  &                  iproc1, 'T' )
+               iobsi2(:) = iobsi1(:)
+               iobsj2(:) = iobsj1(:)
+               iproc2(:) = iproc1(:)
+            ELSEIF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, &
+                  &                  iproc1, 'U' )
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, &
+                  &                  iproc2, 'V' )
+            ENDIF
 
             inowin = 0
@@ -344,10 +357,18 @@
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
                   inowin = inowin + 1
-                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
-                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
-                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
+                  inpfiles(jj)%iproc(ji,1) = iproc1(inowin)
+                  inpfiles(jj)%iobsi(ji,1) = iobsi1(inowin)
+                  inpfiles(jj)%iobsj(ji,1) = iobsj1(inowin)
+                  inpfiles(jj)%iproc(ji,2) = iproc2(inowin)
+                  inpfiles(jj)%iobsi(ji,2) = iobsi2(inowin)
+                  inpfiles(jj)%iobsj(ji,2) = iobsj2(inowin)
+                  IF ( inpfiles(jj)%iproc(ji,1) /= &
+                     & inpfiles(jj)%iproc(ji,2) ) THEN
+                     CALL ctl_stop( 'Error in obs_read_prof:', &
+                        & 'var1 and var2 observation on different processors')
+                  ENDIF
                ENDIF
             END DO
-            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
+            DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 )
 
             DO ji = 1, inpfiles(jj)%nobs
@@ -363,5 +384,5 @@
                   ENDIF
                   llvalprof = .FALSE.
-                  IF ( ldt3d ) THEN
+                  IF ( ldvar1 ) THEN
                      loop_t_count : DO ij = 1,inpfiles(jj)%nlev
                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
@@ -369,9 +390,9 @@
                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                           it3dt0 = it3dt0 + 1
+                           ivar1t0 = ivar1t0 + 1
                         ENDIF
                      END DO loop_t_count
                   ENDIF
-                  IF ( lds3d ) THEN
+                  IF ( ldvar2 ) THEN
                      loop_s_count : DO ij = 1,inpfiles(jj)%nlev
                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
@@ -379,5 +400,5 @@
                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                           is3dt0 = is3dt0 + 1
+                           ivar2t0 = ivar2t0 + 1
                         ENDIF
                      END DO loop_s_count
@@ -388,8 +409,8 @@
                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     ldt3d ) .OR. &
+                        &     ldvar1 ) .OR. &
                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     lds3d ) ) THEN
+                        &     ldvar2 ) ) THEN
                         ip3dt = ip3dt + 1
                         llvalprof = .TRUE.
@@ -405,5 +426,5 @@
 
       END DO prof_files
-      
+
       !-----------------------------------------------------------------------
       ! Get the time ordered indices of the input data
@@ -446,5 +467,5 @@
          &               zdat,     &
          &               iindx   )
-      
+
       iv3dt(:) = -1
       IF (ldsatt) THEN
@@ -452,29 +473,29 @@
          iv3dt(2) = ip3dt
       ELSE
-         iv3dt(1) = it3dt0
-         iv3dt(2) = is3dt0
+         iv3dt(1) = ivar1t0
+         iv3dt(2) = ivar2t0
       ENDIF
       CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
          &                 kstp, jpi, jpj, jpk )
-      
+
       ! * Read obs/positions, QC, all variable and assign to profdata
 
       profdata%nprof     = 0
       profdata%nvprot(:) = 0
-
+      profdata%cvars(:)  = clvars(:)
       iprof = 0
 
       ip3dt = 0
-      it3dt = 0
-      is3dt = 0
-      itypt   (:) = 0
-      ityptmpp(:) = 0
-      
-      ityps   (:) = 0
-      itypsmpp(:) = 0
-      
-      ioserrcount = 0      
+      ivar1t = 0
+      ivar2t = 0
+      itypvar1   (:) = 0
+      itypvar1mpp(:) = 0
+
+      itypvar2   (:) = 0
+      itypvar2mpp(:) = 0
+
+      ioserrcount = 0
       DO jk = 1, iproftot
-         
+
          jj = ifileidx(iindx(jk))
          ji = iprofidx(iindx(jk))
@@ -486,5 +507,5 @@
          IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
             & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
-            
+
             IF ( nproc == 0 ) THEN
                IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
@@ -492,5 +513,5 @@
                IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
             ENDIF
-            
+
             llvalprof = .FALSE.
 
@@ -501,30 +522,30 @@
 
             loop_prof : DO ij = 1, inpfiles(jj)%nlev
-               
+
                IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                   & CYCLE
-               
+
                IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                  
+
                   llvalprof = .TRUE. 
                   EXIT loop_prof
-                  
+
                ENDIF
-               
+
                IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                  
+
                   llvalprof = .TRUE. 
                   EXIT loop_prof
-                  
+
                ENDIF
-               
+
             END DO loop_prof
-            
+
             ! Set profile information
-            
+
             IF ( llvalprof ) THEN
-               
+
                iprof = iprof + 1
 
@@ -545,5 +566,5 @@
                profdata%nhou(iprof) = ihou
                profdata%nmin(iprof) = imin
-               
+
                ! Profile space coordinates
                profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
@@ -551,10 +572,12 @@
 
                ! Coordinate search parameters
-               profdata%mi  (iprof,:) = inpfiles(jj)%iobsi(ji,1)
-               profdata%mj  (iprof,:) = inpfiles(jj)%iobsj(ji,1)
-               
+               profdata%mi  (iprof,1) = inpfiles(jj)%iobsi(ji,1)
+               profdata%mj  (iprof,1) = inpfiles(jj)%iobsj(ji,1)
+               profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2)
+               profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2)
+
                ! Profile WMO number
                profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
-               
+
                ! Instrument type
                READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
@@ -564,7 +587,7 @@
                   itype = 0
                ENDIF
-               
+
                profdata%ntyp(iprof) = itype
-               
+
                ! QC stuff
 
@@ -585,6 +608,6 @@
                profdata%nqc(iprof)  = 0 !TODO
 
-               loop_p : DO ij = 1, inpfiles(jj)%nlev            
-                  
+               loop_p : DO ij = 1, inpfiles(jj)%nlev
+
                   IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                      & CYCLE
@@ -594,135 +617,137 @@
                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     ldt3d ) .OR. &
+                        &     ldvar1 ) .OR. &
                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     lds3d ) ) THEN
+                        &     ldvar2 ) ) THEN
                         ip3dt = ip3dt + 1
                      ELSE
                         CYCLE
                      ENDIF
-                     
+
                   ENDIF
 
                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                      &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                     &       ldt3d ) .OR. ldsatt ) THEN
-                     
+                     &       ldvar1 ) .OR. ldsatt ) THEN
+
                      IF (ldsatt) THEN
 
-                        it3dt = ip3dt
+                        ivar1t = ip3dt
 
                      ELSE
 
-                        it3dt = it3dt + 1
-                        
+                        ivar1t = ivar1t + 1
+
                      ENDIF
 
-                     ! Depth of T observation
-                     profdata%var(1)%vdep(it3dt) = &
+                     ! Depth of var1 observation
+                     profdata%var(1)%vdep(ivar1t) = &
                         &                inpfiles(jj)%pdep(ij,ji)
-                     
-                     ! Depth of T observation QC
-                     profdata%var(1)%idqc(it3dt) = &
+
+                     ! Depth of var1 observation QC
+                     profdata%var(1)%idqc(ivar1t) = &
                         &                inpfiles(jj)%idqc(ij,ji)
-                     
-                     ! Depth of T observation QC flags
-                     profdata%var(1)%idqcf(:,it3dt) = &
+
+                     ! Depth of var1 observation QC flags
+                     profdata%var(1)%idqcf(:,ivar1t) = &
                         &                inpfiles(jj)%idqcf(:,ij,ji)
-                     
+
                      ! Profile index
-                     profdata%var(1)%nvpidx(it3dt) = iprof
-                     
+                     profdata%var(1)%nvpidx(ivar1t) = iprof
+
                      ! Vertical index in original profile
-                     profdata%var(1)%nvlidx(it3dt) = ij
-
-                     ! Profile potential T value
+                     profdata%var(1)%nvlidx(ivar1t) = ij
+
+                     ! Profile var1 value
                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                        profdata%var(1)%vobs(it3dt) = &
+                        profdata%var(1)%vobs(ivar1t) = &
                            &                inpfiles(jj)%pob(ij,ji,1)
                         IF ( ldmod ) THEN
-                           profdata%var(1)%vmod(it3dt) = &
+                           profdata%var(1)%vmod(ivar1t) = &
                               &                inpfiles(jj)%padd(ij,ji,1,1)
                         ENDIF
-                        ! Count number of profile T data as function of type
-                        itypt( profdata%ntyp(iprof) + 1 ) = &
-                           & itypt( profdata%ntyp(iprof) + 1 ) + 1
+                        ! Count number of profile var1 data as function of type
+                        itypvar1( profdata%ntyp(iprof) + 1 ) = &
+                           & itypvar1( profdata%ntyp(iprof) + 1 ) + 1
                      ELSE
-                        profdata%var(1)%vobs(it3dt) = fbrmdi
+                        profdata%var(1)%vobs(ivar1t) = fbrmdi
                      ENDIF
 
-                     ! Profile T qc
-                     profdata%var(1)%nvqc(it3dt) = &
+                     ! Profile var1 qc
+                     profdata%var(1)%nvqc(ivar1t) = &
                         & inpfiles(jj)%ivlqc(ij,ji,1)
 
-                     ! Profile T qc flags
-                     profdata%var(1)%nvqcf(:,it3dt) = &
+                     ! Profile var1 qc flags
+                     profdata%var(1)%nvqcf(:,ivar1t) = &
                         & inpfiles(jj)%ivlqcf(:,ij,ji,1)
 
                      ! Profile insitu T value
-                     profdata%var(1)%vext(it3dt,1) = &
-                        &                inpfiles(jj)%pext(ij,ji,1)
-                     
-                  ENDIF
-                  
+                     IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN
+                        profdata%var(1)%vext(ivar1t,1) = &
+                           &                inpfiles(jj)%pext(ij,ji,1)
+                     ENDIF
+
+                  ENDIF
+
                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                      &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                     &   lds3d ) .OR. ldsatt ) THEN
-                     
+                     &   ldvar2 ) .OR. ldsatt ) THEN
+
                      IF (ldsatt) THEN
 
-                        is3dt = ip3dt
+                        ivar2t = ip3dt
 
                      ELSE
 
-                        is3dt = is3dt + 1
-                        
+                        ivar2t = ivar2t + 1
+
                      ENDIF
 
-                     ! Depth of S observation
-                     profdata%var(2)%vdep(is3dt) = &
+                     ! Depth of var2 observation
+                     profdata%var(2)%vdep(ivar2t) = &
                         &                inpfiles(jj)%pdep(ij,ji)
-                     
-                     ! Depth of S observation QC
-                     profdata%var(2)%idqc(is3dt) = &
+
+                     ! Depth of var2 observation QC
+                     profdata%var(2)%idqc(ivar2t) = &
                         &                inpfiles(jj)%idqc(ij,ji)
-                     
-                     ! Depth of S observation QC flags
-                     profdata%var(2)%idqcf(:,is3dt) = &
+
+                     ! Depth of var2 observation QC flags
+                     profdata%var(2)%idqcf(:,ivar2t) = &
                         &                inpfiles(jj)%idqcf(:,ij,ji)
-                     
+
                      ! Profile index
-                     profdata%var(2)%nvpidx(is3dt) = iprof
-                     
+                     profdata%var(2)%nvpidx(ivar2t) = iprof
+
                      ! Vertical index in original profile
-                     profdata%var(2)%nvlidx(is3dt) = ij
-
-                     ! Profile S value
+                     profdata%var(2)%nvlidx(ivar2t) = ij
+
+                     ! Profile var2 value
                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                        profdata%var(2)%vobs(is3dt) = &
+                        profdata%var(2)%vobs(ivar2t) = &
                            &                inpfiles(jj)%pob(ij,ji,2)
                         IF ( ldmod ) THEN
-                           profdata%var(2)%vmod(is3dt) = &
+                           profdata%var(2)%vmod(ivar2t) = &
                               &                inpfiles(jj)%padd(ij,ji,1,2)
                         ENDIF
-                        ! Count number of profile S data as function of type
-                        ityps( profdata%ntyp(iprof) + 1 ) = &
-                           & ityps( profdata%ntyp(iprof) + 1 ) + 1
+                        ! Count number of profile var2 data as function of type
+                        itypvar2( profdata%ntyp(iprof) + 1 ) = &
+                           & itypvar2( profdata%ntyp(iprof) + 1 ) + 1
                      ELSE
-                        profdata%var(2)%vobs(is3dt) = fbrmdi
+                        profdata%var(2)%vobs(ivar2t) = fbrmdi
                      ENDIF
-                     
-                     ! Profile S qc
-                     profdata%var(2)%nvqc(is3dt) = &
+
+                     ! Profile var2 qc
+                     profdata%var(2)%nvqc(ivar2t) = &
                         & inpfiles(jj)%ivlqc(ij,ji,2)
 
-                     ! Profile S qc flags
-                     profdata%var(2)%nvqcf(:,is3dt) = &
+                     ! Profile var2 qc flags
+                     profdata%var(2)%nvqcf(:,ivar2t) = &
                         & inpfiles(jj)%ivlqcf(:,ij,ji,2)
 
                   ENDIF
-            
+
                END DO loop_p
 
@@ -736,12 +761,12 @@
       ! Sum up over processors
       !-----------------------------------------------------------------------
-      
-      CALL obs_mpp_sum_integer ( it3dt0, it3dtmpp )
-      CALL obs_mpp_sum_integer ( is3dt0, is3dtmpp )
-      CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp )
-      
-      CALL obs_mpp_sum_integers( itypt, ityptmpp, ntyp1770 + 1 )
-      CALL obs_mpp_sum_integers( ityps, itypsmpp, ntyp1770 + 1 )
-      
+
+      CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp )
+      CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp )
+      CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  )
+
+      CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 )
+      CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 )
+
       !-----------------------------------------------------------------------
       ! Output number of observations.
@@ -749,14 +774,14 @@
       IF(lwp) THEN
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile data'
+         WRITE(numout,'(A)') ' Profile data'
          WRITE(numout,'(1X,A)') '------------'
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile T data'
-         WRITE(numout,'(1X,A)') '--------------'
+         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) )
+         WRITE(numout,'(1X,A)') '------------------------'
          DO ji = 0, ntyp1770
-            IF ( ityptmpp(ji+1) > 0 ) THEN
+            IF ( itypvar1mpp(ji+1) > 0 ) THEN
                WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
                   & cwmonam1770(ji)(1:52),' = ', &
-                  & ityptmpp(ji+1)
+                  & itypvar1mpp(ji+1)
             ENDIF
          END DO
@@ -764,16 +789,16 @@
             & '---------------------------------------------------------------'
          WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile T data                                 = ',&
-            & it3dtmpp
+            & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// &
+            & '             = ', ivar1tmpp
          WRITE(numout,'(1X,A)') &
             & '---------------------------------------------------------------'
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile S data'
-         WRITE(numout,'(1X,A)') '--------------'
+         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) )
+         WRITE(numout,'(1X,A)') '------------------------'
          DO ji = 0, ntyp1770
-            IF ( itypsmpp(ji+1) > 0 ) THEN
+            IF ( itypvar2mpp(ji+1) > 0 ) THEN
                WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
                   & cwmonam1770(ji)(1:52),' = ', &
-                  & itypsmpp(ji+1)
+                  & itypvar2mpp(ji+1)
             ENDIF
          END DO
@@ -781,11 +806,11 @@
             & '---------------------------------------------------------------'
          WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile S data                                 = ',&
-            & is3dtmpp
+            & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// &
+            & '             = ', ivar2tmpp
          WRITE(numout,'(1X,A)') &
             & '---------------------------------------------------------------'
          WRITE(numout,*) 
       ENDIF
-      
+
       IF (ldsatt) THEN
          profdata%nvprot(1)    = ip3dt
@@ -794,8 +819,8 @@
          profdata%nvprotmpp(2) = ip3dtmpp
       ELSE
-         profdata%nvprot(1)    = it3dt
-         profdata%nvprot(2)    = is3dt
-         profdata%nvprotmpp(1) = it3dtmpp
-         profdata%nvprotmpp(2) = is3dtmpp
+         profdata%nvprot(1)    = ivar1t
+         profdata%nvprot(2)    = ivar2t
+         profdata%nvprotmpp(1) = ivar1tmpp
+         profdata%nvprotmpp(2) = ivar2tmpp
       ENDIF
       profdata%nprof        = iprof
@@ -804,15 +829,15 @@
       ! Model level search
       !-----------------------------------------------------------------------
-      IF ( ldt3d ) THEN
+      IF ( ldvar1 ) THEN
          CALL obs_level_search( jpk, gdept_1d, &
             & profdata%nvprot(1), profdata%var(1)%vdep, &
             & profdata%var(1)%mvk )
       ENDIF
-      IF ( lds3d ) THEN
+      IF ( ldvar2 ) THEN
          CALL obs_level_search( jpk, gdept_1d, &
             & profdata%nvprot(2), profdata%var(2)%vdep, &
             & profdata%var(2)%mvk )
       ENDIF
-      
+
       !-----------------------------------------------------------------------
       ! Set model equivalent to 99999
@@ -826,5 +851,5 @@
       ! Deallocate temporary data
       !-----------------------------------------------------------------------
-      DEALLOCATE( ifileidx, iprofidx, zdat )
+      DEALLOCATE( ifileidx, iprofidx, zdat, clvars )
 
       !-----------------------------------------------------------------------
@@ -836,5 +861,5 @@
       DEALLOCATE( inpfiles )
 
-   END SUBROUTINE obs_rea_pro_dri
+   END SUBROUTINE obs_rea_prof
 
 END MODULE obs_read_prof
