Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7083)
@@ -246,6 +246,4 @@
                ztest_1 = 1
             ELSE 
-              ! this write is useful
-              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 
                ztest_1 = 0
             ENDIF
@@ -258,7 +256,4 @@
                ztest_2 = 1
             ELSE
-              ! this write is useful
-              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
-                            ' zvt_i_ini = ', zvt_i_ini(i_hemis)
                ztest_2 = 0
             ENDIF
@@ -268,7 +263,4 @@
                ztest_3 = 1
             ELSE
-               ! this write is useful
-               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', &
-               zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
                ztest_3 = 0
             ENDIF
@@ -278,6 +270,4 @@
             DO jl = 1, jpl
                IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
-                  ! this write is useful
-                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
                   ztest_4 = 0
                ENDIF
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90	(revision 7083)
@@ -392,5 +392,4 @@
       INTEGER :: ji,jj,jn
       REAL(wp) :: zalpha
-      REAL(wp), DIMENSION(jpi,jpj,7) :: tabice_agr 
       !!-----------------------------------------------------------------------      
       !
@@ -529,4 +528,10 @@
             END DO
          END DO
+      ELSE
+         DO jj=MAX(j1,2),j2
+            DO ji=MAX(i1,2),i2
+               uice_agr(ji,jj) = tabres(ji,jj)
+            END DO
+         END DO
       ENDIF
 #else
@@ -541,4 +546,10 @@
             END DO
          END DO
+      ELSE
+         DO jj= j1, j2
+            DO ji= i1, i2
+               uice_agr(ji,jj) = tabres(ji,jj)
+            END DO
+         END DO
       ENDIF
 #endif
@@ -566,4 +577,10 @@
                   tabres(ji,jj) = e1f(ji-1,jj-1) * v_ice(ji,jj)
                ENDIF
+            END DO
+         END DO
+      ELSE
+         DO jj=MAX(j1,2),j2
+            DO ji=MAX(i1,2),i2
+               vice_agr(ji,jj) = tabres(ji,jj)
             END DO
          END DO
@@ -580,4 +597,10 @@
 	         END DO
 	      END DO
+      ELSE
+	      DO jj= j1 ,j2
+	         DO ji = i1, i2
+	            vice_agr(ji,jj) = tabres(ji,jj)
+	         END DO
+	      END DO
       ENDIF
 #endif
@@ -585,5 +608,5 @@
 
 
-   SUBROUTINE interp_adv_ice( tabres, i1, i2, j1, j2, before )
+   SUBROUTINE interp_adv_ice( tabres, i1, i2, j1, j2, k1, k2, before )
       !!-----------------------------------------------------------------------
       !!                    *** ROUTINE interp_adv_ice ***                           
@@ -593,6 +616,6 @@
       !!              put -9999 where no ice for correct extrapolation             
       !!-----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,7), INTENT(inout) :: tabres
+      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
       LOGICAL, INTENT(in) :: before
       !!
@@ -601,19 +624,27 @@
       !
       IF( before ) THEN
-	      DO jj=j1,j2
-	         DO ji=i1,i2
-	            IF( tms(ji,jj) == 0. ) THEN
-	               tabres(ji,jj,:) = -9999. 
-	            ELSE
-	               tabres(ji,jj, 1) = frld  (ji,jj)
-	               tabres(ji,jj, 2) = hicif (ji,jj)
-	               tabres(ji,jj, 3) = hsnif (ji,jj)
-	               tabres(ji,jj, 4) = tbif  (ji,jj,1)
-	               tabres(ji,jj, 5) = tbif  (ji,jj,2)
-	               tabres(ji,jj, 6) = tbif  (ji,jj,3)
-	               tabres(ji,jj, 7) = qstoif(ji,jj)
-	            ENDIF
-	         END DO
-	      END DO
+         DO jj=j1,j2
+	    DO ji=i1,i2
+	       IF( tms(ji,jj) == 0. ) THEN
+	          tabres(ji,jj,:) = -9999 
+	       ELSE
+	          tabres(ji,jj, 1) = frld  (ji,jj)
+	          tabres(ji,jj, 2) = hicif (ji,jj)
+	          tabres(ji,jj, 3) = hsnif (ji,jj)
+	          tabres(ji,jj, 4) = tbif  (ji,jj,1)
+	          tabres(ji,jj, 5) = tbif  (ji,jj,2)
+	          tabres(ji,jj, 6) = tbif  (ji,jj,3)
+	          tabres(ji,jj, 7) = qstoif(ji,jj)
+	       ENDIF
+	    END DO
+         END DO
+      ELSE
+	 DO jj=j1,j2
+	    DO ji=i1,i2
+               DO jk=k1, k2
+	          tabice_agr(ji,jj,jk) = tabres(ji,jj,jk)
+               END DO
+	    END DO
+	 END DO
       ENDIF
       !
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 7083)
@@ -212,21 +212,4 @@
       REAL(wp) ::   zztmp  
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
-      ! reading initial file
-      LOGICAL  ::   ln_tsd_init      !: T & S data flag
-      LOGICAL  ::   ln_tsd_tradmp    !: internal damping toward input data flag
-      CHARACTER(len=100)            ::   cn_dir
-      TYPE(FLD_N)                   ::  sn_tem,sn_sal
-      INTEGER  ::   ios=0
-
-      NAMELIST/namtsd/ ln_tsd_init,ln_tsd_tradmp,cn_dir,sn_tem,sn_sal
-      !
-
-      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :
-      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in reference namelist for dia_ar5', lwp )
-      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
-      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in configuration namelist for dia_ar5', lwp )
-      IF(lwm) WRITE ( numond, namtsd )
       !
       !!----------------------------------------------------------------------
@@ -234,5 +217,5 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
       !
-      CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_alloc( jpi, jpj, jpk, 2, zsaldta )
       !                                      ! allocate dia_ar5 arrays
       IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
@@ -250,8 +233,9 @@
       IF( lk_mpp )   CALL mpp_sum( vol0 )
 
-      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_sal%clname), inum )
-      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,1), 1  )
-      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 )
+      CALL iom_open ( 'sali_ref_clim_monthly', inum )
+      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
+      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
       CALL iom_close( inum )
+
       sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )        
       sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
@@ -268,5 +252,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_dealloc( jpi, jpj, jpk, 2, zsaldta )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 7083)
@@ -323,9 +323,9 @@
             ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
             va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                      / ( ze3va * rau0 ) 
+               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)
 #else
             va(ji,jj,1) = vb(ji,jj,1) &
                &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
+               &                                      / ( fse3v(ji,jj,1) * rau0     ) * vmask(ji,jj,1) )
 #endif
          END DO
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 7083)
@@ -120,5 +120,5 @@
       ! first entry with narea for this processor is left hand interior index
       ! last  entry                               is right hand interior index
-      jj = jpj/2
+      jj = nlcj/2
       nicbdi = -1
       nicbei = -1
@@ -136,5 +136,5 @@
       !
       ! repeat for j direction
