Changeset 14056

Timestamp:

2020-12-03T15:08:29+01:00 (3 years ago)

Author:

ayoung

Message:

Adding branch for ticket #2567 to trunk.

Location:

NEMO/trunk

Files:

• cfgs/SHARED/namelist_ref (modified) (3 diffs)
• doc/namelists/namobs (modified) (3 diffs)
• src/OCE/OBS/diaobs.F90 (modified) (21 diffs)
• src/OCE/OBS/obs_oper.F90 (modified) (19 diffs)
• src/OCE/OBS/obs_prep.F90 (modified) (9 diffs)
• src/OCE/OBS/obs_read_prof.F90 (modified) (27 diffs)
• src/OCE/OBS/obs_read_surf.F90 (modified) (6 diffs)
• src/OCE/OBS/obs_write.F90 (modified) (8 diffs)

NEMO/trunk/cfgs/SHARED/namelist_ref

@@ -1400 +1400 @@
    ln_sstnight = .false.             ! Logical switch for calculating night-time average for SST obs
    ln_bound_reject  = .false.        ! Logical to remove obs near boundaries in LAMs.
+   ln_default_fp_indegs = .true.     ! Logical: T=> averaging footprint is in degrees, F=> in metres
    ln_sla_fp_indegs = .true.         ! Logical for SLA: T=> averaging footprint is in degrees, F=> in metres
    ln_sst_fp_indegs = .true.         ! Logical for SST: T=> averaging footprint is in degrees, F=> in metres
@@ -1415 +1416 @@
    cn_gridsearchfile ='gridsearch.nc' ! Grid search file name
    rn_gridsearchres = 0.5            ! Grid search resolution
+   rn_default_avglamscl = 0.         ! Default E/W diameter of observation footprint (metres/degrees)
+   rn_default_avgphiscl = 0.         ! Default N/S diameter of observation footprint (metres/degrees)
    rn_mdtcorr  = 1.61                ! MDT  correction
    rn_mdtcutoff = 65.0               ! MDT cutoff for computed correction
@@ -1428 +1431 @@
    rn_sic_avgphiscl = 0.             ! N/S diameter of SIC observation footprint (metres/degrees)
    nn_1dint = 0                      ! Type of vertical interpolation method
-   nn_2dint = 0                      ! Default horizontal interpolation method
+   nn_2dint_default = 0              ! Default horizontal interpolation method
    nn_2dint_sla = 0                  ! Horizontal interpolation method for SLA
    nn_2dint_sst = 0                  ! Horizontal interpolation method for SST

NEMO/trunk/doc/namelists/namobs

@@ -20 +20 @@
    ln_sstnight = .false.             ! Logical switch for calculating night-time average for SST obs
    ln_bound_reject  = .false.        ! Logical to remove obs near boundaries in LAMs.
+   ln_default_fp_indegs = .true.     ! Logical: T=> averaging footprint is in degrees, F=> in metres
    ln_sla_fp_indegs = .true.         ! Logical for SLA: T=> averaging footprint is in degrees, F=> in metres
    ln_sst_fp_indegs = .true.         ! Logical for SST: T=> averaging footprint is in degrees, F=> in metres
@@ -39 +40 @@
    rn_dobsini  = 00010101.000000     ! Initial date in window YYYYMMDD.HHMMSS
    rn_dobsend  = 00010102.000000     ! Final date in window YYYYMMDD.HHMMSS
+   rn_default_avglamscl = 0.         ! Default E/W diameter of observation footprint (metres/degrees)
+   rn_default_avgphiscl = 0.         ! Default N/S diameter of observation footprint (metres/degrees)
    rn_sla_avglamscl = 0.             ! E/W diameter of SLA observation footprint (metres/degrees)
    rn_sla_avgphiscl = 0.             ! N/S diameter of SLA observation footprint (metres/degrees)
@@ -48 +51 @@
    rn_sic_avgphiscl = 0.             ! N/S diameter of SIC observation footprint (metres/degrees)
    nn_1dint = 0                      ! Type of vertical interpolation method
-   nn_2dint = 0                      ! Default horizontal interpolation method
+   nn_2dint_default = 0              ! Default horizontal interpolation method
    nn_2dint_sla = 0                  ! Horizontal interpolation method for SLA
    nn_2dint_sst = 0                  ! Horizontal interpolation method for SST

NEMO/trunk/src/OCE/OBS/diaobs.F90

@@ -57 +57 @@
    PUBLIC calc_date        ! Compute the date of a timestep
 
-   LOGICAL, PUBLIC :: ln_diaobs          !: Logical switch for the obs operator
-   LOGICAL         :: ln_sstnight        !  Logical switch for night mean SST obs
-   LOGICAL         :: ln_sla_fp_indegs   !  T=> SLA obs footprint size specified in degrees, F=> in metres
-   LOGICAL         :: ln_sst_fp_indegs   !  T=> SST obs footprint size specified in degrees, F=> in metres
-   LOGICAL         :: ln_sss_fp_indegs   !  T=> SSS obs footprint size specified in degrees, F=> in metres
-   LOGICAL         :: ln_sic_fp_indegs   !  T=> sea-ice obs footprint size specified in degrees, F=> in metres
-
-   REAL(wp) ::   rn_sla_avglamscl   ! E/W diameter of SLA observation footprint (metres)
-   REAL(wp) ::   rn_sla_avgphiscl   ! N/S diameter of SLA observation footprint (metres)
-   REAL(wp) ::   rn_sst_avglamscl   ! E/W diameter of SST observation footprint (metres)
-   REAL(wp) ::   rn_sst_avgphiscl   ! N/S diameter of SST observation footprint (metres)
-   REAL(wp) ::   rn_sss_avglamscl   ! E/W diameter of SSS observation footprint (metres)
-   REAL(wp) ::   rn_sss_avgphiscl   ! N/S diameter of SSS observation footprint (metres)
-   REAL(wp) ::   rn_sic_avglamscl   ! E/W diameter of sea-ice observation footprint (metres)
-   REAL(wp) ::   rn_sic_avgphiscl   ! N/S diameter of sea-ice observation footprint (metres)
-
-   INTEGER :: nn_1dint       ! Vertical interpolation method
-   INTEGER :: nn_2dint       ! Default horizontal interpolation method
-   INTEGER :: nn_2dint_sla   ! SLA horizontal interpolation method 
-   INTEGER :: nn_2dint_sst   ! SST horizontal interpolation method 
-   INTEGER :: nn_2dint_sss   ! SSS horizontal interpolation method 
-   INTEGER :: nn_2dint_sic   ! Seaice horizontal interpolation method 
+   LOGICAL, PUBLIC :: ln_diaobs            !: Logical switch for the obs operator
+   LOGICAL         :: ln_sstnight          !  Logical switch for night mean SST obs
+   LOGICAL         :: ln_default_fp_indegs !  T=> Default obs footprint size specified in degrees, F=> in metres
+   LOGICAL         :: ln_sla_fp_indegs     !  T=> SLA obs footprint size specified in degrees, F=> in metres
+   LOGICAL         :: ln_sst_fp_indegs     !  T=> SST obs footprint size specified in degrees, F=> in metres
+   LOGICAL         :: ln_sss_fp_indegs     !  T=> SSS obs footprint size specified in degrees, F=> in metres
+   LOGICAL         :: ln_sic_fp_indegs     !  T=> sea-ice obs footprint size specified in degrees, F=> in metres
+
+   REAL(wp) ::   rn_default_avglamscl      ! E/W diameter of SLA observation footprint (metres)
+   REAL(wp) ::   rn_default_avgphiscl      ! N/S diameter of SLA observation footprint (metre
+   REAL(wp) ::   rn_sla_avglamscl          ! E/W diameter of SLA observation footprint (metres)
+   REAL(wp) ::   rn_sla_avgphiscl          ! N/S diameter of SLA observation footprint (metres)
+   REAL(wp) ::   rn_sst_avglamscl          ! E/W diameter of SST observation footprint (metres)
+   REAL(wp) ::   rn_sst_avgphiscl          ! N/S diameter of SST observation footprint (metres)
+   REAL(wp) ::   rn_sss_avglamscl          ! E/W diameter of SSS observation footprint (metres)
+   REAL(wp) ::   rn_sss_avgphiscl          ! N/S diameter of SSS observation footprint (metres)
+   REAL(wp) ::   rn_sic_avglamscl          ! E/W diameter of sea-ice observation footprint (metres)
+   REAL(wp) ::   rn_sic_avgphiscl          ! N/S diameter of sea-ice observation footprint (metres)
+
+   INTEGER :: nn_1dint                     ! Vertical interpolation method
+   INTEGER :: nn_2dint_default             ! Default horizontal interpolation method
+   INTEGER :: nn_2dint_sla                 ! SLA horizontal interpolation method 
+   INTEGER :: nn_2dint_sst                 ! SST horizontal interpolation method 
+   INTEGER :: nn_2dint_sss                 ! SSS horizontal interpolation method 
+   INTEGER :: nn_2dint_sic                 ! Seaice horizontal interpolation method 
    INTEGER, DIMENSION(imaxavtypes) ::   nn_profdavtypes   ! Profile data types representing a daily average
    INTEGER :: nproftypes     ! Number of profile obs types
@@ -94 +97 @@
    TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) ::   profdataqc   !: Profile data after quality control
 
-   CHARACTER(len=lca), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types
+   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types
 
    !!----------------------------------------------------------------------
@@ -121 +124 @@
       INTEGER :: jvar            ! Counter for variables
       INTEGER :: jfile           ! Counter for files
-      INTEGER :: jnumsstbias
+      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply
+      INTEGER :: n2dint_type     ! Local version of nn_2dint*
       !
       CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
@@ -130 +134 @@
          & cn_sicfbfiles, &      ! Seaice concentration input filenames
          & cn_velfbfiles, &      ! Velocity profile input filenames
-         & cn_sstbiasfiles      ! SST bias input filenames
+         & cn_sstbiasfiles       ! SST bias input filenames
       CHARACTER(LEN=128) :: &
          & cn_altbiasfile        ! Altimeter bias input filename
@@ -136 +140 @@
          & clproffiles, &        ! Profile filenames
          & clsurffiles           ! Surface filenames
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: &
+         & clvars                ! Expected variable names
          !
       LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
@@ -150 +156 @@
       LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
       LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs.
-      LOGICAL :: llvar1          ! Logical for profile variable 1
-      LOGICAL :: llvar2          ! Logical for profile variable 1
+      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
+      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables)
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read
       LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files
       !
-      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
-      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zglam1, zglam2   ! Model longitudes for profile variable 1 & 2
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zgphi1, zgphi2   ! Model latitudes  for profile variable 1 & 2
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask1, zmask2   ! Model land/sea mask associated with variable 1 & 2
+      REAL(dp) :: rn_dobsini      ! Obs window start date YYYYMMDD.HHMMSS
+      REAL(dp) :: rn_dobsend      ! Obs window end date   YYYYMMDD.HHMMSS
+      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl
+      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl
+      REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: zglam   ! Model longitudes for profile variables
+      REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: zgphi   ! Model latitudes  for profile variables
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zmask   ! Model land/sea mask associated with variables
       !!
       NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
@@ -165 +174 @@
          &            ln_grid_global, ln_grid_search_lookup,          &
          &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
-         &            ln_sstnight,                                    &
+         &            ln_sstnight, ln_default_fp_indegs,              &
          &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
          &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
