Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 7352)
@@ -6,8 +6,4 @@
    !! History :  3.2  !  2009-11  (S. Masson)  Original code
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
-   !!----------------------------------------------------------------------
-#if defined key_diaar5
-   !!----------------------------------------------------------------------
-   !!   'key_diaar5'  :                           activate ar5 diagnotics
    !!----------------------------------------------------------------------
    !!   dia_ar5       : AR5 diagnostics
@@ -31,8 +27,6 @@
 
    PUBLIC   dia_ar5        ! routine called in step.F90 module
-   PUBLIC   dia_ar5_init   ! routine called in opa.F90 module
    PUBLIC   dia_ar5_alloc  ! routine called in nemogcm.F90 module
-
-   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .TRUE.   ! coupled flag
+   PUBLIC   dia_ar5_hst    ! heat/salt transport
 
    REAL(wp)                         ::   vol0         ! ocean volume (interior domain)
@@ -41,7 +35,10 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
+
+   LOGICAL  :: l_ar5
       
    !! * Substitutions
 #  include "zdfddm_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -80,5 +77,5 @@
       !
       REAL(wp), POINTER, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
-      REAL(wp), POINTER, DIMENSION(:,:)     :: pe                         ! 2D workspace 
+      REAL(wp), POINTER, DIMENSION(:,:)     :: zpe                         ! 2D workspace 
       REAL(wp), POINTER, DIMENSION(:,:,:)   :: zrhd , zrhop               ! 3D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
@@ -86,116 +83,128 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_ar5')
  
-      CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres, pe )
-      CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
-      CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
-
-      zarea_ssh(:,:) = area(:,:) * sshn(:,:)
-
-      !                                         ! total volume of liquid seawater
-      zvolssh = SUM( zarea_ssh(:,:) ) 
-      IF( lk_mpp )   CALL mpp_sum( zvolssh )
-      zvol = vol0 + zvolssh
+      IF( kt == nit000 )     CALL dia_ar5_init
+
+      IF( l_ar5 ) THEN
+         CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres )
+         CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
+         CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
+         zarea_ssh(:,:) = area(:,:) * sshn(:,:)
+      ENDIF
+      !
+      IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' )  .OR. iom_use( 'sshdyn' )  ) THEN    
+         !                                         ! total volume of liquid seawater
+         zvolssh = SUM( zarea_ssh(:,:) ) 
+         IF( lk_mpp )   CALL mpp_sum( zvolssh )
+         zvol = vol0 + zvolssh
       
-      CALL iom_put( 'voltot', zvol               )
-      CALL iom_put( 'sshtot', zvolssh / area_tot )
-      CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
-
-      !                     
-      ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
-      ztsn(:,:,:,jp_sal) = sn0(:,:,:)
-      CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
-      !
-      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
-      DO jk = 1, jpkm1
-         zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
-      END DO
-      IF( ln_linssh ) THEN
-         IF( ln_isfcav ) THEN
-            DO ji=1,jpi
-               DO jj=1,jpj
-                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
-               END DO
-            END DO
-         ELSE
-            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
-         END IF
+         CALL iom_put( 'voltot', zvol               )
+         CALL iom_put( 'sshtot', zvolssh / area_tot )
+         CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
+         !
+      ENDIF
+
+      IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN    
+         !                     
+         ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
+         ztsn(:,:,:,jp_sal) = sn0(:,:,:)
+         CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
+         !
+         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+         DO jk = 1, jpkm1
+            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+         END DO
+         IF( ln_linssh ) THEN
+            IF( ln_isfcav ) THEN
+               DO ji = 1, jpi
+                  DO jj = 1, jpj
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                  END DO
+               END DO
+            ELSE
+               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+            END IF
 !!gm
 !!gm   riceload should be added in both ln_linssh=T or F, no?
 !!gm
-      END IF
-      !                                         
-      zarho = SUM( area(:,:) * zbotpres(:,:) ) 
-      IF( lk_mpp )   CALL mpp_sum( zarho )
-      zssh_steric = - zarho / area_tot
-      CALL iom_put( 'sshthster', zssh_steric )
+         END IF
+         !                                         
+         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
+         IF( lk_mpp )   CALL mpp_sum( zarho )
+         zssh_steric = - zarho / area_tot
+         CALL iom_put( 'sshthster', zssh_steric )
       
-      !                                         ! steric sea surface height
-      CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
-      zrhop(:,:,jpk) = 0._wp
-      CALL iom_put( 'rhop', zrhop )
-      !
-      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
-      DO jk = 1, jpkm1
-         zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
-      END DO
-      IF( ln_linssh ) THEN
-         IF ( ln_isfcav ) THEN
-            DO ji=1,jpi
-               DO jj=1,jpj
-                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+         !                                         ! steric sea surface height
+         CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
+         zrhop(:,:,jpk) = 0._wp
+         CALL iom_put( 'rhop', zrhop )
+         !
+         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+         DO jk = 1, jpkm1
+            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+         END DO
+         IF( ln_linssh ) THEN
+            IF ( ln_isfcav ) THEN
+               DO ji = 1,jpi
+                  DO jj = 1,jpj
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                  END DO
+               END DO
+            ELSE
+               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+            END IF
+         END IF
+         !    
+         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
+         IF( lk_mpp )   CALL mpp_sum( zarho )
+         zssh_steric = - zarho / area_tot
+         CALL iom_put( 'sshsteric', zssh_steric )
+      
+         !                                         ! ocean bottom pressure
+         zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
+         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
+         CALL iom_put( 'botpres', zbotpres )
+         !
+      ENDIF
+
+      IF( iom_use( 'masstot' ) .OR. iom_use( 'temptot' )  .OR. iom_use( 'saltot' )  ) THEN    
+         !                                         ! Mean density anomalie, temperature and salinity
+         ztemp = 0._wp
+         zsal  = 0._wp
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
+                  ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
+                  zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
                END DO
             END DO
-         ELSE
-            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+         END DO
+         IF( ln_linssh ) THEN
+            IF( ln_isfcav ) THEN
+               DO ji = 1, jpi
+                  DO jj = 1, jpj
+                     ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
+                     zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
+                  END DO
+               END DO
+            ELSE
+               ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
+               zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
+            END IF
+         ENDIF
+         IF( lk_mpp ) THEN  
+            CALL mpp_sum( ztemp )
+            CALL mpp_sum( zsal  )
          END IF
