Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diaar5.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diaar5.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diaar5.F90	(revision 12109)
@@ -71,11 +71,11 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
-      REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass
+      INTEGER  ::   ji, jj, jk, iks, ikb                      ! dummy loop arguments
+      REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass, zsst
       REAL(wp) ::   zaw, zbw, zrw
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zpe                         ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zrhd , zrhop               ! 3D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zpe, z2d                   ! 2D workspace 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zrhd , zrhop, ztpot   ! 3D workspace
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
 
@@ -86,5 +86,5 @@
 
       IF( l_ar5 ) THEN 
-         ALLOCATE( zarea_ssh(jpi,jpj) , zbotpres(jpi,jpj) )
+         ALLOCATE( zarea_ssh(jpi,jpj), zbotpres(jpi,jpj), z2d(jpi,jpj) )
          ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) )
          ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
@@ -92,9 +92,31 @@
       ENDIF
       !
+      CALL iom_put( 'e2u'      , e2u (:,:) )
+      CALL iom_put( 'e1v'      , e1v (:,:) )
+      CALL iom_put( 'areacello', area(:,:) )
+      !
+      IF( iom_use( 'volcello' ) .OR. iom_use( 'masscello' )  ) THEN  
+         zrhd(:,:,jpk) = 0._wp        ! ocean volume ; rhd is used as workspace
+         DO jk = 1, jpkm1
+            zrhd(:,:,jk) = area(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000
+         CALL iom_put( 'masscello' , rau0 * e3t_n(:,:,:) * tmask(:,:,:) )  ! ocean mass
+      ENDIF 
+      !
+      IF( iom_use( 'e3tb' ) )  THEN    ! bottom layer thickness
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikb = mbkt(ji,jj)
+               z2d(ji,jj) = e3t_n(ji,jj,ikb)
+            END DO
+         END DO
+         CALL iom_put( 'e3tb', z2d )
+      ENDIF 
+      !
       IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' )  .OR. iom_use( 'sshdyn' )  ) THEN    
          !                                         ! total volume of liquid seawater
-         zvolssh = SUM( zarea_ssh(:,:) ) 
-         CALL mpp_sum( 'diaar5', zvolssh )
-         zvol = vol0 + zvolssh
+         zvolssh = glob_sum( 'diaar5', zarea_ssh(:,:) ) 
+         zvol    = vol0 + zvolssh
       
          CALL iom_put( 'voltot', zvol               )
@@ -118,5 +140,6 @@
                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)
+                     iks = mikt(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,iks) + riceload(ji,jj)
                   END DO
                END DO
@@ -129,6 +152,5 @@
          END IF
          !                                         
-         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
-         CALL mpp_sum( 'diaar5', zarho )
+         zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) ) 
          zssh_steric = - zarho / area_tot
          CALL iom_put( 'sshthster', zssh_steric )
@@ -147,5 +169,6 @@
                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)
+                     iks = mikt(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,iks) + riceload(ji,jj)
                   END DO
                END DO
@@ -155,9 +178,7 @@
          END IF
          !    
-         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
-         CALL mpp_sum( 'diaar5', zarho )
+         zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) ) 
          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
@@ -168,48 +189,92 @@
 
       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
-         END DO
-         IF( ln_linssh ) THEN
+          !                                         ! Mean density anomalie, temperature and salinity
+          ztsn(:,:,:,:) = 0._wp                    ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity
+          DO jk = 1, jpkm1
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
+                   ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * tsn(ji,jj,jk,jp_tem)
+                   ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * tsn(ji,jj,jk,jp_sal)
+                ENDDO
+             ENDDO
+          ENDDO
+
+          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) 
+                     iks = mikt(ji,jj)
+                     ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * tsn(ji,jj,iks,jp_tem) 
+                     ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * tsn(ji,jj,iks,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) )
+               ztsn(:,:,1,jp_tem) = ztsn(:,:,1,jp_tem) + zarea_ssh(:,:) * tsn(:,:,1,jp_tem) 
+               ztsn(:,:,1,jp_sal) = ztsn(:,:,1,jp_sal) + zarea_ssh(:,:) * tsn(:,:,1,jp_sal) 
             END IF
          ENDIF
