Index: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 11468)
@@ -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/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 11468)
@@ -33,4 +33,5 @@
    USE mpp_map                  ! MPP mapping
    USE lib_mpp                  ! For ctl_warn/stop
+   USE tradmp                   ! For climatological temperature & salinity
 
    IMPLICIT NONE
@@ -53,4 +54,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
@@ -232,4 +234,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(:,:,:) :: &
@@ -254,5 +257,6 @@
          &            ln_grid_global, ln_grid_search_lookup,          &
          &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
-         &            ln_sstnight, ln_default_fp_indegs,              &
+         &            ln_sstnight,  ln_output_clim,                   &
+         &            ln_default_fp_indegs,                           &
          &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
          &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
@@ -415,4 +419,5 @@
          WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
          WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
+         WRITE(numout,*) '             Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim
       ENDIF
       !-----------------------------------------------------------------------
@@ -438,4 +443,13 @@
          RETURN
       ENDIF
+
+      IF ( ln_output_clim .AND. ( .NOT. ln_tradmp ) ) THEN
+         IF(lwp) WRITE(numout,cform_war)
+         IF(lwp) WRITE(numout,*) ' ln_output_clim is true, but ln_tradmp is false', &
+            &                    ' so climatological T/S not available and will not be output'
+         nwarn = nwarn + 1
+         ln_output_clim = .FALSE.
+      ENDIF
+     
 
       IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes
@@ -714,8 +728,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 )
@@ -763,5 +780,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 )
 
@@ -799,5 +816,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
@@ -805,5 +827,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 )
@@ -926,9 +949,11 @@
       REAL(wp) :: tiny             ! small number
       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
-         & zprofvar                ! Model values for variables in a prof ob
+         & zprofvar, &             ! Model values for variables in a prof ob
+         & zprofclim               ! Climatology 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.
+         & zsurfclim, &            ! Climatology values for variables in a surface ob.
          & zsurfmask               ! Mask associated with surface variable
       REAL(wp), POINTER, DIMENSION(:,:,:) :: &
@@ -940,5 +965,5 @@
          & pco2_3d                 ! 3D pCO2 from FABM
 #endif
-
+      
       IF(lwp) THEN
          WRITE(numout,*)
@@ -963,5 +988,6 @@
             CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
             CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
-            
+            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim )    
+                              
             ! Defaults which might change
             DO jvar = 1, profdataqc(jtype)%nvar
@@ -969,4 +995,5 @@
                zglam(:,:,jvar)       = glamt(:,:)
                zgphi(:,:,jvar)       = gphit(:,:)
+               zprofclim(:,:,:,jvar) = 0._wp
             END DO
 
@@ -976,5 +1003,9 @@
                zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
                zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
-
+               IF ( ln_output_clim ) THEN          
+                  zprofclim(:,:,:,1) = tclim(:,:,:)
+                  zprofclim(:,:,:,2) = sclim(:,:,:)
+               ENDIF
+               
             CASE('vel')
                zprofvar(:,:,:,1) = un(:,:,:)
@@ -1155,4 +1186,5 @@
                   &               nit000, idaystp, jvar,                   &
                   &               zprofvar(:,:,:,jvar),                    &
+                  &               zprofclim(:,:,:,jvar),                   &
                   &               fsdept(:,:,:), fsdepw(:,:,:),            & 
                   &               zprofmask(:,:,:,jvar),                   &
@@ -1166,4 +1198,5 @@
             CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
             CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
+            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim  )            
 
          END DO
@@ -1175,4 +1208,5 @@
          !Allocate local work arrays
          CALL wrk_alloc( jpi, jpj, zsurfvar )
+         CALL wrk_alloc( jpi, jpj, zsurfclim )         
          CALL wrk_alloc( jpi, jpj, zsurfmask )
 #if defined key_fabm
@@ -1184,4 +1218,5 @@
             !Defaults which might be changed
             zsurfmask(:,:) = tmask(:,:,1)
+            zsurfclim(:,:) = 0._wp          
             llog10 = .FALSE.
 