-      END IF
-      !    
-      zarho = SUM( area(:,:) * zbotpres(:,:) ) 
-      IF( lk_mpp )   CALL mpp_sum( zarho )
-      zssh_steric = - zarho / area_tot
-      CALL iom_put( 'sshsteric', zssh_steric )
-      
-      !                                         ! ocean bottom pressure
-      zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
-      zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
-      CALL iom_put( 'botpres', zbotpres )
-
-      !                                         ! Mean density anomalie, temperature and salinity
-      ztemp = 0._wp
-      zsal  = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
-               ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
-               zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
-            END DO
-         END DO
-      END DO
-      IF( ln_linssh ) THEN
-         IF( ln_isfcav ) THEN
-            DO ji=1,jpi
-               DO jj=1,jpj
-                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
-                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
-               END DO
-            END DO
-         ELSE
-            ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
-            zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
-         END IF
-      ENDIF
-      IF( lk_mpp ) THEN  
-         CALL mpp_sum( ztemp )
-         CALL mpp_sum( zsal  )
-      END IF
-      !
-      zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater
-      ztemp = ztemp / zvol                            ! potential temperature in liquid seawater
-      zsal  = zsal  / zvol                            ! Salinity of liquid seawater
-      !
-      CALL iom_put( 'masstot', zmass )
-      CALL iom_put( 'temptot', ztemp )
-      CALL iom_put( 'saltot' , zsal  )
+         !
+         zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater
+         ztemp = ztemp / zvol                            ! potential temperature in liquid seawater
+         zsal  = zsal  / zvol                            ! Salinity of liquid seawater
+         !
+         CALL iom_put( 'masstot', zmass )
+         CALL iom_put( 'temptot', ztemp )
+         CALL iom_put( 'saltot' , zsal  )
+         !
+      ENDIF
 
       IF( iom_use( 'tnpeo' )) THEN    
@@ -203,5 +212,6 @@
       ! Exclude points where rn2 is negative as convection kicks in here and
       ! work is not being done against stratification
-          pe(:,:) = 0._wp
+          CALL wrk_alloc( jpi, jpj, zpe )
+          zpe(:,:) = 0._wp
           IF( lk_zdfddm ) THEN
              DO ji=1,jpi
@@ -214,5 +224,5 @@
                       zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
                       !
-                      pe(ji, jj) = pe(ji, jj) - MIN(0._wp, rn2(ji,jj,jk)) * &
+                      zpe(ji, jj) = zpe(ji, jj) - MIN(0._wp, rn2(ji,jj,jk)) * &
                            &       grav * (avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
                            &       - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
@@ -222,23 +232,84 @@
              ENDDO
           ELSE
-             DO ji=1,jpi
-                DO jj=1,jpj
-                   DO jk=1,jpk
-                       pe(ji,jj) = pe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * fse3w(ji, jj, jk)
+             DO ji = 1, jpi
+                DO jj = 1, jpj
+                   DO jk = 1, jpk
+                       zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
                    ENDDO
                 ENDDO
              ENDDO
           ENDIF
-          CALL lbc_lnk(pe, 'T', 1._wp)         
-          CALL iom_put( 'tnpeo', pe )
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres, pe )
-      CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    )
-      CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 )
+          CALL lbc_lnk( zpe, 'T', 1._wp)         
+          CALL iom_put( 'tnpeo', zpe )
+          CALL wrk_dealloc( jpi, jpj, zpe )
+      ENDIF
+      !
+      IF( l_ar5 ) THEN
+        CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres )
+        CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    )
+        CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 )
+      ENDIF
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5')
       !
    END SUBROUTINE dia_ar5
+
+   SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva ) 
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dia_ar5_htr ***
+      !!----------------------------------------------------------------------
+      !! Wrapper for heat transport calculations
+      !! Called from all advection and/or diffusion routines
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
+      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pua   ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      !
+      INTEGER    ::  ji, jj, jk
+      REAL(wp), POINTER, DIMENSION(:,:)  :: z2d
+
+    
+
+      CALL wrk_alloc( jpi, jpj, z2d )
+      z2d(:,:) = pua(:,:,1) 
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk) 
+            END DO
+         END DO
+       END DO
+       CALL lbc_lnk( z2d, 'U', -1. )
+       IF( cptr == 'adv' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in i-direction
+          IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rau0     * z2d )  ! advective salt transport in i-direction
+       ENDIF
+       IF( cptr == 'ldf' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in i-direction
+          IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rau0     * z2d ) ! diffusive salt transport in i-direction
+       ENDIF
+       !
+       z2d(:,:) = pva(:,:,1) 
+       DO jk = 1, jpkm1
+          DO jj = 2, jpjm1
+             DO ji = fs_2, fs_jpim1   ! vector opt.
+                z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk) 
+             END DO
+          END DO
+       END DO
+       CALL lbc_lnk( z2d, 'V', -1. )
+       IF( cptr == 'adv' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in j-direction
+          IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rau0     * z2d )  ! advective salt transport in j-direction
+       ENDIF
+       IF( cptr == 'ldf' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in j-direction
+          IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rau0     * z2d ) ! diffusive salt transport in j-direction
+       ENDIF
+          
+       CALL wrk_dealloc( jpi, jpj, z2d )
+
+   END SUBROUTINE dia_ar5_hst
 
 
@@ -259,59 +330,54 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
       !
-      CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
-      !                                      ! allocate dia_ar5 arrays
-      IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
-
-      area(:,:) = e1e2t(:,:) * tmask_i(:,:)
-
-      area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot )
-
-      vol0        = 0._wp
-      thick0(:,:) = 0._wp
-      DO jk = 1, jpkm1
-         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
-         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( vol0 )
-
-
-      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(:,:,:)
-      IF( ln_zps ) THEN               ! z-coord. partial steps
-         DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
-            DO ji = 1, jpi
-               ik = mbkt(ji,jj)
-               IF( ik > 1 ) THEN
-                  zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                  sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
-               ENDIF
+      l_ar5 = .FALSE.
+      IF(   iom_use( 'voltot'  ) .OR. iom_use( 'sshtot'    )  .OR. iom_use( 'sshdyn' )  .OR.  & 
+         &  iom_use( 'masstot' ) .OR. iom_use( 'temptot'   )  .OR. iom_use( 'saltot' ) .OR.  &    
+         &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) L_ar5 = .TRUE.
+  
+      IF( l_ar5 ) THEN
+         !
+         CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
+         !                                      ! allocate dia_ar5 arrays
+         IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
+
+         area(:,:) = e1e2t(:,:) * tmask_i(:,:)
+
+         area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot )
+
+         vol0        = 0._wp
+         thick0(:,:) = 0._wp
+         DO jk = 1, jpkm1
+            vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
+            thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( vol0 )
+
+
+         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(:,:,:)
+         IF( ln_zps ) THEN               ! z-coord. partial steps
+            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
+               DO ji = 1, jpi
+                  ik = mbkt(ji,jj)
+                  IF( ik > 1 ) THEN
+                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
+                  ENDIF
+               END DO
             END DO
-         END DO
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+         ENDIF
+         !
+         CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+         !
+      ENDIF
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
       !
    END SUBROUTINE dia_ar5_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                         NO diaar5
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .FALSE.   ! coupled flag
-CONTAINS
-   SUBROUTINE dia_ar5_init    ! Dummy routine
-   END SUBROUTINE dia_ar5_init
-   SUBROUTINE dia_ar5( kt )   ! Empty routine
-      INTEGER ::   kt
-      WRITE(*,*) 'dia_ar5: You should not have seen this print! error?', kt
-   END SUBROUTINE dia_ar5
-#endif
 
    !!======================================================================
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 7352)
@@ -39,9 +39,9 @@
    PUBLIC   dia_ptr_init   ! call in step module
    PUBLIC   dia_ptr        ! call in step module