@@ -174 +183 @@
          &            cn_gridsearchfile, rn_gridsearchres,            &
          &            rn_dobsini, rn_dobsend,                         &
+         &            rn_default_avglamscl, rn_default_avgphiscl,     &
          &            rn_sla_avglamscl, rn_sla_avgphiscl,             &
          &            rn_sst_avglamscl, rn_sst_avgphiscl,             &
          &            rn_sss_avglamscl, rn_sss_avgphiscl,             &
          &            rn_sic_avglamscl, rn_sic_avgphiscl,             &
-         &            nn_1dint, nn_2dint,                             &
+         &            nn_1dint, nn_2dint_default,                     &
          &            nn_2dint_sla, nn_2dint_sst,                     &
          &            nn_2dint_sss, nn_2dint_sic,                     &
@@ -234 +244 @@
          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,*) '      Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default
+         WRITE(numout,*) '      Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla
+         WRITE(numout,*) '      Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst
+         WRITE(numout,*) '      Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss        
+         WRITE(numout,*) '      Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic
+         WRITE(numout,*) '      Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl
+         WRITE(numout,*) '      Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl
+         WRITE(numout,*) '      Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs
+         WRITE(numout,*) '      SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl
+         WRITE(numout,*) '      SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl
+         WRITE(numout,*) '      SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs
+         WRITE(numout,*) '      SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl
+         WRITE(numout,*) '      SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl
+         WRITE(numout,*) '      SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs
+         WRITE(numout,*) '      SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl
+         WRITE(numout,*) '      SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl
+         WRITE(numout,*) '      SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs
          WRITE(numout,*) '      Rejection of observations near land switch               ln_nea = ', ln_nea
          WRITE(numout,*) '      Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject
@@ -278 +304 @@
          IF( ln_t3d .OR. ln_s3d ) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', &
-               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
+            cobstypesprof(jtype) = 'prof'
+            clproffiles(jtype,:) = cn_profbfiles
          ENDIF
          IF( ln_vel3d ) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', &
-               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
+            cobstypesprof(jtype) = 'vel'
+            clproffiles(jtype,:) = cn_velfbfiles
          ENDIF
+         !
+         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles )
          !
       ENDIF
@@ -303 +331 @@
          IF( ln_sla ) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', &
-               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      &
-               &                  nn_2dint, nn_2dint_sla,             &
-               &                  rn_sla_avglamscl, rn_sla_avgphiscl, &
-               &                  ln_sla_fp_indegs, .FALSE.,          &
-               &                  n2dintsurf, zavglamscl, zavgphiscl, &
-               &                  lfpindegs, llnightav )
+            cobstypessurf(jtype) = 'sla'
+            clsurffiles(jtype,:) = cn_slafbfiles
          ENDIF
          IF( ln_sst ) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', &
-               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      &
-               &                  nn_2dint, nn_2dint_sst,             &
-               &                  rn_sst_avglamscl, rn_sst_avgphiscl, &
-               &                  ln_sst_fp_indegs, ln_sstnight,      &
-               &                  n2dintsurf, zavglamscl, zavgphiscl, &
-               &                  lfpindegs, llnightav )
+            cobstypessurf(jtype) = 'sst'
+            clsurffiles(jtype,:) = cn_sstfbfiles
          ENDIF
 #if defined key_si3 || defined key_cice
          IF( ln_sic ) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', &
-               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      &
-               &                  nn_2dint, nn_2dint_sic,             &
-               &                  rn_sic_avglamscl, rn_sic_avgphiscl, &
-               &                  ln_sic_fp_indegs, .FALSE.,          &
-               &                  n2dintsurf, zavglamscl, zavgphiscl, &
-               &                  lfpindegs, llnightav )
+            cobstypessurf(jtype) = 'sic'
+            clsurffiles(jtype,:) = cn_sicfbfiles
          ENDIF
 #endif
          IF( ln_sss ) THEN
             jtype = jtype + 1
-            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', &
-               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
-            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      &
-               &                  nn_2dint, nn_2dint_sss,             &
-               &                  rn_sss_avglamscl, rn_sss_avgphiscl, &
-               &                  ln_sss_fp_indegs, .FALSE.,          &
-               &                  n2dintsurf, zavglamscl, zavgphiscl, &
-               &                  lfpindegs, llnightav )
+            cobstypessurf(jtype) = 'sss'
+            clsurffiles(jtype,:) = cn_sssfbfiles
          ENDIF
+         !
+         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles )
+
+         DO jtype = 1, nsurftypes
+
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
+               IF ( nn_2dint_sla == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
+               ELSE
+                  n2dint_type  = nn_2dint_sla
+               ENDIF
+               ztype_avglamscl = rn_sla_avglamscl
+               ztype_avgphiscl = rn_sla_avgphiscl
+               ltype_fp_indegs = ln_sla_fp_indegs
+               ltype_night     = .FALSE.
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
+               IF ( nn_2dint_sst == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
+               ELSE
+                  n2dint_type  = nn_2dint_sst
+               ENDIF
+               ztype_avglamscl = rn_sst_avglamscl
+               ztype_avgphiscl = rn_sst_avgphiscl
+               ltype_fp_indegs = ln_sst_fp_indegs
+               ltype_night     = ln_sstnight
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
+               IF ( nn_2dint_sic == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
+               ELSE
+                  n2dint_type  = nn_2dint_sic
+               ENDIF
+               ztype_avglamscl = rn_sic_avglamscl
+               ztype_avgphiscl = rn_sic_avgphiscl
+               ltype_fp_indegs = ln_sic_fp_indegs
+               ltype_night     = .FALSE.
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
+               IF ( nn_2dint_sss == -1 ) THEN
+                  n2dint_type  = nn_2dint_default
+               ELSE
+                  n2dint_type  = nn_2dint_sss
+               ENDIF
+               ztype_avglamscl = rn_sss_avglamscl
+               ztype_avgphiscl = rn_sss_avgphiscl
+               ltype_fp_indegs = ln_sss_fp_indegs
+               ltype_night     = .FALSE.
+            ELSE
+               n2dint_type     = nn_2dint_default
+               ztype_avglamscl = rn_default_avglamscl
+               ztype_avgphiscl = rn_default_avgphiscl
+               ltype_fp_indegs = ln_default_fp_indegs
+               ltype_night     = .FALSE.
+            ENDIF
+            
+            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), &
+               &                    nn_2dint_default, n2dint_type,                 &
+               &                    ztype_avglamscl, ztype_avgphiscl,              &
+               &                    ltype_fp_indegs, ltype_night,                  &
+               &                    n2dintsurf, zavglamscl, zavgphiscl,            &
+               &                    lfpindegs, llnightav )
+
+         END DO
          !
       ENDIF
@@ -368 +433 @@
       ENDIF
       !
-      IF( nn_2dint < 0  .OR.  nn_2dint > 6  ) THEN
-         CALL ctl_stop('dia_obs_init: Choice of horizontal (2D) interpolation method is not available')
+      IF( nn_2dint_default < 0  .OR.  nn_2dint_default > 6  ) THEN
+         CALL ctl_stop('dia_obs_init: Choice of default horizontal (2D) interpolation method is not available')
       ENDIF
       !
@@ -388 +453 @@
          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
+               nvarsprof(jtype) = 2
+               nextrprof(jtype) = 1             
+               ALLOCATE( llvar (nvarsprof(jtype)) )
+               ALLOCATE( clvars(nvarsprof(jtype)) )
+               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) )
+               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) )
+               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) )
+               llvar(1)       = ln_t3d
+               llvar(2)       = ln_s3d
+               clvars(1)      = 'POTM'
+               clvars(2)      = 'PSAL'
+               zglam(:,:,1)   = glamt(:,:)
+               zglam(:,:,2)   = glamt(:,:)
+               zgphi(:,:,1)   = gphit(:,:)
+               zgphi(:,:,2)   = gphit(:,:)
+               zmask(:,:,:,1) = tmask(:,:,:)
+               zmask(:,:,:,2) = tmask(:,:,:)
+            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
+               nvarsprof(jtype) = 2
                nextrprof(jtype) = 2
-               llvar1 = ln_vel3d
-               llvar2 = ln_vel3d
-               zglam1 = glamu
-               zgphi1 = gphiu
-               zmask1 = umask
-               zglam2 = glamv
-               zgphi2 = gphiv
-               zmask2 = vmask
+               ALLOCATE( llvar (nvarsprof(jtype)) )
+               ALLOCATE( clvars(nvarsprof(jtype)) )
+               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) )
+               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) )
+               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) )
+               llvar(1)       = ln_vel3d
+               llvar(2)       = ln_vel3d
+               clvars(1)      = 'UVEL'
+               clvars(2)      = 'VVEL'
+               zglam(:,:,1)   = glamu(:,:)
+               zglam(:,:,2)   = glamv(:,:)
+               zgphi(:,:,1)   = gphiu(:,:)
+               zgphi(:,:,2)   = gphiv(:,:)
+               zmask(:,:,:,1) = umask(:,:,:)
+               zmask(:,:,:,2) = vmask(:,:,:)
+            ELSE
+               nvarsprof(jtype) = 1
+               nextrprof(jtype) = 0
+               ALLOCATE( llvar (nvarsprof(jtype)) )
+               ALLOCATE( clvars(nvarsprof(jtype)) )
+               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) )
+               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) )
+               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) )
+               llvar(1)       = .TRUE.
+               zglam(:,:,1)   = glamt(:,:)
+               zgphi(:,:,1)   = gphit(:,:)
+               zmask(:,:,:,1) = tmask(:,:,:)
             ENDIF
             !
@@ -416 +507 @@
                &               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., &
+               &               rn_dobsini, rn_dobsend, llvar, &
+               &               ln_ignmis, ln_s_at_t, .FALSE., clvars, &
                &               kdailyavtypes = nn_profdavtypes )
                !
@@ -425 +516 @@
             !
             CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
-               &               llvar1, llvar2, &
+               &               llvar, &
                &               jpi, jpj, jpk, &
-               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
+               &               zmask, zglam, zgphi, &
                &               ln_nea, ln_bound_reject, Kmm, &
                &               kdailyavtypes = nn_profdavtypes )
+            !
+            DEALLOCATE( llvar, clvars, zglam, zgphi, zmask )
+            !
          END DO
          !
@@ -449 +543 @@
             IF( TRIM(cobstypessurf(jtype)) == 'sst' )   llnightav(jtype) = ln_sstnight
             !
+            ALLOCATE( clvars( nvarssurf(jtype) ) )
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
+               clvars(1) = 'SLA'
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
+               clvars(1) = 'SST'
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
+               clvars(1) = 'ICECONC'
+            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
+               clvars(1) = 'SSS'
+            ENDIF
+            !
             ! Read in surface obs types
             CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
                &               clsurffiles(jtype,1:ifilessurf(jtype)), &
                &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
-               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
+               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype), &
+               &               clvars )
                !
             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
@@ -473 +579 @@
                   &                  jnumsstbias      , cn_sstbiasfiles(1:jnumsstbias) ) 
             ENDIF
+            !
+            DEALLOCATE( clvars )
          END DO
          !
@@ -516 +624 @@
       INTEGER :: jvar              ! Variable number
       INTEGER :: ji, jj            ! Loop counters
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
-         & zprofvar1, &            ! Model values for 1st variable in a prof ob
-         & zprofvar2               ! Model values for 2nd variable in a prof ob
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
-         & zprofmask1, &           ! Mask associated with zprofvar1
-         & zprofmask2              ! Mask associated with zprofvar2
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
+         & zprofvar                ! Model values for variables in a prof ob
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
+         & zprofmask               ! Mask associated with zprofvar
       REAL(wp), DIMENSION(jpi,jpj) :: &
          & zsurfvar, &             ! Model values equivalent to surface ob.
          & zsurfmask               ! Mask associated with surface variable