@@ -1189,8 +1224,10 @@
             CASE('sst')
                zsurfvar(:,:) = tsn(:,:,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(:,:) = sclim(:,:,1)              
             CASE('sic')
                IF ( kstp == 0 ) THEN
@@ -1485,7 +1522,8 @@
 
             CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
-               &               nit000, idaystp, zsurfvar, zsurfmask,    &
+               &               nit000, idaystp, zsurfvar,               &
+               &               zsurfclim, zsurfmask,                    &
                &               n2dintsurf(jtype), llnightav(jtype),     &
-               &               ravglamscl(jtype), ravgphiscl(jtype),     &
+               &               ravglamscl(jtype), ravgphiscl(jtype),    &
                &               lfpindegs(jtype) )
 
Index: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 11468)
@@ -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,   &
@@ -187,5 +190,6 @@
       REAL(KIND=wp), DIMENSION(1) :: zmsk
       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
-
+      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim
+      
       LOGICAL :: ld_dailyav
 
@@ -262,4 +266,6 @@
          & )
 
+      IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) )
+
       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
          iobs = jobs - prodatqc%nprofup
@@ -286,4 +292,8 @@
       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 
 
+      IF ( prodatqc%lclim ) 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 +359,6 @@
                   inum_obs = iend - ista + 1 
                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
-
+                  IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
+                  
                   DO iin=1,2 
                      DO ijn=1,2 
@@ -358,4 +369,12 @@
                               &     zobs2k, zgdept(iin,ijn,:,iobs), & 
                               &     zmask(iin,ijn,:,iobs)) 
+
+                           IF ( prodatqc%lclim ) THEN
+                              CALL obs_int_z1d_spl( kpk, & 
+                                 &     zclim(iin,ijn,:,iobs), & 
+                                 &     zclm2k, zgdept(iin,ijn,:,iobs), & 
+                                 &     zmask(iin,ijn,:,iobs)) 
+                           ENDIF
+
                         ENDIF 
        
@@ -371,5 +390,14 @@
                            &    zgdept(iin,ijn,:,iobs), & 
                            &    zmask(iin,ijn,:,iobs)) 
-       
+
+                        IF ( prodatqc%lclim ) 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 
@@ -386,4 +414,5 @@
                inum_obs = iend - ista + 1 
                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
+               IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )                  
                DO iin=1,2  
                   DO ijn=1,2 
@@ -394,4 +423,11 @@
                            &    zobs2k, zgdept(iin,ijn,:,iobs), & 
                            &    zmask(iin,ijn,:,iobs)) 
+
+                        IF ( prodatqc%lclim ) THEN
+                           CALL obs_int_z1d_spl( kpk, & 
+                              &    zclim(iin,ijn,:,iobs),& 
+                              &    zclm2k, zgdept(iin,ijn,:,iobs), & 
+                              &    zmask(iin,ijn,:,iobs)) 
+                        ENDIF
   
                      ENDIF 
@@ -408,4 +444,13 @@
                          &          zgdept(iin,ijn,:,iobs),         & 
                          &          zmask(iin,ijn,:,iobs) )      
+
+                     IF ( prodatqc%lclim ) 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 +496,9 @@
                   &              prodatqc%var(kvar)%vmod(iend:iend) ) 
 
+               IF ( prodatqc%lclim ) 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 +508,6 @@
  
             DEALLOCATE(interp_corner,iv_indic) 
-          
+            IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim )         
+             
          ENDIF
 
@@ -475,4 +526,6 @@
          & )
 
+      IF ( prodatqc%lclim ) DEALLOCATE( zclim )
+      
       ! At the end of the day also get interpolated means
       IF ( ld_dailyav .AND. idayend == 0 ) THEN
@@ -487,5 +540,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 +594,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 +623,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 +631,5 @@
          & zsurfm, &
          & zsurftmp, &
+         & zclim,  &
          & zglam,  &
          & zgphi,  &
@@ -586,5 +641,5 @@
          & zouttmp, &
          & zmeanday    ! to compute model sst in region of 24h daylight (pole)