-         IF( lk_mpp ) THEN  
-            CALL mpp_sum( 'diaar5', ztemp )
-            CALL mpp_sum( 'diaar5', 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
+         !
+         ztemp = glob_sum( 'diaar5', ztsn(:,:,1,jp_tem) )
+         zsal  = glob_sum( 'diaar5', ztsn(:,:,1,jp_sal) )
+         zmass = rau0 * ( zarho + zvol )      
          !
          CALL iom_put( 'masstot', zmass )
-         CALL iom_put( 'temptot', ztemp )
-         CALL iom_put( 'saltot' , zsal  )
-         !
+         CALL iom_put( 'temptot', ztemp / zvol )
+         CALL iom_put( 'saltot' , zsal  / zvol )
+         !
+      ENDIF     
+
+      IF( ln_teos10 ) THEN        ! ! potential temperature (TEOS-10 case)
+         IF( iom_use( 'toce_pot') .OR. iom_use( 'temptot_pot' ) .OR. iom_use( 'sst_pot' )  &
+                                  .OR. iom_use( 'ssttot' ) .OR.  iom_use( 'tosmint_pot' ) ) THEN
+            !
+            ALLOCATE( ztpot(jpi,jpj,jpk) )
+            ztpot(:,:,jpk) = 0._wp
+            DO jk = 1, jpkm1
+               ztpot(:,:,jk) = eos_pt_from_ct( tsn(:,:,jk,jp_tem), tsn(:,:,jk,jp_sal) )
+            END DO
+            !
+            CALL iom_put( 'toce_pot', ztpot(:,:,:) )  ! potential temperature (TEOS-10 case)
+            CALL iom_put( 'sst_pot' , ztpot(:,:,1) )  ! surface temperature
+            !
+            IF( iom_use( 'temptot_pot' ) ) THEN   ! Output potential temperature in case we use TEOS-10
+               z2d(:,:) = 0._wp
+               DO jk = 1, jpkm1
+                 z2d(:,:) = z2d(:,:) + area(:,:) * e3t_n(:,:,jk) * ztpot(:,:,jk)
+               END DO
+               ztemp = glob_sum( 'diaar5', z2d(:,:)  ) 
+               CALL iom_put( 'temptot_pot', ztemp / zvol )
+             ENDIF
+             !
+             IF( iom_use( 'ssttot' ) ) THEN   ! Output potential temperature in case we use TEOS-10
+               zsst = glob_sum( 'diaar5',  area(:,:) * ztpot(:,:,1)  ) 
+               CALL iom_put( 'ssttot', zsst / area_tot )
+             ENDIF
+             ! Vertical integral of temperature
+             IF( iom_use( 'tosmint_pot') ) THEN
+               z2d(:,:) = 0._wp
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi   ! vector opt.
+                        z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  ztpot(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               CALL iom_put( 'tosmint_pot', z2d ) 
+            ENDIF
+            DEALLOCATE( ztpot )
+        ENDIF
+      ELSE       
+         IF( iom_use('ssttot') ) THEN   ! Output sst in case we use EOS-80
+            zsst  = glob_sum( 'diaar5', area(:,:) * tsn(:,:,1,jp_tem) )
+            CALL iom_put('ssttot', zsst / area_tot )
+         ENDIF
       ENDIF
 
       IF( iom_use( 'tnpeo' )) THEN    
-      ! Work done against stratification by vertical mixing
-      ! Exclude points where rn2 is negative as convection kicks in here and
-      ! work is not being done against stratification
+        ! Work done against stratification by vertical mixing
+        ! Exclude points where rn2 is negative as convection kicks in here and
+        ! work is not being done against stratification
          ALLOCATE( zpe(jpi,jpj) )
          zpe(:,:) = 0._wp
@@ -219,14 +284,10 @@
                   DO ji = 1, jpi
                      IF( rn2(ji,jj,jk) > 0._wp ) THEN
-                        zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
-                           &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
-!!gm  this can be reduced to :  (depw-dept) / e3w   (NB idem dans bn2 !)
-!                        zrw =   ( gdept_n(ji,jj,jk) - gdepw_n(ji,jj,jk) ) / e3w_n(ji,jj,jk)
-!!gm end
+                        zrw = ( gdept_n(ji,jj,jk) - gdepw_n(ji,jj,jk) ) / e3w_n(ji,jj,jk)
                         !
                         zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
                         zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
                         !
-                        zpe(ji, jj) = zpe(ji, jj)            &
+                        zpe(ji, jj) = zpe(ji,jj)   &
                            &        -  grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
                            &                   - avs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
@@ -239,11 +300,9 @@
                DO ji = 1, jpi
                   DO jj = 1, jpj
-                     zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
+                     zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rau0 * e3w_n(ji,jj,jk)
                   END DO
                END DO
             END DO
          ENDIF
-!!gm useless lbc_lnk since the computation above is performed over 1:jpi & 1:jpj
-!!gm           CALL lbc_lnk( 'diaar5', zpe, 'T', 1._wp)         
           CALL iom_put( 'tnpeo', zpe )
           DEALLOCATE( zpe )
@@ -251,5 +310,5 @@
 
       IF( l_ar5 ) THEN
-        DEALLOCATE( zarea_ssh , zbotpres )
+        DEALLOCATE( zarea_ssh , zbotpres, z2d )
         DEALLOCATE( zrhd      , zrhop    )
         DEALLOCATE( ztsn                 )
@@ -287,10 +346,10 @@
        CALL lbc_lnk( 'diaar5', 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
+          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
+          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
        !
@@ -305,10 +364,10 @@
        CALL lbc_lnk( 'diaar5', 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
+          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
+          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
           
@@ -323,8 +382,9 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   inum
-      INTEGER  ::   ik
+      INTEGER  ::   ik, idep
       INTEGER  ::   ji, jj, jk  ! dummy loop indices
       REAL(wp) ::   zztmp  
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zvol0     
       !
       !!----------------------------------------------------------------------
@@ -340,18 +400,24 @@
          IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
 
-         area(:,:) = e1e2t(:,:) * tmask_i(:,:)
-
-         area_tot = SUM( area(:,:) )   ;   CALL mpp_sum( 'diaar5', area_tot )
-
-         vol0        = 0._wp
+         area(:,:) = e1e2t(:,:)
+         area_tot  = glob_sum( 'diaar5', area(:,:) )
+
+         ALLOCATE( zvol0(jpi,jpj) )
+         zvol0 (:,:) = 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
-         CALL mpp_sum( 'diaar5', vol0 )
+            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
+               DO ji = 1, jpi
+                  idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
+                  zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * area(ji,jj)
+                  thick0(ji,jj) = thick0(ji,jj) +  idep    
+               END DO
+            END DO
+         END DO
+         vol0 = glob_sum( 'diaar5', zvol0 )
+         DEALLOCATE( zvol0 )
 
          IF( iom_use( 'sshthster' ) ) THEN
-            ALLOCATE( zsaldta(jpi,jpj,jpj,jpts) )
+            ALLOCATE( zsaldta(jpi,jpj,jpk,jpts) )
             CALL iom_open ( 'sali_ref_clim_monthly', inum )
             CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diahth.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diahth.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diahth.F90	(revision 12109)
@@ -11,8 +11,4 @@
    !!            3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag
    !!----------------------------------------------------------------------
-#if   defined key_diahth
-   !!----------------------------------------------------------------------
-   !!   'key_diahth' :                              thermocline depth diag.
-   !!----------------------------------------------------------------------
    !!   dia_hth      : Compute varius diagnostics associated with the mixed layer
    !!----------------------------------------------------------------------
@@ -32,11 +28,15 @@
    PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90
 
-   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag
+   LOGICAL, SAVE  ::   l_hth     !: thermocline-20d depths flag
    
    ! note: following variables should move to local variables once iom_put is always used 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd26   !: depth of 26 C isotherm                         [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc7   !: heat content of first 700 m                    [W]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc20  !: heat content of first 2000 m                   [W]
+
 
    !!----------------------------------------------------------------------
@@ -52,5 +52,6 @@
       !!---------------------------------------------------------------------
       !
-      ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc )
+      ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd26(jpi,jpj), hd28(jpi,jpj), &
+         &      htc3(jpi,jpj), htc7(jpi,jpj), htc20(jpi,jpj), STAT=dia_hth_alloc )
       !
       CALL mpp_sum ( 'diahth', dia_hth_alloc )
@@ -82,270 +83,309 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !!
-      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
-      INTEGER                          ::   iid, ilevel           ! temporary integers
-      INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels
-      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth
-      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
-      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
-      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
-      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars
-      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop
-      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3      
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness 
-      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
+      INTEGER                      ::   ji, jj, jk            ! dummy loop arguments
+      REAL(wp)                     ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
+      REAL(wp)                     ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
+      REAL(wp)                     ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
+      REAL(wp)                     ::   zztmp, zzdep          ! temporary scalars inside do loop
+      REAL(wp)                     ::   zu, zv, zw, zut, zvt  ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3      
+      REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
       !!----------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_hth')
 
       IF( kt == nit000 ) THEN
+         l_hth = .FALSE.
+         IF(   iom_use( 'mlddzt'   ) .OR. iom_use( 'mldr0_3'  ) .OR. iom_use( 'mldr0_1'  )    .OR.  & 
+            &  iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' )    .OR.  &    
+            &  iom_use( '20d'      ) .OR. iom_use( '26d'      ) .OR. iom_use( '28d'      )    .OR.  &    
+            &  iom_use( 'hc300'    ) .OR. iom_use( 'hc700'    ) .OR. iom_use( 'hc2000'   )    .OR.  &    
+            &  iom_use( 'pycndep'  ) .OR. iom_use( 'tinv'     ) .OR. iom_use( 'depti'    )  ) l_hth = .TRUE.
          !                                      ! allocate dia_hth array
-         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' )
-
-         IF(.NOT. ALLOCATED(ik20) ) THEN
-            ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), &
-               &      zabs2(jpi,jpj),   &
-               &      ztm2(jpi,jpj),    &
-               &      zrho10_3(jpi,jpj),&
-               &      zpycn(jpi,jpj),   &
-               &      ztinv(jpi,jpj),   &
-               &      zdepinv(jpi,jpj), &
-               &      zrho0_3(jpi,jpj), &
-               &      zrho0_1(jpi,jpj), &
-               &      zmaxdzT(jpi,jpj), &
-               &      zthick(jpi,jpj),  &
-               &      zdelr(jpi,jpj), STAT=ji)
-            CALL mpp_sum('diahth', ji)
-            IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' )
-         END IF
-
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
-         IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         IF(lwp) WRITE(numout,*)
+         IF( l_hth ) THEN 
+            IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' )
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
+            IF(lwp) WRITE(numout,*) '~~~~~~~ '
+            IF(lwp) WRITE(numout,*)
+         ENDIF
       ENDIF
 
-      ! initialization
-      ztinv  (:,:) = 0._wp  
-      zdepinv(:,:) = 0._wp  
-      zmaxdzT(:,:) = 0._wp  
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
-            hth     (ji,jj) = zztmp
-            zabs2   (ji,jj) = zztmp
-            ztm2    (ji,jj) = zztmp
-            zrho10_3(ji,jj) = zztmp
-            zpycn   (ji,jj) = zztmp
-        END DO
-      END DO
-      IF( nla10 > 1 ) THEN 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
-               zrho0_3(ji,jj) = zztmp
-               zrho0_1(ji,jj) = zztmp
-            END DO
-         END DO
+      IF( l_hth ) THEN
+         !
+         IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
+            ! initialization
+            ztinv  (:,:) = 0._wp  
+            zdepinv(:,:) = 0._wp  
+            zmaxdzT(:,:) = 0._wp  
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
+                  hth     (ji,jj) = zztmp
+                  zabs2   (ji,jj) = zztmp
+                  ztm2    (ji,jj) = zztmp
+                  zrho10_3(ji,jj) = zztmp
+                  zpycn   (ji,jj) = zztmp
+                 END DO
+            END DO
+            IF( nla10 > 1 ) THEN 
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
+                     zrho0_3(ji,jj) = zztmp
+                     zrho0_1(ji,jj) = zztmp
+                  END DO
+               END DO
+            ENDIF
+      
+            ! Preliminary computation
+            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( tmask(ji,jj,nla10) == 1. ) THEN
+                     zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)  &
+                        &           - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
+                        &           - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
+                     zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)  &
+                        &           - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
+                     zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
+                     zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
+                     zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
+                     zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
+                  ELSE
+                     zdelr(ji,jj) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+
+            ! ------------------------------------------------------------- !
+            ! thermocline depth: strongest vertical gradient of temperature !
+            ! turbocline depth (mixing layer depth): avt = zavt5            !
+            ! MLD: rho = rho(1) + zrho3                                     !
+            ! MLD: rho = rho(1) + zrho1                                     !
+            ! ------------------------------------------------------------- !
+            DO jk = jpkm1, 2, -1   ! loop from bottom to 2
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     !
+                     zzdep = gdepw_n(ji,jj,jk)
+                     zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) &
+                            & / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
+                     zzdep = zzdep * tmask(ji,jj,1)
+
+                     IF( zztmp > zmaxdzT(ji,jj) ) THEN                        
+                         zmaxdzT(ji,jj) = zztmp   
+                         hth    (ji,jj) = zzdep                ! max and depth of dT/dz
+                     ENDIF
+               
+                     IF( nla10 > 1 ) THEN 
+                        zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                       ! delta rho(1)
+                        IF( zztmp > zrho3 )   zrho0_3(ji,jj) = zzdep                ! > 0.03
+                        IF( zztmp > zrho1 )   zrho0_1(ji,jj) = zzdep                ! > 0.01
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+         
+            CALL iom_put( 'mlddzt', hth )            ! depth of the thermocline
+            IF( nla10 > 1 ) THEN 
+               CALL iom_put( 'mldr0_3', zrho0_3 )   ! MLD delta rho(surf) = 0.03
+               CALL iom_put( 'mldr0_1', zrho0_1 )   ! MLD delta rho(surf) = 0.01
+            ENDIF
+            !
+         ENDIF
+         !
+         IF(  iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR.  &    
+            &  iom_use( 'pycndep' ) .OR. iom_use( 'tinv'     ) .OR. iom_use( 'depti'    )  ) THEN
+            ! ------------------------------------------------------------- !
+            ! MLD: abs( tn - tn(10m) ) = ztem2                              !
+            ! Top of thermocline: tn = tn(10m) - ztem2                      !
+            ! MLD: rho = rho10m + zrho3                                     !
+            ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
+            ! temperature inversion: max( 0, max of tn - tn(10m) )          !
+            ! depth of temperature inversion                                !
+            ! ------------------------------------------------------------- !
+            DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     !
+                     zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
+                     !
+                     zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
+                     IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
+                     IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
+                     zztmp = -zztmp                                          ! delta T(10m)
+                     IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
+                        ztinv(ji,jj) = zztmp   
+                        zdepinv (ji,jj) = zzdep   ! max value and depth
+                     ENDIF
+
+                     zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
+                     IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
+                     IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
+                     !
+                  END DO
+               END DO
+            END DO
+
+            CALL iom_put( 'mld_dt02', zabs2    )   ! MLD abs(delta t) - 0.2
+            CALL iom_put( 'topthdep', ztm2     )   ! T(10) - 0.2
+            CALL iom_put( 'mldr10_3', zrho10_3 )   ! MLD delta rho(10m) = 0.03
+            CALL iom_put( 'pycndep' , zpycn    )   ! MLD delta rho equi. delta T(10m) = 0.2
+            CALL iom_put( 'tinv'    , ztinv    )   ! max. temp. inv. (t10 ref) 
+            CALL iom_put( 'depti'   , zdepinv  )   ! depth of max. temp. inv. (t10 ref) 
+            !
+         ENDIF
+ 
+         ! ------------------------------- !
+         !  Depth of 20C/26C/28C isotherm  !
+         ! ------------------------------- !
+         IF( iom_use ('20d') ) THEN  ! depth of the 20 isotherm
+            ztem2 = 20.
+            CALL dia_hth_dep( ztem2, hd20 )  
+            CALL iom_put( '20d', hd20 )    
+         ENDIF
+         !
+         IF( iom_use ('26d') ) THEN  ! depth of the 26 isotherm
+            ztem2 = 26.
+            CALL dia_hth_dep( ztem2, hd26 )  
+            CALL iom_put( '26d', hd26 )    
+         ENDIF
+         !
+         IF( iom_use ('28d') ) THEN  ! depth of the 28 isotherm
+            ztem2 = 28.
+            CALL dia_hth_dep( ztem2, hd28 )  
+            CALL iom_put( '28d', hd28 )    
+         ENDIF
+        
+         ! ----------------------------- !
+         !  Heat content of first 300 m  !
+         ! ----------------------------- !
+         IF( iom_use ('hc300') ) THEN  
+            zzdep = 300.
+            CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc3 )
+            CALL iom_put( 'hc300', rau0_rcp * htc3 )  ! vertically integrated heat content (J/m2)
+         ENDIF
+         !
+         ! ----------------------------- !
+         !  Heat content of first 700 m  !
+         ! ----------------------------- !
+         IF( iom_use ('hc700') ) THEN  
+            zzdep = 700.
+            CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc7 )
+            CALL iom_put( 'hc700', rau0_rcp * htc7 )  ! vertically integrated heat content (J/m2)
+  
+         ENDIF
+         !
+         ! ----------------------------- !
+         !  Heat content of first 2000 m  !
+         ! ----------------------------- !
+         IF( iom_use ('hc2000') ) THEN  
+            zzdep = 2000.
+            CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc20 )
+            CALL iom_put( 'hc2000', rau0_rcp * htc20 )  ! vertically integrated heat content (J/m2)  
+         ENDIF
+         !
       ENDIF