-      REAL(wp), DIMENSION(jpi,jpj) :: &
-         & 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
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
+         & zglam,    &             ! Model longitudes for prof variables
+         & zgphi                   ! Model latitudes for prof variables
 
       !-----------------------------------------------------------------------
@@ -549 +653 @@
          DO jtype = 1, nproftypes
 
+            ! Allocate local work arrays
+            ALLOCATE( zprofvar (jpi, jpj, jpk, profdataqc(jtype)%nvar) )
+            ALLOCATE( zprofmask(jpi, jpj, jpk, profdataqc(jtype)%nvar) )
+            ALLOCATE( zglam    (jpi, jpj,      profdataqc(jtype)%nvar) )
+            ALLOCATE( zgphi    (jpi, jpj,      profdataqc(jtype)%nvar) )  
+                              
+            ! Defaults which might change
+            DO jvar = 1, profdataqc(jtype)%nvar
+               zprofmask(:,:,:,jvar) = tmask(:,:,:)
+               zglam(:,:,jvar)       = glamt(:,:)
+               zgphi(:,:,jvar)       = gphit(:,:)
+            END DO
+
             SELECT CASE ( TRIM(cobstypesprof(jtype)) )
             CASE('prof')
-               zprofvar1(:,:,:) = ts(:,:,:,jp_tem,Kmm)
-               zprofvar2(:,:,:) = ts(:,:,:,jp_sal,Kmm)
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zprofmask2(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zglam2(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
-               zgphi2(:,:) = gphit(:,:)
+               zprofvar(:,:,:,1) = ts(:,:,:,jp_tem,Kmm)
+               zprofvar(:,:,:,2) = ts(:,:,:,jp_sal,Kmm)
             CASE('vel')
-               zprofvar1(:,:,:) = uu(:,:,:,Kmm)
-               zprofvar2(:,:,:) = vv(:,:,:,Kmm)
-               zprofmask1(:,:,:) = umask(:,:,:)
-               zprofmask2(:,:,:) = vmask(:,:,:)
-               zglam1(:,:) = glamu(:,:)
-               zglam2(:,:) = glamv(:,:)
-               zgphi1(:,:) = gphiu(:,:)
-               zgphi2(:,:) = gphiv(:,:)
+               zprofvar(:,:,:,1) = uu(:,:,:,Kmm)
+               zprofvar(:,:,:,2) = vv(:,:,:,Kmm)
+               zprofmask(:,:,:,1) = umask(:,:,:)
+               zprofmask(:,:,:,2) = vmask(:,:,:)
+               zglam(:,:,1) = glamu(:,:)
+               zglam(:,:,2) = glamv(:,:)
+               zgphi(:,:,1) = gphiu(:,:)
+               zgphi(:,:,2) = gphiv(:,:)
             CASE DEFAULT
                CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
             END SELECT
 
-            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
-               &               nit000, idaystp,                         &
-               &               zprofvar1, zprofvar2,                    &
-               &               gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm),      & 
-               &               zprofmask1, zprofmask2,                  &
-               &               zglam1, zglam2, zgphi1, zgphi2,          &
-               &               nn_1dint, nn_2dint,                      &
-               &               kdailyavtypes = nn_profdavtypes )
+            DO jvar = 1, profdataqc(jtype)%nvar
+               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
+                  &               nit000, idaystp, jvar,                   &
+                  &               zprofvar(:,:,:,jvar),                    &
+                  &               gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm),      & 
+                  &               zprofmask(:,:,:,jvar),                   &
+                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &
+                  &               nn_1dint, nn_2dint_default,              &
+                  &               kdailyavtypes = nn_profdavtypes )
+            END DO
+            
+            DEALLOCATE( zprofvar, zprofmask, zglam, zgphi )
 
          END DO
@@ -680 +795 @@
                   & )
 
-               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
+               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
 
                DO jo = 1, profdataqc(jtype)%nprof
@@ -896 +1011 @@
    END SUBROUTINE fin_date
    
-    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
-       &                         cfilestype, ifiles, cobstypes, cfiles )
-
-    INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
-    INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
-    INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs
-    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                   ifiles      ! Out appended number of files for this type
-
-    CHARACTER(len=6), INTENT(IN) :: ctypein 
-    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
-       &                   cfilestype  ! In list of files for this obs type
-    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                   cobstypes   ! Out appended list of obs types
-    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
-       &                   cfiles      ! Out appended list of files for all types
-
-    !Local variables
-    INTEGER :: jfile
-
-    cfiles(jtype,:) = cfilestype(:)
-    cobstypes(jtype) = ctypein
-    ifiles(jtype) = 0
-    DO jfile = 1, jpmaxnfiles
-       IF ( trim(cfiles(jtype,jfile)) /= '' ) &
-                 ifiles(jtype) = ifiles(jtype) + 1
-    END DO
-
-    IF ( ifiles(jtype) == 0 ) THEN
-         CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   &
-            &           ' set to true but no files available to read' )
-    ENDIF
-
-    IF(lwp) THEN    
-       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
-       DO jfile = 1, ifiles(jtype)
-          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
-       END DO
-    ENDIF
-
-    END SUBROUTINE obs_settypefiles
-
-    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
-               &                  n2dint_default, n2dint_type,        &
-               &                  zavglamscl_type, zavgphiscl_type,   &
-               &                  lfp_indegs_type, lavnight_type,     &
-               &                  n2dint, zavglamscl, zavgphiscl,     &
-               &                  lfpindegs, lavnight )
-
-    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
-    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
-    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
-    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
-    REAL(wp), INTENT(IN) :: &
-       &                    zavglamscl_type, & !E/W diameter of obs footprint for this type
-       &                    zavgphiscl_type    !N/S diameter of obs footprint for this type
-    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
-    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
-    CHARACTER(len=6), INTENT(IN) :: ctypein 
-
-    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                    n2dint 
-    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                    zavglamscl, zavgphiscl
-    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
-       &                    lfpindegs, lavnight
-
-    lavnight(jtype) = lavnight_type
-
-    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
-       n2dint(jtype) = n2dint_type
-    ELSE
-       n2dint(jtype) = n2dint_default
-    ENDIF
-
-    ! For averaging observation footprints set options for size of footprint 
-    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
-       IF ( zavglamscl_type > 0._wp ) THEN
-          zavglamscl(jtype) = zavglamscl_type
-       ELSE
-          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
-                         'scale (zavglamscl) for observation type '//TRIM(ctypein) )      
-       ENDIF
-
-       IF ( zavgphiscl_type > 0._wp ) THEN
-          zavgphiscl(jtype) = zavgphiscl_type
-       ELSE
-          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
-                         'scale (zavgphiscl) for observation type '//TRIM(ctypein) )      
-       ENDIF
-
-       lfpindegs(jtype) = lfp_indegs_type 
-
-    ENDIF
-
-    ! Write out info 
-    IF(lwp) THEN
-       IF ( n2dint(jtype) <= 4 ) THEN
-          WRITE(numout,*) '             '//TRIM(ctypein)// &
-             &            ' model counterparts will be interpolated horizontally'
-       ELSE IF ( n2dint(jtype) <= 6 ) THEN
-          WRITE(numout,*) '             '//TRIM(ctypein)// &
-             &            ' model counterparts will be averaged horizontally'
-          WRITE(numout,*) '             '//'    with E/W scale: ',zavglamscl(jtype)
-          WRITE(numout,*) '             '//'    with N/S scale: ',zavgphiscl(jtype)
-          IF ( lfpindegs(jtype) ) THEN
-              WRITE(numout,*) '             '//'    (in degrees)'
-          ELSE
-              WRITE(numout,*) '             '//'    (in metres)'
-          ENDIF
-       ENDIF
-    ENDIF
-
-    END SUBROUTINE obs_setinterpopts
+   SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles )
+
+      INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
+      INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
+      INTEGER, DIMENSION(ntypes), INTENT(OUT) :: &
+         &                   ifiles      ! Out number of files for each type
+      CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: &
+         &                   cobstypes   ! List of obs types
+      CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: &
+         &                   cfiles      ! List of files for all types
+
+      !Local variables
+      INTEGER :: jfile
+      INTEGER :: jtype
+
+      DO jtype = 1, ntypes
+
+         ifiles(jtype) = 0
+         DO jfile = 1, jpmaxnfiles
+            IF ( trim(cfiles(jtype,jfile)) /= '' ) &
+                      ifiles(jtype) = ifiles(jtype) + 1
+         END DO
+
+         IF ( ifiles(jtype) == 0 ) THEN
+              CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//   &
+                 &           ' set to true but no files available to read' )
+         ENDIF
+
+         IF(lwp) THEN    
+            WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
+            DO jfile = 1, ifiles(jtype)
+               WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
+            END DO
+         ENDIF
+
+      END DO
+
+   END SUBROUTINE obs_settypefiles
+
+   SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
+              &                  n2dint_default, n2dint_type,        &
+              &                  ravglamscl_type, ravgphiscl_type,   &
+              &                  lfp_indegs_type, lavnight_type,     &
+              &                  n2dint, ravglamscl, ravgphiscl,     &
+              &                  lfpindegs, lavnight )
+
+      INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
+      INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
+      INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
+      INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
+      REAL(wp), INTENT(IN) :: &
+         &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
+         &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
+      LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
+      LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
+      CHARACTER(len=8), INTENT(IN) :: ctypein 
+
+      INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
+         &                    n2dint 
+      REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
+         &                    ravglamscl, ravgphiscl
+      LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
+         &                    lfpindegs, lavnight
+
+      lavnight(jtype) = lavnight_type
+
+      IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN
+         n2dint(jtype) = n2dint_type
+      ELSE IF ( n2dint_type == -1 ) THEN
+         n2dint(jtype) = n2dint_default
+      ELSE
+         CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
+           &                    ' is not available')
+      ENDIF
+
+      ! For averaging observation footprints set options for size of footprint 
+      IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
+         IF ( ravglamscl_type > 0._wp ) THEN
+            ravglamscl(jtype) = ravglamscl_type
+         ELSE
+            CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
+                           'scale (ravglamscl) for observation type '//TRIM(ctypein) )      
+         ENDIF
+
+         IF ( ravgphiscl_type > 0._wp ) THEN
+            ravgphiscl(jtype) = ravgphiscl_type
+         ELSE
+            CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
+                           'scale (ravgphiscl) for observation type '//TRIM(ctypein) )      
+         ENDIF
+
+         lfpindegs(jtype) = lfp_indegs_type 
+
+      ENDIF
+
+      ! Write out info 
+      IF(lwp) THEN
+         IF ( n2dint(jtype) <= 4 ) THEN
+            WRITE(numout,*) '             '//TRIM(ctypein)// &
+               &            ' model counterparts will be interpolated horizontally'
+         ELSE IF ( n2dint(jtype) <= 6 ) THEN
+            WRITE(numout,*) '             '//TRIM(ctypein)// &
+               &            ' model counterparts will be averaged horizontally'
+            WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype)
+            WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype)
+            IF ( lfpindegs(jtype) ) THEN
+                WRITE(numout,*) '             '//'    (in degrees)'
+            ELSE
+                WRITE(numout,*) '             '//'    (in metres)'
+            ENDIF
+         ENDIF
+      ENDIF
+
+   END SUBROUTINE obs_setinterpopts
 
 END MODULE diaobs