-
+         
       !------------------------------------------------------------------------
       ! Local initialization 
@@ -679,4 +734,6 @@
          & )
 
+      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
+
       DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
          iobs = jobs - surfdataqc%nsurfup
@@ -715,8 +772,16 @@
       CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
          &                  igrdi, igrdj, psurf, zsurf )
-      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
-         &                  igrdip1, igrdjp1, glamf, zglamf )
-      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
-         &                  igrdip1, igrdjp1, gphif, zgphif )
+
+      IF ( k2dint > 4 ) THEN         
+         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
+            &                  igrdip1, igrdjp1, glamf, zglamf )
+         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
+            &                  igrdip1, igrdjp1, gphif, zgphif )
+      ENDIF
+      
+      IF ( surfdataqc%lclim ) 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
@@ -775,4 +840,9 @@
             CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
 
+            IF ( surfdataqc%lclim ) THEN  
+               CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
+            ENDIF
+
+
          ELSE
 
@@ -788,4 +858,9 @@
                &              zweig, zsurftmp(:,:,iobs),  zext )
 
+            IF ( surfdataqc%lclim ) THEN  
+               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
+                  &              zweig, zclim(:,:,iobs),  zclm )
+            ENDIF
+
          ENDIF
 
@@ -797,4 +872,6 @@
             surfdataqc%rmod(jobs,1) = zext(1)
          ENDIF
+         
+         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
          
          IF ( zext(1) == obfillflt ) THEN
@@ -821,4 +898,6 @@
          & )
 
+      IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
+
       ! At the end of the day also deallocate night-time mean array
       IF ( idayend == 0 .AND. ldnightav ) THEN
Index: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90	(revision 11468)
@@ -72,6 +72,7 @@
          & vdep,  &       !: Depth coordinate of profile data
          & vobs,  &       !: Profile data
-         & vmod           !: Model counterpart of the profile data vector
-
+         & vmod,  &       !: Model counterpart of the profile data vector
+         & vclm           !: Climatological counterpart of the profile data vector
+         
       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
          & vext           !: Extra variables
@@ -102,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
 
@@ -198,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  ***
@@ -221,4 +224,5 @@
       INTEGER, INTENT(IN) :: kpj
       INTEGER, INTENT(IN) :: kpk
+      LOGICAL, INTENT(IN) :: ldclim
 
       !!* Local variables
@@ -236,4 +240,6 @@
       prof%npj       = kpj
       prof%npk       = kpk
+      
+      prof%lclim     = ldclim
 
       ! Allocate arrays of size number of variables
@@ -503,4 +509,9 @@
             & )
       ENDIF
+      IF (prof%lclim) THEN
+         ALLOCATE( & 
+            & prof%var(kvar)%vclm(kobs) &
+            & )
+      ENDIF
 
    END SUBROUTINE obs_prof_alloc_var
@@ -537,4 +548,9 @@
          DEALLOCATE( & 
             & prof%var(kvar)%vext  &
+            & )
+      ENDIF
+      IF (prof%lclim) THEN
+         DEALLOCATE( & 
+            & prof%var(kvar)%vclm  &
             & )
       ENDIF
@@ -630,5 +646,6 @@
             &                 inprof,    invpro,    &
             &                 prof%nstp, prof%npi,  &
-            &                 prof%npj,  prof%npk )
+            &                 prof%npj,  prof%npk,  &
+            &                 prof%lclim )
       ENDIF
 
@@ -745,5 +762,9 @@
                            &                      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
                      
@@ -870,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/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90	(revision 11468)
@@ -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
@@ -697,4 +698,7 @@
                                  &                inpfiles(jj)%padd(ij,ji,1,jvar)
                            ENDIF
+                           IF ( profdata%lclim ) THEN
+                               profdata%var(jvar)%vclm(ivart(jvar)) = fbrmdi
+                           ENDIF                          
                            ! Count number of profile var1 data as function of type
                            itypvar( profdata%ntyp(iprof) + 1, jvar ) = &