-   PUBLIC   dia_ptr_ohst_components        ! called from tra_ldf/tra_adv routines
+   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines
 
    !                                  !!** namelist  namptr  **
-   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_adv, htr_ldf, htr_eiv, htr_vt   !: Heat TRansports (adv, diff, Bolus.)
-   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   str_adv, str_ldf, str_eiv, str_vs   !: Salt TRansports (adv, diff, Bolus.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_adv, htr_ldf, htr_eiv   !: Heat TRansports (adv, diff, Bolus.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   str_adv, str_ldf, str_eiv   !: Salt TRansports (adv, diff, Bolus.)
    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_ove, str_ove   !: heat Salt TRansports ( overturn.)
    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_btr, str_btr   !: heat Salt TRansports ( barotropic )
@@ -340,35 +340,4 @@
          ENDIF
 
-         IF( iom_use("sopht_vt") .OR. iom_use("sopst_vs") ) THEN   
-            z2d(1,:) = htr_vt(:,1) * rc_pwatt        !  (conversion in PW) 
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sopht_vt'
-            CALL iom_put( TRIM(cl1), z2d )
-            z2d(1,:) = str_vs(:,1) * rc_ggram        !  (conversion in Gg)
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sopst_vs'
-            CALL iom_put( TRIM(cl1), z2d )
-            IF( ln_subbas ) THEN
-              DO jn=2,nptr
-               z2d(1,:) = htr_vt(:,jn) * rc_pwatt        !  (conversion in PW)
-               DO ji = 1, jpi
-                 z2d(ji,:) = z2d(1,:)
-               ENDDO
-               cl1 = TRIM('sopht_vt_'//clsubb(jn))                 
-               CALL iom_put( cl1, z2d )
-               z2d(1,:) = str_vs(:,jn) * rc_ggram        ! (conversion in Gg)
-               DO ji = 1, jpi
-                  z2d(ji,:) = z2d(1,:)
-               ENDDO
-               cl1 = TRIM('sopst_vs_'//clsubb(jn))                 
-               CALL iom_put( cl1, z2d )              
-              ENDDO
-            ENDIF
-         ENDIF
-
          IF( iom_use("sophteiv") .OR. iom_use("sopsteiv") ) THEN 
             z2d(1,:) = htr_eiv(:,1) * rc_pwatt        !  (conversion in PW) 
@@ -482,5 +451,4 @@
          htr_ldf(:,:) = 0._wp  ;  str_ldf(:,:) =  0._wp 
          htr_eiv(:,:) = 0._wp  ;  str_eiv(:,:) =  0._wp 
-         htr_vt(:,:) = 0._wp  ;   str_vs(:,:) =  0._wp
          htr_ove(:,:) = 0._wp  ;   str_ove(:,:) =  0._wp
          htr_btr(:,:) = 0._wp  ;   str_btr(:,:) =  0._wp
@@ -490,7 +458,7 @@
    END SUBROUTINE dia_ptr_init
 
-   SUBROUTINE dia_ptr_ohst_components( ktra, cptr, pva ) 
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE dia_ptr_ohst_components  ***
+   SUBROUTINE dia_ptr_hst( ktra, cptr, pva ) 
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dia_ptr_hst ***
       !!----------------------------------------------------------------------
       !! Wrapper for heat and salt transport calculations to calculate them for each basin
@@ -514,8 +482,4 @@
          IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) )
       ENDIF
-      IF( cptr == 'vts' ) THEN
-         IF( ktra == jp_tem )  htr_vt(:,1) = ptr_sj( pva(:,:,:) )
-         IF( ktra == jp_sal )  str_vs(:,1) = ptr_sj( pva(:,:,:) )
-      ENDIF
       !
       IF( ln_subbas ) THEN
@@ -557,19 +521,7 @@
              ENDIF
          ENDIF
-         IF( cptr == 'vts' ) THEN
-             IF( ktra == jp_tem ) THEN 
-                DO jn = 2, nptr
-                    htr_vt(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                 END DO
-             ENDIF
-             IF( ktra == jp_sal ) THEN 
-                DO jn = 2, nptr
-                   str_vs(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                END DO
-             ENDIF
-         ENDIF
          !
       ENDIF
-   END SUBROUTINE dia_ptr_ohst_components
+   END SUBROUTINE dia_ptr_hst
 
 
@@ -586,5 +538,4 @@
          &      htr_adv(jpj,nptr) , str_adv(jpj,nptr) ,   &
          &      htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) ,   &
-         &      htr_vt(jpj,nptr)  , str_vs(jpj,nptr)  ,   &
          &      htr_ove(jpj,nptr) , str_ove(jpj,nptr) ,   &
          &      htr_btr(jpj,nptr) , str_btr(jpj,nptr) ,   &
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7352)
@@ -17,6 +17,4 @@
    !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
    !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
-   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
-   !!                 !                     change name of output variables in dia_wri_state
    !!----------------------------------------------------------------------
 
@@ -27,12 +25,7 @@
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain
-   USE dynadv, ONLY: ln_dynadv_vec
    USE zdf_oce         ! ocean vertical physics
-   USE ldftra          ! lateral physics: eddy diffusivity coef.
-   USE ldfdyn          ! lateral physics: eddy viscosity   coef.
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
-   USE icb_oce         ! Icebergs
-   USE icbdia          ! Iceberg budgets
    USE sbcssr          ! restoring term toward SST/SSS climatology
    USE phycst          ! physical constants
@@ -41,21 +34,13 @@
    USE zdfddm          ! vertical  physics: double diffusion
    USE diahth          ! thermocline diagnostics
-   !
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
-   USE diatmb          ! Top,middle,bottom output
-   USE dia25h          ! 25h Mean output
    USE iom
    USE ioipsl
-
 #if defined key_lim2
    USE limwri_2 
-#elif defined key_lim3
-   USE limwri 
 #endif
    USE lib_mpp         ! MPP library
    USE timing          ! preformance summary
-   USE diurnal_bulk    ! diurnal warm layer
-   USE cool_skin       ! Cool skin
    USE wrk_nemo        ! working array
 
@@ -68,15 +53,10 @@
 
    INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
-   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
    INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
    INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
-   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
    INTEGER ::   ndex(1)                              ! ???
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
 
    !! * Substitutions
-#  include "zdfddm_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -89,12 +69,8 @@
    INTEGER FUNCTION dia_wri_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(2) :: ierr
-      !!----------------------------------------------------------------------
-      ierr = 0
-      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
-         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
-         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
-         !
-      dia_wri_alloc = MAXVAL(ierr)
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( ndex_hT(jpi*jpj), ndex_hU(jpi*jpj), ndex_hV(jpi*jpj), STAT=dia_wri_alloc )
       IF( lk_mpp )   CALL mpp_sum( dia_wri_alloc )
       !
@@ -115,306 +91,12 @@
       !! ** Purpose :   Standard output of opa: dynamics and tracer fields 
       !!      NETCDF format is used by default 
+      !!      Standalone surface scheme 
       !!
       !! ** Method  :  use iom_put
       !!----------------------------------------------------------------------
-      !!
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      !!
-      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
-      INTEGER                      ::   jkbot                   !
-      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
-      !!
-      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!----------------------------------------------------------------------
       ! 
-      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
-      ! 
-      CALL wrk_alloc( jpi , jpj      , z2d )
-      CALL wrk_alloc( jpi , jpj, jpk , z3d )
-      !
-      ! Output the initial state and forcings
-      IF( ninist == 1 ) THEN                       
-         CALL dia_wri_state( 'output.init', kt )
-         ninist = 0
-      ENDIF
-
-      ! Output of initial vertical scale factor
-      CALL iom_put("e3t_0", e3t_0(:,:,:) )
-      CALL iom_put("e3u_0", e3t_0(:,:,:) )
-      CALL iom_put("e3v_0", e3t_0(:,:,:) )
-      !
-      CALL iom_put( "e3t" , e3t_n(:,:,:) )
-      CALL iom_put( "e3u" , e3u_n(:,:,:) )
-      CALL iom_put( "e3v" , e3v_n(:,:,:) )
-      CALL iom_put( "e3w" , e3w_n(:,:,:) )
-      IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
-
-      CALL iom_put( "ssh" , sshn )                 ! sea surface height
-      
-      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
-      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
-      IF ( iom_use("sbt") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               jkbot = mbkt(ji,jj)
-               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
-            END DO
-         END DO
-         CALL iom_put( "sbt", z2d )                ! bottom temperature
-      ENDIF
-      
-      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
-      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
-      IF ( iom_use("sbs") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               jkbot = mbkt(ji,jj)
-               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
-            END DO
-         END DO
-         CALL iom_put( "sbs", z2d )                ! bottom salinity
-      ENDIF
-
-      IF ( iom_use("taubot") ) THEN                ! bottom stress
-         z2d(:,:) = 0._wp
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
-                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )      
-               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
-                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
-               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
-               !
-            ENDDO
-         ENDDO
-         CALL lbc_lnk( z2d, 'T', 1. )
-         CALL iom_put( "taubot", z2d )           
-      ENDIF
-         
-      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
-      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
-      IF ( iom_use("sbu") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               jkbot = mbku(ji,jj)
-               z2d(ji,jj) = un(ji,jj,jkbot)
-            END DO
-         END DO
-         CALL iom_put( "sbu", z2d )                ! bottom i-current
-      ENDIF
-      
-      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
-      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
-      IF ( iom_use("sbv") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               jkbot = mbkv(ji,jj)
-               z2d(ji,jj) = vn(ji,jj,jkbot)
-            END DO
-         END DO
-         CALL iom_put( "sbv", z2d )                ! bottom j-current
-      ENDIF
-
-      CALL iom_put( "woce", wn )                   ! vertical velocity
-      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
-         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rau0 * e1e2t(:,:)
-         DO jk = 1, jpk
-            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
-         END DO
-         CALL iom_put( "w_masstr" , z3d )  
-         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-      ENDIF
-
-      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
-      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
-      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
-
-      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) )
-      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) )
-
-      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
-         DO jj = 2, jpjm1                                    ! sst gradient
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmp  = tsn(ji,jj,1,jp_tem)
-               zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj)
-               zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
-               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
-                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
-            END DO
-         END DO
-         CALL lbc_lnk( z2d, 'T', 1. )
-         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
-         z2d(:,:) = SQRT( z2d(:,:) )
-         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
-      ENDIF
-         
-      ! clem: heat and salt content
-      IF( iom_use("heatc") ) THEN
-         z2d(:,:)  = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
-      ENDIF
-
-      IF( iom_use("saltc") ) THEN
-         z2d(:,:)  = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
-      ENDIF
-      !
-      IF ( iom_use("eken") ) THEN
-         rke(:,:,jk) = 0._wp                               !      kinetic energy 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zztmp   = 1._wp / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
-                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)    &
-                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) )  &
-                     &          *  zztmp 
-                  !
-                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)    &
-                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) )  &
-                     &          *  zztmp 
-                  !
-                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
-                  !
-               ENDDO
-            ENDDO
-         ENDDO
-         CALL lbc_lnk( rke, 'T', 1. )
-         CALL iom_put( "eken", rke )           
-      ENDIF
-      !
-      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
-      !
-      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
-         z3d(:,:,jpk) = 0.e0
-         z2d(:,:) = 0.e0
-         DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
-            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
-         END DO
-         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
-         CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum
-      ENDIF
-      
-      IF( iom_use("u_heattr") ) THEN
-         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( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
-      ENDIF
-
-      IF( iom_use("u_salttr") ) THEN
-         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_sal) + tsn(ji+1,jj,jk,jp_sal) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( z2d, 'U', -1. )
-         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
-      ENDIF
-
-      
-      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
-         z3d(:,:,jpk) = 0.e0
-         DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
-         END DO
-         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
-      ENDIF
-      
-      IF( iom_use("v_heattr") ) THEN
-         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( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
-      ENDIF
-
-      IF( iom_use("v_salttr") ) THEN
-         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_sal) + tsn(ji,jj+1,jk,jp_sal) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( z2d, 'V', -1. )
-         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
-      ENDIF
-
-      ! Vertical integral of temperature
-      IF( iom_use("tosmint") ) THEN
-         z2d(:,:)=0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( z2d, 'T', -1. )
-         CALL iom_put( "tosmint", z2d ) 
-      ENDIF
-
-      ! Vertical integral of salinity
-      IF( iom_use("somint") ) THEN
-         z2d(:,:)=0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( z2d, 'T', -1. )
-         CALL iom_put( "somint", z2d ) 
-      ENDIF
-
-      CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2)
-      !
-      CALL wrk_dealloc( jpi , jpj      , z2d )
-      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
-      !
-      ! If we want tmb values 
-
-      IF (ln_diatmb) THEN
-         CALL dia_tmb 
-      ENDIF 
-      IF (ln_dia25h) THEN
-         CALL dia_25h( kt )
-      ENDIF 
-
-      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
+      !! no relevant 2D arrays to write in iomput case
       !
    END SUBROUTINE dia_wri