+
+      !
+      IF( ln_timing )   CALL timing_stop('dia_hth')
+      !
+   END SUBROUTINE dia_hth
+
+   SUBROUTINE dia_hth_dep( ptem, pdept )
+      !
+      REAL(wp), INTENT(in) :: ptem
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept     
+      !
+      INTEGER  :: ji, jj, jk, iid
+      REAL(wp) :: zztmp, zzdep
+      INTEGER, DIMENSION(jpi,jpj) :: iktem
       
-      ! Preliminary computation
-      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( tmask(ji,jj,nla10) == 1. ) THEN
-               zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
-                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
-                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
-               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
-                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
-               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
-               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
-               zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
-               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
-            ELSE
-               zdelr(ji,jj) = 0._wp
-            ENDIF
-         END DO
-      END DO
-
-      ! ------------------------------------------------------------- !
-      ! thermocline depth: strongest vertical gradient of temperature !
-      ! turbocline depth (mixing layer depth): avt = zavt5            !
-      ! MLD: rho = rho(1) + zrho3                                     !
-      ! MLD: rho = rho(1) + zrho1                                     !
-      ! ------------------------------------------------------------- !
-      DO jk = jpkm1, 2, -1   ! loop from bottom to 2
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               zzdep = gdepw_n(ji,jj,jk)
-               zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
-               zzdep = zzdep * tmask(ji,jj,1)
-
-               IF( zztmp > zmaxdzT(ji,jj) ) THEN                        
-                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
-               ENDIF
-               
-               IF( nla10 > 1 ) THEN 
-                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
-                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
-                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
-               ENDIF
-
-            END DO
-         END DO
-      END DO
-      
-      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline
-      IF( nla10 > 1 ) THEN 
-         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03
-         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01
-      ENDIF
-
-      ! ------------------------------------------------------------- !
-      ! MLD: abs( tn - tn(10m) ) = ztem2                              !
-      ! Top of thermocline: tn = tn(10m) - ztem2                      !
-      ! MLD: rho = rho10m + zrho3                                     !
-      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
-      ! temperature inversion: max( 0, max of tn - tn(10m) )          !
-      ! depth of temperature inversion                                !
-      ! ------------------------------------------------------------- !
-      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
-               !
-               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
-               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
-               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
-               zztmp = -zztmp                                          ! delta T(10m)
-               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
-                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
-               ENDIF
-
-               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
-               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
-               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
-               !
-            END DO
-         END DO
-      END DO
-
-      CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2
-      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2
-      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
-      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
-      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref) 
-      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref) 
-
-
-      ! ----------------------------------- !
-      ! search deepest level above 20C/28C  !
-      ! ----------------------------------- !
-      ik20(:,:) = 1
-      ik28(:,:) = 1
+      ! --------------------------------------- !
+      ! search deepest level above ptem         !
+      ! --------------------------------------- !
+      iktem(:,:) = 1
       DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
          DO jj = 1, jpj
             DO ji = 1, jpi
                zztmp = tsn(ji,jj,jk,jp_tem)
-               IF( zztmp >= 20. )   ik20(ji,jj) = jk
-               IF( zztmp >= 28. )   ik28(ji,jj) = jk
+               IF( zztmp >= ptem )   iktem(ji,jj) = jk
             END DO
          END DO
       END DO
 
-      ! --------------------------- !
-      !  Depth of 20C/28C isotherm  !
-      ! --------------------------- !
+      ! ------------------------------- !
+      !  Depth of ptem isotherm         !
+      ! ------------------------------- !
       DO jj = 1, jpj
          DO ji = 1, jpi
             !
-            zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
+            zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean bottom
             !
-            iid = ik20(ji,jj)
+            iid = iktem(ji,jj)
             IF( iid /= 1 ) THEN 
-               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
+                zztmp =     gdept_n(ji,jj,iid  )   &                     ! linear interpolation
                   &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
                   &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
                   &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
-               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
+               pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
             ELSE 
-               hd20(ji,jj) = 0._wp
+               pdept(ji,jj) = 0._wp
             ENDIF
-            !
-            iid = ik28(ji,jj)
-            IF( iid /= 1 ) THEN 
-               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
-                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
-                  &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
-                  &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
-               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
-            ELSE 
-               hd28(ji,jj) = 0._wp
-            ENDIF
-
          END DO
       END DO
-      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm
-      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm
-
-      ! ----------------------------- !
-      !  Heat content of first 300 m  !
-      ! ----------------------------- !
-
-      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...)
-      ilevel   = 0
-      zthick_0 = 0._wp
-      DO jk = 1, jpkm1                      
-         zthick_0 = zthick_0 + e3t_1d(jk)
-         IF( zthick_0 < 300. )   ilevel = jk
-      END DO
+      !
+   END SUBROUTINE dia_hth_dep
+
+
+   SUBROUTINE dia_hth_htc( pdep, ptn, phtc )
+      !
+      REAL(wp), INTENT(in) :: pdep     ! depth over the heat content
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptn   
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc  
+      !
+      INTEGER  :: ji, jj, jk, ik
+      REAL(wp), DIMENSION(jpi,jpj) :: zthick
+      INTEGER , DIMENSION(jpi,jpj) :: ilevel
+
+
       ! surface boundary condition
-      IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)  
-      ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
+      
+      IF( .NOT. ln_linssh ) THEN   ;   zthick(:,:) = 0._wp       ;   phtc(:,:) = 0._wp                                   
+      ELSE                         ;   zthick(:,:) = sshn(:,:)   ;   phtc(:,:) = ptn(:,:,1) * sshn(:,:) * tmask(:,:,1)   
       ENDIF
-      ! integration down to ilevel
-      DO jk = 1, ilevel
-         zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk)
-         htc3  (:,:) = htc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
-      END DO
-      ! deepest layer
-      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
+      !
+      ilevel(:,:) = 1
+      DO jk = 2, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( ( gdept_n(ji,jj,jk) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
+                   ilevel(ji,jj) = jk
+                   zthick(ji,jj) = zthick(ji,jj) + e3t_n(ji,jj,jk)
+                   phtc  (ji,jj) = phtc  (ji,jj) + e3t_n(ji,jj,jk) * ptn(ji,jj,jk)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+      !
       DO jj = 1, jpj
          DO ji = 1, jpi
-            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  &
-               &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
+            ik = ilevel(ji,jj)
+            zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep
+            phtc(ji,jj)   = phtc(ji,jj) + ptn(ji,jj,ik+1) * MIN( e3t_n(ji,jj,ik+1), zthick(ji,jj) ) &
+                                                          * tmask(ji,jj,ik+1)
          END DO
-      END DO
-      ! from temperature to heat contain
-      zcoef = rau0 * rcp
-      htc3(:,:) = zcoef * htc3(:,:)
-      CALL iom_put( "hc300", htc3 )      ! first 300m heat content
-      !
-      IF( ln_timing )   CALL timing_stop('dia_hth')
-      !
-   END SUBROUTINE dia_hth
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Empty module
-   !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
-CONTAINS
-   SUBROUTINE dia_hth( kt )         ! Empty routine
-      IMPLICIT NONE
-      INTEGER, INTENT( in ) :: kt
-      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
-   END SUBROUTINE dia_hth
-#endif
+      ENDDO
+      !
+      !
+   END SUBROUTINE dia_hth_htc
 
    !!======================================================================
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diaptr.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diaptr.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/DIA/diaptr.F90	(revision 12109)
@@ -10,4 +10,5 @@
    !!            3.6  ! 2014-12  (C. Ethe) use of IOM
    !!            3.6  ! 2016-06  (T. Graham) Addition of diagnostics for CMIP6
+   !!            4.0  ! 2010-08  ( C. Ethe, J. Deshayes ) Improvment
    !!----------------------------------------------------------------------
 
@@ -42,24 +43,21 @@
 
    !                                  !!** namelist  namptr  **
-   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 )
-
-   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
-   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
-   INTEGER, PUBLIC ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T) 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
+
+   LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini)
+   INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc)
 
    REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
    REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