NEMO/trunk/src/OCE/OBS/obs_oper.F90

@@ -40 +40 @@
 CONTAINS
 
-   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          &
-      &                     kit000, kdaystp,                      &
-      &                     pvar1, pvar2, pgdept, pgdepw,         &
-      &                     pmask1, pmask2,                       &  
-      &                     plam1, plam2, pphi1, pphi2,           &
+   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
+      &                     kit000, kdaystp, kvar,       &
+      &                     pvar, pgdept, pgdepw,        &
+      &                     pmask,                       &  
+      &                     plam, pphi,                  &
       &                     k1dint, k2dint, kdailyavtypes )
       !!-----------------------------------------------------------------------
@@ -105 +105 @@
       INTEGER       , INTENT(in   ) ::   k2dint          ! Horizontal interpolation type (see header)
       INTEGER       , INTENT(in   ) ::   kdaystp         ! Number of time steps per day
-      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pvar1 , pvar2    ! Model field     1 and 2
-      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pmask1, pmask2   ! Land-sea mask   1 and 2
-      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   plam1 , plam2    ! Model longitude 1 and 2
-      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   pphi1 , pphi2    ! Model latitudes 1 and 2
+      INTEGER       , INTENT(in   ) ::   kvar            ! Number of variables in prodatqc
+      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pvar             ! Model field
+      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pmask            ! Land-sea mask
+      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   plam             ! Model longitude
+      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   pphi             ! Model latitudes
       REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pgdept, pgdepw   ! depth of T and W levels 
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL ::   kdailyavtypes             ! Types for daily averages
@@ -128 +129 @@
          & idailyavtypes
       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi1, &
-         & igrdi2, &
-         & igrdj1, &
-         & igrdj2
+         & igrdi, &
+         & igrdj
       INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
 
@@ -138 +137 @@
       REAL(KIND=wp) :: zdaystp
       REAL(KIND=wp), DIMENSION(kpk) :: &
-         & zobsmask1, &
-         & zobsmask2, &
-         & zobsk,    &
+         & zobsk,  &
          & zobs2k
       REAL(KIND=wp), DIMENSION(2,2,1) :: &
          & zweig1, &
-         & zweig2, &
          & zweig
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
-         & zmask1, &
-         & zmask2, &
-         & zint1,  &
-         & zint2,  &
-         & zinm1,  &
-         & zinm2,  &
+         & zmask,  &
+         & zint,   &
+         & zinm,   &
          & zgdept, & 
          & zgdepw
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zglam1, &
-         & zglam2, &
-         & zgphi1, &
-         & zgphi2
-      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2   
+         & zglam,  &
+         & zgphi
+      REAL(KIND=wp), DIMENSION(1) :: zmsk
       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
 
@@ -190 +181 @@
          IF ( idayend == 1 .OR. kt == 0 ) THEN
             DO_3D( 1, 1, 1, 1, 1, jpk )
-               prodatqc%vdmean(ji,jj,jk,1) = 0.0
-               prodatqc%vdmean(ji,jj,jk,2) = 0.0
+               prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
             END_3D
          ENDIF
@@ -197 +187 @@
          DO_3D( 1, 1, 1, 1, 1, jpk )
             ! Increment field 1 for computing daily mean
-            prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-               &                        + pvar1(ji,jj,jk)
-            ! Increment field 2 for computing daily mean
-            prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-               &                        + pvar2(ji,jj,jk)
+            prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
+               &                           + pvar(ji,jj,jk)
          END_3D
 
@@ -210 +197 @@
             CALL FLUSH(numout)
             DO_3D( 1, 1, 1, 1, 1, jpk )
-               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                  &                        * zdaystp
-               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                  &                        * zdaystp
+               prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
+                  &                           * zdaystp
             END_3D
          ENDIF
@@ -221 +206 @@
       ! Get the data for interpolation
       ALLOCATE( &
-         & igrdi1(2,2,ipro),      &
-         & igrdi2(2,2,ipro),      &
-         & igrdj1(2,2,ipro),      &
-         & igrdj2(2,2,ipro),      &
-         & zglam1(2,2,ipro),      &
-         & zglam2(2,2,ipro),      &
-         & zgphi1(2,2,ipro),      &
-         & zgphi2(2,2,ipro),      &
-         & zmask1(2,2,kpk,ipro),  &
-         & zmask2(2,2,kpk,ipro),  &
-         & zint1(2,2,kpk,ipro),   &
-         & zint2(2,2,kpk,ipro),   &
-         & zgdept(2,2,kpk,ipro),  & 
-         & zgdepw(2,2,kpk,ipro)   & 
+         & igrdi(2,2,ipro),      &
+         & igrdj(2,2,ipro),      &
+         & zglam(2,2,ipro),      &
+         & zgphi(2,2,ipro),      &
+         & zmask(2,2,kpk,ipro),  &
+         & zint(2,2,kpk,ipro),   &
+         & zgdept(2,2,kpk,ipro), & 
+         & zgdepw(2,2,kpk,ipro)  & 
          & )
 
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
          iobs = jobs - prodatqc%nprofup
-         igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1
-         igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1
-         igrdj1(1,2,iobs) = prodatqc%mj(jobs,1)
-         igrdi1(2,1,iobs) = prodatqc%mi(jobs,1)
-         igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1
-         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1)
-         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1)
-         igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1
-         igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1
-         igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1
-         igrdj2(1,2,iobs) = prodatqc%mj(jobs,2)
-         igrdi2(2,1,iobs) = prodatqc%mi(jobs,2)
-         igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1
-         igrdi2(2,2,iobs) = prodatqc%mi(jobs,2)
-         igrdj2(2,2,iobs) = prodatqc%mj(jobs,2)
+         igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1
+         igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1
+         igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1
+         igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar)
+         igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar)
+         igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1
+         igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar)
+         igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar)
       END DO
 
@@ -261 +232 @@
       zgdepw(:,:,:,:) = 0.0
 
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 )
-      
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 )
-      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 )
-
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 
-      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam )
+      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar,   zint )
+
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 
+      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 
 
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
 
-         ALLOCATE( &
-            & zinm1(2,2,kpk,ipro),  &
-            & zinm2(2,2,kpk,ipro)   &
-            & )
-
-         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, &
-            &                  prodatqc%vdmean(:,:,:,1), zinm1 )
-         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
-            &                  prodatqc%vdmean(:,:,:,2), zinm2 )
+         ALLOCATE( zinm(2,2,kpk,ipro) )
+
+         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &
+            &                  prodatqc%vdmean(:,:,:,kvar), zinm )
 
       ENDIF
@@ -320 +281 @@
          ! Horizontal weights 
          ! Masked values are calculated later.  
-         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
+         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
 
             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
-               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), &
-               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 )
-
-         ENDIF
-
-         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
-
-            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
-               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), &
-               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2)
- 
-         ENDIF
-
-         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
+               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
+               &                   zmask(:,:,1,iobs), zweig1, zmsk )
+
+         ENDIF
+
+         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
 
             zobsk(:) = obfillflt
@@ -346 +299 @@
 
                   ! vertically interpolate all 4 corners 
-                  ista = prodatqc%npvsta(jobs,1) 
-                  iend = prodatqc%npvend(jobs,1) 
+                  ista = prodatqc%npvsta(jobs,kvar) 
+                  iend = prodatqc%npvend(jobs,kvar) 
                   inum_obs = iend - ista + 1 
                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
@@ -356 +309 @@
                         IF ( k1dint == 1 ) THEN 
                            CALL obs_int_z1d_spl( kpk, & 
-                              &     zinm1(iin,ijn,:,iobs), & 
+                              &     zinm(iin,ijn,:,iobs), & 
                               &     zobs2k, zgdept(iin,ijn,:,iobs), & 
-                              &     zmask1(iin,ijn,:,iobs)) 
+                              &     zmask(iin,ijn,:,iobs)) 
                         ENDIF 
        
                         CALL obs_level_search(kpk, & 
                            &    zgdept(iin,ijn,:,iobs), & 
-                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                           &    inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 
                            &    iv_indic) 
 
                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
-                           &    prodatqc%var(1)%vdep(ista:iend), & 
-                           &    zinm1(iin,ijn,:,iobs), & 
+                           &    prodatqc%var(kvar)%vdep(ista:iend), & 
+                           &    zinm(iin,ijn,:,iobs), & 
                            &    zobs2k, interp_corner(iin,ijn,:), & 
                            &    zgdept(iin,ijn,:,iobs), & 
-                           &    zmask1(iin,ijn,:,iobs)) 
+                           &    zmask(iin,ijn,:,iobs)) 
        
                      ENDDO 
@@ -383 +336 @@
      
                ! vertically interpolate all 4 corners 
-               ista = prodatqc%npvsta(jobs,1) 
-               iend = prodatqc%npvend(jobs,1) 
+               ista = prodatqc%npvsta(jobs,kvar) 
+               iend = prodatqc%npvend(jobs,kvar) 
                inum_obs = iend - ista + 1 
                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
@@ -392 +345 @@
                      IF ( k1dint == 1 ) THEN 
                         CALL obs_int_z1d_spl( kpk, & 
-                           &    zint1(iin,ijn,:,iobs),& 
+                           &    zint(iin,ijn,:,iobs),& 
                            &    zobs2k, zgdept(iin,ijn,:,iobs), & 
-                           &    zmask1(iin,ijn,:,iobs)) 
+                           &    zmask(iin,ijn,:,iobs)) 
   
                      ENDIF 
@@ -400 +353 @@
                      CALL obs_level_search(kpk, & 
                          &        zgdept(iin,ijn,:,iobs),& 
-                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
+                         &        inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 
                          &        iv_indic) 
 
                      CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
-                         &          prodatqc%var(1)%vdep(ista:iend),     & 
-                         &          zint1(iin,ijn,:,iobs),            & 
+                         &          prodatqc%var(kvar)%vdep(ista:iend),     & 
+                         &          zint(iin,ijn,:,iobs),            & 
                          &          zobs2k,interp_corner(iin,ijn,:), & 
                          &          zgdept(iin,ijn,:,iobs),         & 
-                         &          zmask1(iin,ijn,:,iobs) )      
+                         &          zmask(iin,ijn,:,iobs) )      
          
                   ENDDO 
@@ -432 +385 @@
                   DO ijn=1,2 
      
-                     depth_loop1: DO ik=kpk,2,-1 
-                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
+                     depth_loop: DO ik=kpk,2,-1 
+                        IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
                             
                            zweig(iin,ijn,1) = &  
                               & zweig1(iin,ijn,1) * & 
                               & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
-                              &  - prodatqc%var(1)%vdep(iend)),0._wp) 
+                              &  - prodatqc%var(kvar)%vdep(iend)),0._wp) 
                             
-                           EXIT depth_loop1 
+                           EXIT depth_loop 
 
                         ENDIF 
 
-                     ENDDO depth_loop1 
+                     ENDDO depth_loop
      
                   ENDDO 
@@ -450 +403 @@
    
                CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
-                  &              prodatqc%var(1)%vmod(iend:iend) ) 
+                  &              prodatqc%var(kvar)%vmod(iend:iend) ) 
 
                   ! Set QC flag for any observations found below the bottom
                   ! needed as the check here is more strict than that in obs_prep
-               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
+               IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
  
             ENDDO 
