Index: /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 9210)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 9211)
@@ -747,5 +747,4 @@
                &               kdailyavtypes = nn_profdavtypes )
             
-            ! Is allocating and deallocating repeatedly in a loop good practice?
             DEALLOCATE( llvar )
             CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zglam )
@@ -889,32 +888,16 @@
       INTEGER :: ji, jj            ! Loop counters
       REAL(wp) :: tiny             ! small number
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & zprofvar1, &            ! Model values for 1st variable in a prof ob
-         & zprofvar2               ! Model values for 2nd variable in a prof ob
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
-         & zprofmask1, &           ! Mask associated with zprofvar1
-         & zprofmask2              ! Mask associated with zprofvar2
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
+         & zprofvar                ! Model values for variables in a prof ob
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
+         & zprofmask               ! Mask associated with zprofvar
       REAL(wp), POINTER, DIMENSION(:,:) :: &
          & zsurfvar, &             ! Model values equivalent to surface ob.
          & zsurfmask               ! Mask associated with surface variable
-      REAL(wp), POINTER, DIMENSION(:,:) :: &
-         & zglam1,    &            ! Model longitudes for prof variable 1
-         & zglam2,    &            ! Model longitudes for prof variable 2
-         & zgphi1,    &            ! Model latitudes for prof variable 1
-         & zgphi2                  ! Model latitudes for prof variable 2
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
+         & zglam,    &             ! Model longitudes for prof variables
+         & zgphi                   ! Model latitudes for prof variables
       LOGICAL :: llog10            ! Perform log10 transform of variable
 
-
-      !Allocate local work arrays
-      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
-      CALL wrk_alloc( jpi, jpj, zsurfvar )
-      CALL wrk_alloc( jpi, jpj, zsurfmask )
-      CALL wrk_alloc( jpi, jpj, zglam1 )
-      CALL wrk_alloc( jpi, jpj, zglam2 )
-      CALL wrk_alloc( jpi, jpj, zgphi1 )
-      CALL wrk_alloc( jpi, jpj, zgphi2 )
 
       IF(lwp) THEN
@@ -935,87 +918,87 @@
          DO jtype = 1, nproftypes
 
+            ! Allocate local work arrays
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
+            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
+            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
+            
+            ! 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(:,:,:) = tsn(:,:,:,jp_tem)
-               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zprofmask2(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zglam2(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
-               zgphi2(:,:) = gphit(:,:)
+               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
+               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
+
             CASE('vel')
-               zprofvar1(:,:,:) = un(:,:,:)
-               zprofvar2(:,:,:) = vn(:,:,:)
-               zprofmask1(:,:,:) = umask(:,:,:)
-               zprofmask2(:,:,:) = vmask(:,:,:)
-               zglam1(:,:) = glamu(:,:)
-               zglam2(:,:) = glamv(:,:)
-               zgphi1(:,:) = gphiu(:,:)
-               zgphi2(:,:) = gphiv(:,:)
+               zprofvar(:,:,:,1) = un(:,:,:)
+               zprofvar(:,:,:,2) = vn(:,:,:)
+               zprofmask(:,:,:,1) = umask(:,:,:)
+               zprofmask(:,:,:,2) = vmask(:,:,:)
+               zglam(:,:,1) = glamu(:,:)
+               zglam(:,:,2) = glamv(:,:)
+               zgphi(:,:,1) = gphiu(:,:)
+               zgphi(:,:,2) = gphiv(:,:)
 
             CASE('plchltot')
 #if defined key_hadocc
                ! Chlorophyll from HadOCC
-               zprofvar1(:,:,:) = HADOCC_CHL(:,:,:)
+               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
 #elif defined key_medusa && defined key_foam_medusa
                ! Add non-diatom and diatom chlorophyll from MEDUSA
-               zprofvar1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
+               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
 #elif defined key_fabm
                ! Add all chlorophyll groups from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &
-                  &               trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &
+                  &                trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)
 #else
                CALL ctl_stop( ' Trying to run plchltot observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
                ! Take the log10 where we can, otherwise exclude
                tiny = 1.0e-20
-               WHERE(zprofvar1(:,:,:) > tiny .AND. zprofvar1(:,:,:) /= obfillflt )
-                  zprofvar1(:,:,:)  = LOG10(zprofvar1(:,:,:))
+               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt )
+                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:))
                ELSEWHERE
