Index: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 5967)
+++ branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 5972)
@@ -6,8 +6,6 @@
    !! History :  OPA  ! 1991-03  (O. Marti, G. Madec)  Original code
    !!                 ! 1992-06  (M. Imbard)  doctor norme
-   !!                 ! 1996-01  (G. Madec)  statement function for e3
-   !!                 ! 1997-05  (G. Madec)  macro-tasked on jk-slab
    !!                 ! 1998-07  (M. Imbard, G. Madec) ORCA version
-   !!            7.0  ! 2001-02  (M. Imbard)  cofdis, Original code
+   !!            7.0  ! 2001-02  (M. Imbard)  add distance to coast, Original code
    !!            8.1  ! 2001-02  (G. Madec, E. Durand)  cleaning
    !!  NEMO      1.0  ! 2002-08  (G. Madec, E. Durand)  free form + modules
@@ -15,4 +13,6 @@
    !!            3.3  ! 2010-06  (C. Ethe, G. Madec) merge TRA-TRC 
    !!            3.4  ! 2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys
+   !!            3.6  ! 2015-06  (T. Graham)  read restoring coefficient in a file
+   !!            3.7  ! 2015-10  (G. Madec)  remove useless trends arrays
    !!----------------------------------------------------------------------
 
@@ -42,16 +42,11 @@
 
    PUBLIC   tra_dmp      ! routine called by step.F90
-   PUBLIC   tra_dmp_init ! routine called by opa.F90
-
-   !                               !!* Namelist namtra_dmp : T & S newtonian damping *
-   ! nn_zdmp and cn_resto are public as they are used by C1D/dyndmp.F90
-   LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag
-   INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer
-   CHARACTER(LEN=200) , PUBLIC :: cn_resto      ! name of netcdf file containing restoration coefficient field
+   PUBLIC   tra_dmp_init ! routine called by nemogcm.F90
+
+   !                                           !!* Namelist namtra_dmp : T & S newtonian damping *
+   LOGICAL            , PUBLIC ::   ln_tradmp   !: internal damping flag
+   INTEGER            , PUBLIC ::   nn_zdmp     !: = 0/1/2 flag for damping in the mixed layer
+   CHARACTER(LEN=200) , PUBLIC ::   cn_resto    !: name of netcdf file containing restoration coefficient field
    !
-
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   strdmp   !: damping salinity trend (psu/s)
-   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)
 
@@ -70,5 +65,5 @@
       !!                ***  FUNCTION tra_dmp_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), STAT= tra_dmp_alloc )
+      ALLOCATE( resto(jpi,jpj,jpk), STAT= tra_dmp_alloc )
       !
       IF( lk_mpp            )   CALL mpp_sum ( tra_dmp_alloc )
@@ -96,69 +91,58 @@
       !! ** Action  : - (ta,sa)   tracer trends updated with the damping trend
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zta, zsa             ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta 
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start( 'tra_dmp')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, jpts,  zts_dta )
-      !
-      !                           !==   input T-S data at kt   ==!
+      !
+      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta, ztrdts
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('tra_dmp')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,jpts,   zts_dta )
+      !
+      IF( l_trdtra )   THEN                    !* Save ta and sa trends
+         CALL wrk_alloc( jpi,jpj,jpk,jpts,   ztrdts ) 
+         ztrdts(:,:,:,:) = tsa(:,:,:,:) 
+      ENDIF
+      !                           !==  input T-S data at kt  ==!
       CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt
       !
-      SELECT CASE ( nn_zdmp )     !==    type of damping   ==!
-      !
-      CASE( 0 )                   !==  newtonian damping throughout the water column  ==!
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zta = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                  zsa = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
-                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
-                  strdmp(ji,jj,jk) = zsa           ! save the trend (used in asmtrj)
-                  ttrdmp(ji,jj,jk) = zta      
+      SELECT CASE ( nn_zdmp )     !==  type of damping  ==!
+      !
+      CASE( 0 )                        !*  newtonian damping throughout the water column  *!
+         DO jn = 1, jpts
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) )
+                  END DO
                END DO
             END DO
          END DO
          !
-      CASE ( 1 )                  !==  no damping in the turbocline (avt > 5 cm2/s)  ==!
+      CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *!
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( avt(ji,jj,jk) <= 5.e-4_wp ) THEN
-                     zta = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     zsa = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
-                  ELSE
-                     zta = 0._wp
-                     zsa = 0._wp  
+                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
+                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
                   ENDIF
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
-                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
-                  strdmp(ji,jj,jk) = zsa           ! save the salinity trend (used in asmtrj)
-                  ttrdmp(ji,jj,jk) = zta
                END DO
             END DO
          END DO
          !
-      CASE ( 2 )                  !==  no damping in the mixed layer   ==!
+      CASE ( 2 )                       !*  no damping in the mixed layer   *!
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
-                     zta = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     zsa = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
-                  ELSE
-                     zta = 0._wp
-                     zsa = 0._wp  
+                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
+                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
+                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
                   ENDIF
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
-                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
-                  strdmp(ji,jj,jk) = zsa           ! save the salinity trend (used in asmtrj)
-                  ttrdmp(ji,jj,jk) = zta
                END DO
             END DO
@@ -168,6 +152,8 @@
       !
       IF( l_trdtra )   THEN       ! trend diagnostic
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ttrdmp )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, strdmp )
+         ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:)
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
+         CALL wrk_dealloc( jpi,jpj,jpk,jpts,   ztrdts ) 
       ENDIF
       !                           ! Control print
@@ -175,7 +161,7 @@
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, jpts,  zts_dta )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop( 'tra_dmp')
+      CALL wrk_dealloc( jpi,jpj,jpk,jpts,   zts_dta )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('tra_dmp')
       !
    END SUBROUTINE tra_dmp
@@ -190,8 +176,7 @@
       !! ** Method  :   read the namtra_dmp namelist and check the parameters
       !!----------------------------------------------------------------------
+      INTEGER ::   ios, imask   ! local integers 
+      !!
       NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto
-      INTEGER ::  ios         ! Local integer for output status of namelist read
-      INTEGER :: imask        ! File handle 
-      !!
       !!----------------------------------------------------------------------
       !
@@ -204,9 +189,9 @@
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namtra_dmp )
-
-      IF(lwp) THEN                 !Namelist print
+      !
+      IF(lwp) THEN                  ! Namelist print
          WRITE(numout,*)
          WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation'
-         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_dmp : set relaxation parameters'
          WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
@@ -215,34 +200,31 @@
          WRITE(numout,*)
       ENDIF
-
+      !
       IF( ln_tradmp) THEN
-         !
-         !Allocate arrays
+         !                          ! Allocate arrays
          IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' )
-
-         !Check values of nn_zdmp
-         SELECT CASE (nn_zdmp)
-         CASE ( 0 )  ; IF(lwp) WRITE(numout,*) '   tracer damping as specified by mask'
-         CASE ( 1 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline'
-         CASE ( 2 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer'
+         !
+         SELECT CASE (nn_zdmp)      ! Check values of nn_zdmp
+         CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping as specified by mask'
+         CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixing layer (kz > 5 cm2/s)'
+         CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed  layer'
+         CASE DEFAULT
+            CALL ctl_stop('tra_dmp_init : wrong value of nn_zdmp')
          END SELECT
-
-         !TG: Initialisation of dtatsd - Would it be better to have dmpdta routine
-         !so can damp to something other than intitial conditions files?
+         !
+         !!TG: Initialisation of dtatsd - Would it be better to have dmpdta routine
+         !    so can damp to something other than intitial conditions files?
+         !!gm: In principle yes. Nevertheless, we can't anticipate demands that have never been formulated.
          IF( .NOT.ln_tsd_tradmp ) THEN
-            CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' )
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout, *)  '   read T-S data not initialized, we force ln_tsd_tradmp=T'
             CALL dta_tsd_init( ld_tradmp=ln_tradmp )        ! forces the initialisation of T-S data
          ENDIF
-
-         !initialise arrays - Are these actually used anywhere else?
-         strdmp(:,:,:) = 0._wp
-         ttrdmp(:,:,:) = 0._wp
-
-         !Read in mask from file
+         !                          ! Read in mask from file
          CALL iom_open ( cn_resto, imask)
-         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
+         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto )
          CALL iom_close( imask )
-       ENDIF
-
+      ENDIF
+      !
    END SUBROUTINE tra_dmp_init
 
