Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 11455)
@@ -1232,4 +1232,5 @@
    ln_s_at_t  = .false.             ! Logical switch for computing model S at T obs if not there
    ln_sstnight = .false.            ! Logical switch for calculating night-time average for SST obs
+   ln_output_clim = .false.         ! Logical switch for writing climatological values to fdbk files
    ln_default_fp_indegs = .true.    ! Logical: T=> averaging footprint is in degrees, F=> in metres
    ln_sla_fp_indegs = .true.        ! Logical: T=> averaging footprint is in degrees, F=> in metres
Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 11455)
@@ -50,5 +50,4 @@
    LOGICAL :: ln_diaobs            !: Logical switch for the obs operator
    LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs
-   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and outputting T/S climatology
    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
@@ -56,4 +55,5 @@
    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
+   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and writing T/S climatology
 
    REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint
@@ -235,4 +235,5 @@
       LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
       LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables)
+      LOGICAL :: ltype_clim      ! Local version of ln_output_clim
 
       REAL(wp), POINTER, DIMENSION(:,:,:) :: &
@@ -728,8 +729,11 @@
 
          DO jtype = 1, nproftypes
-
+            
+            ltype_clim = .FALSE. 
+            
             IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
                nvarsprof(jtype) = 2
                nextrprof(jtype) = 1
+               IF ( ln_output_clim ) ltype_clim = .TRUE.              
                ALLOCATE(llvar(nvarsprof(jtype)))
                CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
@@ -777,5 +781,5 @@
                &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
                &               rn_dobsini, rn_dobsend, llvar, &
-               &               ln_ignmis, ln_s_at_t, .FALSE., &
+               &               ln_ignmis, ln_s_at_t, .FALSE., ltype_clim, &
                &               kdailyavtypes = nn_profdavtypes )
 
@@ -813,5 +817,10 @@
             nvarssurf(jtype) = 1
             nextrsurf(jtype) = 0
+            ltype_clim = .FALSE.
             IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
+            IF ( ln_output_clim .AND. &
+               & ( ( TRIM(cobstypessurf(jtype)) == 'sst' ) .OR. &
+               &   ( TRIM(cobstypessurf(jtype)) == 'sss' ) ) ) &
+               & ltype_clim = .TRUE.
 
             !Read in surface obs types
@@ -819,5 +828,6 @@
                &               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), ltype_clim )
 
             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
@@ -928,4 +938,7 @@
          & jp_spm
 #endif
+      USE tradmp, ONLY: &
+         & tclim, &
+         & sclim
 
       IMPLICIT NONE
@@ -956,5 +969,4 @@
          & pco2_3d                 ! 3D pCO2 from FABM
 #endif
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta 
       
       IF(lwp) THEN
@@ -966,7 +978,4 @@
 
       idaystp = NINT( rday / rdt )
-
-      ! Get the climatological T & S fields on this time step
-      IF ( ln_output_clim ) CALL dta_tsd( kstp, zts_dta )
 
       !-----------------------------------------------------------------------
@@ -998,7 +1007,7 @@
                zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
                zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
-               IF ( ln_output_clim ) THEN               
-                  zprofclim(:,:,:,1) = zts_dta(:,:,:,jp_tem)
-                  zprofclim(:,:,:,2) = zts_dta(:,:,:,jp_sal)
+               IF ( ln_output_clim ) THEN          
+                  zprofclim(:,:,:,1) = tclim(:,:,:)
+                  zprofclim(:,:,:,2) = sclim(:,:,:)
                ENDIF
                
@@ -1219,10 +1228,10 @@
             CASE('sst')
                zsurfvar(:,:) = tsn(:,:,1,jp_tem)
-               IF ( ln_output_clim ) zsurfclim(:,:) = zts_dta(:,:,1,jp_tem)
+               IF ( ln_output_clim ) zsurfclim(:,:) = tclim(:,:,1)
             CASE('sla')
                zsurfvar(:,:)  = sshn(:,:)
             CASE('sss')
                zsurfvar(:,:) = tsn(:,:,1,jp_sal)