-   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg
-
-   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb
-   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
-   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   :: btm30   ! mask out Southern Ocean (=0 south of 30°S)
-
-   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d
-   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d
-
+   REAL(wp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rau0)
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk34 ! mask out Southern Ocean (=0 south of 34°S)
+
+   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)   :: p_fval1d
+   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d
+
+   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini)
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -71,8 +69,9 @@
 CONTAINS
 
-   SUBROUTINE dia_ptr( pvtr )
+   SUBROUTINE dia_ptr( kt, pvtr )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dia_ptr  ***
       !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index     
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
       !
@@ -80,51 +79,41 @@
       REAL(wp) ::   zsfc,zvfc               ! local scalar
       REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d    ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
-      REAL(wp), DIMENSION(jpj)     ::  vsum   ! 1D workspace
-      REAL(wp), DIMENSION(jpj,jpts)     ::  tssum   ! 1D workspace
- 
+      REAL(wp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
       !
       !overturning calculation
-      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse
-      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace
-
-
-      CHARACTER( len = 12 )  :: cl1
+      REAL(wp), DIMENSION(jpj,jpk,nptr) :: sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
+      REAL(wp), DIMENSION(jpj,jpk,nptr) :: zt_jk, zs_jk ! i-mean T and S, j-Stream-Function
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk,nptr)  :: z4d1, z4d2
+      REAL(wp), DIMENSION(jpi,jpj,nptr)      :: z3dtr ! i-mean T and S, j-Stream-Function
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('dia_ptr')
 
-      !
+      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init
+      !
+      IF( .NOT. l_diaptr )   RETURN
+
       IF( PRESENT( pvtr ) ) THEN
-         IF( iom_use("zomsfglo") ) THEN    ! effective MSF
-            z3d(1,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport
-            DO jk = 2, jpkm1 
-              z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)   ! effective j-Stream-Function (MSF)
+         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
+            DO jn = 1, nptr                                    ! by sub-basins
+               z4d1(1,:,:,jn) =  ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )  ! zonal cumulative effective transport excluding closed seas
+               DO jk = jpkm1, 1, -1 
+                  z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn)    ! effective j-Stream-Function (MSF)
+               END DO
+               DO ji = 1, jpi
+                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
+               ENDDO
             END DO
-            DO ji = 1, jpi
-               z3d(ji,:,:) = z3d(1,:,:)
-            ENDDO
-            cl1 = TRIM('zomsf'//clsubb(1) )
-            CALL iom_put( cl1, z3d * rc_sv )
-            DO jn = 2, nptr                                    ! by sub-basins
-               z3d(1,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) 
-               DO jk = 2, jpkm1 
-                  z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)    ! effective j-Stream-Function (MSF)
-               END DO
-               DO ji = 1, jpi
-                  z3d(ji,:,:) = z3d(1,:,:)
-               ENDDO
-               cl1 = TRIM('zomsf'//clsubb(jn) )
-               CALL iom_put( cl1, z3d * rc_sv )
-            END DO
-         ENDIF
-         IF( iom_use("sopstove") .OR. iom_use("sophtove") .OR. iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
+            CALL iom_put( 'zomsf', z4d1 * rc_sv )
+         ENDIF
+         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.   &
+            & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
             ! define fields multiplied by scalar
             zmask(:,:,:) = 0._wp
             zts(:,:,:,:) = 0._wp
-            zvn(:,:,:) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpjm1
@@ -134,114 +123,72 @@
                      zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
                      zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
-                     zvn(ji,jj,jk)        = vn(ji,jj,jk)         * zvfc
                   ENDDO
                ENDDO
              ENDDO
          ENDIF
-         IF( iom_use("sopstove") .OR. iom_use("sophtove") ) THEN
-             sjk(:,:,1) = ptr_sjk( zmask(:,:,:), btmsk(:,:,1) )
-             r1_sjk(:,:,1) = 0._wp
-             WHERE( sjk(:,:,1) /= 0._wp )   r1_sjk(:,:,1) = 1._wp / sjk(:,:,1)
-
-             ! i-mean T and S, j-Stream-Function, global
-             tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1)
-             sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1)
-             v_msf(:,:,1) = ptr_sjk( zvn(:,:,:) )
-
-             htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 )
-             str_ove(:,1) = SUM( v_msf(:,:,1)*sn_jk(:,:,1) ,2 )
-
-             z2d(1,:) = htr_ove(:,1) * rc_pwatt        !  (conversion in PW)
-             DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-             ENDDO
-             cl1 = 'sophtove'
-             CALL iom_put( TRIM(cl1), z2d )
-             z2d(1,:) = str_ove(:,1) * rc_ggram        !  (conversion in Gg)
-             DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-             ENDDO
-             cl1 = 'sopstove'
-             CALL iom_put( TRIM(cl1), z2d )
-             IF( ln_subbas ) THEN
-                DO jn = 2, nptr
-                    sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
-                    r1_sjk(:,:,jn) = 0._wp
-                    WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
-
-                    ! i-mean T and S, j-Stream-Function, basin
-                    tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
-                    sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
-                    v_msf(:,:,jn) = ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) ) 
-                    htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 )
-                    str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 )
-
-                    z2d(1,:) = htr_ove(:,jn) * rc_pwatt !  (conversion in PW)
-                    DO ji = 1, jpi
-                        z2d(ji,:) = z2d(1,:)
-                    ENDDO
-                    cl1 = TRIM('sophtove_'//clsubb(jn))
-                    CALL iom_put( cl1, z2d )
-                    z2d(1,:) = str_ove(:,jn) * rc_ggram        ! (conversion in Gg)
-                    DO ji = 1, jpi
-                        z2d(ji,:) = z2d(1,:)
-                    ENDDO
-                    cl1 = TRIM('sopstove_'//clsubb(jn))
-                    CALL iom_put( cl1, z2d )
-                END DO
-             ENDIF
-         ENDIF
-         IF( iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
-         ! Calculate barotropic heat and salt transport here 
-             sjk(:,1,1) = ptr_sj( zmask(:,:,:), btmsk(:,:,1) )
-             r1_sjk(:,1,1) = 0._wp
-             WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1)
-            
-            vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,1))
-            tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) )
-            tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) )
-            htr_btr(:,1) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,1)
-            str_btr(:,1) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,1)
-            z2d(1,:) = htr_btr(:,1) * rc_pwatt        !  (conversion in PW)
-            DO ji = 2, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sophtbtr'
-            CALL iom_put( TRIM(cl1), z2d )
-            z2d(1,:) = str_btr(:,1) * rc_ggram        !  (conversion in Gg)
-            DO ji = 2, jpi
-              z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sopstbtr'
-            CALL iom_put( TRIM(cl1), z2d )
-            IF( ln_subbas ) THEN
-                DO jn = 2, nptr
-                    sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) )
-                    r1_sjk(:,1,jn) = 0._wp
-                    WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
-                    vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,jn))
-                    tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
-                    tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
-                    htr_btr(:,jn) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,jn)
-                    str_btr(:,jn) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,jn)
-                    z2d(1,:) = htr_btr(:,jn) * rc_pwatt !  (conversion in PW)
-                    DO ji = 1, jpi
-                        z2d(ji,:) = z2d(1,:)
-                    ENDDO
-                    cl1 = TRIM('sophtbtr_'//clsubb(jn))
-                    CALL iom_put( cl1, z2d )
-                    z2d(1,:) = str_btr(:,jn) * rc_ggram        ! (conversion in Gg)
-                    DO ji = 1, jpi
-                        z2d(ji,:) = z2d(1,:)
-                    ENDDO
-                    cl1 = TRIM('sopstbtr_'//clsubb(jn))
-                    CALL iom_put( cl1, z2d )
-               ENDDO
-            ENDIF !ln_subbas
-         ENDIF !iom_use("sopstbtr....)
+         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
+            DO jn = 1, nptr
+               sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
+               r1_sjk(:,:,jn) = 0._wp
+               WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
+               ! i-mean T and S, j-Stream-Function, basin
+               zt_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
+               zs_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
+               v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) ) 
+               hstr_ove(:,jp_tem,jn) = SUM( v_msf(:,:,jn)*zt_jk(:,:,jn), 2 )
+               hstr_ove(:,jp_sal,jn) = SUM( v_msf(:,:,jn)*zs_jk(:,:,jn), 2 )
+               !
+            ENDDO
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sophtove', z3dtr )
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sopstove', z3dtr )
+         ENDIF
+
+         IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
+            ! Calculate barotropic heat and salt transport here 
+            DO jn = 1, nptr
+               sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) )
+               r1_sjk(:,1,jn) = 0._wp
+               WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
+               !
+               zvsum(:) = ptr_sj( pvtr(:,:,:), btmsk34(:,:,jn) )
+               ztsum(:) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
+               zssum(:) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
+               hstr_btr(:,jp_tem,jn) = zvsum(:) * ztsum(:) * r1_sjk(:,1,jn)
+               hstr_btr(:,jp_sal,jn) = zvsum(:) * zssum(:) * r1_sjk(:,1,jn)
+               !
+            ENDDO
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sophtbtr', z3dtr )
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sopstbtr', z3dtr )
+         ENDIF 
          !
       ELSE
          !
-         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface 
+         zmask(:,:,:) = 0._wp
+         zts(:,:,:,:) = 0._wp
+         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface 
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -254,119 +201,119 @@
                END DO
             END DO
+            !
             DO jn = 1, nptr
                zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
-               cl1 = TRIM('zosrf'//clsubb(jn) )
-               CALL iom_put( cl1, zmask )
-               !
-               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
-                  &            / MAX( zmask(1,:,:), 10.e-15 )
-               DO ji = 1, jpi
-                  z3d(ji,:,:) = z3d(1,:,:)
-               ENDDO
-               cl1 = TRIM('zotem'//clsubb(jn) )
-               CALL iom_put( cl1, z3d )
-               !
-               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
-                  &            / MAX( zmask(1,:,:), 10.e-15 )
-               DO ji = 1, jpi
-                  z3d(ji,:,:) = z3d(1,:,:)
-               ENDDO
-               cl1 = TRIM('zosal'//clsubb(jn) )
-               CALL iom_put( cl1, z3d )
-            END DO
+               z4d1(:,:,:,jn) = zmask(:,:,:)
+            ENDDO
+            CALL iom_put( 'zosrf', z4d1 )
+            !
+            DO jn = 1, nptr
+               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
+                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 )
+               DO ji = 1, jpi
+                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'zotem', z4d2 )
+            !
+            DO jn = 1, nptr
+               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
+                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 )
+               DO ji = 1, jpi
+                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'zosal', z4d2 )
+            !
          ENDIF
          !
          !                                ! Advective and diffusive heat and salt transport