-      ji = jpi/2
+      ji = nlci/2
       nicbdj = -1
       nicbej = -1
@@ -153,9 +153,9 @@
       ! special for east-west boundary exchange we save the destination index
       i1 = MAX( nicbdi-1, 1)
-      i3 = INT( src_calving(i1,jpj/2) )
+      i3 = INT( src_calving(i1,nlcj/2) )
       jj = INT( i3/nicbpack )
       ricb_left = REAL( i3 - nicbpack*jj, wp )
       i1 = MIN( nicbei+1, jpi )
-      i3 = INT( src_calving(i1,jpj/2) )
+      i3 = INT( src_calving(i1,nlcj/2) )
       jj = INT( i3/nicbpack )
       ricb_right = REAL( i3 - nicbpack*jj, wp )
@@ -196,5 +196,5 @@
          WRITE(numicb,*) 'berg left       ', ricb_left
          WRITE(numicb,*) 'berg right      ', ricb_right
-         jj = jpj/2
+         jj = nlcj/2
          WRITE(numicb,*) "central j line:"
          WRITE(numicb,*) "i processor"
@@ -202,5 +202,5 @@
          WRITE(numicb,*) "i point"
          WRITE(numicb,*) (INT(src_calving(ji,jj)), ji=1,jpi)
-         ji = jpi/2
+         ji = nlci/2
          WRITE(numicb,*) "central i line:"
          WRITE(numicb,*) "j processor"
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 7083)
@@ -157,4 +157,14 @@
          END DO
       ENDIF
+
+      ! ORCA R1: Take the minimum between aeiw  and aeiv0
+      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 )
+            END DO
+         END DO
+      ENDIF
+
       CALL lbc_lnk( aeiw, 'W', 1. )      ! lateral boundary condition on aeiw 
 
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7083)
@@ -206,5 +206,7 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
          qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
-         qsr_ice(:,:,1)   = sf(jp_qsr)%fnow(:,:,1)
+         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
+         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
+         ENDIF
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)         
          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7083)
@@ -1618,39 +1618,48 @@
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
       !!
-      !! ** Purpose :   provide the heat and freshwater fluxes of the 
-      !!              ocean-ice system.
+      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system
       !!
       !! ** Method  :   transform the fields received from the atmosphere into
       !!             surface heat and fresh water boundary condition for the 
       !!             ice-ocean system. The following fields are provided:
-      !!              * total non solar, solar and freshwater fluxes (qns_tot, 
+      !!               * total non solar, solar and freshwater fluxes (qns_tot, 
       !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
       !!             NB: emp_tot include runoffs and calving.
-      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
+      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
       !!             emp_ice = sublimation - solid precipitation as liquid
       !!             precipitation are re-routed directly to the ocean and 
-      !!             runoffs and calving directly enter the ocean.
-      !!              * solid precipitation (sprecip), used to add to qns_tot 
+      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90)
+      !!               * solid precipitation (sprecip), used to add to qns_tot 
       !!             the heat lost associated to melting solid precipitation
       !!             over the ocean fraction.
-      !!       ===>> CAUTION here this changes the net heat flux received from
-      !!             the atmosphere
-      !!
-      !!                  - the fluxes have been separated from the stress as
-      !!                 (a) they are updated at each ice time step compare to
-      !!                 an update at each coupled time step for the stress, and
-      !!                 (b) the conservative computation of the fluxes over the
-      !!                 sea-ice area requires the knowledge of the ice fraction
-      !!                 after the ice advection and before the ice thermodynamics,
-      !!                 so that the stress is updated before the ice dynamics
-      !!                 while the fluxes are updated after it.
+      !!               * heat content of rain, snow and evap can also be provided,
+      !!             otherwise heat flux associated with these mass flux are
+      !!             guessed (qemp_oce, qemp_ice)
+      !!
+      !!             - the fluxes have been separated from the stress as
+      !!               (a) they are updated at each ice time step compare to
+      !!               an update at each coupled time step for the stress, and
+      !!               (b) the conservative computation of the fluxes over the
+      !!               sea-ice area requires the knowledge of the ice fraction
+      !!               after the ice advection and before the ice thermodynamics,
+      !!               so that the stress is updated before the ice dynamics
+      !!               while the fluxes are updated after it.
+      !!
+      !! ** Details
+      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided
+      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
+      !!
+      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided
+      !!
+      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce)
+      !!                                                                      river runoff (rnf) is provided but not included here
       !!
       !! ** Action  :   update at each nf_ice time step:
       !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
       !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
-      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
-      !!                   emp_ice            ice sublimation - solid precipitation over the ice
-      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
-      !!                   sprecip             solid precipitation over the ocean  
+      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving)
+      !!                   emp_ice           ice sublimation - solid precipitation over the ice
+      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice
+      !!                   sprecip           solid precipitation over the ocean  
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
@@ -1662,5 +1671,5 @@
       INTEGER ::   jl         ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
@@ -1670,5 +1679,5 @@
       !
       CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
-      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
@@ -1679,11 +1688,11 @@
       !
       !                                                      ! ========================= !
-      !                                                      !    freshwater budget      !   (emp)
+      !                                                      !    freshwater budget      !   (emp_tot)
       !                                                      ! ========================= !
       !
-      !                                                           ! total Precipitation - total Evaporation (emp_tot)
-      !                                                           ! solid precipitation - sublimation       (emp_ice)
-      !                                                           ! solid Precipitation                     (sprecip)
-      !                                                           ! liquid + solid Precipitation            (tprecip)
+      !                                                           ! solid Precipitation                                (sprecip)
+      !                                                           ! liquid + solid Precipitation                       (tprecip)
+      !                                                           ! total Evaporation - total Precipitation            (emp_tot)
+      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice)
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
       CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
@@ -1717,5 +1726,5 @@
          ENDIF
 #else         
-         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
+         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:)
 #endif                  
          CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1) * tmask(:,:,1)      )   ! liquid precipitation 
@@ -1733,5 +1742,5 @@
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
-         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
+         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:)
          zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
@@ -1739,9 +1748,15 @@
 
 #if defined key_lim3
-      ! zsnw = snow percentage over ice after wind blowing
-      zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( p_frld, zsnw )
+      ! zsnw = snow fraction over ice after wind blowing
+      zsnw(:,:) = 0._wp  ;  CALL lim_thd_snwblow( p_frld, zsnw )
       
-      ! --- evaporation (kg/m2/s) --- !
+      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
+      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( zicefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
+      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
+
+      ! --- evaporation over ocean (used later for qemp) --- !
+      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
+
+      ! --- evaporation over ice (kg/m2/s) --- !
       zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
       ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
@@ -1749,14 +1764,11 @@
       zdevap_ice(:,:) = 0._wp
       
-      ! --- evaporation minus precipitation corrected for the effect of wind blowing on snow --- !
-      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) - zsprecip * (1._wp - zsnw)
-      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip * (1._wp - zsnw)          
-
-      ! Sublimation over sea-ice (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:) )
-      ! runoffs and calving (put in emp_tot)
+      ! --- runoffs (included in emp later on) --- !
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+
+      ! --- calving (put in emp_tot and emp_oce) --- !
       IF( srcv(jpr_cal)%laction ) THEN 
          zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