@@ -437,6 +119,7 @@
       !!      Each nwrite time step, output the instantaneous or mean fields
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !
+      !!
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!
       LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
       CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
@@ -445,15 +128,8 @@
       INTEGER  ::   ierr                                     ! error code return from allocation
       INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
-      INTEGER  ::   jn, ierror                               ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
-      !
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
       IF( nn_timing == 1 )   CALL timing_start('dia_wri')
-      !
-                             CALL wrk_alloc( jpi,jpj      , zw2d )
-      IF( .NOT.ln_linssh )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d )
       !
       ! Output the initial state and forcings
@@ -471,5 +147,7 @@
 
       ! Define frequency of output and means
-      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
+      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
+      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
+      ENDIF
 #if defined key_diainstant
       zsto = nwrite * rdt
@@ -526,29 +204,5 @@
             &           "m", ipk, gdept_1d, nz_T, "down" )
          !                                                            ! Index of ocean points
-         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
          CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
-         !
-         IF( ln_icebergs ) THEN
-            !
-            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after 
-            !! that routine is called from nemogcm, so do it here immediately before its needed
-            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
-            IF( lk_mpp )   CALL mpp_sum( ierror )
-            IF( ierror /= 0 ) THEN
-               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
-               RETURN
-            ENDIF
-            !
-            !! iceberg vertical coordinate is class number
-            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
-               &           "number", nclasses, class_num, nb_T )
-            !
-            !! each class just needs the surface index pattern
-            ndim_bT = 3
-            DO jn = 1,nclasses
-               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
-            ENDDO
-            !
-         ENDIF
 
          ! Define the U grid FILE ( nid_U )