-         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN   
-            z2d(1,:) = htr_adv(:,1) * rc_pwatt        !  (conversion in PW)
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sophtadv'                 
-            CALL iom_put( TRIM(cl1), z2d )
-            z2d(1,:) = str_adv(:,1) * rc_ggram        ! (conversion in Gg)
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sopstadv'
-            CALL iom_put( TRIM(cl1), z2d )
-            IF( ln_subbas ) THEN
-              DO jn=2,nptr
-               z2d(1,:) = htr_adv(:,jn) * rc_pwatt        !  (conversion in PW)
-               DO ji = 1, jpi
-                 z2d(ji,:) = z2d(1,:)
-               ENDDO
-               cl1 = TRIM('sophtadv_'//clsubb(jn))                 
-               CALL iom_put( cl1, z2d )
-               z2d(1,:) = str_adv(:,jn) * rc_ggram        ! (conversion in Gg)
-               DO ji = 1, jpi
-                  z2d(ji,:) = z2d(1,:)
-               ENDDO
-               cl1 = TRIM('sopstadv_'//clsubb(jn))                 
-               CALL iom_put( cl1, z2d )              
-              ENDDO
-            ENDIF
-         ENDIF
-         !
-         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN   
-            z2d(1,:) = htr_ldf(:,1) * rc_pwatt        !  (conversion in PW) 
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sophtldf'
-            CALL iom_put( TRIM(cl1), z2d )
-            z2d(1,:) = str_ldf(:,1) * rc_ggram        !  (conversion in Gg)
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sopstldf'
-            CALL iom_put( TRIM(cl1), z2d )
-            IF( ln_subbas ) THEN
-              DO jn=2,nptr
-               z2d(1,:) = htr_ldf(:,jn) * rc_pwatt        !  (conversion in PW)
-               DO ji = 1, jpi
-                 z2d(ji,:) = z2d(1,:)
-               ENDDO
-               cl1 = TRIM('sophtldf_'//clsubb(jn))                 
-               CALL iom_put( cl1, z2d )
-               z2d(1,:) = str_ldf(:,jn) * rc_ggram        ! (conversion in Gg)
-               DO ji = 1, jpi
-                  z2d(ji,:) = z2d(1,:)
-               ENDDO
-               cl1 = TRIM('sopstldf_'//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) 
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sophteiv'
-            CALL iom_put( TRIM(cl1), z2d )
-            z2d(1,:) = str_eiv(:,1) * rc_ggram        !  (conversion in Gg)
-            DO ji = 1, jpi
-               z2d(ji,:) = z2d(1,:)
-            ENDDO
-            cl1 = 'sopsteiv'
-            CALL iom_put( TRIM(cl1), z2d )
-            IF( ln_subbas ) THEN
-               DO jn=2,nptr
-                  z2d(1,:) = htr_eiv(:,jn) * rc_pwatt        !  (conversion in PW)
+         IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN  
+            ! 
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sophtadv', z3dtr )
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sopstadv', z3dtr )
+         ENDIF
+         !
+         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN  
+            ! 
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sophtldf', z3dtr )
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sopstldf', z3dtr )
+         ENDIF
+         !
+         IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN  
+            ! 
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sophteiv', z3dtr )
+            DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sopsteiv', z3dtr )
+         ENDIF
+         !
+         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
+            zts(:,:,:,:) = 0._wp
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1
                   DO ji = 1, jpi
-                     z2d(ji,:) = z2d(1,:)
+                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
+                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
+                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
                   ENDDO
-                  cl1 = TRIM('sophteiv_'//clsubb(jn))                 
-                  CALL iom_put( cl1, z2d )
-                  z2d(1,:) = str_eiv(:,jn) * rc_ggram        ! (conversion in Gg)
-                  DO ji = 1, jpi
-                     z2d(ji,:) = z2d(1,:)
-                  ENDDO
-                  cl1 = TRIM('sopsteiv_'//clsubb(jn)) 
-                  CALL iom_put( cl1, z2d )              
-               ENDDO
-            ENDIF
+               ENDDO
+             ENDDO
+             CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) )
+             CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) )
+             DO jn = 1, nptr
+                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+                DO ji = 1, jpi
+                   z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+                ENDDO
+             ENDDO
+             CALL iom_put( 'sophtvtr', z3dtr )
+             DO jn = 1, nptr
+               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = 1, jpi
+                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               ENDDO
+            ENDDO
+            CALL iom_put( 'sopstvtr', z3dtr )
+         ENDIF
+         !
+         IF( iom_use( 'uocetr_vsum_cumul' ) ) THEN
+            CALL iom_get_var(  'uocetr_vsum_op', z2d ) ! get uocetr_vsum_op from xml
+            z2d(:,:) = ptr_ci_2d( z2d(:,:) )  
+            CALL iom_put( 'uocetr_vsum_cumul', z2d )
          ENDIF
          !
@@ -384,20 +331,18 @@
       !! ** Purpose :   Initialization, namelist read
       !!----------------------------------------------------------------------
-      INTEGER ::  jn           ! local integers
-      INTEGER ::  inum, ierr   ! local integers
-      INTEGER ::  ios          ! Local integer output status for namelist read
-      !!
-      NAMELIST/namptr/ ln_diaptr, ln_subbas
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
-      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist' )
-
-      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
-      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist' )
-      IF(lwm) WRITE ( numond, namptr )
-
+      INTEGER ::  inum, jn           ! local integers
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) :: zmsk
+      !!----------------------------------------------------------------------
+
+      l_diaptr = .FALSE.
+      IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  &
+         &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  &
+         &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  &
+         &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  & 
+         &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  &
+         &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE.
+
+ 
       IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
@@ -405,52 +350,47 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namptr : set ptr parameters'
-         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
-         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas
+         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr
       ENDIF
 
-      IF( ln_diaptr ) THEN  
-         !
-         IF( ln_subbas ) THEN 
-            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific
-            ALLOCATE( clsubb(nptr) )
-            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc'
-         ELSE               
-            nptr = 1       ! Global only
-            ALLOCATE( clsubb(nptr) )
-            clsubb(1) = 'glo' 
-         ENDIF
-
-         !                                      ! allocate dia_ptr arrays
+      IF( l_diaptr ) THEN  
+         !
          IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
 
          rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt
+         rc_ggram = rc_ggram * rau0              ! conversion from m3/s to Gg/s
 
          IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
 
-         IF( ln_subbas ) THEN                ! load sub-basin mask
-            CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
-            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
-            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
-            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
-            CALL iom_close( inum )
-            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
-            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
-            ELSE WHERE                     ;   btm30(:,:) = ssmask(:,:)
-            END WHERE
-         ENDIF
-   
-         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
-      
-         DO jn = 1, nptr
+         btmsk(:,:,1) = tmask_i(:,:)                 
+         CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
+         CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
+         CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
+         CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
+         CALL iom_close( inum )
+         btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
+         DO jn = 2, nptr
             btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
          END DO
+         ! JD : modification so that overturning streamfunction is available in Atlantic at 34S to compare with observations
+         WHERE( gphit(:,:)*tmask_i(:,:) < -34._wp)
+           zmsk(:,:) = 0._wp      ! mask out Southern Ocean
+         ELSE WHERE                  
+           zmsk(:,:) = ssmask(:,:)
+         END WHERE
+         btmsk34(:,:,1) = btmsk(:,:,1)                 
+         DO jn = 2, nptr
+            btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)               ! interior domain only
+         ENDDO
 
          ! Initialise arrays to zero because diatpr is called before they are first calculated
          ! Note that this means diagnostics will not be exactly correct when model run is restarted.
-         htr_adv(:,:) = 0._wp  ;  str_adv(:,:) =  0._wp 
-         htr_ldf(:,:) = 0._wp  ;  str_ldf(:,:) =  0._wp 
-         htr_eiv(:,:) = 0._wp  ;  str_eiv(:,:) =  0._wp 
-         htr_ove(:,:) = 0._wp  ;   str_ove(:,:) =  0._wp
-         htr_btr(:,:) = 0._wp  ;   str_btr(:,:) =  0._wp
+         hstr_adv(:,:,:) = 0._wp           
+         hstr_ldf(:,:,:) = 0._wp           
+         hstr_eiv(:,:,:) = 0._wp           
+         hstr_ove(:,:,:) = 0._wp           
+         hstr_btr(:,:,:) = 0._wp           !
+         hstr_vtr(:,:,:) = 0._wp           !
+         !
+         ll_init = .FALSE.
          !
       ENDIF 
@@ -471,57 +411,57 @@
       INTEGER                                        :: jn    !
 
+      !
       IF( cptr == 'adv' ) THEN
-         IF( ktra == jp_tem )  htr_adv(:,1) = ptr_sj( pva(:,:,:) )
-         IF( ktra == jp_sal )  str_adv(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_tem )  THEN
+             DO jn = 1, nptr
+                hstr_adv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
+         IF( ktra == jp_sal )  THEN
+             DO jn = 1, nptr
+                hstr_adv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
       ENDIF
+      !
       IF( cptr == 'ldf' ) THEN