Index: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90	(revision 11468)
@@ -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
@@ -462,5 +463,8 @@
                   IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi
                ENDIF
-
+               
+               ! Initialise climatology if set
+               IF ( surfdata%lclim ) surfdata%rclm(iobs,1) = fbrmdi
+               
                ! STD (obs error standard deviation) read from file and passed through obs operator
                IF ( iadd_std(jj) /= -1 ) THEN
Index: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90	(revision 11468)
@@ -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
 
@@ -84,6 +86,7 @@
       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
          & robs, &        !: Surface observation 
-         & rmod           !: Model counterpart of the surface observation vector
-
+         & rmod, &        !: Model counterpart of the surface observation vector
+         & rclm           !: Climatological counterpart of the surface observation vector
+         
       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
          & rext           !: Extra fields interpolated to observation points
@@ -124,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  ***
@@ -143,4 +146,5 @@
       INTEGER, INTENT(IN) :: kpi     ! Number of 3D grid points
       INTEGER, INTENT(IN) :: kpj
+      LOGICAL, INTENT(IN) :: ldclim  
 
       !!* Local variables
@@ -157,4 +161,5 @@
       surf%npi      = kpi
       surf%npj      = kpj
+      surf%lclim    = ldclim
 
       ! Allocate arrays of size number of variables
@@ -197,7 +202,9 @@
       ALLOCATE( & 
          & surf%robs(ksurf,kvar), &
-         & surf%rmod(ksurf,kvar)  &
+         & surf%rmod(ksurf,kvar) &
          & )   
 
+      IF (surf%lclim) ALLOCATE( surf%rclm(ksurf,kvar) )
+      
       ! Allocate arrays of number of extra fields at observation points
 
@@ -293,4 +300,5 @@
          & )
 
+      IF (surf%lclim) DEALLOCATE( surf%rclm )
       ! Deallocate arrays of number of extra fields at observation points
 
@@ -371,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
 
@@ -418,4 +426,5 @@
                newsurf%robs(insurf,jk)  = surf%robs(ji,jk)
                newsurf%rmod(insurf,jk)  = surf%rmod(ji,jk)
+               IF (newsurf%lclim) newsurf%rclm(insurf,jk) = surf%rclm(ji,jk) 
                
             END DO
@@ -514,4 +523,5 @@
             oldsurf%robs(jj,jk)  = surf%robs(ji,jk)
             oldsurf%rmod(jj,jk)  = surf%rmod(ji,jk)
+            IF (surf%lclim) oldsurf%rclm(jj,jk)  = surf%rclm(ji,jk)            
 
          END DO
Index: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 11468)
@@ -95,7 +95,12 @@
       INTEGER :: je
       INTEGER :: iadd
+      INTEGER :: iadd_clm ! 1 if climatology present
       INTEGER :: iext
       REAL(wp) :: zpres
 
+
+      iadd_clm = 0 
+      IF ( profdata%lclim ) iadd_clm = 1
+      
       IF ( PRESENT( padd ) ) THEN
          iadd = padd%inum
@@ -123,5 +128,5 @@
          clfiletype='profb'
          CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
-            &                 1 + iadd, 1 + iext, .TRUE. )
+            &                 1 + iadd_clm + iadd, 1 + iext, .TRUE. )
          fbdata%cname(1)      = profdata%cvars(1)
          fbdata%cname(2)      = profdata%cvars(2)
@@ -137,4 +142,10 @@
          fbdata%caddunit(1,1) = 'Degrees centigrade'
          fbdata%caddunit(1,2) = 'PSU'
+         IF ( profdata%lclim ) THEN
+            fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature'
+            fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity'
+            fbdata%caddunit(2,1) = 'Degrees centigrade'
+            fbdata%caddunit(2,2) = 'PSU'
+         ENDIF
          fbdata%cgrid(:)      = 'T'
          DO je = 1, iext
@@ -144,8 +155,8 @@
          END DO
          DO ja = 1, iadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