+         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
          CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
       ENDIF
@@ -1784,10 +1796,9 @@
       ENDIF
 
-                                     CALL iom_put( 'snowpre'    , sprecip                         )  ! Snow
-      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw ) )  ! Snow over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw   )  ! Snow over sea-ice         (cell average)    
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:)         )  ! Sublimation over sea-ice (cell average)
+                                     CALL iom_put( 'snowpre'    , sprecip(:,:)                         )  ! Snow
+      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw(:,:) ) )  ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw(:,:)   )  ! Snow over sea-ice         (cell average)
 #else
-      ! Sublimation over sea-ice (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )
       ! runoffs and calving (put in emp_tot)
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
@@ -1821,18 +1832,17 @@
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
       !                                                      ! ========================= !
-      CASE( 'oce only' )                                     ! the required field is directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
-      CASE( 'conservative' )                                      ! the required fields are directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
+      CASE( 'oce only' )         ! the required field is directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'conservative' )     ! the required fields are directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
-            ! Set all category values equal for the moment
             DO jl=1,jpl
-               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
             ENDDO
          ENDIF
-      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
-         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
+         zqns_tot(:,:) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             DO jl=1,jpl
@@ -1841,5 +1851,5 @@
             ENDDO
          ELSE
-            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+            qns_tot(:,:) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl=1,jpl
                zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
@@ -1847,10 +1857,10 @@
             ENDDO
          ENDIF
-      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
+      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
 ! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
          zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
             &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
-            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
+            &                                           + pist(:,:,1) * zicefr(:,:) ) )
       END SELECT
 !!gm
@@ -1862,18 +1872,11 @@
 !! similar job should be done for snow and precipitation temperature
       !                                     
-      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting 
-         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting 
-         zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:)
-         IF( iom_use('hflx_cal_cea') )   &
-            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
-      ENDIF
-
-      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus
-      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
+      IF( srcv(jpr_cal)%laction ) THEN   ! Iceberg melting 
+         zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! add the latent heat of iceberg melting
+                                                                         ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
+         IF( iom_use('hflx_cal_cea') )   CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus )   ! heat flux from calving
+      ENDIF
 
 #if defined key_lim3      
-      ! --- evaporation --- !
-      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
-
       ! --- non solar flux over ocean --- !
       !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
@@ -1882,18 +1885,18 @@
 
       ! --- heat flux associated with emp (W/m2) --- !
-      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
-         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
-         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
+      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn(:,:)   &       ! evap
+         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &       ! liquid precip
+         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus )  ! solid precip over ocean + snow melting
 !      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
 !         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
       zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only)
-                                                                                                       ! qevap_ice=0 since we consider Tice=0°C
+                                                                                                       ! qevap_ice=0 since we consider Tice=0degC
       
-      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
+      ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
       zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
       DO jl = 1, jpl
-         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0°C
+         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0degC
       END DO
 
@@ -1921,4 +1924,12 @@
          qemp_ice (:,:  ) = zqemp_ice (:,:  )
       ENDIF
+
+      !! clem: we should output qemp_oce and qemp_ice (at least)
+      IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average)
+      !! these diags are not outputed yet
+!!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average)
+!!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average)
+!!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average)
+
 #else
       ! clem: this formulation is certainly wrong... but better than it was...
@@ -1927,5 +1938,5 @@
          &          - (p_frld(:,:) * zsprecip(:,:) * lfus)  &          ! remove the latent heat flux of solid precip. melting
          &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
-         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
+         &             - zemp_ice(:,:) ) * zcptn(:,:) 
 
      IF( ln_mixcpl ) THEN
@@ -2047,5 +2058,5 @@
 
       CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
-      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
       CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7083)
@@ -1018,5 +1018,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zs= SQRT( ABS( psal(ji,jj) ) * r1_S0 )           ! square root salinity
+               zs= SQRT( ABS( psal(ji,jj) ) / 35.16504_wp )           ! square root salinity
                ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
                   &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
@@ -1066,5 +1066,5 @@
       CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==!
          !
-         zs  = SQRT( ABS( psal ) * r1_S0 )           ! square root salinity
+         zs  = SQRT( ABS( psal ) / 35.16504_wp )           ! square root salinity
          ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
                   &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7083)
@@ -26,4 +26,6 @@
    USE cla             ! cross land advection      (cla_traadv     routine)
    USE ldftra_oce      ! lateral diffusion coefficient on tracers
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
    !
    USE in_out_manager  ! I/O manager
@@ -79,4 +81,5 @@
       INTEGER ::   jk   ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -120,5 +123,10 @@
       IF( ln_diaptr )   CALL dia_ptr( zvn )                                     ! diagnose the effective MSF 
       !
-   
+      IF( l_trdtra )   THEN                    !* Save ta and sa trends
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+      ENDIF
+      !
       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
       CASE ( 1 )   ;    CALL tra_adv_cen2   ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered
@@ -151,4 +159,13 @@
       END SELECT
       !
+      IF( l_trdtra )   THEN                      ! save the advective trends for further diagnostics
+         DO jk = 1, jpkm1
+            ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk)
+         END DO
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_totad, ztrds )
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
+      ENDIF
       !                                              ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv  - Ta: ', mask1=tmask,               &
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 7083)
@@ -79,5 +79,5 @@
 # endif  
       REAL(wp), POINTER, DIMENSION(:,:) :: zu_eiv, zv_eiv, zw_eiv, z2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, z3d_T
       !!----------------------------------------------------------------------
       !
@@ -86,5 +86,5 @@
 # if defined key_diaeiv 
       CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
-      IF( ln_diaptr ) CALL wrk_alloc( jpi, jpj, jpk, z3d )
+      CALL wrk_alloc( jpi, jpj, jpk, z3d, z3d_T )
 # else
       CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
@@ -170,31 +170,60 @@
            CALL iom_put( "weiv_masstr" , z3d )  
          ENDIF
-         IF( iom_use("ueiv_masstr") .OR. iom_use("ueiv_heattr") ) THEN
+         IF( iom_use("ueiv_masstr") .OR. iom_use("ueiv_heattr") .OR. iom_use('ueiv_heattr3d')        &
+                                    .OR. iom_use("ueiv_salttr") .OR. iom_use('ueiv_salttr3d') ) THEN
             z3d(:,:,jpk) = 0.e0