@@ -562,5 +216,4 @@
             &           "m", ipk, gdept_1d, nz_U, "down" )
          !                                                            ! Index of ocean points
-         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
          CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
 
@@ -575,145 +228,33 @@
             &          "m", ipk, gdept_1d, nz_V, "down" )
          !                                                            ! Index of ocean points
-         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
          CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
 
-         ! Define the W grid FILE ( nid_W )
-
-         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
-         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
-         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
-            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
-            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
-         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
-            &          "m", ipk, gdepw_1d, nz_W, "down" )
-
+         ! No W grid FILE
 
          ! Declare all the output fields as NETCDF variables
 
          !                                                                                      !!! nid_T : 3D
-         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
-            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
-            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-         IF(  .NOT.ln_linssh  ) THEN
-            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
-            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
-            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
-            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-         ENDIF
-         !                                                                                      !!! nid_T : 2D
-         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
+         CALL histdef( nid_T, "sst_m", "Sea Surface temperature"            , "C"      ,   &  ! sst
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sss_m", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         IF(  ln_linssh  ) THEN
-            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
-            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
-            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
-            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
-            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ENDIF
+         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! (sfx)
+             &         jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!
-         IF( ln_icebergs ) THEN
-            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
-               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
-            IF( ln_bergdia ) THEN
-               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
-                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
-                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
-            ENDIF
-         ENDIF
-
-         IF( .NOT. ln_cpl ) THEN
-            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ENDIF
-
-         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
-            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ENDIF
-         
-         clmx ="l_max(only(x))"    ! max index on a period
-!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX 
-!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
-#if defined key_diahth
-         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#endif
-
-         IF( ln_cpl .AND. nn_ice == 2 ) THEN
-            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ENDIF
 
          CALL histend( nid_T, snc4chunks=snc4set )
 
          !                                                                                      !!! nid_U : 3D
-         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
-            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
-         !                                                                                      !!! nid_U : 2D
+         CALL histdef( nid_U, "ssu_m", "Velocity component in x-direction", "m/s"   ,         &  ! ssu
+            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
          CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
             &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
@@ -722,26 +263,10 @@
 
          !                                                                                      !!! nid_V : 3D
-         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
-            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
-         !                                                                                      !!! nid_V : 2D
+         CALL histdef( nid_V, "ssv_m", "Velocity component in y-direction", "m/s",            &  ! ssv_m
+            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
          CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
             &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
 
          CALL histend( nid_V, snc4chunks=snc4set )
-
-         !                                                                                      !!! nid_W : 3D
-         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
-            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
-         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
-            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
-         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
-            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
-
-         IF( lk_zdfddm ) THEN
-            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
-               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
-         ENDIF
-         !                                                                                      !!! nid_W : 2D
-         CALL histend( nid_W, snc4chunks=snc4set )
 
          IF(lwp) WRITE(numout,*)
@@ -754,5 +279,5 @@
       ! ---------------------
 
-      ! ndex(1) est utilise ssi l'avant dernier argument est different de 
+      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de 
       ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
       ! donne le nombre d'elements, et ndex la liste des indices a sortir
@@ -763,102 +288,24 @@
       ENDIF
 
-      IF( .NOT.ln_linssh ) THEN
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
-         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
-         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
-      ELSE
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
-         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
-         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
-      ENDIF
-      IF( .NOT.ln_linssh ) THEN
-         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
-         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
-         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
-         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
-      ENDIF
-      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
-      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
-      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
+      ! Write fields on T grid
+      CALL histwrite( nid_T, "sst_m", it, sst_m, ndim_hT, ndex_hT )   ! sea surface temperature
+      CALL histwrite( nid_T, "sss_m", it, sss_m, ndim_hT, ndex_hT )   ! sea surface salinity
+      CALL histwrite( nid_T, "sowaflup", it, (emp - rnf )  , ndim_hT, ndex_hT )   ! upward water flux
       CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux 
                                                                                   ! (includes virtual salt flux beneath ice 
                                                                                   ! in linear free surface case)
-      IF( ln_linssh ) THEN
-         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
-         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
-         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
-         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
-      ENDIF
+
       CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
       CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
-      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
-      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
-!
-      IF( ln_icebergs ) THEN
-         !
-         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT )  
-         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
-         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT )  
-         !
-         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
-         !
-         IF( ln_bergdia ) THEN
-            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   )  
-            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   )  
-            !
-            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
-         ENDIF
-      ENDIF
-
-      IF( .NOT. ln_cpl ) THEN
-         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
-         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
-         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
-      ENDIF
-      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
-         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
-         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
-         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
-      ENDIF
-!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
-!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
-
-#if defined key_diahth
-      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
-      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
-      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
-      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
-#endif
-
-      IF( ln_cpl .AND. nn_ice == 2 ) THEN
-         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
-         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
-      ENDIF
-
-      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
+
+         ! Write fields on U grid
+      CALL histwrite( nid_U, "ssu_m"   , it, ssu_m         , ndim_hU, ndex_hU )   ! i-current speed
       CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
-      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
+         ! Write fields on V grid
+      CALL histwrite( nid_V, "ssv_m"   , it, ssv_m         , ndim_hV, ndex_hV )   ! j-current speed
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
-
-      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
-      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
-      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
-      IF( lk_zdfddm ) THEN
-         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
-      ENDIF
 
       ! 3. Close all files
@@ -868,9 +315,5 @@
          CALL histclo( nid_U )
          CALL histclo( nid_V )
-         CALL histclo( nid_W )
-      ENDIF
-      !
-                             CALL wrk_dealloc( jpi , jpj        , zw2d )
-      IF( .NOT.ln_linssh )   CALL wrk_dealloc( jpi , jpj , jpk  , zw3d )
+      ENDIF
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
@@ -902,4 +345,6 @@
       !!----------------------------------------------------------------------
       ! 