-               IF ( ln_output_clim ) zsurfclim(:,:) = zts_dta(:,:,1,jp_sal)               
+               IF ( ln_output_clim ) zsurfclim(:,:) = sclim(:,:,1)              
             CASE('sic')
                IF ( kstp == 0 ) THEN
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 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 11455)
@@ -190,7 +190,7 @@
       REAL(KIND=wp), DIMENSION(1) :: zmsk
       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
-
+      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim
+      
       LOGICAL :: ld_dailyav
-      LOGICAL :: ld_clim
 
       !------------------------------------------------------------------------
@@ -200,11 +200,4 @@
       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
@@ -273,5 +266,5 @@
          & )
 
-      IF ( ld_clim ) ALLOCATE( zclim(2,2,kpk,ipro) )
+      IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) )
 
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
@@ -299,5 +292,5 @@
       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 
 
-      IF ( ld_clim ) THEN
+      IF ( prodatqc%lclim ) THEN
          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim )            
       ENDIF 
@@ -366,5 +359,5 @@
                   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) )
+                  IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
                   
                   DO iin=1,2 
@@ -377,5 +370,5 @@
                               &     zmask(iin,ijn,:,iobs)) 
 
-                           IF ( ld_clim ) THEN
+                           IF ( prodatqc%lclim ) THEN
                               CALL obs_int_z1d_spl( kpk, & 
                                  &     zclim(iin,ijn,:,iobs), & 
@@ -398,5 +391,5 @@
                            &    zmask(iin,ijn,:,iobs)) 
 
-                        IF ( ld_clim ) THEN
+                        IF ( prodatqc%lclim ) THEN
                            CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
                               &    prodatqc%var(kvar)%vdep(ista:iend), & 
@@ -422,5 +415,5 @@
                ALLOCATE( interp_corner(2,2,inum_obs),      &
                   &      iv_indic(inum_obs) ) 
-               IF ( ld_clim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )                  
+               IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )                  
                DO iin=1,2  
                   DO ijn=1,2 
@@ -432,5 +425,5 @@
                            &    zmask(iin,ijn,:,iobs)) 
 
-                        IF ( ld_clim ) THEN
+                        IF ( prodatqc%lclim ) THEN
                            CALL obs_int_z1d_spl( kpk, & 
                               &    zclim(iin,ijn,:,iobs),& 
@@ -453,5 +446,5 @@
                          &          zmask(iin,ijn,:,iobs) )      
 
-                     IF ( ld_clim ) THEN
+                     IF ( prodatqc%lclim ) THEN
                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
                             &          prodatqc%var(kvar)%vdep(ista:iend),     & 
@@ -504,5 +497,5 @@
                   &              prodatqc%var(kvar)%vmod(iend:iend) ) 
 
-               IF ( ld_clim ) THEN
+               IF ( prodatqc%lclim ) THEN
                   CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), & 
                      &              prodatqc%var(kvar)%vclm(iend:iend) )
@@ -516,5 +509,5 @@
  
             DEALLOCATE(interp_corner,iv_indic) 
-            IF ( ld_clim ) DEALLOCATE( interp_corner_clim )         
+            IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim )         
              
          ENDIF
@@ -534,5 +527,5 @@
          & )
 
-      IF ( ld_clim ) DEALLOCATE( zclim )
+      IF ( prodatqc%lclim ) DEALLOCATE( zclim )
       
       ! At the end of the day also get interpolated means
@@ -650,5 +643,4 @@
          & zmeanday    ! to compute model sst in region of 24h daylight (pole)
          
-      LOGICAL :: ld_clim ! T => climatological data is available
       !------------------------------------------------------------------------
       ! Local initialization 
@@ -657,11 +649,4 @@
       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 