-            z2d(:,:) = 0.e0
             DO jk = 1, jpkm1
                z3d(:,:,jk) = rau0 * u_eiv(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
-               z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
             END DO
             CALL iom_put( "ueiv_masstr", z3d )                  ! mass transport in i-direction
          ENDIF
 
-         IF( iom_use('ueiv_heattr') ) THEN
+         IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
             zztmp = 0.5 * rcp 
             z2d(:,:) = 0.e0 
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( z2d, 'U', -1. )
-            CALL iom_put( "ueiv_heattr", zztmp * z2d )                  ! heat transport in i-direction
-         ENDIF
-
-         IF( iom_use("veiv_masstr") .OR. iom_use("veiv_heattr") ) THEN
+            z3d_T(:,:,:) = 0.e0 
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
+                  END DO
+               END DO
+            END DO
+            IF (iom_use('ueiv_heattr') ) THEN
+               CALL lbc_lnk( z2d, 'U', -1. )
+               CALL iom_put( "ueiv_heattr", zztmp * z2d )                  ! 2D heat transport in i-direction
+            ENDIF
+            IF (iom_use('ueiv_heattr3d') ) THEN
+               CALL lbc_lnk( z3d_T, 'U', -1. )
+               CALL iom_put( "ueiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in i-direction
+            ENDIF
+         ENDIF
+
+         IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d') ) THEN
+            zztmp = 0.5 * 0.001
+            z2d(:,:) = 0.e0 
+            z3d_T(:,:,:) = 0.e0 
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
+                  END DO
+               END DO
+            END DO
+            IF (iom_use('ueiv_salttr') ) THEN
+               CALL lbc_lnk( z2d, 'U', -1. )
+               CALL iom_put( "ueiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
+            ENDIF
+            IF (iom_use('ueiv_salttr3d') ) THEN
+               CALL lbc_lnk( z3d_T, 'U', -1. )
+               CALL iom_put( "ueiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
+            ENDIF
+         ENDIF
+
+         IF( iom_use("veiv_masstr") .OR. iom_use("veiv_heattr") .OR. iom_use('veiv_heattr3d')       &
+                                    .OR. iom_use("veiv_salttr") .OR. iom_use('veiv_salttr3d') ) THEN
             z3d(:,:,jpk) = 0.e0
-            z2d(:,:) = 0.e0
             DO jk = 1, jpkm1
                z3d(:,:,jk) = rau0 * v_eiv(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
@@ -202,18 +231,82 @@
             CALL iom_put( "veiv_masstr", z3d )                  ! mass transport in j-direction
          ENDIF
-            
-         IF( iom_use('veiv_heattr') ) THEN
+         IF( iom_use('veiv_heattr') .OR. iom_use('veiv_heattr3d') ) THEN
             zztmp = 0.5 * rcp 
             z2d(:,:) = 0.e0 
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( z2d, 'V', -1. )
-            CALL iom_put( "veiv_heattr", zztmp * z2d )                  !  heat transport in j-direction
-         ENDIF
+            z3d_T(:,:,:) = 0.e0 
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
+                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
+                  END DO
+               END DO
+            END DO
+            IF (iom_use('veiv_heattr') ) THEN
+               CALL lbc_lnk( z2d, 'V', -1. )
+               CALL iom_put( "veiv_heattr", zztmp * z2d )                  ! 2D heat transport in j-direction
+            ENDIF
+            IF (iom_use('veiv_heattr3d') ) THEN
+               CALL lbc_lnk( z3d_T, 'V', -1. )
+               CALL iom_put( "veiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in j-direction
+            ENDIF
+         ENDIF
+
+         IF( iom_use('veiv_salttr') .OR. iom_use('veiv_salttr3d') ) THEN
+            zztmp = 0.5 * 0.001
+            z2d(:,:) = 0.e0 
+            z3d_T(:,:,:) = 0.e0 
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
+                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            IF (iom_use('veiv_salttr') ) THEN
+               CALL lbc_lnk( z2d, 'V', -1. )
+               CALL iom_put( "veiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
+            ENDIF
+            IF (iom_use('veiv_salttr3d') ) THEN
+               CALL lbc_lnk( z3d_T, 'V', -1. )
+               CALL iom_put( "veiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
+            ENDIF
+         ENDIF
+
+         IF( iom_use('weiv_masstr') .OR. iom_use('weiv_heattr3d') .OR. iom_use('weiv_salttr3d')) THEN   ! vertical mass transport & its square value
+           z2d(:,:) = rau0 * e12t(:,:)
+           DO jk = 1, jpk
+              z3d(:,:,jk) = w_eiv(:,:,jk) * z2d(:,:)
+           END DO
+           CALL iom_put( "weiv_masstr" , z3d )                  ! mass transport in k-direction
+         ENDIF
+
+         IF( iom_use('weiv_heattr3d') ) THEN
+            zztmp = 0.5 * rcp 
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj,jk+1,jp_tem) )
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( z3d_T, 'T', 1. )
+            CALL iom_put( "weiv_heattr3d", zztmp * z3d_T )                 ! 3D heat transport in k-direction
+         ENDIF
+
+         IF( iom_use('weiv_salttr3d') ) THEN
+            zztmp = 0.5 * 0.001 
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj,jk+1,jp_sal) )
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( z3d_T, 'T', 1. )
+            CALL iom_put( "weiv_salttr3d", zztmp * z3d_T )                 ! 3D salt transport in k-direction
+         ENDIF
+
     END IF
 !
@@ -244,5 +337,5 @@
 # if defined key_diaeiv 
       CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
-      IF( ln_diaptr ) CALL wrk_dealloc( jpi, jpj, jpk, z3d )
+      CALL wrk_dealloc( jpi, jpj, jpk, z3d, z3d_T )
 # else
       CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 7083)
@@ -184,12 +184,11 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
                   ! total intermediate advective trends
-                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj)
                   ! update and guess with monotonic sheme
-                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra   * tmask(ji,jj,jk)
-                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
+                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -454,12 +453,11 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
                   ! total intermediate advective trends
-                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj)
                   ! update and guess with monotonic sheme
-                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra
-                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
+                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -479,8 +477,7 @@
          ! --------------------------------------------------
          ! antidiffusive flux on i and j
-
-
+         !
          DO jk = 1, jpkm1
-
+            !
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -513,4 +510,6 @@
          !
          ztrs(:,:,:,1) = ptb(:,:,:,jn)
+         ztrs(:,:,1,2) = ptb(:,:,1,jn)
+         ztrs(:,:,1,3) = ptb(:,:,1,jn)
          zwzts(:,:,:) = 0._wp
 
@@ -614,4 +613,5 @@
    END SUBROUTINE tra_adv_tvd_zts
 
+
    SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7083)
@@ -26,4 +26,6 @@
    USE ldfslp          ! iso-neutral slopes
    USE diaptr          ! poleward transport diagnostics
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O library
@@ -105,9 +107,11 @@
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::  ikt
-      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars
-      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
-      REAL(wp) ::  zcoef0, zbtr, ztra            !   -      -
+      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3       ! local scalars
+      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4       !   -      -
+      REAL(wp) ::  zcoef0, zbtr                      !   -      -
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ztrax, ztray, ztraz 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ztrax_T, ztray_T, ztraz_T
       !!----------------------------------------------------------------------
       !
@@ -115,5 +119,7 @@
       !
       CALL wrk_alloc( jpi, jpj,      z2d ) 