+      IF( nn_timing == 1 )   CALL timing_start('dia_wri_state')
+
       ! 0. Initialisation
       ! -----------------
@@ -932,26 +377,4 @@
       ! Declare all the output fields as NetCDF variables
 
-      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
-      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
-      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
-         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
-      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
-      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
-      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
-         !
-      CALL histdef( id_i, "ahtu"    , "u-eddy diffusivity"    , "m2/s"    ,   &   ! zonal current
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
-      CALL histdef( id_i, "ahtv"    , "v-eddy diffusivity"    , "m2/s"    ,   &   ! meridonal current
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
-      CALL histdef( id_i, "ahmt"    , "t-eddy viscosity"      , "m2/s"    ,   &   ! zonal current
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
-      CALL histdef( id_i, "ahmf"    , "f-eddy viscosity"      , "m2/s"    ,   &   ! meridonal current
-         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
-         !
       CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater 
          &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -966,15 +389,7 @@
       CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
          &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-      IF( .NOT.ln_linssh ) THEN
-         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      , &   ! t-point depth
-            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
-         CALL histdef( id_i, "vovvle3t", "T point thickness"     , "m"      , &   ! t-point depth
-            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
-      ENDIF
 
 #if defined key_lim2
       CALL lim_wri_state_2( kt, id_i, nh_i )
-#elif defined key_lim3
-      CALL lim_wri_state( kt, id_i, nh_i )
 #else
       CALL histend( id_i, snc4chunks=snc4set )
@@ -989,17 +404,5 @@
 
       ! Write all fields on T grid
-      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
-      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
-      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
-      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
-      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
-      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
-      !
-      CALL histwrite( id_i, "ahtu"    , kt, ahtu             , jpi*jpj*jpk, idex )    ! aht at u-point
-      CALL histwrite( id_i, "ahtv"    , kt, ahtv             , jpi*jpj*jpk, idex )    !  -  at v-point
-      CALL histwrite( id_i, "ahmt"    , kt, ahmt             , jpi*jpj*jpk, idex )    ! ahm at t-point
-      CALL histwrite( id_i, "ahmf"    , kt, ahmf             , jpi*jpj*jpk, idex )    !  -  at f-point
-      !
-      CALL histwrite( id_i, "sowaflup", kt, emp-rnf          , jpi*jpj    , idex )    ! freshwater budget
+      CALL histwrite( id_i, "sowaflup", kt, emp              , jpi*jpj    , idex )    ! freshwater budget
       CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
       CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
@@ -1008,8 +411,4 @@
       CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
 
-      IF(  .NOT.ln_linssh  ) THEN             
-         CALL histwrite( id_i, "vovvldep", kt, gdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth 
-         CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:) , jpi*jpj*jpk, idex )!  T-cell thickness  
-      END IF 
       ! 3. Close the file
       ! -----------------
@@ -1020,10 +419,11 @@
          CALL histclo( nid_U )
          CALL histclo( nid_V )
-         CALL histclo( nid_W )
       ENDIF
 #endif
+       
+      IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state')
       ! 
+
    END SUBROUTINE dia_wri_state
-
    !!======================================================================
 END MODULE diawri
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 7352)
@@ -24,5 +24,4 @@
    USE ldfslp          ! lateral diffusion: slope of iso-neutral surfaces
    USE ldfc1d_c2d      ! lateral diffusion: 1D & 2D cases 
-   USE diaar5, ONLY:   lk_diaar5
    USE diaptr
    !
@@ -733,77 +732,77 @@
       !      
       !
-      IF( lk_diaar5 .OR. ln_diaptr ) THEN                     !==  eiv heat transport: calculate and output  ==!
-         CALL wrk_alloc( jpi,jpj,   zw2d )
-         !
-         zztmp = 0.5_wp * rau0 * rcp 
-         IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d')) THEN
-           zw2d(:,:) = 0._wp 
-           zw3d(:,:,:) = 0._wp 
-           DO jk = 1, jpkm1
-              DO jj = 2, jpjm1
-                 DO ji = fs_2, fs_jpim1   ! vector opt.
-                    zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
-                       &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
-                    zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
-                 END DO
+      CALL wrk_alloc( jpi,jpj,   zw2d )
+      !
+      zztmp = 0.5_wp * rau0 * rcp 
+      IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
+        zw2d(:,:)   = 0._wp 
+        zw3d(:,:,:) = 0._wp 
+        DO jk = 1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
+                    &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
+                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
               END DO
            END DO
-           CALL lbc_lnk( zw2d, 'U', -1. )
-           CALL lbc_lnk( zw3d, 'U', -1. )
-           CALL iom_put( "ueiv_heattr", zztmp * zw2d )                  ! heat transport in i-direction
-           CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction
-         ENDIF
-         zw2d(:,:) = 0._wp 
-         zw3d(:,:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
-                     &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
-                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( zw2d, 'V', -1. )
-         CALL iom_put( "veiv_heattr", zztmp * zw2d )                  !  heat transport in i-direction
-         CALL iom_put( "veiv_heattr", zztmp * zw3d )                  !  heat transport in i-direction
-         CALL dia_ptr_ohst_components( jp_tem, 'eiv', 0.5 * zw3d )
-         !
-         zztmp = 0.5_wp * 0.5
-         IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
-           zw2d(:,:) = 0._wp 
-           zw3d(:,:,:) = 0._wp 
-           DO jk = 1, jpkm1
-              DO jj = 2, jpjm1
-                 DO ji = fs_2, fs_jpim1   ! vector opt.
-                    zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
-                       &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) ) 
-                    zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
-                 END DO
+        END DO
+        CALL lbc_lnk( zw2d, 'U', -1. )
+        CALL lbc_lnk( zw3d, 'U', -1. )
+        CALL iom_put( "ueiv_heattr"  , zztmp * zw2d )                  ! heat transport in i-direction
+        CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction
+      ENDIF
+      zw2d(:,:)   = 0._wp 
+      zw3d(:,:,:) = 0._wp 
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
+                  &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
+               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zw2d, 'V', -1. )
+      CALL iom_put( "veiv_heattr", zztmp * zw2d )                  !  heat transport in j-direction
+      CALL iom_put( "veiv_heattr", zztmp * zw3d )                  !  heat transport in j-direction
+      !
+      IF( ln_diaptr )  CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
+      !
+      zztmp = 0.5_wp * 0.5
+      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
+        zw2d(:,:) = 0._wp 
+        zw3d(:,:,:) = 0._wp 
+        DO jk = 1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
+                    &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) ) 
+                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
               END DO
            END DO
