Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 11437)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 11449)
@@ -62,6 +62,6 @@
    SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
       &                     kit000, kdaystp, kvar,       &
-      &                     pvar, pgdept, pgdepw,        &
-      &                     pmask,                       &  
+      &                     pvar, pclim,                 &
+      &                     pgdept, pgdepw, pmask,       &  
       &                     plam, pphi,                  &
       &                     k1dint, k2dint, kdailyavtypes )
@@ -137,4 +137,5 @@
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
          & pvar,   &                 ! Model field for variable
+         & pclim,  &                 ! Climatology field for variable         
          & pmask                     ! Land-sea mask for variable
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
@@ -172,5 +173,6 @@
       REAL(KIND=wp), DIMENSION(kpk) :: &
          & zobsk,    &
-         & zobs2k
+         & zobs2k,   &
+         & zclm2k         
       REAL(KIND=wp), DIMENSION(2,2,1) :: &
          & zweig1, &
@@ -178,4 +180,5 @@
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
          & zmask,  &
+         & zclim,  &         
          & zint,   &
          & zinm,   &
@@ -189,4 +192,5 @@
 
       LOGICAL :: ld_dailyav
+      LOGICAL :: ld_clim
 
       !------------------------------------------------------------------------
@@ -196,4 +200,11 @@
       inrc = kt - kit000 + 2
       ipro = prodatqc%npstp(inrc)
+
+      ! Check if climatology is available and set flag
+      IF ( SUM( pclim(:,:,:) ) == 0. ) THEN
+         ld_clim = .FALSE.
+      ELSE
+         ld_clim = .TRUE.      
+      ENDIF
 
       ! Daily average types
@@ -262,4 +273,6 @@
          & )
 
+      IF ( ld_clim ) ALLOCATE( zclim(2,2,kpk,ipro) )
+
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
          iobs = jobs - prodatqc%nprofup
@@ -286,4 +299,8 @@
       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 
 
+      IF ( ld_clim ) THEN
+         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim )            
+      ENDIF 
+      
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
@@ -349,5 +366,6 @@
                   inum_obs = iend - ista + 1 
                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
-
+                  IF ( ld_clim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
+                  
                   DO iin=1,2 
                      DO ijn=1,2 
@@ -358,4 +376,12 @@
                               &     zobs2k, zgdept(iin,ijn,:,iobs), & 
                               &     zmask(iin,ijn,:,iobs)) 
+
+                           IF ( ld_clim ) THEN
+                              CALL obs_int_z1d_spl( kpk, & 
+                                 &     zclim(iin,ijn,:,iobs), & 
+                                 &     zclm2k, zgdept(iin,ijn,:,iobs), & 
+                                 &     zmask(iin,ijn,:,iobs)) 
+                           ENDIF
+
                         ENDIF 
        
@@ -371,5 +397,14 @@
                            &    zgdept(iin,ijn,:,iobs), & 
                            &    zmask(iin,ijn,:,iobs)) 
-       
+
+                        IF ( ld_clim ) THEN
+                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
+                              &    prodatqc%var(kvar)%vdep(ista:iend), & 
+                              &    zclim(iin,ijn,:,iobs), & 
+                              &    zclm2k, interp_corner_clim(iin,ijn,:), & 
+                              &    zgdept(iin,ijn,:,iobs), & 
+                              &    zmask(iin,ijn,:,iobs)) 
+                        ENDIF
+                        
                      ENDDO 
                   ENDDO 
@@ -385,5 +420,7 @@
                iend = prodatqc%npvend(jobs,kvar) 
                inum_obs = iend - ista + 1 
-               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
+               ALLOCATE( interp_corner(2,2,inum_obs),      &
+                  &      iv_indic(inum_obs) ) 
+               IF ( ld_clim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )                  
                DO iin=1,2  
                   DO ijn=1,2 
@@ -394,4 +431,11 @@
                            &    zobs2k, zgdept(iin,ijn,:,iobs), & 
                            &    zmask(iin,ijn,:,iobs)) 
+
+                        IF ( ld_clim ) THEN
+                           CALL obs_int_z1d_spl( kpk, & 
+                              &    zclim(iin,ijn,:,iobs),& 
+                              &    zclm2k, zgdept(iin,ijn,:,iobs), & 
+                              &    zmask(iin,ijn,:,iobs)) 
+                        ENDIF
   
                      ENDIF 