-      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t) 
+      ALLOCATE( ztrax(jpi,jpj,jpk), ztray(jpi,jpj,jpk), ztraz(jpi,jpj,jpk) ) 
+      IF( l_trdtra .and. cdtype == 'TRA' ) ALLOCATE( ztrax_T(jpi,jpj,jpk), ztray_T(jpi,jpj,jpk), ztraz_T(jpi,jpj,jpk) ) 
       !
 
@@ -127,4 +133,5 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
+         ztrax(:,:,:) = 0._wp ; ztray(:,:,:) = 0._wp ; ztraz(:,:,:) = 0._wp ; 
          !                                               
          !!----------------------------------------------------------------------
@@ -226,6 +233,6 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  ztrax(ji,jj,jk) = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) )
+                  ztray(ji,jj,jk) = zbtr * ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk) )
                END DO
             END DO
@@ -233,4 +240,6 @@
          END DO                                        !   End of slab  
          !                                             ! ===============
+         !
+         pta(:,:,:,jn) = pta(:,:,:,jn) + ztrax(:,:,:) + ztray(:,:,:)
          !
          ! "Poleward" diffusive heat or salt transports (T-S case only)
@@ -311,14 +320,35 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  ztra = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  ztraz(ji,jj,jk) = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
                END DO
             END DO
          END DO
+         pta(:,:,:,jn) = pta(:,:,:,jn) + ztraz(:,:,:)
          !
+         IF( l_trdtra .AND. cdtype == "TRA" .AND. jn .eq. 1 )  THEN      ! save the temperature trends
+            ztrax_T(:,:,:) = ztrax(:,:,:)
+            ztray_T(:,:,:) = ztray(:,:,:)
+            ztraz_T(:,:,:) = ztraz(:,:,:)
+         ENDIF
+         IF( l_trdtra .AND. cdtype == "TRC" )   THEN      ! save the horizontal component of diffusive trends for further diagnostics
+            CALL trd_tra( kt, cdtype, jn, jptra_iso_x, ztrax )
+            CALL trd_tra( kt, cdtype, jn, jptra_iso_y, ztray ) 
+            CALL trd_tra( kt, cdtype, jn, jptra_iso_z1, ztraz )  ! This is the first part of the vertical component.
+         ENDIF
       END DO
+      !
+      IF( l_trdtra .AND. cdtype == "TRA" )   THEN      ! save the horizontal component of diffusive trends for further diagnostics
+         CALL trd_tra( kt, cdtype, jp_tem, jptra_iso_x, ztrax_T )
+         CALL trd_tra( kt, cdtype, jp_sal, jptra_iso_x, ztrax )
+         CALL trd_tra( kt, cdtype, jp_tem, jptra_iso_y, ztray_T )
+         CALL trd_tra( kt, cdtype, jp_sal, jptra_iso_y, ztray )
+         CALL trd_tra( kt, cdtype, jp_tem, jptra_iso_z1, ztraz_T )  ! This is the first part of the vertical component
+         CALL trd_tra( kt, cdtype, jp_sal, jptra_iso_z1, ztraz )    !
+      ENDIF
       !
       CALL wrk_dealloc( jpi, jpj, z2d ) 
       CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
+      DEALLOCATE( ztrax, ztray, ztraz ) 
+      IF( l_trdtra ) DEALLOCATE( ztrax_T, ztray_T, ztraz_T ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso')
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 7083)
@@ -129,4 +129,5 @@
       IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+         ! Asselin filter trend
          ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
          ztrds(:,:,:) = tsn(:,:,:,jp_sal)
@@ -135,4 +136,13 @@
             CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
          ENDIF
+         ! total trend (note slightly inaccurate because tsb is time-filtered and tsa isn't)
+         ztrdt(:,:,:) = 0._wp
+         ztrds(:,:,:) = 0._wp
+         DO jk = 1, jpkm1
+            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / r2dtra(jk) 
+            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) 
+         END DO
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
       ENDIF
 
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 7083)
@@ -158,4 +158,5 @@
          ELSE                                         ! No restart or restart not found: Euler forward time stepping
             zfact = 1._wp
+            sbc_tsc(:,:,:) = 0._wp
             sbc_tsc_b(:,:,:) = 0._wp
          ENDIF
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90	(revision 7083)
@@ -33,5 +33,5 @@
 # endif
    !                                                  !!!* Active tracers trends indexes
-   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 14     !: Total trend nb: change it when adding/removing one indice below
+   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 20     !: Total trend nb: change it when adding/removing one indice below
    !                               ===============     !  
    INTEGER, PUBLIC, PARAMETER ::   jptra_xad  =  1     !: x- horizontal advection
@@ -39,19 +39,25 @@
    INTEGER, PUBLIC, PARAMETER ::   jptra_zad  =  3     !: z- vertical   advection
    INTEGER, PUBLIC, PARAMETER ::   jptra_sad  =  4     !: z- vertical   advection
-   INTEGER, PUBLIC, PARAMETER ::   jptra_ldf  =  5     !: lateral       diffusion
-   INTEGER, PUBLIC, PARAMETER ::   jptra_zdf  =  6     !: vertical      diffusion
-   INTEGER, PUBLIC, PARAMETER ::   jptra_zdfp =  7     !: "PURE" vert.  diffusion (ln_traldf_iso=T)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  =  8     !: Bottom Boundary Condition (geoth. heating) 
-   INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  =  9     !: Bottom Boundary Layer (diffusive and/or advective)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 10     !: non-penetrative convection treatment
-   INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 11     !: internal restoring (damping)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_qsr  = 12     !: penetrative solar radiation
-   INTEGER, PUBLIC, PARAMETER ::   jptra_nsr  = 13     !: non solar radiation / C/D on salinity  (+runoff if ln_rnf=T)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_atf  = 14     !: Asselin time filter
+   INTEGER, PUBLIC, PARAMETER ::   jptra_totad  =  5   !: total         advection
+   INTEGER, PUBLIC, PARAMETER ::   jptra_ldf  =  6     !: lateral       diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_iso_x  =  7   !: x-component of isopycnal diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_iso_y  =  8   !: y-component of isopycnal diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_iso_z1 =  9   !: z-component of isopycnal diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_zdf  = 10     !: vertical      diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_zdfp = 11     !: "PURE" vert.  diffusion (ln_traldf_iso=T)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_evd  = 12     !: EVD term (convection)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  = 13     !: Bottom Boundary Condition (geoth. heating) 
+   INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  = 14     !: Bottom Boundary Layer (diffusive and/or advective)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 15     !: non-penetrative convection treatment
+   INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 16     !: internal restoring (damping)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_qsr  = 17     !: penetrative solar radiation
+   INTEGER, PUBLIC, PARAMETER ::   jptra_nsr  = 18     !: non solar radiation / C/D on salinity  (+runoff if ln_rnf=T)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_atf  = 19     !: Asselin time filter
+   INTEGER, PUBLIC, PARAMETER ::   jptra_tot  = 20     !: Model total trend
    !
    !                                                  !!!* Passive tracers trends indices (use if "key_top" defined)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_sms  = 15     !: sources m. sinks