-                  zprofvar1(:,:,:)  = obfillflt
-                  zprofmask1(:,:,:) = 0
+                  zprofvar(:,:,:,:)  = obfillflt
+                  zprofmask(:,:,:,:) = 0
                END WHERE
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('pchltot')
 #if defined key_hadocc
                ! Chlorophyll from HadOCC
-               zprofvar1(:,:,:) = HADOCC_CHL(:,:,:)
+               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
 #elif defined key_medusa && defined key_foam_medusa
                ! Add non-diatom and diatom chlorophyll from MEDUSA
-               zprofvar1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
+               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
 #elif defined key_fabm
                ! Add all chlorophyll groups from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &
-                  &               trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &
+                  &                trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)
 #else
                CALL ctl_stop( ' Trying to run pchltot observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('pno3')
 #if defined key_hadocc
                ! Dissolved inorganic nitrogen from HadOCC
-               zprofvar1(:,:,:) = trn(:,:,:,jp_had_nut)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut)
 #elif defined key_medusa && defined key_foam_medusa
                ! Dissolved inorganic nitrogen from MEDUSA
-               zprofvar1(:,:,:) = trn(:,:,:,jpdin)
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdin)
 #elif defined key_fabm
                ! Nitrate from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n3n)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n3n)
 #else
                CALL ctl_stop( ' Trying to run pno3 observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('psi4')
@@ -1025,15 +1008,12 @@
 #elif defined key_medusa && defined key_foam_medusa
                ! Silicate from MEDUSA
-               zprofvar1(:,:,:) = trn(:,:,:,jpsil)
+               zprofvar(:,:,:,1) = trn(:,:,:,jpsil)
 #elif defined key_fabm
                ! Silicate from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n5s)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n5s)
 #else
                CALL ctl_stop( ' Trying to run psi4 observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('ppo4')
@@ -1046,48 +1026,39 @@
 #elif defined key_fabm
                ! Phosphate from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n1p)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n1p)
 #else
                CALL ctl_stop( ' Trying to run ppo4 observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('pdic')
 #if defined key_hadocc
                ! Dissolved inorganic carbon from HadOCC
-               zprofvar1(:,:,:) = trn(:,:,:,jp_had_dic)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic)
 #elif defined key_medusa && defined key_foam_medusa
                ! Dissolved inorganic carbon from MEDUSA
-               zprofvar1(:,:,:) = trn(:,:,:,jpdic)
+               zprofvar(:,:,:,1) = trn(:,:,:,jpdic)
 #elif defined key_fabm
                ! Dissolved inorganic carbon from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3c)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3c)
 #else
                CALL ctl_stop( ' Trying to run pdic observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('palk')
 #if defined key_hadocc
                ! Alkalinity from HadOCC
-               zprofvar1(:,:,:) = trn(:,:,:,jp_had_alk)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk)
 #elif defined key_medusa && defined key_foam_medusa
                ! Alkalinity from MEDUSA
-               zprofvar1(:,:,:) = trn(:,:,:,jpalk)
+               zprofvar(:,:,:,1) = trn(:,:,:,jpalk)
 #elif defined key_fabm
                ! Alkalinity from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3a)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3a)
 #else
                CALL ctl_stop( ' Trying to run palk observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('pph')
@@ -1097,15 +1068,12 @@
 #elif defined key_medusa && defined key_foam_medusa
                ! pH from MEDUSA
-               zprofvar1(:,:,:) = f3_pH(:,:,:)
+               zprofvar(:,:,:,1) = f3_pH(:,:,:)
 #elif defined key_fabm
                ! pH from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3ph)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3ph)
 #else
                CALL ctl_stop( ' Trying to run pph observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE('po2')
@@ -1115,35 +1083,33 @@
 #elif defined key_medusa && defined key_foam_medusa
                ! Oxygen from MEDUSA
-               zprofvar1(:,:,:) = trn(:,:,:,jpoxy)
+               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)
 #elif defined key_fabm
                ! Oxygen from ERSEM