@@ -750,5 +735,5 @@
          & )
 
-      IF ( ld_clim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
+      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
 
       DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
@@ -793,5 +778,5 @@
          &                  igrdip1, igrdjp1, gphif, zgphif )
 
-      IF ( ld_clim ) THEN 
+      IF ( surfdataqc%lclim ) THEN 
          CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
             &                  igrdi, igrdj, pclim, zclim )
@@ -853,6 +838,9 @@
             CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
 
-            IF ( ld_clim ) THEN  
+            IF ( surfdataqc%lclim ) THEN  
                CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
+               IF (lwp) THEN
+                  WRITE(numout,*)'zclim: ', iobs, zclim(:,:,iobs), zclm
+               ENDIF
             ENDIF
 
@@ -871,5 +859,5 @@
                &              zweig, zsurftmp(:,:,iobs),  zext )
 
-            IF ( ld_clim ) THEN  
+            IF ( surfdataqc%lclim ) THEN  
                CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
                   &              zweig, zsurftmp(:,:,iobs),  zclm )
@@ -886,5 +874,5 @@
          ENDIF
          
-         IF ( ld_clim ) surfdataqc%rclm(jobs,1) = zclm(1)
+         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
          
          IF ( zext(1) == obfillflt ) THEN
@@ -911,5 +899,5 @@
          & )
 
-      IF ( ld_clim ) DEALLOCATE( zclim )
+      IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
 
       ! At the end of the day also deallocate night-time mean array
Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90	(revision 11455)
@@ -73,5 +73,5 @@
          & vobs,  &       !: Profile data
          & vmod,  &       !: Model counterpart of the profile data vector
-         & vclm,  &       !: Climatological counterpart of the profile data vector
+         & vclm           !: Climatological counterpart of the profile data vector
          
       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
@@ -103,4 +103,6 @@
       INTEGER :: nprofup  !: Observation counter used in obs_oper
 
+      LOGICAL :: lclim    !: Climatology will be calculated for this structure
+      
       ! Bookkeeping arrays with sizes equal to number of variables
 
@@ -199,5 +201,5 @@
    
    SUBROUTINE obs_prof_alloc( prof,  kvar, kext, kprof,  &
-      &                       ko3dt, kstp, kpi, kpj, kpk )
+      &                       ko3dt, kstp, kpi, kpj, kpk, ldclim )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE obs_prof_alloc  ***
@@ -222,4 +224,5 @@
       INTEGER, INTENT(IN) :: kpj
       INTEGER, INTENT(IN) :: kpk
+      LOGICAL, INTENT(IN) :: ldclim
 
       !!* Local variables
@@ -237,4 +240,6 @@
       prof%npj       = kpj
       prof%npk       = kpk
+      
+      prof%lclim     = ldclim
 
       ! Allocate arrays of size number of variables
@@ -493,5 +498,4 @@
          & prof%var(kvar)%vobs(kobs),      &
          & prof%var(kvar)%vmod(kobs),      &
-         & prof%var(kvar)%vclm(kobs),      &         
          & prof%var(kvar)%nvind(kobs)      &
          & )
@@ -503,4 +507,9 @@
          ALLOCATE( & 
             & prof%var(kvar)%vext(kobs,kext) &
+            & )
+      ENDIF
+      IF (prof%lclim) THEN
+         ALLOCATE( & 
+            & prof%var(kvar)%vclm(kobs) &
             & )
       ENDIF
@@ -532,5 +541,4 @@
          & prof%var(kvar)%vobs,   &
          & prof%var(kvar)%vmod,   &
-         & prof%var(kvar)%vclm,   &         
          & prof%var(kvar)%nvind,  &
          & prof%var(kvar)%idqcf,  &
@@ -540,4 +548,9 @@
          DEALLOCATE( & 
             & prof%var(kvar)%vext  &
+            & )
+      ENDIF
+      IF (prof%lclim) THEN
+         DEALLOCATE( & 
+            & prof%var(kvar)%vclm  &
             & )
       ENDIF