-         IF( ktra == jp_tem )  htr_ldf(:,1) = ptr_sj( pva(:,:,:) )
-         IF( ktra == jp_sal )  str_ldf(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_tem )  THEN
+             DO jn = 1, nptr
+                hstr_ldf(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
+         IF( ktra == jp_sal )  THEN
+             DO jn = 1, nptr
+                hstr_ldf(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
       ENDIF
+      !
       IF( cptr == 'eiv' ) THEN
-         IF( ktra == jp_tem )  htr_eiv(:,1) = ptr_sj( pva(:,:,:) )
-         IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_tem )  THEN
+             DO jn = 1, nptr
+                hstr_eiv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
+         IF( ktra == jp_sal )  THEN
+             DO jn = 1, nptr
+                hstr_eiv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
       ENDIF
       !
-      IF( ln_subbas ) THEN
-         !
-         IF( cptr == 'adv' ) THEN
-             IF( ktra == jp_tem ) THEN 
-                DO jn = 2, nptr
-                   htr_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                END DO
-             ENDIF
-             IF( ktra == jp_sal ) THEN 
-                DO jn = 2, nptr
-                   str_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                END DO
-             ENDIF
-         ENDIF
-         IF( cptr == 'ldf' ) THEN
-             IF( ktra == jp_tem ) THEN 
-                DO jn = 2, nptr
-                    htr_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                 END DO
-             ENDIF
-             IF( ktra == jp_sal ) THEN 
-                DO jn = 2, nptr
-                   str_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                END DO
-             ENDIF
-         ENDIF
-         IF( cptr == 'eiv' ) THEN
-             IF( ktra == jp_tem ) THEN 
-                DO jn = 2, nptr
-                    htr_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                 END DO
-             ENDIF
-             IF( ktra == jp_sal ) THEN 
-                DO jn = 2, nptr
-                   str_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
-                END DO
-             ENDIF
-         ENDIF
-         !
+      IF( cptr == 'vtr' ) THEN
+         IF( ktra == jp_tem )  THEN
+             DO jn = 1, nptr
+                hstr_vtr(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
+         IF( ktra == jp_sal )  THEN
+             DO jn = 1, nptr
+                hstr_vtr(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+             ENDDO
+         ENDIF
       ENDIF
+      !
    END SUBROUTINE dia_ptr_hst
 
@@ -536,18 +476,15 @@
       ierr(:) = 0
       !
-      ALLOCATE( btmsk(jpi,jpj,nptr) ,              &
-         &      htr_adv(jpj,nptr) , str_adv(jpj,nptr) ,   &
-         &      htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) ,   &
-         &      htr_ove(jpj,nptr) , str_ove(jpj,nptr) ,   &
-         &      htr_btr(jpj,nptr) , str_btr(jpj,nptr) ,   &
-         &      htr_ldf(jpj,nptr) , str_ldf(jpj,nptr) , STAT=ierr(1)  )
-         !
-      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
-      !
-      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  )
-
-         !
-      dia_ptr_alloc = MAXVAL( ierr )
-      CALL mpp_sum( 'diaptr', dia_ptr_alloc )
+      IF( .NOT. ALLOCATED( btmsk ) ) THEN
+         ALLOCATE( btmsk(jpi,jpj,nptr)    , btmsk34(jpi,jpj,nptr),   &
+            &      hstr_adv(jpj,jpts,nptr), hstr_eiv(jpj,jpts,nptr), &
+            &      hstr_ove(jpj,jpts,nptr), hstr_btr(jpj,jpts,nptr), &
+            &      hstr_ldf(jpj,jpts,nptr), hstr_vtr(jpj,jpts,nptr), STAT=ierr(1)  )
+            !
+         ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
+         !
+         dia_ptr_alloc = MAXVAL( ierr )
+         CALL mpp_sum( 'diaptr', dia_ptr_alloc )
+      ENDIF
       !
    END FUNCTION dia_ptr_alloc
@@ -565,6 +502,6 @@
       !! ** Action  : - p_fval: i-k-mean poleward flux of pva
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pva   ! mask flux array at V-point
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)      ::   pmsk   ! Optional 2D basin mask
       !
       INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
@@ -577,21 +514,11 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      IF( PRESENT( pmsk ) ) THEN 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! Vector opt.
-                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj)
-               END DO
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! Vector opt.
+               p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
             END DO
          END DO
-      ELSE
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! Vector opt.
-                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
-               END DO
-            END DO
-         END DO
-      ENDIF
+      END DO
 #if defined key_mpp_mpi
       CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl)
@@ -612,6 +539,6 @@
       !! ** Action  : - p_fval: i-k-mean poleward flux of pva
       !!----------------------------------------------------------------------
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
       !
       INTEGER                  ::   ji,jj       ! dummy loop arguments
@@ -624,17 +551,9 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      IF( PRESENT( pmsk ) ) THEN 
-         DO jj = 2, jpjm1
-            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
-               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj)
-            END DO
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! Vector opt.
+            p_fval(jj) = p_fval(jj) + pva(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
          END DO
-      ELSE
-         DO jj = 2, jpjm1
-            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
-               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj)
-            END DO
-         END DO
-      ENDIF
+      END DO
 #if defined key_mpp_mpi
       CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl )
@@ -643,4 +562,35 @@
    END FUNCTION ptr_sj_2d
 
+   FUNCTION ptr_ci_2d( pva )   RESULT ( p_fval )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ptr_ci_2d  ***
+      !!
+      !! ** Purpose :   "meridional" cumulated sum computation of a j-flux array
+      !!
+      !! ** Method  : - j cumulated sum of pva using the interior 2D vmask (umask_i).
+      !!
+      !! ** Action  : - p_fval: j-cumulated sum of pva
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)  ::   pva   ! mask flux array at V-point
+      !
+      INTEGER                  ::   ji,jj,jc       ! dummy loop arguments
+      INTEGER                  ::   ijpj        ! ??? 
+      REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
+      !!--------------------------------------------------------------------
+      ! 
+      ijpj = jpj  ! ???
+      p_fval(:,:) = 0._wp
+      DO jc = 1, jpnj ! looping over all processors in j axis
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! Vector opt.
+               p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
+            END DO
+         END DO
+         CALL lbc_lnk( 'diaptr', p_fval, 'U', -1. )
+      END DO
+      ! 
+   END FUNCTION ptr_ci_2d
+
+
 
    FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
@@ -656,6 +606,6 @@
       !!
       IMPLICIT none
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)     ::   pmsk   ! Optional 2D basin mask
       !!
       INTEGER                           :: ji, jj, jk ! dummy loop arguments
@@ -673,22 +623,11 @@
       p_fval(:,:) = 0._wp
       !
-      IF( PRESENT( pmsk ) ) THEN 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
-               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
-                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj)
-               END DO
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! Vector opt.
+               p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
             END DO
          END DO
-      ELSE 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
-                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
-      END IF
+      END DO
       !
 #if defined key_mpp_mpi
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/IOM/iom.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/IOM/iom.F90	(revision 12109)
@@ -56,5 +56,5 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_iomput = .FALSE.       !: iom_put flag
 #endif
-   PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get
+   PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_get_var
    PUBLIC iom_chkatt, iom_getatt, iom_putatt, iom_getszuld, iom_rstput, iom_delay_rst, iom_put
    PUBLIC iom_use, iom_context_finalize, iom_miss_val
@@ -62,5 +62,5 @@
    PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d
    PRIVATE iom_g0d, iom_g1d, iom_g2d, iom_g3d, iom_get_123d
-   PRIVATE iom_p1d, iom_p2d, iom_p3d
+   PRIVATE iom_p1d, iom_p2d, iom_p3d, iom_p4d
 #if defined key_iomput
    PRIVATE iom_set_domain_attr, iom_set_axis_attr, iom_set_field_attr, iom_set_file_attr, iom_get_file_attr, iom_set_grid_attr
@@ -83,5 +83,5 @@
    END INTERFACE
    INTERFACE iom_put
-      MODULE PROCEDURE iom_p0d, iom_p1d, iom_p2d, iom_p3d
+      MODULE PROCEDURE iom_p0d, iom_p1d, iom_p2d, iom_p3d, iom_p4d
    END INTERFACE iom_put
   
@@ -108,4 +108,5 @@
       TYPE(xios_date)     :: start_date
       CHARACTER(len=lc) :: clname
+      INTEGER             :: irefyear, irefmonth, irefday
       INTEGER           :: ji, jkmin
       LOGICAL :: llrst_context              ! is context related to restart
@@ -139,10 +140,14 @@
 
       ! Calendar type is now defined in xml file 
+      IF (.NOT.(xios_getvar('ref_year' ,irefyear ))) irefyear  = 1900
+      IF (.NOT.(xios_getvar('ref_month',irefmonth))) irefmonth = 01
+      IF (.NOT.(xios_getvar('ref_day'  ,irefday  ))) irefday   = 01
+
       SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