-   INTEGER, PUBLIC, PARAMETER ::   jptra_radn = 16     !: corr. trn<0 in trcrad
-   INTEGER, PUBLIC, PARAMETER ::   jptra_radb = 17     !: corr. trb<0 in trcrad (like atf)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_sms  = 19     !: sources m. sinks
+   INTEGER, PUBLIC, PARAMETER ::   jptra_radn = 20     !: corr. trn<0 in trcrad
+   INTEGER, PUBLIC, PARAMETER ::   jptra_radb = 21     !: corr. trb<0 in trcrad (like atf)
    !
    !                                                  !!!* Momentum trends indices
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90	(revision 7083)
@@ -91,5 +91,4 @@
 !!gm end
       !
-      IF( lk_vvl .AND. ( l_trdtra .OR. l_trddyn ) )  CALL ctl_stop( 'trend diagnostics with variable volume not validated' )
       
 !!gm  : Potential BUG : 3D output only for vector invariant form!  add a ctl_stop or code the flux form case
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 7083)
@@ -38,5 +38,6 @@
    REAL(wp) ::   r2dt   ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt   ! use to store the temperature trends
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt  ! use to store the temperature trends
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_evd  ! store avt_evd to calculate EVD trend
 
    !! * Substitutions
@@ -55,5 +56,5 @@
       !!                  ***  FUNCTION trd_tra_alloc  ***
       !!---------------------------------------------------------------------
-      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc )
+      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , avt_evd(jpi,jpj,jpk), STAT= trd_tra_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( trd_tra_alloc )
@@ -104,4 +105,5 @@
                                  ztrds(:,:,:) = 0._wp
                                  CALL trd_tra_mng( trdt, ztrds, ktrd, kt )
+         CASE( jptra_evd )   ;   avt_evd(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
          CASE DEFAULT                 ! other trends: masked trends
             trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)              ! mask & store
@@ -138,4 +140,18 @@
             END DO
             CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )  
+            !
+            !                         ! Also calculate EVD trend at this point. 
+            zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes
+            DO jk = 2, jpk
+               zwt(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
+               zws(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk)
+            END DO
+            !
+            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
+            DO jk = 1, jpkm1
+               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk)
+               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk) 
+            END DO
+            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt )  
             !
             CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt )
@@ -312,4 +328,6 @@
                                   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 )
@@ -318,4 +336,12 @@
       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 )
@@ -330,4 +356,6 @@
       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
       !
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 7083)
@@ -19,4 +19,6 @@
    USE zdf_oce         ! ocean vertical physics variables
    USE zdfkpp          ! KPP vertical mixing
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
    USE in_out_manager  ! I/O manager
    USE iom             ! for iom_put
@@ -122,4 +124,5 @@
       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd
       CALL iom_put( "avt_evd", zavt_evd )              ! output this change
+      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_evd')
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 7083)
@@ -140,8 +140,4 @@
                    PCO2a_in_cpl(jpi,jpj), Dust_in_cpl(jpi,jpj),     STAT=ierr(5) )
 
-         ! RSRH Temporarily initialise output coupling fields while we await clarification
-         ! of exactly how these will be initialised at model startup!
-         DMS_out_cpl(:,:) = 0.0
-         CO2Flux_out_cpl(:,:) = 0.0
       ENDIF
 #endif
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 7083)
@@ -76,5 +76,6 @@
       REAL(wp) ::   zchl
       REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zqsr100, zqsr_corr
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3
       !!---------------------------------------------------------------------
@@ -83,5 +84,6 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_alloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_alloc( jpi, jpj,      zqsr100, zqsr_corr )
       CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 )
 
@@ -112,8 +114,10 @@
       !                                        !  --------------------------------------
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
-         ! 1% of qsr to compute euphotic layer
+         !                                       ! 1% of qsr to compute euphotic layer
          zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr
          !
-         CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         !
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
          !
          DO jk = 1, nksrp      
@@ -123,5 +127,7 @@
          END DO
          !
-         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         !
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
          !
          DO jk = 1, nksrp      
@@ -133,5 +139,7 @@
          zqsr100(:,:) = 0.01 * qsr(:,:)
          !
-         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         !
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
          !
          DO jk = 1, nksrp      
@@ -226,5 +234,6 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_dealloc( jpi, jpj,      zqsr100, zqsr_corr )
       CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 )
       !
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7083)
@@ -136,5 +136,5 @@
                   zval = MAX( 1., zstrn(ji,jj) )
                   zval = 1.5 * zval / ( 12. + zval )
-                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
+                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval * ( 1. - fr_i(ji,jj) )
                   zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
                ENDIF
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 7083)
@@ -26,4 +26,6 @@
    USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine)
    USE ldftra_oce      ! lateral diffusion coefficient on tracers
+   USE trd_oce
+   USE trdtra
    USE prtctl_trc      ! Print control
 
@@ -74,4 +76,5 @@
       CHARACTER (len=22) ::   charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn  ! effective velocity
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd
       !!----------------------------------------------------------------------
       !
@@ -111,4 +114,9 @@
       IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )    ! add the mle transport (if necessary)
       !
+      IF( l_trdtrc )  THEN
+         CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd )
+         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+      ENDIF
+      !
       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
       CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
@@ -140,4 +148,12 @@
          !
       END SELECT
+      !
+      IF( l_trdtrc )   THEN                      ! save the advective trends for further diagnostics
+        DO jn = 1, jptra
+           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
+           CALL trd_tra( kt, 'TRC', jn, jptra_totad, ztrtrd(:,:,:,jn) )
+        END DO
+        CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
+      ENDIF
 
       !                                              ! print mean trends (used for debugging)
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 7083)
@@ -35,5 +35,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   restotr   ! restoring coeff. on tracers (s-1)
 
-   INTEGER, PARAMETER           ::   npncts   = 5        ! number of closed sea
+   INTEGER, PARAMETER           ::   npncts   = 8        ! number of closed sea
    INTEGER, DIMENSION(npncts)   ::   nctsi1, nctsj1      ! south-west closed sea limits (i,j)
    INTEGER, DIMENSION(npncts)   ::   nctsi2, nctsj2      ! north-east closed sea limits (i,j)
@@ -107,6 +107,5 @@
                
                jl = n_trc_index(jn) 
-               CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000
-               ztrcdta(:,:,:) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl)
+               CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
 
                SELECT CASE ( nn_zdmp_tr )
@@ -187,4 +186,5 @@
       INTEGER :: ji , jj, jk, jn, jl, jc                     ! dummy loop indicesa
       INTEGER :: isrow                                      ! local index
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta       ! 3D  workspace
 
       !!----------------------------------------------------------------------
@@ -207,6 +207,27 @@
             !
                                                         ! Caspian Sea