@@ -633,5 +646,6 @@
             &                 inprof,    invpro,    &
             &                 prof%nstp, prof%npi,  &
-            &                 prof%npj,  prof%npk )
+            &                 prof%npj,  prof%npk,  &
+            &                 prof%lclim )
       ENDIF
 
@@ -744,11 +758,13 @@
                      newprof%var(jvar)%vmod(invpro(jvar))   = &
                         &                           prof%var(jvar)%vmod(jj)
-                     newprof%var(jvar)%vclm(invpro(jvar))   = &
-                        &                           prof%var(jvar)%vclm(jj)
                      DO jext = 1, prof%next
                         newprof%var(jvar)%vext(invpro(jvar),jext) = &
                            &                      prof%var(jvar)%vext(jj,jext)
                      END DO
-                  
+                     IF (newprof%lclim) THEN
+                        newprof%var(jvar)%vclm(invpro(jvar))   = &
+                           &                           prof%var(jvar)%vclm(jj)
+                     ENDIF
+                                    
                      ! nvind is the index of the original variable data
                      
@@ -869,5 +885,4 @@
                oldprof%var(jvar)%vobs(jl)   = prof%var(jvar)%vobs(jj)
                oldprof%var(jvar)%vmod(jl)   = prof%var(jvar)%vmod(jj)
-               oldprof%var(jvar)%vclm(jl)   = prof%var(jvar)%vclm(jj)               
                oldprof%var(jvar)%idqcf(:,jl) = prof%var(jvar)%idqcf(:,jj)
                oldprof%var(jvar)%nvqcf(:,jl) = prof%var(jvar)%nvqcf(:,jj)
@@ -876,5 +891,7 @@
                      &                        prof%var(jvar)%vext(jj,jext)
                END DO
-               
+               IF (prof%lclim) THEN
+                  oldprof%var(jvar)%vclm(jl)   = prof%var(jvar)%vclm(jj)
+               ENDIF              
             END DO
 
Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 11455)
@@ -46,5 +46,5 @@
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
       &                     ldvar, ldignmis, ldsatt, &
-      &                     ldmod, kdailyavtypes )
+      &                     ldmod, ldclim, kdailyavtypes )
       !!---------------------------------------------------------------------
       !!
@@ -78,4 +78,5 @@
       LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points
       LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
+      LOGICAL, INTENT(IN) :: ldclim     ! Set flag to show climatology will be output
       REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
       REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
@@ -500,5 +501,5 @@
       ENDIF
       CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
-         &                 kstp, jpi, jpj, jpk )
+         &                 kstp, jpi, jpj, jpk, ldclim )
 
       ! * Read obs/positions, QC, all variable and assign to profdata
Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 11455)
@@ -40,5 +40,5 @@
    SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                     ldignmis, ldmod, ldnightav )
+      &                     ldignmis, ldmod, ldnightav, ldclim )
       !!---------------------------------------------------------------------
       !!
@@ -71,4 +71,5 @@
       LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
       LOGICAL, INTENT(IN) :: ldnightav  ! Observations represent a night-time average
+      LOGICAL, INTENT(IN) :: ldclim     ! Will include climatology at obs points.
       REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
       REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
@@ -359,5 +360,5 @@
          &               iindx   )
 
-      CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj )
+      CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj, ldclim )
 
       ! Read obs/positions, QC, all variable and assign to surfdata
Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 11455)
@@ -52,4 +52,6 @@
       INTEGER :: nrec       !: Number of surface observation records in window
 
+      LOGICAL :: lclim      !: Climatology will be calculated for this structure
+      
       ! Arrays with size equal to the number of surface observations
 
@@ -125,5 +127,5 @@
 CONTAINS
    