@@ -460 +413 @@
             DEALLOCATE(interp_corner,iv_indic) 
           
-         ENDIF 
-
-         ! For the second variable
-         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
-
-            zobsk(:) = obfillflt
-
-            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
-
-               IF ( idayend == 0 )  THEN
-                  ! Daily averaged data
-
-                  ! vertically interpolate all 4 corners 
-                  ista = prodatqc%npvsta(jobs,2) 
-                  iend = prodatqc%npvend(jobs,2) 
-                  inum_obs = iend - ista + 1 
-                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
-
-                  DO iin=1,2 
-                     DO ijn=1,2 
-
-                        IF ( k1dint == 1 ) THEN 
-                           CALL obs_int_z1d_spl( kpk, & 
-                              &     zinm2(iin,ijn,:,iobs), & 
-                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
-                              &     zmask2(iin,ijn,:,iobs)) 
-                        ENDIF 
-       
-                        CALL obs_level_search(kpk, & 
-                           &    zgdept(iin,ijn,:,iobs), & 
-                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
-                           &    iv_indic) 
-
-                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
-                           &    prodatqc%var(2)%vdep(ista:iend), & 
-                           &    zinm2(iin,ijn,:,iobs), & 
-                           &    zobs2k, interp_corner(iin,ijn,:), & 
-                           &    zgdept(iin,ijn,:,iobs), & 
-                           &    zmask2(iin,ijn,:,iobs)) 
-       
-                     ENDDO 
-                  ENDDO 
-
-               ENDIF !idayend
-
-            ELSE   
-
-               ! Point data 
-     
-               ! vertically interpolate all 4 corners 
-               ista = prodatqc%npvsta(jobs,2) 
-               iend = prodatqc%npvend(jobs,2) 
-               inum_obs = iend - ista + 1 
-               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
-               DO iin=1,2  
-                  DO ijn=1,2 
-                    
-                     IF ( k1dint == 1 ) THEN 
-                        CALL obs_int_z1d_spl( kpk, & 
-                           &    zint2(iin,ijn,:,iobs),& 
-                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
-                           &    zmask2(iin,ijn,:,iobs)) 
-  
-                     ENDIF 
-       
-                     CALL obs_level_search(kpk, & 
-                         &        zgdept(iin,ijn,:,iobs),& 
-                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
-                         &        iv_indic) 
-
-                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
-                         &          prodatqc%var(2)%vdep(ista:iend),     & 
-                         &          zint2(iin,ijn,:,iobs),            & 
-                         &          zobs2k,interp_corner(iin,ijn,:), & 
-                         &          zgdept(iin,ijn,:,iobs),         & 
-                         &          zmask2(iin,ijn,:,iobs) )      
-         
-                  ENDDO 
-               ENDDO 
-             
-            ENDIF 
-
-            !------------------------------------------------------------- 
-            ! Compute the horizontal interpolation for every profile level 
-            !------------------------------------------------------------- 
-             
-            DO ikn=1,inum_obs 
-               iend=ista+ikn-1
-                  
-               zweig(:,:,1) = 0._wp 
-   
-               ! This code forces the horizontal weights to be  
-               ! zero IF the observation is below the bottom of the  
-               ! corners of the interpolation nodes, Or if it is in  
-               ! the mask. This is important for observations near  
-               ! steep bathymetry 
-               DO iin=1,2 
-                  DO ijn=1,2 
-     
-                     depth_loop2: DO ik=kpk,2,-1 
-                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
-                            
-                           zweig(iin,ijn,1) = &  
-                              & zweig2(iin,ijn,1) * & 
-                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
-                              &  - prodatqc%var(2)%vdep(iend)),0._wp) 
-                            
-                           EXIT depth_loop2 
-
-                        ENDIF 
-
-                     ENDDO depth_loop2 
-     
-                  ENDDO 
-               ENDDO 
-   
-               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
-                  &              prodatqc%var(2)%vmod(iend:iend) ) 
-
-                  ! Set QC flag for any observations found below the bottom
-                  ! needed as the check here is more strict than that in obs_prep
-               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
- 
-            ENDDO 
- 
-            DEALLOCATE(interp_corner,iv_indic) 
-          
-         ENDIF 
+         ENDIF
 
       ENDDO
 
       ! Deallocate the data for interpolation
-      DEALLOCATE( &
-         & igrdi1, &
-         & igrdi2, &
-         & igrdj1, &
-         & igrdj2, &
-         & zglam1, &
-         & zglam2, &
-         & zgphi1, &
-         & zgphi2, &
-         & zmask1, &
-         & zmask2, &
-         & zint1,  &
-         & zint2,  &
+      DEALLOCATE(  &
+         & igrdi,  &
+         & igrdj,  &
+         & zglam,  &
+         & zgphi,  &
+         & zmask,  &
+         & zint,   &
          & zgdept, &
          & zgdepw  &
@@ -611 +431 @@
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
-         DEALLOCATE( &
-            & zinm1,  &
-            & zinm2   &
-            & )
+         DEALLOCATE( zinm )
       ENDIF
 
-      prodatqc%nprofup = prodatqc%nprofup + ipro 
+      IF ( kvar == prodatqc%nvar ) THEN
+         prodatqc%nprofup = prodatqc%nprofup + ipro 
+      ENDIF
 
    END SUBROUTINE obs_prof_opt

NEMO/trunk/src/OCE/OBS/obs_prep.F90

@@ -241 +241 @@
 
 
-   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var1, ld_var2, &
+   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, &
       &                     kpi, kpj, kpk, &
-      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  &
+      &                     zmask, pglam, pgphi,  &
       &                     ld_nea, ld_bound_reject, Kmm, kdailyavtypes,  kqc_cutoff )
 
@@ -269 +269 @@
       TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data
       TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening
-      LOGICAL, INTENT(IN) :: ld_var1              ! Observed variables switches
-      LOGICAL, INTENT(IN) :: ld_var2
+      LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: &
+         & ld_var                                 ! Observed variables switches
       LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land
       LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary
@@ -277 +277 @@
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes                          ! Types for daily averages
-      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
-         & zmask1, &
-         & zmask2
-      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & pglam1, &
-         & pglam2, &
-         & pgphi1, &
-         & pgphi2
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: &
+         & zmask
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: &
+         & pglam, &
+         & pgphi
       INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value
 
@@ -295 +292 @@
       INTEGER :: imin0
       INTEGER :: icycle       ! Current assimilation cycle
-                              ! Counters for observations that are
-      INTEGER :: iotdobs      !  - outside time domain
-      INTEGER :: iosdv1obs    !  - outside space domain (variable 1)
-      INTEGER :: iosdv2obs    !  - outside space domain (variable 2)
-      INTEGER :: ilanv1obs    !  - within a model land cell (variable 1)
-      INTEGER :: ilanv2obs    !  - within a model land cell (variable 2)
-      INTEGER :: inlav1obs    !  - close to land (variable 1)
-      INTEGER :: inlav2obs    !  - close to land (variable 2)
-      INTEGER :: ibdyv1obs    !  - boundary (variable 1) 
-      INTEGER :: ibdyv2obs    !  - boundary (variable 2)      
-      INTEGER :: igrdobs      !  - fail the grid search
-      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa
-      INTEGER :: iuvchkv      !
-                              ! Global counters for observations that are
-      INTEGER :: iotdobsmpp   !  - outside time domain
-      INTEGER :: iosdv1obsmpp !  - outside space domain (variable 1)
-      INTEGER :: iosdv2obsmpp !  - outside space domain (variable 2)
-      INTEGER :: ilanv1obsmpp !  - within a model land cell (variable 1)
-      INTEGER :: ilanv2obsmpp !  - within a model land cell (variable 2)
-      INTEGER :: inlav1obsmpp !  - close to land (variable 1)
-      INTEGER :: inlav2obsmpp !  - close to land (variable 2)
-      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1) 
-      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)      
-      INTEGER :: igrdobsmpp   !  - fail the grid search
-      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa
-      INTEGER :: iuvchkvmpp   !
+                                                       ! Counters for observations that are
+      INTEGER                           :: iotdobs     !  - outside time domain
+      INTEGER, DIMENSION(profdata%nvar) :: iosdvobs    !  - outside space domain
+      INTEGER, DIMENSION(profdata%nvar) :: ilanvobs    !  - within a model land cell
+      INTEGER, DIMENSION(profdata%nvar) :: inlavobs    !  - close to land
+      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs    !  - boundary   
+      INTEGER                           :: igrdobs     !  - fail the grid search
+      INTEGER                           :: iuvchku     !  - reject UVEL if VVEL rejected
+      INTEGER                           :: iuvchkv     !  - reject VVEL if UVEL rejected
+                                                       ! Global counters for observations that are
+      INTEGER                           :: iotdobsmpp  !  - outside time domain
+      INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp !  - outside space domain
+      INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp !  - within a model land cell
+      INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp !  - close to land
+      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp !  - boundary
+      INTEGER :: igrdobsmpp                            !  - fail the grid search
+      INTEGER :: iuvchkumpp                            !  - reject UVEL if VVEL rejected
+      INTEGER :: iuvchkvmpp                            !  - reject VVEL if UVEL rejected
       TYPE(obs_prof_valid) ::  llvalid      ! Profile selection 
       TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
-         & llvvalid           ! var1,var2 selection 
+         & llvvalid           ! var selection 
       INTEGER :: jvar         ! Variable loop variable
       INTEGER :: jobs         ! Obs. loop variable
       INTEGER :: jstp         ! Time loop variable
       INTEGER :: inrc         ! Time index variable
+      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line
+      CHARACTER(LEN=256) :: cout2  ! Diagnostic output line
       !!----------------------------------------------------------------------
 
@@ -342 +333 @@
       icycle = nn_no     ! Assimilation cycle
 
-      ! Diagnotics counters for various failures.
-
-      iotdobs   = 0
-      igrdobs   = 0
-      iosdv1obs = 0
-      iosdv2obs = 0
-      ilanv1obs = 0
-      ilanv2obs = 0
-      inlav1obs = 0
-      inlav2obs = 0
-      ibdyv1obs = 0
-      ibdyv2obs = 0
-      iuvchku   = 0
-      iuvchkv   = 0
+      ! Diagnostic counters for various failures.
+
+      iotdobs     = 0
+      igrdobs     = 0
+      iosdvobs(:) = 0
+      ilanvobs(:) = 0
+      inlavobs(:) = 0
+      ibdyvobs(:) = 0
+      iuvchku     = 0
+      iuvchkv     = 0
 
 
@@ -388 +375 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), &
-         &              profdata%nqc,     igrdobs                         )
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), &
-         &              profdata%nqc,     igrdobs                         )
+      DO jvar = 1, profdata%nvar
+         CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,jvar), profdata%mj(:,jvar), &
+            &              profdata%nqc,     igrdobs                         )
+      END DO
 
       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
@@ -406 +393 @@
       ! -----------------------------------------------------------------------
 