-            nctsi1(1)   = 332  ; nctsj1(1)   = 243 - isrow
-            nctsi2(1)   = 344  ; nctsj2(1)   = 275 - isrow
+            nctsi1(1)   = 333  ; nctsj1(1)   = 243 - isrow
+            nctsi2(1)   = 342  ; nctsj2(1)   = 274 - isrow
+            !                                           ! Lake Superior
+            nctsi1(2)   = 198  ; nctsj1(2)   = 258 - isrow
+            nctsi2(2)   = 204  ; nctsj2(2)   = 262 - isrow
+            !                                           ! Lake Michigan
+            nctsi1(3)   = 201  ; nctsj1(3)   = 250 - isrow
+            nctsi2(3)   = 203  ; nctsj2(3)   = 256 - isrow
+            !                                           ! Lake Huron
+            nctsi1(4)   = 204  ; nctsj1(4)   = 252 - isrow
+            nctsi2(4)   = 209  ; nctsj2(4)   = 256 - isrow
+            !                                           ! Lake Erie
+            nctsi1(5)   = 206  ; nctsj1(5)   = 249 - isrow
+            nctsi2(5)   = 209  ; nctsj2(5)   = 251 - isrow
+            !                                           ! Lake Ontario
+            nctsi1(6)   = 210  ; nctsj1(6)   = 252 - isrow
+            nctsi2(6)   = 212  ; nctsj2(6)   = 252 - isrow
+            !                                           ! Victoria Lake
+            nctsi1(7)   = 321  ; nctsj1(7)   = 180 - isrow
+            nctsi2(7)   = 322  ; nctsj2(7)   = 189 - isrow
+            !                                           ! Baltic Sea
+            nctsi1(8)   = 297  ; nctsj1(8)   = 270 - isrow
+            nctsi2(8)   = 308  ; nctsj2(8)   = 293 - isrow
             !                                        
             !                                           ! =======================
@@ -277,13 +298,15 @@
          IF(lwp)  WRITE(numout,*)
          !
+         CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )   ! Memory allocation
+         !
          DO jn = 1, jptra
             IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
                 jl = n_trc_index(jn)
-                CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000
+                CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
                 DO jc = 1, npncts
                    DO jk = 1, jpkm1
                       DO jj = nctsj1(jc), nctsj2(jc)
                          DO ji = nctsi1(jc), nctsi2(jc)
-                            trn(ji,jj,jk,jn) = sf_trcdta(jl)%fnow(ji,jj,jk) * tmask(ji,jj,jk) * rf_trfac(jl)
+                            trn(ji,jj,jk,jn) = ztrcdta(ji,jj,jk)
                             trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
                          ENDDO
@@ -293,5 +316,5 @@
              ENDIF
           ENDDO
-          !
+          CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta )
       ENDIF
       !
@@ -313,4 +336,6 @@
       IF( nn_timing == 1 )  CALL timing_start('trc_dmp_init')
       !
+      !Allocate arrays
+      IF( trc_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_dmp_init: unable to allocate arrays' )
 
       IF( lzoom )   nn_zdmp_tr = 0           ! restoring to climatology at closed north or south boundaries
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/trcdta.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 7083)
@@ -77,5 +77,5 @@
       ALLOCATE( n_trc_index(ntrc), slf_i(ntrc), STAT=ierr0 )
       IF( ierr0 > 0 ) THEN
-         CALL ctl_stop( 'trc_nam: unable to allocate n_trc_index' )   ;   RETURN
+         CALL ctl_stop( 'trc_dta_init: unable to allocate n_trc_index' )   ;   RETURN
       ENDIF
       nb_trcdta      = 0
@@ -91,4 +91,6 @@
       IF(lwp) THEN
          WRITE(numout,*) ' '
+         WRITE(numout,*) 'trc_dta_init : Passive tracers Initial Conditions '
+         WRITE(numout,*) '~~~~~~~~~~~~~~ '
          WRITE(numout,*) ' number of passive tracers to be initialize by data :', ntra
          WRITE(numout,*) ' '
@@ -107,14 +109,18 @@
          DO jn = 1, ntrc
             IF( ln_trc_ini(jn) )  THEN    ! open input file only if ln_trc_ini(jn) is true
-               clndta = TRIM( sn_trcdta(jn)%clvar ) 
-               clntrc = TRIM( ctrcnm   (jn)       ) 
+               clndta = TRIM( sn_trcdta(jn)%clvar )
+               if (jn > jptra) then
+                  clntrc='Dummy' ! By pass weird formats in ocean.output if ntrc > jptra
+               else
+                  clntrc = TRIM( ctrcnm   (jn)       )
+               endif
                zfact  = rn_trfac(jn)
-               IF( clndta /=  clntrc ) THEN 
-                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation :  ',   &
-                  &              'the variable name in the data file : '//clndta//   & 
-                  &              '  must be the same than the name of the passive tracer : '//clntrc//' ')
+               IF( clndta /=  clntrc ) THEN
+                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation    ',   &
+                  &              'Input name of data file : '//TRIM(clndta)//   &
+                  &              ' differs from that of tracer : '//TRIM(clntrc)//' ')
                ENDIF
-               WRITE(numout,*) ' read an initial file for passive tracer number :', jn, ' name : ', clndta, & 
-               &               ' multiplicative factor : ', zfact
+               WRITE(numout,'(a, i4,3a,e11.3)') ' Read IC file for tracer number :', &
+               &            jn, ', name : ', TRIM(clndta), ', Multiplicative Scaling factor : ', zfact
             ENDIF
          END DO
@@ -124,5 +130,5 @@
          ALLOCATE( sf_trcdta(nb_trcdta), rf_trfac(nb_trcdta), STAT=ierr1 )
          IF( ierr1 > 0 ) THEN
-            CALL ctl_stop( 'trc_dta_ini: unable to allocate  sf_trcdta structure' )   ;   RETURN
+            CALL ctl_stop( 'trc_dta_init: unable to allocate  sf_trcdta structure' )   ;   RETURN
          ENDIF
          !
@@ -135,5 +141,5 @@
                IF( sn_trcdta(jn)%ln_tint )  ALLOCATE( sf_trcdta(jl)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
                IF( ierr2 + ierr3 > 0 ) THEN
-                 CALL ctl_stop( 'trc_dta : unable to allocate passive tracer data arrays' )   ;   RETURN
+                 CALL ctl_stop( 'trc_dta_init : unable to allocate passive tracer data arrays' )   ;   RETURN
                ENDIF
             ENDIF
@@ -141,5 +147,5 @@
          ENDDO
          !                         ! fill sf_trcdta with slf_i and control print
-         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta', 'Passive tracer data', 'namtrc' )
+         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta_init', 'Passive tracer data', 'namtrc' )
          !
       ENDIF
@@ -151,5 +157,5 @@
 
 
-   SUBROUTINE trc_dta( kt, sf_dta )
+   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc)
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_dta  ***
@@ -164,9 +170,12 @@
       !!----------------------------------------------------------------------
       INTEGER                     , INTENT(in   ) ::   kt     ! ocean time-step