-   SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj )
+   SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj, ldclim )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE obs_surf_alloc  ***
@@ -144,4 +146,5 @@
       INTEGER, INTENT(IN) :: kpi     ! Number of 3D grid points
       INTEGER, INTENT(IN) :: kpj
+      LOGICAL, INTENT(IN) :: ldclim  
 
       !!* Local variables
@@ -158,4 +161,5 @@
       surf%npi      = kpi
       surf%npj      = kpj
+      surf%lclim    = ldclim
 
       ! Allocate arrays of size number of variables
@@ -198,8 +202,9 @@
       ALLOCATE( & 
          & surf%robs(ksurf,kvar), &
-         & surf%rmod(ksurf,kvar), &
-         & surf%rclm(ksurf,kvar)  &         
+         & surf%rmod(ksurf,kvar) &
          & )   
 
+      IF (surf%lclim) ALLOCATE( surf%rclm(ksurf,kvar) )
+      
       ! Allocate arrays of number of extra fields at observation points
 
@@ -292,8 +297,8 @@
       DEALLOCATE( & 
          & surf%robs,    &
-         & surf%rmod,    &
-         & surf%rclm     &         
-         & )
-
+         & surf%rmod     &
+         & )
+
+      IF (surf%lclim) DEALLOCATE( surf%rclm )
       ! Deallocate arrays of number of extra fields at observation points
 
@@ -374,5 +379,5 @@
       IF ( lallocate ) THEN
          CALL obs_surf_alloc( newsurf,  insurf, surf%nvar, &
-            & surf%nextra, surf%nstp, surf%npi, surf%npj )
+            & surf%nextra, surf%nstp, surf%npi, surf%npj, surf%lclim )
       ENDIF
 
@@ -421,5 +426,5 @@
                newsurf%robs(insurf,jk)  = surf%robs(ji,jk)
                newsurf%rmod(insurf,jk)  = surf%rmod(ji,jk)
-               newsurf%rclm(insurf,jk)  = surf%rclm(ji,jk)               
+               IF (newsurf%lclim) newsurf%rclm(insurf,jk) = surf%rclm(ji,jk) 
                
             END DO
@@ -518,5 +523,5 @@
             oldsurf%robs(jj,jk)  = surf%robs(ji,jk)
             oldsurf%rmod(jj,jk)  = surf%rmod(ji,jk)
-            oldsurf%rclm(jj,jk)  = surf%rclm(ji,jk)            
+            IF (surf%lclim) oldsurf%rclm(jj,jk)  = surf%rclm(ji,jk)            
 
          END DO
Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 11455)
@@ -29,6 +29,4 @@
    USE obs_mpp              ! MPP support routines for observation diagnostics
    USE lib_mpp		    ! MPP routines
-   USE diaobs, ONLY: &
-      & ln_output_clim
 
    IMPLICIT NONE
@@ -97,5 +95,5 @@
       INTEGER :: je
       INTEGER :: iadd
-      INTEGER :: iadd_exp  ! expected additional variables
+      INTEGER :: iadd_clm ! 1 if climatology present
       INTEGER :: iext
       REAL(wp) :: zpres
@@ -104,6 +102,6 @@
       ! Set up number of additional variables to be ouput:
       ! Hx, CLIM, ...
-      iadd_exp = 1   ! Hx
-      IF ( ln_output_clim ) iadd_exp = iadd_exp + 1
+      iadd_clm = 0 
+      IF ( profdata%lclim ) iadd_clm = 1
       
       IF ( PRESENT( padd ) ) THEN
@@ -132,5 +130,5 @@
          clfiletype='profb'
          CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
-            &                 iadd_exp + iadd, 1 + iext, .TRUE. )
+            &                 1 + iadd_clm + iadd, 1 + iext, .TRUE. )
          fbdata%cname(1)      = profdata%cvars(1)
          fbdata%cname(2)      = profdata%cvars(2)