-      ! Variable 1
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   &
-         &                 profdata%npvsta(:,1),  profdata%npvend(:,1), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          &
-         &                 profdata%var(1)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(1)%vdep,                        &
-         &                 pglam1,                pgphi1,               &
-         &                 gdept_1d,              zmask1,               &
-         &                 profdata%nqc,          profdata%var(1)%nvqc, &
-         &                 iosdv1obs,             ilanv1obs,            &
-         &                 inlav1obs,             ld_nea,               &
-         &                 ibdyv1obs,             ld_bound_reject,      &
-         &                 iqc_cutoff,            Kmm                 )
-
-      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp )
-      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp )
-      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp )
-      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp )
-
-      ! Variable 2
-      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   &
-         &                 profdata%npvsta(:,2),  profdata%npvend(:,2), &
-         &                 jpi,                   jpj,                  &
-         &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          & 
-         &                 profdata%var(2)%mvk,                         &
-         &                 profdata%rlam,         profdata%rphi,        &
-         &                 profdata%var(2)%vdep,                        &
-         &                 pglam2,                pgphi2,               &
-         &                 gdept_1d,              zmask2,               &
-         &                 profdata%nqc,          profdata%var(2)%nvqc, &
-         &                 iosdv2obs,             ilanv2obs,            &
-         &                 inlav2obs,             ld_nea,               &
-         &                 ibdyv2obs,             ld_bound_reject,      &
-         &                 iqc_cutoff,            Kmm                 )
-
-      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp )
-      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp )
-      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp )
-      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp )
+      DO jvar = 1, profdata%nvar
+         CALL obs_coo_spc_3d( profdata%nprof,          profdata%nvprot(jvar),   &
+            &                 profdata%npvsta(:,jvar), profdata%npvend(:,jvar), &
+            &                 jpi,                     jpj,                     &
+            &                 jpk,                                              &
+            &                 profdata%mi,             profdata%mj,             &
+            &                 profdata%var(jvar)%mvk,                           &
+            &                 profdata%rlam,           profdata%rphi,           &
+            &                 profdata%var(jvar)%vdep,                          &
+            &                 pglam(:,:,jvar),         pgphi(:,:,jvar),         &
+            &                 gdept_1d,                zmask(:,:,:,jvar),       &
+            &                 profdata%nqc,            profdata%var(jvar)%nvqc, &
+            &                 iosdvobs(jvar),          ilanvobs(jvar),          &
+            &                 inlavobs(jvar),          ld_nea,                  &
+            &                 ibdyvobs(jvar),          ld_bound_reject,         &
+            &                 iqc_cutoff,              Kmm       )
+
+         CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) )
+         CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) )
+      END DO
 
       ! -----------------------------------------------------------------------
@@ -499 +465 @@
       
          WRITE(numout,*)
-         WRITE(numout,*) ' Profiles outside time domain                     = ', &
+         WRITE(numout,*) ' Profiles outside time domain                       = ', &
             &            iotdobsmpp
-         WRITE(numout,*) ' Remaining profiles that failed grid search       = ', &
+         WRITE(numout,*) ' Remaining profiles that failed grid search         = ', &
             &            igrdobsmpp
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain       = ', &
-            &            iosdv1obsmpp
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points             = ', &
-            &            ilanv1obsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',&
-               &            inlav1obsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept)    = ',&
-               &            inlav1obsmpp
-         ENDIF
-         IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
-            WRITE(numout,*) ' U observation rejected since V rejected     = ', &
-               &            iuvchku
-         ENDIF
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',&
-               &            ibdyv1obsmpp
-         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', &
-            &            prodatqc%nvprotmpp(1)
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain       = ', &
-            &            iosdv2obsmpp
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points             = ', &
-            &            ilanv2obsmpp
-         IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',&
-               &            inlav2obsmpp
-         ELSE
-            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept)    = ',&
-               &            inlav2obsmpp
-         ENDIF
-         IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
-            WRITE(numout,*) ' V observation rejected since U rejected     = ', &
-               &            iuvchkv
-         ENDIF
-         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',&
-               &            ibdyv2obsmpp
-         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', &
-            &            prodatqc%nvprotmpp(2)
+         DO jvar = 1, profdata%nvar
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain       = ', &
+               &            iosdvobsmpp(jvar)
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points             = ', &
+               &            ilanvobsmpp(jvar)
+            IF (ld_nea) THEN
+               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ',&
+                  &            inlavobsmpp(jvar)
+            ELSE
+               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept)    = ',&
+                  &            inlavobsmpp(jvar)
+            ENDIF
+            IF ( TRIM(profdata%cvars(jvar)) == 'UVEL' ) THEN
+               WRITE(numout,*) ' U observation rejected since V rejected     = ', &
+                  &            iuvchku
+            ELSE IF ( TRIM(profdata%cvars(jvar)) == 'VVEL' ) THEN
+               WRITE(numout,*) ' V observation rejected since U rejected     = ', &
+                  &            iuvchkv
+            ENDIF
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ',&
+                  &            ibdyvobsmpp(jvar)
+            WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted                             = ', &
+               &            prodatqc%nvprotmpp(jvar)
+         END DO
 
          WRITE(numout,*)
          WRITE(numout,*) ' Number of observations per time step :'
          WRITE(numout,*)
-         WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', &
-            &                               '     '//prodatqc%cvars(1)//'     ', &
-            &                               '     '//prodatqc%cvars(2)//'     '
-         WRITE(numout,998)
+         WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles'
+         WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------'
+         DO jvar = 1, prodatqc%nvar
+            WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar))
+            WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------'
+         END DO
+         WRITE(numout,*) cout1
+         WRITE(numout,*) cout2
       ENDIF
       
@@ -575 +531 @@
          DO jstp = nit000 - 1, nitend
             inrc = jstp - nit000 + 2
-            WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
-               &                    prodatqc%nvstpmpp(inrc,1), &
-               &                    prodatqc%nvstpmpp(inrc,2)
+            WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc)
+            DO jvar = 1, prodatqc%nvar
+               WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar)
+            END DO
+            WRITE(numout,*) cout1
          END DO
       ENDIF
-
-998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
-999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
 
    END SUBROUTINE obs_pre_prof

NEMO/trunk/src/OCE/OBS/obs_read_prof.F90

@@ -45 +45 @@
    SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, &
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                     ldvar1, ldvar2, ldignmis, ldsatt, &
-      &                     ldmod, kdailyavtypes )
+      &                     ldvar, ldignmis, ldsatt, &
+      &                     ldmod, cdvars, kdailyavtypes )
       !!---------------------------------------------------------------------
       !!
@@ -74 +74 @@
       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, DIMENSION(kvars), INTENT(IN) :: ldvar     ! Observed variables switches
       LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
       LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points
@@ -81 +80 @@
       REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
       REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
+      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes                ! Types of daily average observations
@@ -87 +87 @@
       CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
       CHARACTER(len=8) :: clrefdate
-      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
       INTEGER :: jvar
       INTEGER :: ji
@@ -105 +105 @@
       INTEGER :: iprof
       INTEGER :: iproftot
-      INTEGER :: ivar1t0
-      INTEGER :: ivar2t0
-      INTEGER :: ivar1t
-      INTEGER :: ivar2t
+      INTEGER, DIMENSION(kvars) :: ivart0
+      INTEGER, DIMENSION(kvars) :: ivart
       INTEGER :: ip3dt
       INTEGER :: ios
       INTEGER :: ioserrcount
-      INTEGER :: ivar1tmpp
-      INTEGER :: ivar2tmpp
+      INTEGER, DIMENSION(kvars) :: ivartmpp
       INTEGER :: ip3dtmpp
       INTEGER :: itype
       INTEGER, DIMENSION(knumfiles) :: &
          & irefdate
-      INTEGER, DIMENSION(ntyp1770+1) :: &
-         & itypvar1,    &
-         & itypvar1mpp, &
-         & itypvar2,    &
-         & itypvar2mpp 
+      INTEGER, DIMENSION(ntyp1770+1,kvars) :: &
+         & itypvar,    &
+         & itypvarmpp
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
+         & iobsi,    &
+         & iobsj,    &
+         & iproc
       INTEGER, DIMENSION(:), ALLOCATABLE :: &
-         & iobsi1,    &
-         & iobsj1,    &
-         & iproc1,    &
-         & iobsi2,    &
-         & iobsj2,    &
-         & iproc2,    &
          & iindx,    &
          & ifileidx, &
@@ -147 +140 @@
       LOGICAL :: llvalprof
       LOGICAL :: lldavtimset
+      LOGICAL :: llcycle
       TYPE(obfbdata), POINTER, DIMENSION(:) :: &
          & inpfiles
@@ -152 +146 @@
       ! Local initialization
       iprof = 0
-      ivar1t0 = 0
-      ivar2t0 = 0
+      ivart0(:) = 0
       ip3dt = 0
 
@@ -219 +212 @@
                &                ldgrid = .TRUE. )
 
-            IF ( inpfiles(jj)%nvar < 2 ) THEN
+            IF ( inpfiles(jj)%nvar /= kvars ) THEN
                CALL ctl_stop( 'Feedback format error: ', &
-                  &           ' less than 2 vars in profile file' )
+                  &           ' unexpected number of vars in profile file' )
             ENDIF
 
@@ -229 +222 @@
 
             IF ( jj == 1 ) THEN
-               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
+               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
                DO ji = 1, inpfiles(jj)%nvar
-                 clvars(ji) = inpfiles(jj)%cname(ji)
+                 clvarsin(ji) = inpfiles(jj)%cname(ji)
+                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN
+                    CALL ctl_stop( 'Feedback file variables do not match', &
+                        &           ' expected variable names for this type' )
+                 ENDIF
                END DO
             ELSE
                DO ji = 1, inpfiles(jj)%nvar
-                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
+                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
                      CALL ctl_stop( 'Feedback file variables not consistent', &
                         &           ' with previous files for this type' )
@@ -308 +305 @@
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -317 +320 @@
             ALLOCATE( zlam(inowin)  )
             ALLOCATE( zphi(inowin)  )
-            ALLOCATE( iobsi1(inowin) )
-            ALLOCATE( iobsj1(inowin) )
-            ALLOCATE( iproc1(inowin) )
-            ALLOCATE( iobsi2(inowin) )
-            ALLOCATE( iobsj2(inowin) )
-            ALLOCATE( iproc2(inowin) )
+            ALLOCATE( iobsi(inowin,kvars) )
+            ALLOCATE( iobsj(inowin,kvars) )
+            ALLOCATE( iproc(inowin,kvars) )
             inowin = 0
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -336 +342 @@
             END DO
 
-            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' )
+            ! Assume anything other than velocity is on T grid
+            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
+                  &                  iproc(:,1), 'U' )
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), &
+                  &                  iproc(:,2), 'V' )
+            ELSE
+               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
+                  &                  iproc(:,1), 'T' )
+               IF ( kvars > 1 ) THEN
+                  DO jvar = 2, kvars
+                     iobsi(:,jvar) = iobsi(:,1)
+                     iobsj(:,jvar) = iobsj(:,1)
+                     iproc(:,jvar) = iproc(:,1)
+                  END DO
+               ENDIF
             ENDIF
 
@@ -352 +363 @@
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
                   inowin = inowin + 1
-                  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')
+                  DO jvar = 1, kvars
+                     inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar)
+                     inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar)
+                     inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar)
+                  END DO
+                  IF ( kvars > 1 ) THEN
+                     DO jvar = 2, kvars
+                        IF ( inpfiles(jj)%iproc(ji,jvar) /= &
+                           & inpfiles(jj)%iproc(ji,1) ) THEN
+                           CALL ctl_stop( 'Error in obs_read_prof:', &
+                              & 'observation on different processors for different vars')
+                        ENDIF
+                     END DO
                   ENDIF
                ENDIF
             END DO
-            DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 )
+            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
 
             DO ji = 1, inpfiles(jj)%nobs
                IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+               llcycle = .TRUE.