-               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o2o)
+               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o2o)
 #else
                CALL ctl_stop( ' Trying to run po2 observation operator', &
                   &           ' but no biogeochemical model appears to have been defined' )
 #endif
-               zprofmask1(:,:,:) = tmask(:,:,:)
-               zglam1(:,:) = glamt(:,:)
-               zgphi1(:,:) = gphit(:,:)
 
             CASE DEFAULT
                CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
+
             END SELECT
-            
-            IF ( ( TRIM(cobstypesprof(jtype)) /= 'prof' ) .AND. ( TRIM(cobstypesprof(jtype)) /= 'vel' ) ) THEN
-               zprofvar2(:,:,:)  = zprofvar1(:,:,:)
-               zprofmask2(:,:,:) = zprofmask1(:,:,:)
-               zglam2(:,:)       = zglam1(:,:)
-               zgphi2(:,:)       = zgphi1(:,:)
-            ENDIF
-
-            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
-               &               nit000, idaystp,                         &
-               &               zprofvar1, zprofvar2,                    &
-               &               fsdept(:,:,:), fsdepw(:,:,:),            & 
-               &               zprofmask1, zprofmask2,                  &
-               &               zglam1, zglam2, zgphi1, zgphi2,          &
-               &               nn_1dint, nn_2dint_default,              &
-               &               kdailyavtypes = nn_profdavtypes )
+
+            DO jvar = 1, profdataqc(jtype)%nvar
+               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
+                  &               nit000, idaystp, jvar,                   &
+                  &               zprofvar(:,:,:,jvar),                    &
+                  &               fsdept(:,:,:), fsdepw(:,:,:),            & 
+                  &               zprofmask(:,:,:,jvar),                   &
+                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &
+                  &               nn_1dint, nn_2dint_default,              &
+                  &               kdailyavtypes = nn_profdavtypes )
+            END DO
+
+            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
+            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
+            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
+            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
 
          END DO
@@ -1152,4 +1118,8 @@
 
       IF ( nsurftypes > 0 ) THEN
+
+         !Allocate local work arrays
+         CALL wrk_alloc( jpi, jpj, zsurfvar )
+         CALL wrk_alloc( jpi, jpj, zsurfmask )
 
          DO jtype = 1, nsurftypes
@@ -1414,16 +1384,8 @@
          END DO
 
+         CALL wrk_dealloc( jpi, jpj, zsurfvar )
+         CALL wrk_dealloc( jpi, jpj, zsurfmask )
+
       ENDIF
-
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
-      CALL wrk_dealloc( jpi, jpj, zsurfvar )
-      CALL wrk_dealloc( jpi, jpj, zsurfmask )
-      CALL wrk_dealloc( jpi, jpj, zglam1 )
-      CALL wrk_dealloc( jpi, jpj, zglam2 )
-      CALL wrk_dealloc( jpi, jpj, zgphi1 )
-      CALL wrk_dealloc( jpi, jpj, zgphi2 )
 
    END SUBROUTINE dia_obs
Index: /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 9210)
+++ /branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 9211)
@@ -36,5 +36,6 @@
       & sbc_dcy, nday_qsr
    USE obs_grid,      ONLY : & 
-      & obs_level_search     
+      & obs_level_search   
+   USE dom_oce                  ! Ocean space and time domain variables  
 
    IMPLICIT NONE
@@ -60,9 +61,9 @@
 
 
-   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 )
 