-           CALL lbc_lnk( zw2d, 'U', -1. )
-           CALL lbc_lnk( zw3d, 'U', -1. )
-           CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
-           CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
-         ENDIF
-         zw2d(:,:) = 0._wp 
-         zw3d(:,:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
-                     &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) ) 
-                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( zw2d, 'V', -1. )
-         CALL iom_put( "veiv_salttr", zztmp * zw2d )                  !  salt transport in j-direction
-         CALL iom_put( "veiv_salttr", zztmp * zw3d )                  !  salt transport in j-direction
-         CALL dia_ptr_ohst_components( jp_sal, 'eiv', 0.5 * zw3d )
-
-         CALL wrk_dealloc( jpi,jpj,   zw2d )
-      ENDIF
+        END DO
+        CALL lbc_lnk( zw2d, 'U', -1. )
+        CALL lbc_lnk( zw3d, 'U', -1. )
+        CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
+        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
+      ENDIF
+      zw2d(:,:) = 0._wp 
+      zw3d(:,:,:) = 0._wp 
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
+                  &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) ) 
+               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zw2d, 'V', -1. )
+      CALL iom_put( "veiv_salttr", zztmp * zw2d )                  !  salt transport in j-direction
+      CALL iom_put( "veiv_salttr", zztmp * zw3d )                  !  salt transport in j-direction
+      !
+      IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
+      !
+      CALL wrk_dealloc( jpi,jpj,   zw2d )
       CALL wrk_dealloc( jpi,jpj,jpk,   zw3d )
       !
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90	(revision 7352)
@@ -18,4 +18,5 @@
    USE trdtra         ! trends manager: tracers 
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    !
    USE in_out_manager ! I/O manager
@@ -32,4 +33,8 @@
    
    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
+
+   LOGICAL :: l_trd   ! flag to compute trends
+   LOGICAL :: l_ptr   ! flag to compute poleward transport
+   LOGICAL :: l_hst   ! flag to compute heat/salt transport
 
    !! * Substitutions
@@ -88,4 +93,12 @@
       ENDIF
       !
+      l_trd = .FALSE.
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )        l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
+      !
       !                    
       zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
@@ -184,11 +197,13 @@
          END DO
          !                             ! trend diagnostics
-         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN
+         IF( l_trd ) THEN
             CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
          END IF
-         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:) )
+         !                                 ! "Poleward" heat and salt transports 
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         !                                 !  heat and salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
       END DO
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7352)
@@ -20,4 +20,5 @@
    USE trdtra         ! tracers trends
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    USE phycst, ONLY: rau0_rcp
    !
@@ -38,5 +39,6 @@
 
    LOGICAL  ::   l_trd   ! flag to compute trends
-   LOGICAL  ::   l_trans   ! flag to output vertically integrated transports
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
 
@@ -98,14 +100,17 @@
       !
       l_trd = .FALSE.
-      l_trans = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
-      IF( cdtype == 'TRA' .AND. (iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") ) ) l_trans = .TRUE.
-      !
-      IF( l_trd .OR. l_trans )  THEN
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA'   .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )     l_trd = .TRUE.
+      IF(   cdtype == 'TRA'   .AND. ln_diaptr )                                              l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
+      !
+      IF( l_trd .OR. l_hst )  THEN
          CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
          ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
       ENDIF
       !
-      IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+      IF( l_ptr ) THEN  
          CALL wrk_alloc( jpi, jpj, jpk, zptry )
          zptry(:,:,:) = 0._wp
@@ -171,9 +176,9 @@
          CALL lbc_lnk( zwi, 'T', 1. )  ! Lateral boundary conditions on zwi  (unchanged sign)
          !                
-         IF( l_trd .OR. l_trans )  THEN             ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
          END IF
          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr )    zptry(:,:,:) = zwy(:,:,:) 
+         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:) 
          !
          !        !==  anti-diffusive flux : high order minus low order  ==!
@@ -299,5 +304,5 @@
          END DO
          !
-         IF( l_trd .OR. l_trans ) THEN     ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
@@ -311,43 +316,18 @@
             !
          END IF
-
-         IF( l_trans .AND. jn==jp_tem ) THEN
-            CALL wrk_alloc( jpi, jpj, z2d )
-            z2d(:,:) = 0._wp 
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + ztrdx(ji,jj,jk) 
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( z2d, 'U', -1. )
-            CALL iom_put( "uadv_heattr", rau0_rcp * z2d )       ! heat transport in i-direction
-              !
-            z2d(:,:) = 0._wp 
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + ztrdy(ji,jj,jk) 
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( z2d, 'V', -1. )
-            CALL iom_put( "vadv_heattr", rau0_rcp * z2d )       ! heat transport in j-direction
-            CALL wrk_dealloc( jpi, jpj, z2d )
+         !                                !  heat/salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
+
+         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( l_ptr ) THEN  
+            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
          ENDIF
-         ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
-            CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
-         ENDIF
          !
       END DO                     ! end of tracer loop
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
-      IF( l_trd .OR. l_trans )  THEN 
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
-      ENDIF
-      IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry )
+                              CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
+      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct')
@@ -412,12 +392,17 @@
       !
       l_trd = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
-      !
-      IF( l_trd )  THEN
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
+      !
+      IF( l_trd .OR. l_hst )  THEN
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
          ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
       ENDIF
       !
-      IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+      IF( l_ptr ) THEN  
          CALL wrk_alloc( jpi, jpj,jpk, zptry )
          zptry(:,:,:) = 0._wp
@@ -488,9 +473,9 @@
          CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign)
          !                
-         IF( l_trd )  THEN                ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd .OR. l_hst )  THEN                ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
          END IF
          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr )  zptry(:,:,:) = zwy(:,:,:)
+         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:)
 
          ! 3. anti-diffusive flux : high order minus low order
@@ -608,28 +593,33 @@
          END DO
 
-         !                                 ! trend diagnostics (contribution of upstream fluxes)
-         IF( l_trd )  THEN 
+        !
+         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
             ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
-            !
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )   
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )  
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
-            !
-            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
+         ENDIF
+            !
+         IF( l_trd ) THEN 
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
+            !
          END IF
-         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) 
-            CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
+         !                                             ! heat/salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
+
+         !                                            ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( l_ptr ) THEN  
+            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
          ENDIF
          !
       END DO
       !
-      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
-      CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
-      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
-      IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry )
+                              CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
+                              CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
+                              CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
+      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts')
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7352)
@@ -23,5 +23,8 @@
    USE sbcrnf         ! river runoffs
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
+
    !
+   USE iom
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
@@ -40,4 +43,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
    
+   LOGICAL  ::   l_trd   ! flag to compute trends
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -116,4 +123,12 @@
       ENDIF 
       !      
+      l_trd = .FALSE.
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
+      !
       DO jn = 1, kjpt            !==  loop over the tracers  ==!
          !
@@ -192,12 +207,12 @@
          END DO        
          !                                ! trend diagnostics
-         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   &
-            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN
+         IF( l_trd )  THEN
             CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
          END IF
-         !                                 ! "Poleward" heat and salt transports
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:)  )
-         !                                ! "Poleward" heat and salt transports
+         !                                 ! "Poleward" heat and salt transports 
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         !                                 !  heat transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
          !                          !* Vertical advective fluxes
@@ -260,7 +275,5 @@
          END DO
          !                                ! send trends for diagnostic