+               DO jvar = 1, kvars
+                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                     llcycle = .FALSE.
+                     EXIT
+                  ENDIF
+               END DO
+               IF ( llcycle ) CYCLE
                IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                   & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -384 +410 @@
                   ENDIF
                   llvalprof = .FALSE.
-                  IF ( ldvar1 ) THEN
-                     loop_t_count : DO ij = 1,inpfiles(jj)%nlev
-                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
-                           & CYCLE
-                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-                           ivar1t0 = ivar1t0 + 1
-                        ENDIF
-                     END DO loop_t_count
-                  ENDIF
-                  IF ( ldvar2 ) THEN
-                     loop_s_count : DO ij = 1,inpfiles(jj)%nlev
-                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
-                           & CYCLE
-                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-                           ivar2t0 = ivar2t0 + 1
-                        ENDIF
-                     END DO loop_s_count
-                  ENDIF
-                  loop_p_count : DO ij = 1,inpfiles(jj)%nlev
+                  DO jvar = 1, kvars
+                     IF ( ldvar(jvar) ) THEN
+                        DO ij = 1,inpfiles(jj)%nlev
+                           IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
+                              & CYCLE
+                           IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                              & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
+                              ivart0(jvar) = ivart0(jvar) + 1
+                           ENDIF
+                        END DO
+                     ENDIF
+                  END DO
+                  DO ij = 1,inpfiles(jj)%nlev
                      IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                         & CYCLE
-                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &    ldvar1 ) .OR. &
-                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &     ldvar2 ) ) THEN
-                        ip3dt = ip3dt + 1
-                        llvalprof = .TRUE.
-                     ENDIF
-                  END DO loop_p_count
+                     DO jvar = 1, kvars
+                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                           &    ldvar(jvar) ) ) THEN
+                           ip3dt = ip3dt + 1
+                           llvalprof = .TRUE.
+                           EXIT
+                        ENDIF
+                     END DO
+                  END DO
 
                   IF ( llvalprof ) iprof = iprof + 1
@@ -438 +456 @@
          DO ji = 1, inpfiles(jj)%nobs
             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            llcycle = .TRUE.
+            DO jvar = 1, kvars
+               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                  llcycle = .FALSE.
+                  EXIT
+               ENDIF
+            END DO
+            IF ( llcycle ) CYCLE
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -453 +477 @@
          DO ji = 1, inpfiles(jj)%nobs
             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            llcycle = .TRUE.
+            DO jvar = 1, kvars
+               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                  llcycle = .FALSE.
+                  EXIT
+               ENDIF
+            END DO
+            IF ( llcycle ) CYCLE
             IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
                & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
@@ -470 +500 @@
       iv3dt(:) = -1
       IF (ldsatt) THEN
-         iv3dt(1) = ip3dt
-         iv3dt(2) = ip3dt
+         iv3dt(:) = ip3dt
       ELSE
-         iv3dt(1) = ivar1t0
-         iv3dt(2) = ivar2t0
+         iv3dt(:) = ivart0(:)
       ENDIF
       CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
@@ -483 +511 @@
       profdata%nprof     = 0
       profdata%nvprot(:) = 0
-      profdata%cvars(:)  = clvars(:)
+      profdata%cvars(:)  = clvarsin(:)
       iprof = 0
 
       ip3dt = 0
-      ivar1t = 0
-      ivar2t = 0
-      itypvar1   (:) = 0
-      itypvar1mpp(:) = 0
-
-      itypvar2   (:) = 0
-      itypvar2mpp(:) = 0
+      ivart(:) = 0
+      itypvar   (:,:) = 0
+      itypvarmpp(:,:) = 0
 
       ioserrcount = 0
@@ -501 +525 @@
          ji = iprofidx(iindx(jk))
 
-            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+         IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
+         llcycle = .TRUE.
+         DO jvar = 1, kvars
+            IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+               llcycle = .FALSE.
+               EXIT
+            ENDIF
+         END DO
+         IF ( llcycle ) CYCLE
 
          IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
@@ -519 +549 @@
 
             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
-            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. &
-               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE
+            llcycle = .TRUE.
+            DO jvar = 1, kvars
+               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
+                  llcycle = .FALSE.
+                  EXIT
+               ENDIF
+            END DO
+            IF ( llcycle ) CYCLE
 
             loop_prof : DO ij = 1, inpfiles(jj)%nlev
@@ -527 +563 @@
                   & CYCLE
 
-               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-
-                  llvalprof = .TRUE. 
-                  EXIT loop_prof
-
-               ENDIF
-
-               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-
-                  llvalprof = .TRUE. 
-                  EXIT loop_prof
-
-               ENDIF
+               DO jvar = 1, kvars
+                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                     & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
+
+                     llvalprof = .TRUE. 
+                     EXIT loop_prof
+
+                  ENDIF
+               END DO
 
             END DO loop_prof
@@ -573 +603 @@
 
                ! Coordinate search parameters
-               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)
+               DO jvar = 1, kvars
+                  profdata%mi  (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar)
+                  profdata%mj  (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar)
+               END DO
 
                ! Profile WMO number
@@ -616 +646 @@
                   IF (ldsatt) THEN
 
-                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &    ldvar1 ) .OR. &
-                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                        &   ldvar2 ) ) THEN
-                        ip3dt = ip3dt + 1
-                     ELSE
-                        CYCLE
+                     DO jvar = 1, kvars
+                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                           &    ldvar(jvar) ) ) THEN
+                           ip3dt = ip3dt + 1
+                           EXIT
+                        ELSE IF ( jvar == kvars ) THEN
+                           CYCLE loop_p
+                        ENDIF
+                     END DO
+
+                  ENDIF
+
+                  DO jvar = 1, kvars
+                  
+                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                       &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
+                       &    ldvar(jvar) ) .OR. ldsatt ) THEN
+
+                        IF (ldsatt) THEN
+
+                           ivart(jvar) = ip3dt
+
+                        ELSE
+
+                           ivart(jvar) = ivart(jvar) + 1
+
+                        ENDIF
+
+                        ! Depth of jvar observation
+                        profdata%var(jvar)%vdep(ivart(jvar)) = &
+                           &                inpfiles(jj)%pdep(ij,ji)
+
+                        ! Depth of jvar observation QC
+                        profdata%var(jvar)%idqc(ivart(jvar)) = &
+                           &                inpfiles(jj)%idqc(ij,ji)
+
+                        ! Depth of jvar observation QC flags
+                        profdata%var(jvar)%idqcf(:,ivart(jvar)) = &
+                           &                inpfiles(jj)%idqcf(:,ij,ji)
+
+                        ! Profile index
+                        profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof
+
+                        ! Vertical index in original profile
+                        profdata%var(jvar)%nvlidx(ivart(jvar)) = ij
+
+                        ! Profile jvar value
+                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
+                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
+                           profdata%var(jvar)%vobs(ivart(jvar)) = &
+                              &                inpfiles(jj)%pob(ij,ji,jvar)
+                           IF ( ldmod ) THEN
+                              profdata%var(jvar)%vmod(ivart(jvar)) = &
+                                 &                inpfiles(jj)%padd(ij,ji,1,jvar)
+                           ENDIF
+                           ! Count number of profile var1 data as function of type
+                           itypvar( profdata%ntyp(iprof) + 1, jvar ) = &
+                              & itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1
+                        ELSE
+                           profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi
+                        ENDIF
+
+                        ! Profile jvar qc
+                        profdata%var(jvar)%nvqc(ivart(jvar)) = &
+                           & inpfiles(jj)%ivlqc(ij,ji,jvar)
+
+                        ! Profile jvar qc flags
+                        profdata%var(jvar)%nvqcf(:,ivart(jvar)) = &
+                           & inpfiles(jj)%ivlqcf(:,ij,ji,jvar)
+
+                        ! Profile insitu T value
+                        IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN
+                           profdata%var(jvar)%vext(ivart(jvar),1) = &
+                              &                inpfiles(jj)%pext(ij,ji,1)
+                        ENDIF
+
                      ENDIF
-
-                  ENDIF
-
-                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                    &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
-                    &    ldvar1 ) .OR. ldsatt ) THEN
-
-                     IF (ldsatt) THEN
-
-                        ivar1t = ip3dt
-
-                     ELSE
-
-                        ivar1t = ivar1t + 1
-
-                     ENDIF
-
-                     ! Depth of var1 observation
-                     profdata%var(1)%vdep(ivar1t) = &
-                        &                inpfiles(jj)%pdep(ij,ji)
-
-                     ! Depth of var1 observation QC
-                     profdata%var(1)%idqc(ivar1t) = &
-                        &                inpfiles(jj)%idqc(ij,ji)
-
-                     ! Depth of var1 observation QC flags
-                     profdata%var(1)%idqcf(:,ivar1t) = &
-                        &                inpfiles(jj)%idqcf(:,ij,ji)
-
-                     ! Profile index
-                     profdata%var(1)%nvpidx(ivar1t) = iprof
-
-                     ! Vertical index in original profile
-                     profdata%var(1)%nvlidx(ivar1t) = ij
-
-                     ! Profile var1 value
-                     IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. &
-                        & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
-                        profdata%var(1)%vobs(ivar1t) = &
-                           &                inpfiles(jj)%pob(ij,ji,1)
-                        IF ( ldmod ) THEN
-                           profdata%var(1)%vmod(ivar1t) = &
-                              &                inpfiles(jj)%padd(ij,ji,1,1)
-                        ENDIF
-                        ! 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(ivar1t) = fbrmdi
-                     ENDIF
-
-                     ! Profile var1 qc
-                     profdata%var(1)%nvqc(ivar1t) = &
-                        & inpfiles(jj)%ivlqc(ij,ji,1)
-
-                     ! Profile var1 qc flags
-                     profdata%var(1)%nvqcf(:,ivar1t) = &
-                        & inpfiles(jj)%ivlqcf(:,ij,ji,1)
-
-                     ! Profile insitu T value
-                     IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN
-                        profdata%var(1)%vext(ivar1t,1) = &
-                           &                inpfiles(jj)%pext(ij,ji,1)
-                     ENDIF
-
-                  ENDIF
-
-                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. &
-                     &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    .AND. &
-                     &   ldvar2 ) .OR. ldsatt ) THEN
-
-                     IF (ldsatt) THEN
-
-                        ivar2t = ip3dt
-
-                     ELSE
-
-                        ivar2t = ivar2t + 1
-
-                     ENDIF
-
-                     ! Depth of var2 observation
-                     profdata%var(2)%vdep(ivar2t) = &
-                        &                inpfiles(jj)%pdep(ij,ji)
-
-                     ! Depth of var2 observation QC
-                     profdata%var(2)%idqc(ivar2t) = &
-                        &                inpfiles(jj)%idqc(ij,ji)
-
-                     ! Depth of var2 observation QC flags
-                     profdata%var(2)%idqcf(:,ivar2t) = &
-                        &                inpfiles(jj)%idqcf(:,ij,ji)
-
-                     ! Profile index
-                     profdata%var(2)%nvpidx(ivar2t) = iprof
-
-                     ! Vertical index in original profile
-                     profdata%var(2)%nvlidx(ivar2t) = ij
-
-                     ! Profile var2 value
-                  IF (  ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. &
-                    &   ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    )  ) THEN
-                        profdata%var(2)%vobs(ivar2t) = &
-                           &                inpfiles(jj)%pob(ij,ji,2)
-                        IF ( ldmod ) THEN
-                           profdata%var(2)%vmod(ivar2t) = &
-                              &                inpfiles(jj)%padd(ij,ji,1,2)
-                        ENDIF
-                        ! 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(ivar2t) = fbrmdi
-                     ENDIF
-
-                     ! Profile var2 qc
-                     profdata%var(2)%nvqc(ivar2t) = &
-                        & inpfiles(jj)%ivlqc(ij,ji,2)
-
-                     ! Profile var2 qc flags
-                     profdata%var(2)%nvqcf(:,ivar2t) = &
-                        & inpfiles(jj)%ivlqcf(:,ij,ji,2)
-
-                  ENDIF
+                  
+                  END DO
 
                END DO loop_p