@@ -148,5 +146,5 @@
          fbdata%caddunit(1,1) = 'Degrees centigrade'
          fbdata%caddunit(1,2) = 'PSU'
-         IF ( ln_output_clim ) THEN
+         IF ( profdata%lclim ) THEN
             fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature'
             fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity'
@@ -161,8 +159,8 @@
          END DO
          DO ja = 1, iadd
-            fbdata%caddname(iadd_exp+ja) = padd%cdname(ja)
+            fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja)
             DO jvar = 1, 2
-               fbdata%caddlong(iadd_exp+ja,jvar) = padd%cdlong(ja,jvar)
-               fbdata%caddunit(iadd_exp+ja,jvar) = padd%cdunit(ja,jvar)
+               fbdata%caddlong(1+iadd_clm+ja,jvar) = padd%cdlong(ja,jvar)
+               fbdata%caddunit(1+iadd_clm+ja,jvar) = padd%cdunit(ja,jvar)
             END DO
          END DO
@@ -172,5 +170,5 @@
          clfiletype='velfb'
          CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
-            &                 iadd_exp + iadd, 0, .TRUE. )
+            &                 1 + iadd_clm + iadd, 0, .TRUE. )
          fbdata%cname(1)      = profdata%cvars(1)
          fbdata%cname(2)      = profdata%cvars(2)
@@ -188,5 +186,5 @@
          fbdata%caddunit(1,1) = 'm/s'
          fbdata%caddunit(1,2) = 'm/s'
-         IF ( ln_output_clim ) THEN
+         IF ( profdata%lclim ) THEN
             fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity'
             fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity'
@@ -197,7 +195,7 @@
          fbdata%cgrid(2)      = 'V'
          DO ja = 1, iadd
-            fbdata%caddname(iadd_exp+ja) = padd%cdname(ja)
-            fbdata%caddlong(iadd_exp+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(iadd_exp+ja,1) = padd%cdunit(ja,1)
+            fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja)
+            fbdata%caddlong(1+iadd_clm+ja,1) = padd%cdlong(ja,1)
+            fbdata%caddunit(1+iadd_clm+ja,1) = padd%cdunit(ja,1)
          END DO
 
@@ -270,5 +268,5 @@
          & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
          CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
-            &                 iadd_expt + iadd, iext, .TRUE. )
+            &                 1 + iadd_clm + iadd, iext, .TRUE. )
          fbdata%cname(1)      = profdata%cvars(1)
          fbdata%coblong(1)    = cllongname
@@ -276,5 +274,5 @@
          fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
          fbdata%caddunit(1,1) = clunits
-         IF ( ln_output_clim ) THEN
+         IF ( profdata%lclim ) THEN
             fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname)
             fbdata%caddunit(2,1) = clunits
@@ -287,12 +285,12 @@
          END DO
          DO ja = 1, iadd
-            fbdata%caddname(iadd_expt+ja) = padd%cdname(ja)
-            fbdata%caddlong(iadd_expt+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(iadd_expt+ja,1) = padd%cdunit(ja,1)
+            fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja)
+            fbdata%caddlong(1+iadd_clm+ja,1) = padd%cdlong(ja,1)
+            fbdata%caddunit(1+iadd_clm+ja,1) = padd%cdunit(ja,1)
          END DO
       ENDIF
 
       fbdata%caddname(1)   = 'Hx'
-      IF ( ln_output_clim ) fbdata%caddname(2)   = 'CLM'
+      IF ( profdata%lclim ) fbdata%caddname(1+iadd_clm)   = 'CLM'
       
       WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
@@ -348,8 +346,4 @@
             DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
                ik = profdata%var(jvar)%nvlidx(jk)
-               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
-               IF ( ln_output_clim ) THEN           
-                  fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vclm(jk)     
-               ENDIF              
                fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
                fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
@@ -365,6 +359,11 @@
                ENDIF
                fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