-      CASE ( 1)   ; CALL xios_define_calendar( TYPE = "Gregorian", time_origin = xios_date(1900,01,01,00,00,00), &
+      CASE ( 1)   ; CALL xios_define_calendar( TYPE = "Gregorian", time_origin = xios_date(irefyear,irefmonth,irefday,00,00,00), &
           &                                    start_date = xios_date(nyear,nmonth,nday,0,0,0) )
-      CASE ( 0)   ; CALL xios_define_calendar( TYPE = "NoLeap"   , time_origin = xios_date(1900,01,01,00,00,00), &
+      CASE ( 0)   ; CALL xios_define_calendar( TYPE = "NoLeap"   , time_origin = xios_date(irefyear,irefmonth,irefday,00,00,00), &
           &                                    start_date = xios_date(nyear,nmonth,nday,0,0,0) )
-      CASE (30)   ; CALL xios_define_calendar( TYPE = "D360"     , time_origin = xios_date(1900,01,01,00,00,00), &
+      CASE (30)   ; CALL xios_define_calendar( TYPE = "D360"     , time_origin = xios_date(irefyear,irefmonth,irefday,00,00,00), &
           &                                    start_date = xios_date(nyear,nmonth,nday,0,0,0) )
       END SELECT
@@ -223,5 +228,7 @@
           CALL iom_set_axis_attr( "icbcla", class_num )
           CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )   ! strange syntaxe and idea...
+          CALL iom_set_axis_attr( "iax_26C", (/ REAL(26,wp) /) )   ! strange syntaxe and idea...
           CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )   ! strange syntaxe and idea...
+          CALL iom_set_axis_attr( "basin"  , (/ (REAL(ji,wp), ji=1,5) /) )
       ENDIF
       !
@@ -1340,4 +1347,17 @@
    END SUBROUTINE iom_get_123d
 
+   SUBROUTINE iom_get_var( cdname, z2d)
+      CHARACTER(LEN=*), INTENT(in ) ::   cdname
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2d 
+#if defined key_iomput
+      IF( xios_field_is_active( cdname, at_current_timestep_arg = .TRUE. ) ) THEN
+         z2d(:,:) = 0._wp
+         CALL xios_recv_field( cdname, z2d)
+      ENDIF
+#else
+      IF( .FALSE. )   WRITE(numout,*) cdname, z2d ! useless test to avoid compilation warnings
+#endif
+   END SUBROUTINE iom_get_var
+
 
    FUNCTION iom_getszuld ( kiomid )  
@@ -1709,4 +1729,15 @@
    END SUBROUTINE iom_p3d
 
+   SUBROUTINE iom_p4d( cdname, pfield4d )
+      CHARACTER(LEN=*)                , INTENT(in) ::   cdname
+      REAL(wp),       DIMENSION(:,:,:,:), INTENT(in) ::   pfield4d
+#if defined key_iomput
+      CALL xios_send_field(cdname, pfield4d)
+#else
+      IF( .FALSE. )   WRITE(numout,*) cdname, pfield4d   ! useless test to avoid compilation warnings
+#endif
+   END SUBROUTINE iom_p4d
+
+
 #if defined key_iomput
    !!----------------------------------------------------------------------
@@ -2051,12 +2082,11 @@
       ALLOCATE( zlon(ni*nj) )       ;       zlon(:) = 0._wp
       !
-      CALL dom_ngb( -168.53, 65.03, ix, iy, 'T' ) !  i-line that passes through Bering Strait: Reference latitude (used in plots)
-!      CALL dom_ngb( 180., 90., ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
+!      CALL dom_ngb( -168.53, 65.03, ix, iy, 'T' ) !  i-line that passes through Bering Strait: Reference latitude (used in plots)
+      CALL dom_ngb( 180., 90., ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
       CALL iom_set_domain_attr("gznl", ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj)
       CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj)
       CALL iom_set_domain_attr("gznl", lonvalue = zlon,   &
          &                             latvalue = RESHAPE(plat(nldi:nlei, nldj:nlej),(/ ni*nj /)))  
-      CALL iom_set_zoom_domain_attr("znl_T", ibegin=ix-1, jbegin=0, ni=1, nj=jpjglo)
-      CALL iom_set_zoom_domain_attr("znl_W", ibegin=ix-1, jbegin=0, ni=1, nj=jpjglo)
+      CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=jpjglo)
       !
       CALL iom_update_file_name('ptr')
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/LDF/ldftra.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/LDF/ldftra.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/LDF/ldftra.F90	(revision 12109)
@@ -851,4 +851,19 @@
       CALL iom_put( "woce_eiv", zw3d )
       !
+      IF( iom_use('weiv_masstr') ) THEN   ! vertical mass transport & its square value
+         zw2d(:,:) = rau0 * e1e2t(:,:)
+         DO jk = 1, jpk
+            zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:)
+         END DO
+         CALL iom_put( "weiv_masstr" , zw3d )  
+      ENDIF
+      !
+      IF( iom_use('ueiv_masstr') ) THEN
+         zw3d(:,:,:) = 0.e0
+         DO jk = 1, jpkm1
+            zw3d(:,:,jk) = rau0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 
+         END DO
+         CALL iom_put( "ueiv_masstr", zw3d )                  ! mass transport in i-direction
+      ENDIF
       !
       zztmp = 0.5_wp * rau0 * rcp 
@@ -870,4 +885,13 @@
         CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction
       ENDIF
+      !
+      IF( iom_use('veiv_masstr') ) THEN
+         zw3d(:,:,:) = 0.e0
+         DO jk = 1, jpkm1
+            zw3d(:,:,jk) = rau0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 
+         END DO
+         CALL iom_put( "veiv_masstr", zw3d )                  ! mass transport in i-direction
+      ENDIF
+      !
       zw2d(:,:)   = 0._wp 
       zw3d(:,:,:) = 0._wp 
@@ -885,5 +909,5 @@
       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 )
+      IF( iom_use( 'sophteiv' ) )   CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
       !
       zztmp = 0.5_wp * 0.5
@@ -920,5 +944,5 @@
       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 )
+      IF( iom_use( 'sopsteiv' ) ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
       !
       !
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcblk.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcblk.F90	(revision 12109)
@@ -801,4 +801,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
       !!---------------------------------------------------------------------
       !
@@ -913,4 +914,24 @@
          qtr_ice_top(:,:,:) = 0._wp 
       END WHERE
+      !
+
+      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
+         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
+         IF( iom_use('evap_ao_cea'  ) )  CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average)
+         IF( iom_use('hflx_evap_cea') )  CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average)
+      ENDIF
+      IF( iom_use('hflx_rain_cea') ) THEN
+         ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
+         IF( iom_use('hflx_rain_cea') )  CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average)
+      ENDIF
+      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
+          WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
+          ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)    
+          ENDWHERE
+          ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 
+          IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average)
+          IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
+          IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
+      ENDIF
       !
       IF(ln_ctl) THEN
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbccpl.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbccpl.F90	(revision 12109)
@@ -1774,4 +1774,5 @@
       IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
       IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
+      IF( iom_use('rain_ao_cea') )  CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average)
       IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average)
       IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
@@ -1897,13 +1898,15 @@
 #endif
       ! outputs
-      IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving
-      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting
-      IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
-      IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) &
+      IF( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving
+      IF( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting
+      IF( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
+      IF( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) &
            &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
-      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average)
-      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
+      IF( iom_use('hflx_prec_cea')    ) CALL iom_put('hflx_prec_cea'   ,  sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) +  &                    ! heat flux from all precip (cell avg)
+         &                                                               ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
+      IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average)
+      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
            &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean)
-      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
+      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
            &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice)
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
@@ -2301,5 +2304,8 @@
       !                                                      !  CO2 flux from PISCES     ! 
       !                                                      ! ------------------------- !
-      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
+      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   THEN 
+         ztmp1(:,:) = oce_co2(:,:) * 1000.  ! conversion in molC/m2/s
+         CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info )
+      ENDIF
       !
       !                                                      ! ------------------------- !
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcmod.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcmod.F90	(revision 12109)
@@ -244,4 +244,12 @@
          fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp
       END IF
+      !
+      IF( sbc_ssr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_ssr arrays' )
+      IF( .NOT.ln_ssr ) THEN               !* Initialize qrp and erp if no restoring 
+         qrp(:,:) = 0._wp
+         erp(:,:) = 0._wp
+      ENDIF
+      !
+
       IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
          IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
@@ -552,4 +560,6 @@
          CALL iom_put( "taum"  , taum       )                   ! wind stress module
          CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice
+         CALL iom_put( "qrp", qrp )                             ! heat flux damping
+         CALL iom_put( "erp", erp )                             ! freshwater flux damping
       ENDIF
       !
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcrnf.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcrnf.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcrnf.F90	(revision 12109)
@@ -43,8 +43,10 @@
    REAL(wp)                   ::      rn_dep_max        !: depth over which runoffs is spread       (ln_rnf_depth_ini =T)
    INTEGER                    ::      nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
+   LOGICAL                    ::   ln_rnf_icb        !: iceberg flux is specified in a file
    LOGICAL                    ::   ln_rnf_tem        !: temperature river runoffs attribute specified in a file
    LOGICAL           , PUBLIC ::   ln_rnf_sal        !: salinity    river runoffs attribute specified in a file
    TYPE(FLD_N)       , PUBLIC ::   sn_rnf            !: information about the runoff file to be read
    TYPE(FLD_N)                ::   sn_cnf            !: information about the runoff mouth file to be read
+   TYPE(FLD_N)                ::   sn_i_rnf        !: information about the iceberg flux file to be read
    TYPE(FLD_N)                ::   sn_s_rnf          !: information about the salinities of runoff file to be read
    TYPE(FLD_N)                ::   sn_t_rnf          !: information about the temperatures of runoff file to be read
@@ -65,4 +67,5 @@
 
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read)
+   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_i_rnf     ! structure: iceberg flux (file information, fields read)
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)  
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)  
@@ -112,5 +115,9 @@
       !                                            !-------------------!
       !
-      IF( .NOT. l_rnfcpl )   CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt
+      !
+      IF( .NOT. l_rnfcpl )  THEN
+                            CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt ( runoffs + iceberg )
+         IF( ln_rnf_icb )   CALL fld_read ( kt, nn_fsbc, sf_i_rnf )    ! idem for iceberg flux if required
+      ENDIF
       IF(   ln_rnf_tem   )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required
       IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
@@ -118,5 +125,12 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
          !
-         IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)       ! updated runoff value at time step kt
+         IF( .NOT. l_rnfcpl ) THEN
+             rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt
+             IF( ln_rnf_icb ) THEN
+                fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt
+                CALL iom_put( 'iceberg_cea'  , fwficb(:,:)  )         ! output iceberg flux
+                CALL iom_put( 'hflx_icb_cea' , fwficb(:,:) * rLfus )   ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
+             ENDIF
+         ENDIF
          !
          !                                                           ! set temperature & salinity content of runoffs
@@ -132,10 +146,10 @@
          ELSE                                                        ! use SST as runoffs temperature
             !CEOD River is fresh water so must at least be 0 unless we consider ice
-            rnf_tsc(:,:,jp_tem) = MAX(sst_m(:,:),0.0_wp) * rnf(:,:) * r1_rau0
+            rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rau0
          ENDIF
          !                                                           ! use runoffs salinity data
          IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
          !                                                           ! else use S=0 for runoffs (done one for all in the init)
-         IF( iom_use('runoffs') )        CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux
+                                         CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux
          IF( iom_use('hflx_rnf_cea') )   CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rau0 * rcp )   ! output runoff sensible heat (W/m2)
       ENDIF
@@ -242,6 +256,6 @@
       REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl    
       !!