@@ -763 +739 @@
       !-----------------------------------------------------------------------
 
-      CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp )
-      CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp )
+      DO jvar = 1, kvars
+         CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) )
+      END DO
       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 )
+      DO jvar = 1, kvars
+         CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 )
+      END DO
 
       !-----------------------------------------------------------------------
@@ -778 +756 @@
          WRITE(numout,'(1X,A)') '------------'
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) )
-         WRITE(numout,'(1X,A)') '------------------------'
-         DO ji = 0, ntyp1770
-            IF ( itypvar1mpp(ji+1) > 0 ) THEN
-               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
-                  & cwmonam1770(ji)(1:52),' = ', &
-                  & itypvar1mpp(ji+1)
-            ENDIF
+         DO jvar = 1, kvars
+            WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) )
+            WRITE(numout,'(1X,A)') '------------------------'
+            DO ji = 0, ntyp1770
+               IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN
+                  WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
+                     & cwmonam1770(ji)(1:52),' = ', &
+                     & itypvarmpp(ji+1,jvar)
+               ENDIF
+            END DO
+            WRITE(numout,'(1X,A)') &
+               & '---------------------------------------------------------------'
+            WRITE(numout,'(1X,A55,I8)') &
+               & 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// &
+               & '             = ', ivartmpp(jvar)
+            WRITE(numout,'(1X,A)') &
+               & '---------------------------------------------------------------'
+            WRITE(numout,*) 
          END DO
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// &
-            & '             = ', ivar1tmpp
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) )
-         WRITE(numout,'(1X,A)') '------------------------'
-         DO ji = 0, ntyp1770
-            IF ( itypvar2mpp(ji+1) > 0 ) THEN
-               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
-                  & cwmonam1770(ji)(1:52),' = ', &
-                  & itypvar2mpp(ji+1)
-            ENDIF
+      ENDIF
+
+      IF (ldsatt) THEN
+         profdata%nvprot(:)    = ip3dt
+         profdata%nvprotmpp(:) = ip3dtmpp
+      ELSE
+         DO jvar = 1, kvars
+            profdata%nvprot(jvar)    = ivart(jvar)
+            profdata%nvprotmpp(jvar) = ivartmpp(jvar)
          END DO
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// &
-            & '             = ', ivar2tmpp
-         WRITE(numout,'(1X,A)') &
-            & '---------------------------------------------------------------'
-         WRITE(numout,*) 
-      ENDIF
-
-      IF (ldsatt) THEN
-         profdata%nvprot(1)    = ip3dt
-         profdata%nvprot(2)    = ip3dt
-         profdata%nvprotmpp(1) = ip3dtmpp
-         profdata%nvprotmpp(2) = ip3dtmpp
-      ELSE
-         profdata%nvprot(1)    = ivar1t
-         profdata%nvprot(2)    = ivar2t
-         profdata%nvprotmpp(1) = ivar1tmpp
-         profdata%nvprotmpp(2) = ivar2tmpp
       ENDIF
       profdata%nprof        = iprof
@@ -830 +791 @@
       ! Model level search
       !-----------------------------------------------------------------------
-      IF ( ldvar1 ) THEN
-         CALL obs_level_search( jpk, gdept_1d, &
-            & profdata%nvprot(1), profdata%var(1)%vdep, &
-            & profdata%var(1)%mvk )
-      ENDIF
-      IF ( ldvar2 ) THEN
-         CALL obs_level_search( jpk, gdept_1d, &
-            & profdata%nvprot(2), profdata%var(2)%vdep, &
-            & profdata%var(2)%mvk )
-      ENDIF
+      DO jvar = 1, kvars
+         IF ( ldvar(jvar) ) THEN
+            CALL obs_level_search( jpk, gdept_1d, &
+               & profdata%nvprot(jvar), profdata%var(jvar)%vdep, &
+               & profdata%var(jvar)%mvk )
+         ENDIF
+      END DO
 
       !-----------------------------------------------------------------------
@@ -852 +810 @@
       ! Deallocate temporary data
       !-----------------------------------------------------------------------
-      DEALLOCATE( ifileidx, iprofidx, zdat, clvars )
+      DEALLOCATE( ifileidx, iprofidx, zdat, clvarsin )
 
       !-----------------------------------------------------------------------

NEMO/trunk/src/OCE/OBS/obs_read_surf.F90

@@ -40 +40 @@
    SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                     ldignmis, ldmod, ldnightav )
+      &                     ldignmis, ldmod, ldnightav, cdvars )
       !!---------------------------------------------------------------------
       !!
@@ -73 +73 @@
       REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
       REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
+      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
 
       !! * Local declarations
       CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
       CHARACTER(len=8) :: clrefdate
-      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars
+      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
       INTEGER :: ji
       INTEGER :: jj
@@ -178 +179 @@
                &                ldgrid = .TRUE. )
 
+            IF ( inpfiles(jj)%nvar /= kvars ) THEN
+               CALL ctl_stop( 'Feedback format error: ', &
+                  &           ' unexpected number of vars in feedback file' )
+            ENDIF
+
             IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
                CALL ctl_stop( 'Model not in input data' )
@@ -184 +190 @@
 
             IF ( jj == 1 ) THEN
-               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
+               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
                DO ji = 1, inpfiles(jj)%nvar
-                 clvars(ji) = inpfiles(jj)%cname(ji)
+                 clvarsin(ji) = inpfiles(jj)%cname(ji)
+                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN
+                    CALL ctl_stop( 'Feedback file variables do not match', &
+                        &           ' expected variable names for this type' )
+                 ENDIF
                END DO
             ELSE
                DO ji = 1, inpfiles(jj)%nvar
-                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
+                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
                      CALL ctl_stop( 'Feedback file variables not consistent', &
                         &           ' with previous files for this type' )
@@ -347 +357 @@
       iobs = 0
 
-      surfdata%cvars(:)  = clvars(:)
+      surfdata%cvars(:)  = clvarsin(:)
 
       ityp   (:) = 0
@@ -480 +490 @@
       ! Deallocate temporary data
       !-----------------------------------------------------------------------
-      DEALLOCATE( ifileidx, isurfidx, zdat, clvars )
+      DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin )
 
       !-----------------------------------------------------------------------

NEMO/trunk/src/OCE/OBS/obs_write.F90

@@ -86 +86 @@
       CHARACTER(LEN=40) :: clfname
       CHARACTER(LEN=10) :: clfiletype
+      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
+      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
+      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
       CHARACTER(LEN=12) :: clfmt            ! writing format
       INTEGER :: idg                        ! number of digits
@@ -115 +118 @@
       ! Find maximum level
       ilevel = 0
-      DO jvar = 1, 2
+      DO jvar = 1, profdata%nvar
          ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
       END DO
@@ -180 +183 @@
 
       END SELECT
+      
+      IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. &
+         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
+         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
+            &                 1 + iadd, iext, .TRUE. )
+         fbdata%cname(1)      = profdata%cvars(1)
+         fbdata%coblong(1)    = cllongname
+         fbdata%cobunit(1)    = clunits
+         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
+         fbdata%caddunit(1,1) = clunits
+         fbdata%cgrid(:)      = clgrid
+         DO je = 1, iext
+            fbdata%cextname(je) = pext%cdname(je)
+            fbdata%cextlong(je) = pext%cdlong(je,1)
+            fbdata%cextunit(je) = pext%cdunit(je,1)
+         END DO
+         DO ja = 1, iadd
+            fbdata%caddname(1+ja) = padd%cdname(ja)
+            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
+            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
+         END DO
+      ENDIF
 
       fbdata%caddname(1)   = 'Hx'
@@ -234 +259 @@
             &           krefdate = 19500101 )
          ! Reform the profiles arrays for output
-         DO jvar = 1, 2
+         DO jvar = 1, profdata%nvar
             DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
                ik = profdata%var(jvar)%nvlidx(jk)
@@ -329 +354 @@
       CHARACTER(LEN=40) :: clfname         ! netCDF filename
       CHARACTER(LEN=10) :: clfiletype
+      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
+      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
+      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
       CHARACTER(LEN=12) :: clfmt           ! writing format
@@ -354 +382 @@
       SELECT CASE ( TRIM(surfdata%cvars(1)) )
       CASE('SLA')
+         
+         ! SLA needs special treatment because of MDT, so is all done here
+         ! Other variables are done more generically
+         ! No climatology for SLA, MDT is our best estimate of that and is already output.
 
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
@@ -384 +416 @@
       CASE('SST')
 
+         clfiletype = 'sstfb'
+         cllongname = 'Sea surface temperature'
+         clunits    = 'Degree centigrade'
+         clgrid     = 'T'
+         
+      CASE('ICECONC')
+
+         clfiletype = 'sicfb'
+         cllongname = 'Sea ice concentration'
+         clunits    = 'Fraction'
+         clgrid     = 'T'
+
+      CASE('SSS')
+
+         clfiletype = 'sssfb'
+         cllongname = 'Sea surface salinity'
+         clunits    = 'psu'
+         clgrid     = 'T'
+
+      CASE DEFAULT
+
+         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
+
+      END SELECT
+
+      ! SLA needs special treatment because of MDT, so is done above
+      ! Remaining variables treated more generically
+
+      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
+      
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
             &                 1 + iadd, iext, .TRUE. )
 
-         clfiletype = 'sstfb'
          fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'Sea surface temperature'
-         fbdata%cobunit(1)    = 'Degree centigrade'
+         fbdata%coblong(1)    = cllongname
+         fbdata%cobunit(1)    = clunits
          DO je = 1, iext
             fbdata%cextname(je) = pext%cdname(je)
             fbdata%cextlong(je) = pext%cdlong(je,1)
             fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated SST'
-         fbdata%caddunit(1,1) = 'Degree centigrade'
-         fbdata%cgrid(1)      = 'T'
+         END DO        
+         IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
+            fbdata%caddlong(1,1) = 'Model interpolated ICE'
+         ELSE
+            fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1))
+         ENDIF
+         fbdata%caddunit(1,1) = clunits
+         fbdata%cgrid(1)      = clgrid
          DO ja = 1, iadd
             fbdata%caddname(1+ja) = padd%cdname(ja)
@@ -404 +469 @@
             fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
          END DO
-
-      CASE('ICECONC')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'sicfb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'Sea ice'
-         fbdata%cobunit(1)    = 'Fraction'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated ICE'
-         fbdata%caddunit(1,1) = 'Fraction'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE('SSS')
-
-         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd, iext, .TRUE. )
-
-         clfiletype = 'sssfb'
-         fbdata%cname(1)      = surfdata%cvars(1)
-         fbdata%coblong(1)    = 'Sea surface salinity'
-         fbdata%cobunit(1)    = 'psu'
-         DO je = 1, iext
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated SSS'
-         fbdata%caddunit(1,1) = 'psu'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
-      CASE DEFAULT
-
-         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
-
-      END SELECT
+      ENDIF
 
       fbdata%caddname(1)   = 'Hx'