@@ -408,4 +452,13 @@
                          &          zgdept(iin,ijn,:,iobs),         & 
                          &          zmask(iin,ijn,:,iobs) )      
+
+                     IF ( ld_clim ) THEN
+                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
+                            &          prodatqc%var(kvar)%vdep(ista:iend),     & 
+                            &          zclim(iin,ijn,:,iobs),            & 
+                            &          zclm2k,interp_corner_clim(iin,ijn,:), & 
+                            &          zgdept(iin,ijn,:,iobs),         & 
+                            &          zmask(iin,ijn,:,iobs) )   
+                     ENDIF   
          
                   ENDDO 
@@ -451,4 +504,9 @@
                   &              prodatqc%var(kvar)%vmod(iend:iend) ) 
 
+               IF ( ld_clim ) THEN
+                  CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), & 
+                     &              prodatqc%var(kvar)%vclm(iend:iend) )
+               ENDIF
+
                   ! Set QC flag for any observations found below the bottom
                   ! needed as the check here is more strict than that in obs_prep
@@ -458,5 +516,6 @@
  
             DEALLOCATE(interp_corner,iv_indic) 
-          
+            IF ( ld_clim ) DEALLOCATE( interp_corner_clim )         
+             
          ENDIF
 
@@ -475,4 +534,6 @@
          & )
 
+      IF ( ld_clim ) DEALLOCATE( zclim )
+      
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
@@ -487,5 +548,5 @@
 
    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
-      &                     kit000, kdaystp, psurf, psurfmask,   &
+      &                     kit000, kdaystp, psurf, pclim, psurfmask,   &
       &                     k2dint, ldnightav, plamscl, pphiscl, &
       &                     lindegrees )
@@ -541,4 +602,5 @@
       REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
          & psurf,  &                   ! Model surface field
+         & pclim,  &                   ! Climatological surface field         
          & psurfmask                   ! Land-sea mask
       LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
@@ -569,5 +631,5 @@
       REAL(wp) :: zlam
       REAL(wp) :: zphi
-      REAL(wp), DIMENSION(1) :: zext, zobsmask
+      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm
       REAL(wp) :: zdaystp
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
@@ -577,4 +639,5 @@
          & zsurfm, &
          & zsurftmp, &
+         & zclim,  &
          & zglam,  &
          & zgphi,  &
@@ -586,5 +649,6 @@
          & zouttmp, &
          & zmeanday    ! to compute model sst in region of 24h daylight (pole)
-
+         
+      LOGICAL :: ld_clim ! T => climatological data is available
       !------------------------------------------------------------------------
       ! Local initialization 
@@ -593,4 +657,11 @@
       inrc = kt - kit000 + 2
       isurf = surfdataqc%nsstp(inrc)
+      
+      ! Check if climatological information is available
+      IF ( SUM(pclim(:,:)) == 0._wp ) THEN
+        ld_clim = .FALSE.
+      ELSE
+        ld_clim = .TRUE.      
+      ENDIF
 
       ! Work out the maximum footprint size for the 
@@ -679,4 +750,6 @@
          & )
 
+      IF ( ld_clim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
+
       DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
          iobs = jobs - surfdataqc%nsurfup
@@ -720,4 +793,9 @@
          &                  igrdip1, igrdjp1, gphif, zgphif )
 
+      IF ( ld_clim ) THEN 
+         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
+            &                  igrdi, igrdj, pclim, zclim )
+      ENDIF
+
       ! At the end of the day get interpolated means
       IF ( idayend == 0 .AND. ldnightav ) THEN
@@ -775,4 +853,9 @@
             CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
 
+            IF ( ld_clim ) THEN  
+               CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
+            ENDIF
+
+
          ELSE
 
@@ -788,4 +871,9 @@
                &              zweig, zsurftmp(:,:,iobs),  zext )
 
+            IF ( ld_clim ) THEN  
+               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
+                  &              zweig, zsurftmp(:,:,iobs),  zclm )
+            ENDIF
+
          ENDIF
 
@@ -797,4 +885,6 @@
             surfdataqc%rmod(jobs,1) = zext(1)
          ENDIF
+         
+         IF ( ld_clim ) surfdataqc%rclm(jobs,1) = zclm(1)
          
          IF ( zext(1) == obfillflt ) THEN
@@ -821,4 +911,6 @@
          & )
 
+      IF ( ld_clim ) DEALLOCATE( zclim )
+
       ! At the end of the day also deallocate night-time mean array
       IF ( idayend == 0 .AND. ldnightav ) THEN