+               
+               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
+               IF ( profdata%lclim ) THEN           
+                  fbdata%padd(ik,jo,1+iadd_clm,jvar) = profdata%var(jvar)%vclm(jk)     
+               ENDIF              
                DO ja = 1, iadd
-                  fbdata%padd(ik,jo,iadd_exp+ja,jvar) = &
+                  fbdata%padd(ik,jo,1+iadd_clm+ja,jvar) = &
                      & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
                END DO
@@ -449,15 +448,19 @@
       INTEGER :: je
       INTEGER :: iadd
-      INTEGER :: iadd_exp
       INTEGER :: iext
       INTEGER :: indx_std
       INTEGER :: iadd_std
-      INTEGER :: iadd_clm      
+      INTEGER :: iadd_clm     
+      INTEGER :: iadd_mdt 
+
+      IF ( PRESENT( pext ) ) THEN
+         iext = pext%inum
+      ELSE
+         iext = 0
+      ENDIF
 
 
       ! Set up number of additional variables to be ouput:
-      ! Hx, CLIM, ...
-      iadd_exp = 1   ! Hx
-      IF ( ln_output_clim ) iadd_exp = iadd_exp + 1
+      ! Hx, CLM, STD, MDT...
  
       IF ( PRESENT( padd ) ) THEN
@@ -466,11 +469,5 @@
          iadd = 0
       ENDIF
-
-      IF ( PRESENT( pext ) ) THEN
-         iext = pext%inum
-      ELSE
-         iext = 0
-      ENDIF
-
+      
       iadd_std = 0
       indx_std = -1
@@ -485,4 +482,8 @@
       
       iadd_clm = 0
+      IF ( surfdata%lclim ) iadd_clm = 1
+      
+      iadd_mdt = 0
+      IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_mdt = 1
       
       CALL init_obfbdata( fbdata )
@@ -496,5 +497,6 @@
 
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 2 + iadd + iadd_std, 1 + iext, .TRUE. )
+            &                 1 + iadd_mdt + iadd_std + iadd, &
+            &                 1 + iext, .TRUE. )
 
          clfiletype = 'slafb'
@@ -517,7 +519,7 @@
          fbdata%cgrid(1)      = 'T'
          DO ja = 1, iadd
-            fbdata%caddname(2+iadd_std+ja) = padd%cdname(ja)
-            fbdata%caddlong(2+iadd_std+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(2+iadd_std+ja,1) = padd%cdunit(ja,1)
+            fbdata%caddname(1+iadd_mdt+iadd_std+ja) = padd%cdname(ja)
+            fbdata%caddlong(1+iadd_mdt+iadd_std+ja,1) = padd%cdlong(ja,1)
+            fbdata%caddunit(1+iadd_mdt+iadd_std+ja,1) = padd%cdunit(ja,1)
          END DO
 
@@ -528,5 +530,4 @@
          clunits    = 'Degree centigrade'
          clgrid     = 'T'
-         IF ( ln_output_clim ) iadd_clm = 1
          
       CASE('ICECONC')
@@ -543,5 +544,4 @@
          clunits    = 'psu'
          clgrid     = 'T'
-         IF ( ln_output_clim ) iadd_clm = 1
          
       CASE('SLCHLTOT','LOGCHL','LogChl','logchl')
@@ -655,5 +655,5 @@
       
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd + iadd_std + iadd_clm, iext, .TRUE. )
+            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )
 
          fbdata%cname(1)      = surfdata%cvars(1)
@@ -673,7 +673,7 @@
          fbdata%cgrid(1)      = clgrid
          DO ja = 1, iadd
-            fbdata%caddname(1+iadd_std+iadd_clm+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1)
+            fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja)
+            fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1)
+            fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1)
          END DO
 