-      TYPE(FLD), DIMENSION(1)   , INTENT(inout) ::   sf_dta     ! array of information on the field to read
+      TYPE(FLD), DIMENSION(1)     , INTENT(inout) ::   sf_dta     ! array of information on the field to read
+      REAL(wp)                    , INTENT(in   ) ::   ptrfac  ! multiplication factor
+      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL  , INTENT(out  ) ::   ptrc
       !
       INTEGER ::   ji, jj, jk, jl, jkk, ik    ! dummy loop indices
       REAL(wp)::   zl, zi
       REAL(wp), DIMENSION(jpk) ::  ztp                ! 1D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta   ! 3D  workspace
       CHARACTER(len=100) :: clndta
       !!----------------------------------------------------------------------
@@ -176,5 +185,8 @@
       IF( nb_trcdta > 0 ) THEN
          !
+         CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )    ! Memory allocation
+         !
          CALL fld_read( kt, 1, sf_dta )      !==   read data at kt time step   ==!
+         ztrcdta(:,:,:) = sf_dta(1)%fnow(:,:,:) * tmask(:,:,:)    ! Mask
          !
          IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
@@ -185,51 +197,58 @@
             ENDIF
             !
-               DO jj = 1, jpj                         ! vertical interpolation of T & S
+            DO jj = 1, jpj                         ! vertical interpolation of T & S
+               DO ji = 1, jpi
+                  DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+                     zl = fsdept_n(ji,jj,jk)
+                     IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
+                        ztp(jk) = ztrcdta(ji,jj,1)
+                     ELSEIF( zl > gdept_1d(jpk) ) THEN         ! below the last level of data
+                        ztp(jk) =  ztrcdta(ji,jj,jpkm1)
+                     ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
+                        DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
+                           IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
+                              zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
+                              ztp(jk) = ztrcdta(ji,jj,jkk) + ( ztrcdta(ji,jj,jkk+1) - &
+                                        ztrcdta(ji,jj,jkk) ) * zi 
+                           ENDIF
+                        END DO
+                     ENDIF
+                  END DO
+                  DO jk = 1, jpkm1
+                    ztrcdta(ji,jj,jk) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
+                  END DO
+                  ztrcdta(ji,jj,jpk) = 0._wp
+                END DO
+            END DO
+            ! 
+         ELSE                                !==   z- or zps- coordinate   ==!
+            !
+            IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
+               DO jj = 1, jpj
                   DO ji = 1, jpi
-                     DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                        zl = fsdept_n(ji,jj,jk)
-                        IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
-                           ztp(jk) =  sf_dta(1)%fnow(ji,jj,1)
-                        ELSEIF( zl > gdept_1d(jpk) ) THEN         ! below the last level of data
-                           ztp(jk) =  sf_dta(1)%fnow(ji,jj,jpkm1)
-                        ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                           DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                              IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                                 zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                                 ztp(jk) = sf_dta(1)%fnow(ji,jj,jkk) + ( sf_dta(1)%fnow(ji,jj,jkk+1) - &
-                                           sf_dta(1)%fnow(ji,jj,jkk) ) * zi 
-                              ENDIF
-                           END DO
-                        ENDIF
-                     END DO
-                     DO jk = 1, jpkm1
-                        sf_dta(1)%fnow(ji,jj,jk) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                     END DO
-                     sf_dta(1)%fnow(ji,jj,jpk) = 0._wp
+                     ik = mbkt(ji,jj) 
+                     IF( ik > 1 ) THEN
+                        zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik-1)
+                     ENDIF
+                     ik = mikt(ji,jj)
+                     IF( ik > 1 ) THEN
+                        zl = ( fsdept_n(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
+                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik+1)
+                     ENDIF
                   END DO
                END DO
-            ! 
-         ELSE                                !==   z- or zps- coordinate   ==!
-            !                             
-               sf_dta(1)%fnow(:,:,:) = sf_dta(1)%fnow(:,:,:) * tmask(:,:,:)    ! Mask
-               !
-               IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        ik = mbkt(ji,jj) 
-                        IF( ik > 1 ) THEN
-                           zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                           sf_dta(1)%fnow(ji,jj,ik) = (1.-zl) * sf_dta(1)%fnow(ji,jj,ik) + zl * sf_dta(1)%fnow(ji,jj,ik-1)
-                        ENDIF
-                        ik = mikt(ji,jj)
-                        IF( ik > 1 ) THEN
-                           zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
-                           sf_dta(1)%fnow(ji,jj,ik) = (1.-zl) * sf_dta(1)%fnow(ji,jj,ik) + zl * sf_dta(1)%fnow(ji,jj,ik+1)
-                        ENDIF
-                     END DO
-                  END DO
-               ENDIF
-            !
-         ENDIF
+            ENDIF
+            !
+         ENDIF
+         !
+         ! Add multiplicative factor
+         ztrcdta(:,:,:) = ztrcdta(:,:,:) * ptrfac
+         !
+         ! Data structure for trc_ini (and BFMv5.1 coupling)
+         IF( .NOT. PRESENT(ptrc) ) sf_dta(1)%fnow(:,:,:) = ztrcdta(:,:,:)
+         !
+         ! Data structure for trc_dmp
+         IF( PRESENT(ptrc) )  ptrc(:,:,:) = ztrcdta(:,:,:)
          !
          IF( lwp .AND. kt == nit000 ) THEN
@@ -238,11 +257,14 @@
                WRITE(numout,*)
                WRITE(numout,*)'  level = 1'
-               CALL prihre( sf_dta(1)%fnow(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+               CALL prihre( ztrcdta(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
                WRITE(numout,*)'  level = ', jpk/2
-               CALL prihre( sf_dta(1)%fnow(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+               CALL prihre( ztrcdta(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
                WRITE(numout,*)'  level = ', jpkm1
-               CALL prihre( sf_dta(1)%fnow(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+               CALL prihre( ztrcdta(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
                WRITE(numout,*)
          ENDIF
+         !
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta )
+         !
       ENDIF
       !
@@ -255,5 +277,5 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_dta( kt, sf_dta, zrf_trfac )        ! Empty routine
+   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc)        ! Empty routine
       WRITE(*,*) 'trc_dta: You should not have seen this print! error?', kt
    END SUBROUTINE trc_dta
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 7082)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 7083)
@@ -125,7 +125,6 @@
                IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
                   jl = n_trc_index(jn) 
-                  CALL trc_dta( nit000, sf_trcdta(jl) )   ! read tracer data at nit000
-                  trn(:,:,:,jn) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl)
-                  !
+                  CALL trc_dta( nit000, sf_trcdta(jl), rf_trfac(jl) )   ! read tracer data at nit000
+                  trn(:,:,:,jn) = sf_trcdta(jl)%fnow(:,:,:) 
                   IF( .NOT.ln_trcdmp .AND. .NOT.ln_trcdmp_clo ) THEN      !== deallocate data structure   ==!
                      !                                                    (data used only for initialisation)