+            fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja)
             DO jvar = 1, 2
-               fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)
-               fbdata%caddunit(1+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
@@ -154,5 +165,6 @@
 
          clfiletype='velfb'
-         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 1, 0, .TRUE. )
+         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
+            &                 1 + iadd_clm + iadd, 0, .TRUE. )
          fbdata%cname(1)      = profdata%cvars(1)
          fbdata%cname(2)      = profdata%cvars(2)
@@ -170,10 +182,16 @@
          fbdata%caddunit(1,1) = 'm/s'
          fbdata%caddunit(1,2) = 'm/s'
+         IF ( profdata%lclim ) THEN
+            fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity'
+            fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity'
+            fbdata%caddunit(2,1) = 'm/s'
+            fbdata%caddunit(2,2) = 'm/s'
+         ENDIF         
          fbdata%cgrid(1)      = 'U' 
          fbdata%cgrid(2)      = 'V'
          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)
+            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
 
@@ -246,5 +264,5 @@
          & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
          CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
-            &                 1 + iadd, iext, .TRUE. )
+            &                 1 + iadd_clm + iadd, iext, .TRUE. )
          fbdata%cname(1)      = profdata%cvars(1)
          fbdata%coblong(1)    = cllongname
@@ -252,4 +270,8 @@
          fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
          fbdata%caddunit(1,1) = clunits
+         IF ( profdata%lclim ) THEN
+            fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname)
+            fbdata%caddunit(2,1) = clunits
+         ENDIF         
          fbdata%cgrid(:)      = clgrid
          DO je = 1, iext
@@ -259,12 +281,13 @@
          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)
+            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 ( profdata%lclim ) fbdata%caddname(1+iadd_clm)   = 'CLM'
+      
       WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
 
@@ -319,5 +342,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)
                fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
                fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
@@ -333,6 +355,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,1+ja,jvar) = &
+                  fbdata%padd(ik,jo,1+iadd_clm+ja,jvar) = &
                      & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
                END DO
@@ -420,5 +447,17 @@
       INTEGER :: indx_std
       INTEGER :: iadd_std
-
+      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, CLM, STD, MDT...
+ 
       IF ( PRESENT( padd ) ) THEN
          iadd = padd%inum
@@ -426,11 +465,5 @@
          iadd = 0
       ENDIF
-
-      IF ( PRESENT( pext ) ) THEN
-         iext = pext%inum
-      ELSE
-         iext = 0
-      ENDIF
-
+      
       iadd_std = 0
       indx_std = -1
@@ -444,14 +477,22 @@
       ENDIF
       
+      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 )
 
       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, &
-            &                 2 + iadd + iadd_std, 1 + iext, .TRUE. )
+            &                 1 + iadd_mdt + iadd_std + iadd, &
+            &                 1 + iext, .TRUE. )
 
          clfiletype = 'slafb'
@@ -474,7 +515,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
 
@@ -485,5 +526,5 @@
          clunits    = 'Degree centigrade'
          clgrid     = 'T'
-
+         
       CASE('ICECONC')
 
@@ -499,5 +540,5 @@
          clunits    = 'psu'
          clgrid     = 'T'
-
+         
       CASE('SLCHLTOT','LOGCHL','LogChl','logchl')
 
@@ -610,5 +651,5 @@
       
          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-            &                 1 + iadd + iadd_std, iext, .TRUE. )
+            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )
 
          fbdata%cname(1)      = surfdata%cvars(1)
@@ -628,7 +669,7 @@
          fbdata%cgrid(1)      = clgrid
          DO ja = 1, iadd
-            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+iadd_std+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
 
@@ -637,9 +678,15 @@
       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
+         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
+      
       
       WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
@@ -689,6 +736,5 @@
             &           fbdata%ptim(jo),   &
             &           krefdate = 19500101 )
-         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
-         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
+                    
          fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
          fbdata%pdep(1,jo)     = 0.0
@@ -706,13 +752,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/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 11461)
+++ branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 11468)
@@ -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   ==!