-         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     &
-            &( cdtype == 'TRC' .AND. l_trdtrc )      )   &
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
          !
       END DO                     ! end of tracer loop
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 7352)
@@ -34,6 +34,9 @@
    PUBLIC   tra_adv_qck   ! routine called by step.F90
 
-   LOGICAL  :: l_trd           ! flag to compute trends
    REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
+
+   LOGICAL  ::   l_trd   ! flag to compute trends
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+
 
    !! * Substitutions
@@ -103,5 +106,8 @@
       !
       l_trd = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      !
       !
       !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
@@ -224,5 +230,5 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
+         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
          !
       END DO
@@ -347,7 +353,7 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
+         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr )  CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )                     CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
          !
       END DO
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 7352)
@@ -19,5 +19,8 @@
    USE trdtra         ! trends manager: tracers 
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
+
    !
+   USE iom
    USE lib_mpp        ! I/O library
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
@@ -32,5 +35,8 @@
    PUBLIC   tra_adv_ubs   ! routine called by traadv module
 
-   LOGICAL :: l_trd  ! flag to compute trends or not
+   LOGICAL :: l_trd   ! flag to compute trends
+   LOGICAL :: l_ptr   ! flag to compute poleward transport
+   LOGICAL :: l_hst   ! flag to compute heat transport
+
 
    !! * Substitutions
@@ -109,5 +115,10 @@
       !
       l_trd = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
       !
       ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers
@@ -176,6 +187,10 @@
              CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) )
          END IF
-         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'adv', ztv(:,:,:) )
+         !     
+         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
+         !                                !  heati/salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
+         !
          !
          !                       !== vertical advective trend  ==!
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7352)
@@ -24,4 +24,5 @@
    USE ldfslp         ! iso-neutral slopes
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    !
    USE in_out_manager ! I/O manager
@@ -36,4 +37,7 @@
 
    PUBLIC   tra_ldf_iso   ! routine called by step.F90
+
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat transport
 
    !! * Substitutions
@@ -107,7 +111,4 @@
       REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
-#if defined key_diaar5
-      REAL(wp) ::   zztmp   ! local scalar
-#endif
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zdkt, zdk1t, z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdit, zdjt, zftu, zftv, ztfw 
@@ -127,5 +128,12 @@
          ah_wslp2(:,:,:) = 0._wp
       ENDIF
-      !                                               ! set time step size (Euler/Leapfrog)
+      !   
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
+      !
+      !                                            ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
@@ -369,38 +377,8 @@
             !
             !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
-               ! note sign is reversed to give down-gradient diffusive transports (#1043)
-            IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'ldf', -zftv(:,:,:)  )
-            !
-            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
-              !
-              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
-                  z2d(:,:) = zftu(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-!!gm CAUTION I think there is an error of sign when using BLP operator....
-!!gm         a multiplication by zsign is required (to be checked twice !)
-                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
-                  CALL lbc_lnk( z2d, 'U', -1. )
-                  CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
-                  !
-                  z2d(:,:) = zftv(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
-                  CALL lbc_lnk( z2d, 'V', -1. )
-                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
-               END IF
-               !
-            ENDIF
+               ! note sign is reversed to give down-gradient diffusive transports )
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:)  )
+            !                          ! Diffusive heat transports
+            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
             !
          ENDIF                                                    !== end pass selection  ==!
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90	(revision 7352)
@@ -17,4 +17,5 @@
    USE traldf_triad   ! iso-neutral lateral diffusion (triad    operator)     (tra_ldf_triad routine)
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    USE trc_oce        ! share passive tracers/Ocean variables
    USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
@@ -25,4 +26,5 @@
    USE timing         ! Timing
    USE wrk_nemo       ! Memory allocation
+   USE iom
 
    IMPLICIT NONE
@@ -39,4 +41,7 @@
    INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator
    INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator
+
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat transport
 
    !! * Substitutions
@@ -95,4 +100,10 @@
       CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
       !
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
+      !
       !                                !==  Initialization of metric arrays used for all tracers  ==!
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -150,5 +161,7 @@
          IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
              ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
-            IF( cdtype == 'TRA' .AND. ln_diaptr )    CALL dia_ptr_ohst_components( jn, 'ldf', -ztv(:,:,:) )
+
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:)  )
+            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
          ENDIF
          !                          ! ==================
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 7352)
@@ -20,4 +20,5 @@
    USE traldf_iso     ! lateral diffusion (Madec operator)         (tra_ldf_iso routine)
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
    !
@@ -35,4 +36,8 @@
 
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels
+
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat transport
+
 
    !! * Substitutions
@@ -89,7 +94,4 @@
       REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
       REAL(wp) ::   zah, zah_slp, zaei_slp
-#if defined key_diaar5
-      REAL(wp) ::   zztmp              ! local scalar
-#endif
       REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d                                            ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
@@ -112,5 +114,12 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
-      !                                               ! set time step size (Euler/Leapfrog)
+      !   
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
+      !
+      !                                                        ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
@@ -416,35 +425,7 @@
             !
             !                          ! "Poleward" diffusive heat or salt transports (T-S case only)
-            IF( cdtype == 'TRA' .AND. ln_diaptr )  CALL dia_ptr_ohst_components( jn, 'ldf', zftv(:,:,:) )
-            !
-            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
-              !
-              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
-                  z2d(:,:) = zftu(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-                  z2d(:,:) = rau0_rcp * z2d(:,:) 
-                  CALL lbc_lnk( z2d, 'U', -1. )
-                  CALL iom_put( "udiff_heattr", z2d )                  ! heat i-transport
-                  !
-                  z2d(:,:) = zftv(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-                  z2d(:,:) = rau0_rcp * z2d(:,:)     
-                  CALL lbc_lnk( z2d, 'V', -1. )
-                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat j-transport
-               ENDIF
-               !
-            ENDIF
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  )
+            !                          ! Diffusive heat transports
+            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
             !
          ENDIF                                                    !== end pass selection  ==!
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7352)
@@ -490,5 +490,4 @@
       IF( lk_floats     )   CALL     flo_init   ! drifting Floats
                             CALL dia_cfl_init   ! Initialise CFL diagnostics
-      IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag
                             CALL dia_ptr_init   ! Poleward TRansports initialization
       IF( lk_diadct     )   CALL dia_dct_init   ! Sections tranports
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7352)
@@ -234,5 +234,5 @@
       IF(.NOT.ln_cpl )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
       IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
-      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
+                         CALL dia_ar5( kstp )         ! ar5 diag
       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
                          CALL dia_wri( kstp )         ! ocean model: outputs
Index: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 7351)
+++ branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 7352)
@@ -114,5 +114,4 @@
    USE zdfmxl , ONLY :   hmlpt       =>   hmlpt       !: mixed layer depth at t-points (m)
 
-   USE diaar5 , ONLY :   lk_diaar5  =>   lk_diaar5
 #else
    !!----------------------------------------------------------------------