@@ -682,14 +682,13 @@
       fbdata%caddname(1)   = 'Hx'
       IF ( indx_std /= -1 ) THEN
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_std = iadd_std + 1
-         fbdata%caddname(1+iadd_std)   = surfdata%cext(indx_std)
-         fbdata%caddlong(1+iadd_std,1) = 'Obs error standard deviation'
-         fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1)
-      ENDIF
-      
-      IF ( ln_output_clim .AND. ( iadd_clm > 0 ) ) THEN
-         fbdata%caddname(1+iadd_std+iadd_clm)   = 'CLM'
-         fbdata%caddlong(1+iadd_std+iadd_clm,1) = 'Climatology'
-         fbdata%caddunit(1+iadd_std+iadd_clm,1) = fbdata%cobunit(1)
+         fbdata%caddname(1+iadd_mdt+iadd_std)   = surfdata%cext(indx_std)
+         fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation'
+         fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1)
+      ENDIF
+      
+      IF ( surfdata%lclim ) THEN
+         fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm)   = 'CLM'
+         fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,1) = 'Climatology'
+         fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,1) = fbdata%cobunit(1)
       ENDIF
       
@@ -741,12 +740,5 @@
             &           fbdata%ptim(jo),   &
             &           krefdate = 19500101 )
-         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
-            fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
-         ENDIF    
-         IF ( ln_output_clim .AND. ( iadd_clm > 0 ) ) THEN
-            fbdata%padd(1,jo,2,1) = surfdata%rclm(jo,1)
-         ENDIF
-                     
+                    
          fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
          fbdata%pdep(1,jo)     = 0.0
@@ -764,13 +756,28 @@
          ENDIF
          fbdata%iobsk(1,jo,1)  = 0
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
+ 
+         ! Additional variables.
+         ! Hx is always the first additional variable
+         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
+         ! MDT is output as an additional variable if SLA obs type
+         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
+            fbdata%padd(1,jo,1+iadd_mdt,1) = surfdata%rext(jo,1)
+         ENDIF    
+         ! STD is output as an additional variable if available
          IF ( indx_std /= -1 ) THEN
-            fbdata%padd(1,jo,1+iadd_std,1) = surfdata%rext(jo,indx_std)
+            fbdata%padd(1,jo,1+iadd_mdt+iadd_std,1) = surfdata%rext(jo,indx_std)
          ENDIF
+         ! CLM is output as an additional variable if available
+         IF ( surfdata%lclim ) THEN
+            fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,1) = surfdata%rclm(jo,1)
+         ENDIF
+         ! Then other additional variables are output
+         DO ja = 1, iadd
+            fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,1) = &
+               & surfdata%rext(jo,padd%ipoint(ja))
+         END DO
          
-         DO ja = 1, iadd
-            fbdata%padd(1,jo,2+iadd_std+ja,1) = &
-               & surfdata%rext(jo,padd%ipoint(ja))
-         END DO
+         ! Extra variables
+         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)         
          DO je = 1, iext
             fbdata%pext(1,jo,1+je) = &
Index: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 11449)
+++ branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 11455)
@@ -55,5 +55,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ttrdmp   !: damping temperature trend (Celcius/s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1)
-
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tclim    !: temperature climatology on each time step(Celcius)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sclim    !: salinity climatology on each time step (psu)
+   
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -70,5 +72,6 @@
       !!                ***  FUNCTION tra_dmp_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), STAT= tra_dmp_alloc )
+      ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), &
+         &      tclim(jpi,jpj,jpk) , sclim(jpi,jpj,jpk), STAT= tra_dmp_alloc )
       !
       IF( lk_mpp            )   CALL mpp_sum ( tra_dmp_alloc )
@@ -110,4 +113,7 @@
       !                           !==   input T-S data at kt   ==!
       CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt
+      
+      tclim(:,:,:) = zts_dta(:,:,:,jp_tem)
+      sclim(:,:,:) = zts_dta(:,:,:,jp_sal)
       !
       SELECT CASE ( nn_zdmp )     !==    type of damping   ==!