-      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   &
-         &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &
+      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb,   &
+         &                 sn_rnf, sn_cnf    , sn_i_rnf, sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &
          &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact,     &
          &                 ln_rnf_depth_ini  , rn_dep_max  , rn_rnf_max, nn_rnf_depth_file
@@ -299,4 +313,19 @@
          IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
          CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print )
+         !
+         IF( ln_rnf_icb ) THEN                      ! Create (if required) sf_i_rnf structure
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '          iceberg flux read in a file'
+            ALLOCATE( sf_i_rnf(1), STAT=ierror  )
+            IF( ierror > 0 ) THEN
+               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' )   ;   RETURN
+            ENDIF
+            ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1)   )
+            IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
+         ELSE
+            fwficb(:,:) = 0._wp
+         ENDIF
+
       ENDIF
       !
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcssr.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcssr.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/SBC/sbcssr.F90	(revision 12109)
@@ -30,7 +30,9 @@
    PUBLIC   sbc_ssr        ! routine called in sbcmod
    PUBLIC   sbc_ssr_init   ! routine called in sbcmod
+   PUBLIC   sbc_ssr_alloc  ! routine called in sbcmod
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   erp   !: evaporation damping   [kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qrp   !: heat flux damping        [w/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   coefice   !: under ice relaxation coefficient
 
    !                                   !!* Namelist namsbc_ssr *
@@ -41,4 +43,5 @@
    LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
    REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
+   INTEGER         ::   nn_icedmp       ! Control of restoring under ice
 
    REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange
@@ -97,5 +100,17 @@
                   END DO
                END DO
-               CALL iom_put( "qrp", qrp )                             ! heat flux damping
+            ENDIF
+            !
+            IF( nn_sssr /= 0 .AND. nn_icedmp /= 1 ) THEN
+              ! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_icedmp .ne. 1
+              ! n.b. coefice is initialised and fixed to 1._wp if nn_icedmp = 1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     SELECT CASE ( nn_icedmp )
+                       CASE ( 0 )    ;  coefice(ji,jj) = 1._wp - fr_i(ji,jj)              ! no/reduced damping under ice
+                       CASE  DEFAULT ;  coefice(ji,jj) = 1._wp +(nn_icedmp-1)*fr_i(ji,jj) ! reinforced damping (x nn_icedmp) under ice )
+                     END SELECT
+                  END DO
+               END DO
             ENDIF
             !
@@ -105,4 +120,5 @@
                   DO ji = 1, jpi
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                        &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
                      sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
@@ -110,5 +126,4 @@
                   END DO
                END DO
-               CALL iom_put( "erp", erp )                             ! freshwater flux damping
                !
             ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns)
@@ -118,4 +133,5 @@
                   DO ji = 1, jpi                            
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                        &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
                         &        / MAX(  sss_m(ji,jj), 1.e-20   ) * tmask(ji,jj,1)
@@ -126,5 +142,4 @@
                   END DO
                END DO
-               CALL iom_put( "erp", erp )                             ! freshwater flux damping
             ENDIF
             !
@@ -154,5 +169,6 @@
       CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
       TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
-      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd
+      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, &
+              & sn_sss, ln_sssr_bnd, rn_sssr_bnd, nn_icedmp
       INTEGER     ::  ios
       !!----------------------------------------------------------------------
@@ -182,9 +198,9 @@
          WRITE(numout,*) '         flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd
          WRITE(numout,*) '         ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day'
-      ENDIF
-      !
-      !                            !* Allocate erp and qrp array
-      ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), STAT=ierror )
-      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate erp and qrp array' )
+         WRITE(numout,*) '      Cntrl of surface restoration under ice nn_icedmp      = ', nn_icedmp
+         WRITE(numout,*) '          ( 0 = no restoration under ice)'
+         WRITE(numout,*) '          ( 1 = restoration everywhere  )'
+         WRITE(numout,*) '          (>1 = enhanced restoration under ice  )'
+      ENDIF
       !
       IF( nn_sstr == 1 ) THEN      !* set sf_sst structure & allocate arrays
@@ -216,4 +232,5 @@
       ENDIF
       !
+      coefice(:,:) = 1._wp         !  Initialise coefice to 1._wp ; will not need to be changed if nn_icedmp=1
       !                            !* Initialize qrp and erp if no restoring 
       IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
@@ -221,4 +238,18 @@
       !
    END SUBROUTINE sbc_ssr_init
+         
+   INTEGER FUNCTION sbc_ssr_alloc()
+      !!----------------------------------------------------------------------
+      !!               ***  FUNCTION sbc_ssr_alloc  ***
+      !!----------------------------------------------------------------------
+      sbc_ssr_alloc = 0       ! set to zero if no array to be allocated
+      IF( .NOT. ALLOCATED( erp ) ) THEN
+         ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), coefice(jpi,jpj), STAT= sbc_ssr_alloc )
+         !
+         IF( lk_mpp                  )   CALL mpp_sum ( 'sbcssr', sbc_ssr_alloc )
+         IF( sbc_ssr_alloc /= 0 )   CALL ctl_warn('sbc_ssr_alloc: failed to allocate arrays.')
+         !
+      ENDIF
+   END FUNCTION
       
    !!======================================================================
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv.F90	(revision 12109)
@@ -134,7 +134,8 @@
       !
 !!gm ???
-      IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF 
+      CALL dia_ptr( kt, zvn )                                    ! diagnose the effective MSF 
 !!gm ???
       !
+
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_cen.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_cen.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_cen.F90	(revision 12109)
@@ -61,5 +61,5 @@
       !! ** Action : - update pta  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
-      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
+      !!             - poleward advective heat and salt transport (l_diaptr=T)
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
@@ -89,8 +89,8 @@
       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. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    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.
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
       !
       !                    
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_fct.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_fct.F90	(revision 12109)
@@ -68,5 +68,5 @@
       !! ** Action : - update pta  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostics (l_trdtra=T)
-      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
+      !!             - poleward advective heat and salt transport (ln_diaptr=T)
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
@@ -101,7 +101,7 @@
       l_ptr = .FALSE.
       ll_zAimp = .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.  &
+      IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    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.
       !
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_mus.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_mus.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_mus.F90	(revision 12109)
@@ -68,5 +68,5 @@
       !! ** Action : - update pta  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
-      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
+      !!             - poleward advective heat and salt transport (ln_diaptr=T)
       !!
       !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
@@ -120,5 +120,5 @@
       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( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   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.
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_qck.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_qck.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_qck.F90	(revision 12109)
@@ -21,4 +21,5 @@
    USE trdtra          ! trends manager: tracers 
    USE diaptr          ! poleward transport diagnostics
+   USE iom
    !
    USE in_out_manager  ! I/O manager
@@ -80,5 +81,5 @@
       !! ** Action : - update pta  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
-      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
+      !!             - poleward advective heat and salt transport (ln_diaptr=T)
       !!
       !! ** Reference : Leonard (1979, 1991)
@@ -103,6 +104,6 @@
       l_trd = .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. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 
       !
       !
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_ubs.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_ubs.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traadv_ubs.F90	(revision 12109)
@@ -79,5 +79,5 @@
       !! ** Action : - update pta  with the now advective tracer trends
       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
-      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
+      !!             - poleward advective heat and salt transport (ln_diaptr=T)
       !!
       !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 
@@ -111,5 +111,5 @@
       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( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   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.
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/trabbc.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/trabbc.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/trabbc.F90	(revision 12109)
@@ -100,4 +100,6 @@
       ENDIF
       !
+      CALL iom_put ( "hfgeou" , rau0_rcp * qgh_trd0(:,:) )
+      !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_iso.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_iso.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_iso.F90	(revision 12109)
@@ -124,5 +124,5 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     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.
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_lap_blp.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_lap_blp.F90	(revision 12109)
@@ -89,5 +89,5 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     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.
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_triad.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/TRA/traldf_triad.F90	(revision 12109)
@@ -110,5 +110,5 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      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.
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/nemogcm.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/nemogcm.F90	(revision 12109)
@@ -479,5 +479,5 @@
                            CALL     flo_init    ! drifting Floats
       IF( ln_diacfl    )   CALL dia_cfl_init    ! Initialise CFL diagnostics
-                           CALL dia_ptr_init    ! Poleward TRansports initialization
+!                           CALL dia_ptr_init    ! Poleward TRansports initialization
                            CALL dia_dct_init    ! Sections tranports
                            CALL dia_hsb_init    ! heat content, salt content and volume budgets
Index: NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/step.F90	(revision 12108)
+++ NEMO/branches/2019/dev_r12072_TOP-01_ENHANCE-11_CEthe/src/OCE/step.F90	(revision 12109)
@@ -154,5 +154,4 @@
       IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp )       !       and/or eiv coeff.
       IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp )       ! eddy viscosity coeff. 
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       !  Ocean dynamics : hdiv, ssh, e3, u, v, w
@@ -188,10 +187,10 @@
          IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
       ENDIF
-                         CALL dyn_zdf       ( kstp )  ! vertical diffusion
-
-      IF( ln_dynspg_ts ) THEN                         
+                            CALL dyn_zdf    ( kstp )  ! vertical diffusion
+      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
                             CALL wzv        ( kstp )              ! now cross-level velocity 
          IF( ln_zad_Aimp )  CALL wAimp      ( kstp )  ! Adaptive-implicit vertical advection partitioning
       ENDIF
+      
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -205,7 +204,8 @@
       IF( ln_floats  )   CALL flo_stp ( kstp )        ! drifting Floats
       IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics
-      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth)
+                         CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth)
       IF( ln_diadct  )   CALL dia_dct ( kstp )        ! Transports
                          CALL dia_ar5 ( kstp )        ! ar5 diag
+                         CALL dia_ptr ( kstp )        ! Poleward adv/ldf TRansports diagnostics
       IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
                          CALL dia_wri ( kstp )        ! ocean model: outputs
@@ -243,7 +243,4 @@
                          CALL tra_ldf       ( kstp )  ! lateral mixing
 
-!!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
-      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics
-!!gm
                          CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields
       IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection
