Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 8235)
@@ -66,6 +66,6 @@
 
       IF( lk_oasis) THEN
-      ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
-      IF( ln_coupled_iceshelf_fluxes ) THEN
+      ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+      IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
 
         ! Adjust total calving rates so that sum of iceberg calving and iceshelf melting in the northern
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8235)
@@ -1627,4 +1627,6 @@
 #if ! defined key_xios2
       WRITE(cl1,'(i1)')        1   ;   CALL iom_set_field_attr('field_definition', freq_op=cl1//'ts', freq_offset='0ts')
+      WRITE(cl1,'(i1)')        2   ;   CALL iom_set_field_attr('trendT_even'      , freq_op=cl1//'ts', freq_offset='0ts')
+      WRITE(cl1,'(i1)')        2   ;   CALL iom_set_field_attr('trendT_odd'       , freq_op=cl1//'ts', freq_offset='-1ts')
       WRITE(cl1,'(i1)')  nn_fsbc   ;   CALL iom_set_field_attr('SBC'             , freq_op=cl1//'ts', freq_offset='0ts')
       WRITE(cl1,'(i1)')  nn_fsbc   ;   CALL iom_set_field_attr('SBC_scalar'      , freq_op=cl1//'ts', freq_offset='0ts')
@@ -1633,4 +1635,6 @@
 #else
       f_op%timestep = 1        ;  f_of%timestep = 0  ; CALL iom_set_field_attr('field_definition', freq_op=f_op, freq_offset=f_of)
+      f_op%timestep = 2        ;  f_of%timestep = 0  ; CALL iom_set_field_attr('trendT_even'      , freq_op=f_op, freq_offset=f_of)
+      f_op%timestep = 2        ;  f_of%timestep = -1 ; CALL iom_set_field_attr('trendT_odd'       , freq_op=f_op, freq_offset=f_of)
       f_op%timestep = nn_fsbc  ;  f_of%timestep = 0  ; CALL iom_set_field_attr('SBC'             , freq_op=f_op, freq_offset=f_of)
       f_op%timestep = nn_fsbc  ;  f_of%timestep = 0  ; CALL iom_set_field_attr('SBC_scalar'      , freq_op=f_op, freq_offset=f_of)
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 8235)
@@ -158,6 +158,6 @@
 #endif
                      IF( lk_oasis) THEN
-                     ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
-                     IF( ln_coupled_iceshelf_fluxes ) THEN
+                     ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+                     IF( nn_coupled_iceshelf_fluxes .eq. 1 ) THEN
                         CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass', greenland_icesheet_mass )
                         CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed )
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8235)
@@ -153,9 +153,15 @@
    ! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to 
    ! avoid circular dependencies.
-   LOGICAL, PUBLIC     ::   ln_coupled_iceshelf_fluxes     ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the 
+   INTEGER, PUBLIC     ::   nn_coupled_iceshelf_fluxes     ! =0 : total freshwater input from iceberg calving and ice shelf basal melting 
+                                                           ! taken from climatologies used (no action in coupling routines).
+                                                           ! =1 :  use rate of change of mass of Greenland and Antarctic icesheets to set the 
                                                            ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
+                                                           ! =2 :  specify constant freshwater inputs in this namelist to set the combined
+                                                           ! magnitude of iceberg calving and iceshelf melting freshwater fluxes.
    LOGICAL, PUBLIC     ::   ln_iceshelf_init_atmos         ! If true force ocean to initialise iceshelf masses from atmospheric values rather
-                                                           ! than values in ocean restart.
+                                                           ! than values in ocean restart (applicable if nn_coupled_iceshelf_fluxes=1).
+   REAL(wp), PUBLIC    ::   rn_greenland_total_fw_flux    ! Constant total rate of freshwater input (kg/s) for Greenland (if nn_coupled_iceshelf_fluxes=2) 
    REAL(wp), PUBLIC    ::   rn_greenland_calving_fraction  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   REAL(wp), PUBLIC    ::   rn_antarctica_total_fw_flux   ! Constant total rate of freshwater input (kg/s) for Antarctica (if nn_coupled_iceshelf_fluxes=2) 
    REAL(wp), PUBLIC    ::   rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
    REAL(wp), PUBLIC    ::   rn_iceshelf_fluxes_tolerance   ! Absolute tolerance for detecting differences in icesheet masses. 
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8235)
@@ -242,7 +242,7 @@
          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  &
          &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  &
-         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, ln_iceshelf_init_atmos,        &
-         &                  rn_greenland_calving_fraction, &
-         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance
+         &                  ln_usecplmask, nn_coupled_iceshelf_fluxes, ln_iceshelf_init_atmos,        &
+         &                  rn_greenland_total_fw_flux, rn_greenland_calving_fraction, &
+         &                  rn_antarctica_total_fw_flux, rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance
       !!---------------------------------------------------------------------
 
@@ -314,6 +314,8 @@
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
-         WRITE(numout,*)'  ln_coupled_iceshelf_fluxes          = ', ln_coupled_iceshelf_fluxes
+         WRITE(numout,*)'  nn_coupled_iceshelf_fluxes          = ', nn_coupled_iceshelf_fluxes
          WRITE(numout,*)'  ln_iceshelf_init_atmos              = ', ln_iceshelf_init_atmos
+         WRITE(numout,*)'  rn_greenland_total_fw_flux         = ', rn_greenland_total_fw_flux
+         WRITE(numout,*)'  rn_antarctica_total_fw_flux        = ', rn_antarctica_total_fw_flux
          WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction
          WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction
@@ -957,5 +959,5 @@
       ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
-      IF( ln_coupled_iceshelf_fluxes ) THEN
+      IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
           ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something
           ! more complicated could be done if required.
@@ -1345,5 +1347,5 @@
       ENDIF
 
-      IF( srcv(jpr_grnm)%laction ) THEN
+      IF( srcv(jpr_grnm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
          greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
          ! take average over ocean points of input array to avoid cumulative error over time
@@ -1377,8 +1379,10 @@
          IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
          IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         greenland_icesheet_mass_rate_of_change = rn_greenland_total_fw_flux
       ENDIF
 
       !                                                        ! land ice masses : Antarctica
-      IF( srcv(jpr_antm)%laction ) THEN
+      IF( srcv(jpr_antm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
          antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
          ! take average over ocean points of input array to avoid cumulative error from rounding errors over time
@@ -1412,4 +1416,6 @@
          IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
          IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
       ENDIF
 
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 8235)
@@ -255,6 +255,6 @@
 
             IF( lk_oasis) THEN
-            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
-            IF( ln_coupled_iceshelf_fluxes ) THEN
+            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
 
               ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
@@ -304,6 +304,6 @@
 
             IF( lk_oasis) THEN
-            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
-            IF( ln_coupled_iceshelf_fluxes ) THEN
+            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
 
               ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 8235)
@@ -135,16 +135,27 @@
             CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
          ENDIF
-         ! total trend for the non-time-filtered variables. 
-         DO jk = 1, jpkm1
-            zfact = 1.0 / rdttra(jk)
-            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
-            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
-         END DO
+         ! total trend for the non-time-filtered variables.
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpkm1
+               zfact = 1.0 / rdttra(jk)
+               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact
+               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact
+            END DO
+         ELSE
+            DO jk = 1, jpkm1
+               zfact = 1.0 / rdttra(jk)
+               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
+               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
+            END DO
+         END IF
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
-         ! Store now fields before applying the Asselin filter 
-         ! in order to calculate Asselin filter trend later.
-         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
-         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
+         IF( .NOT.lk_vvl )  THEN
+            ! Store now fields before applying the Asselin filter 
+            ! in order to calculate Asselin filter trend later.
+            ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
+            ztrds(:,:,:) = tsn(:,:,:,jp_sal)
+         END IF
       ENDIF
 
@@ -155,4 +166,11 @@
             END DO
          END DO
+         IF (l_trdtra.AND.lk_vvl) THEN      ! Zero Asselin filter contribution must be explicitly written out since for vvl
+                                            ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
+            ztrdt(:,:,:) = 0._wp
+            ztrds(:,:,:) = 0._wp
+            CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
+            CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
+         END IF
       ELSE                                            ! Leap-Frog + Asselin filter time stepping
          !
@@ -164,5 +182,5 @@
       !
      ! trends computation
-      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
+      IF( l_trdtra.AND..NOT.lk_vvl) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
          DO jk = 1, jpkm1
             zfact = 1._wp / r2dtra(jk)             
@@ -172,6 +190,6 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
       END IF
+      IF( l_trdtra) CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
       !
       !                        ! control print
@@ -289,6 +307,7 @@
       LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
-      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
+      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
       REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrd_atf
       !!----------------------------------------------------------------------
       !
@@ -315,6 +334,11 @@
       ENDIF
       !
+      IF( l_trdtra )   THEN
+         CALL wrk_alloc( jpi, jpj, jpk, kjpt, ztrd_atf )
+         ztrd_atf(:,:,:,:) = 0.0_wp
+      ENDIF
       DO jn = 1, kjpt      
          DO jk = 1, jpkm1
+            zfact = 1._wp / r2dtra(jk)
             zfact1 = atfp * p2dt(jk)
             zfact2 = zfact1 / rau0
@@ -371,4 +395,7 @@
                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
                   ENDIF
+                  IF( l_trdtra ) THEN
+                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
+                  ENDIF
                END DO
             END DO
@@ -377,4 +404,10 @@
       END DO
       !
+      IF( l_trdtra ) THEN
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
+         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
+      ENDIF
+
    END SUBROUTINE tra_nxt_vvl
 
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 8235)
@@ -94,8 +94,18 @@
 
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / r2dtra(jk) ) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk)
-         END DO
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn.
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpkm1
+               ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*fse3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*fse3t_b(:,:,jk) ) &
+                    & / (fse3t_n(:,:,jk)*r2dtra(jk)) ) - ztrdt(:,:,jk)
+               ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*fse3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*fse3t_b(:,:,jk) ) &
+                    & / (fse3t_n(:,:,jk)*r2dtra(jk)) ) - ztrds(:,:,jk)
+            END DO
+         ELSE
+            DO jk = 1, jpkm1
+               ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / r2dtra(jk) ) - ztrdt(:,:,jk)
+               ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk)
+            END DO
+         END IF
          CALL lbc_lnk( ztrdt, 'T', 1. )
          CALL lbc_lnk( ztrds, 'T', 1. )
Index: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 8136)
+++ branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 8235)
@@ -301,63 +301,70 @@
       !! ** Purpose :   output 3D tracer trends using IOM
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend 
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
-      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
-      INTEGER                   , INTENT(in   ) ::   kt      ! time step
-      !!
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
-      INTEGER ::   ikbu, ikbv   ! local integers
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace 
-      !!----------------------------------------------------------------------
-      !
-!!gm Rq: mask the trends already masked in trd_tra, but lbc_lnk should probably be added
-      !
-      SELECT CASE( ktrd )
-      CASE( jptra_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection
-                               CALL iom_put( "strd_xad" , ptrdy )
-      CASE( jptra_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection
-                               CALL iom_put( "strd_yad" , ptrdy )
-      CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection
-                               CALL iom_put( "strd_zad" , ptrdy )
-                               IF( .NOT. lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface
-                               	 CALL wrk_alloc( jpi, jpj, z2dx, z2dy )
-                                  z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1)
-                                  z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1)
-                                  CALL iom_put( "ttrd_sad", z2dx )
-                                  CALL iom_put( "strd_sad", z2dy )
-                                  CALL wrk_dealloc( jpi, jpj, z2dx, z2dy )
-                               ENDIF
-      CASE( jptra_totad  ) ;   CALL iom_put( "ttrd_totad" , ptrdx )        ! total   advection
-                               CALL iom_put( "strd_totad" , ptrdy )
-      CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion
-                               CALL iom_put( "strd_ldf" , ptrdy )
-      CASE( jptra_zdf  )   ;   CALL iom_put( "ttrd_zdf" , ptrdx )        ! vertical diffusion (including Kz contribution)
-                               CALL iom_put( "strd_zdf" , ptrdy )
-      CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution)
-                               CALL iom_put( "strd_zdfp", ptrdy )
-      CASE( jptra_evd )    ;   CALL iom_put( "ttrd_evd", ptrdx )         ! EVD trend (convection)
-                               CALL iom_put( "strd_evd", ptrdy )
-      CASE( jptra_iso_x )  ;   CALL iom_put( "ttrd_iso_x", ptrdx )       ! x-component of isopycnal mixing
-                               CALL iom_put( "strd_iso_x", ptrdy )
-      CASE( jptra_iso_y )  ;   CALL iom_put( "ttrd_iso_y", ptrdx )       ! y-component of isopycnal mixing
-                               CALL iom_put( "strd_iso_y", ptrdy )
-      CASE( jptra_iso_z1 ) ;   CALL iom_put( "ttrd_iso_z1", ptrdx )      ! first part of z-component of isopycnal mixing
-                               CALL iom_put( "strd_iso_z1", ptrdy )
-      CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping)
-                               CALL iom_put( "strd_dmp" , ptrdy )
-      CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl" , ptrdx )        ! bottom boundary layer
-                               CALL iom_put( "strd_bbl" , ptrdy )
-      CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing
-                               CALL iom_put( "strd_npc" , ptrdy )
-      CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx(:,:,1) )        ! surface forcing + runoff (ln_rnf=T)
-                               CALL iom_put( "strd_cdt" , ptrdy(:,:,1) )        ! output as 2D surface fields
-      CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature)
-      CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature)
-      CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter
-                               CALL iom_put( "strd_atf" , ptrdy )
-      CASE( jptra_tot  )   ;   CALL iom_put( "ttrd_tot" , ptrdx )        ! model total trend
-                               CALL iom_put( "strd_tot" , ptrdy )
-      END SELECT
-      !
+     REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend 
+     REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
+     INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
+     INTEGER                   , INTENT(in   ) ::   kt      ! time step
+     !!
+     INTEGER ::   ji, jj, jk   ! dummy loop indices
+     INTEGER ::   ikbu, ikbv   ! local integers
+     REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace 
+     !!----------------------------------------------------------------------
+     !
+     !!gm Rq: mask the trends already masked in trd_tra, but lbc_lnk should probably be added
+     !
+     ! Trends evaluated every time step that could go to the standard T file and can be output every ts into a 1ts file if 1ts output is selected
+     SELECT CASE( ktrd )
+     ! This total trend is done every time step
+     CASE( jptra_tot  )   ;   CALL iom_put( "ttrd_tot" , ptrdx )           ! model total trend
+        CALL iom_put( "strd_tot" , ptrdy )
+     END SELECT
+
+     ! These trends are done every second time step. When 1ts output is selected must go different (2ts) file from standard T-file
+     IF( MOD( kt, 2 ) == 0 ) THEN
+        SELECT CASE( ktrd )
+        CASE( jptra_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection
+           CALL iom_put( "strd_xad" , ptrdy )
+        CASE( jptra_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection
+           CALL iom_put( "strd_yad" , ptrdy )
+        CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection
+           CALL iom_put( "strd_zad" , ptrdy )
+           IF( .NOT. lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface
+              CALL wrk_alloc( jpi, jpj, z2dx, z2dy )
+              z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1)
+              z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1)
+              CALL iom_put( "ttrd_sad", z2dx )
+              CALL iom_put( "strd_sad", z2dy )
+              CALL wrk_dealloc( jpi, jpj, z2dx, z2dy )
+           ENDIF
+        CASE( jptra_totad  ) ;   CALL iom_put( "ttrd_totad" , ptrdx )      ! total   advection
+           CALL iom_put( "strd_totad" , ptrdy )
+        CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion
+           CALL iom_put( "strd_ldf" , ptrdy )
+        CASE( jptra_zdf  )   ;   CALL iom_put( "ttrd_zdf" , ptrdx )        ! vertical diffusion (including Kz contribution)
+           CALL iom_put( "strd_zdf" , ptrdy )
+        CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution)
+           CALL iom_put( "strd_zdfp", ptrdy )
+        CASE( jptra_evd )    ;   CALL iom_put( "ttrd_evd", ptrdx )         ! EVD trend (convection)
+           CALL iom_put( "strd_evd", ptrdy )
+        CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping)
+           CALL iom_put( "strd_dmp" , ptrdy )
+        CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl" , ptrdx )        ! bottom boundary layer
+           CALL iom_put( "strd_bbl" , ptrdy )
+        CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing
+           CALL iom_put( "strd_npc" , ptrdy )
+        CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature)
+        CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx(:,:,1) ) ! surface forcing + runoff (ln_rnf=T)
+           CALL iom_put( "strd_cdt" , ptrdy(:,:,1) )        ! output as 2D surface fields
+        CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature)
+        END SELECT
+        ! the Asselin filter trend  is also every other time step but needs to be lagged one time step
+        ! Even when 1ts output is selected can go to the same (2ts) file as the trends plotted every even time step.
+     ELSE IF( MOD( kt, 2 ) == 1 ) THEN
+        SELECT CASE( ktrd )
+        CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter
+           CALL iom_put( "strd_atf" , ptrdy )
+        END SELECT
+     END IF
+     !
    END SUBROUTINE trd_tra_iom
 