@@ -134,14 +135,11 @@
       INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header)
       INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
+      INTEGER, INTENT(IN) :: kvar    ! Number of variable in prodatqc
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
-         & pvar1,    &               ! Model field 1
-         & pvar2,    &               ! Model field 2
-         & pmask1,   &               ! Land-sea mask 1
-         & pmask2                    ! Land-sea mask 2
+         & pvar,   &                 ! Model field for variable
+         & pmask                     ! Land-sea mask for variable
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & plam1,    &               ! Model longitudes for variable 1
-         & plam2,    &               ! Model longitudes for variable 2
-         & pphi1,    &               ! Model latitudes for variable 1
-         & pphi2                     ! Model latitudes for variable 2
+         & plam,   &                 ! Model longitudes for variable
+         & pphi                      ! Model latitudes for variable
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
          & pgdept, &                 ! Model array of depth T levels 
@@ -166,8 +164,6 @@
          & idailyavtypes
       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi1, &
-         & igrdi2, &
-         & igrdj1, &
-         & igrdj2
+         & igrdi, &
+         & igrdj
       INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
 
@@ -176,27 +172,19 @@
       REAL(KIND=wp) :: zdaystp
       REAL(KIND=wp), DIMENSION(kpk) :: &
-         & zobsmask1, &
-         & zobsmask2, &
          & 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
 
@@ -230,8 +218,5 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     prodatqc%vdmean(ji,jj,jk,1) = 0.0
-                     IF ( prodatqc%nvar == 2 ) THEN
-                        prodatqc%vdmean(ji,jj,jk,2) = 0.0
-                     ENDIF
+                     prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
                   END DO
                END DO
@@ -242,12 +227,7 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! Increment field 1 for computing daily mean
-                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                     &                        + pvar1(ji,jj,jk)
-                  IF ( prodatqc%nvar == 2 ) THEN
-                     ! Increment field 2 for computing daily mean
-                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                        &                        + pvar2(ji,jj,jk)
-                  ENDIF
+                  ! Increment field for computing daily mean
+                  prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
+                     &                           + pvar(ji,jj,jk)
                END DO
             END DO
@@ -262,10 +242,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                        &                        * zdaystp
-                     IF ( prodatqc%nvar == 2 ) THEN
-                        prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                           &                        * zdaystp
-                     ENDIF
+                     prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
+                        &                           * zdaystp
                   END DO
                END DO
@@ -277,16 +253,10 @@
       ! 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),   &
+         & 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)   & 
@@ -295,22 +265,12 @@
       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)
-         IF ( prodatqc%nvar == 2 ) THEN
-            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)
-         ENDIF
+         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
 
@@ -319,33 +279,18 @@
       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 )
-      
-      IF ( prodatqc%nvar == 2 ) THEN
-         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 )
-      ENDIF
-
-      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, 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 )
-         IF ( prodatqc%nvar == 2 ) THEN
-            CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
-               &                  prodatqc%vdmean(:,:,:,2), zinm2 )
-         ENDIF
+         ALLOCATE( zinm(2,2,kpk,ipro) )
+
+         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &
+            &                  prodatqc%vdmean(:,:,:,kvar), zinm )
 
       ENDIF
@@ -382,23 +327,13 @@
          ! 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%nvar == 2 ) THEN
-            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
-         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
@@ -420,20 +355,20 @@
                         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 
@@ -447,6 +382,6 @@
      
                ! 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)) 
@@ -456,7 +391,7 @@
                      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 
@@ -464,13 +399,13 @@
                      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 
@@ -496,17 +431,17 @@
                   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 
@@ -514,9 +449,9 @@
    
                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 
@@ -524,151 +459,16 @@
             DEALLOCATE(interp_corner,iv_indic) 
           
-         ENDIF 
-
-         IF ( prodatqc%nvar == 2 ) THEN
-            ! 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 
+         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  &
@@ -677,11 +477,10 @@
       ! 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
