Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyini.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyini.F90	(revision 12706)
@@ -398,5 +398,6 @@
       jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
       ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) )
-    
+      nbrdta(:,:,:) = 0   ! initialize nbrdta as it may not be completely defined for each bdy
+      
       ! Calculate global boundary index arrays or read in from file
       !------------------------------------------------------------               
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyvol.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyvol.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyvol.F90	(revision 12706)
@@ -143,5 +143,5 @@
       ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
       ! ------------------------------------------------------
-      IF( MOD( kt, nn_write ) == 0 .AND. ( kc == 1 ) ) THEN
+      IF( MOD( kt, MAX(nn_write,1) ) == 0 .AND. ( kc == 1 ) ) THEN
          !
          ! compute residual transport across boundary
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diaar5.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diaar5.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diaar5.F90	(revision 12706)
@@ -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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diahth.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diahth.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diahth.F90	(revision 12706)
@@ -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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diaptr.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diaptr.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diaptr.F90	(revision 12706)
@@ -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,23 +43,20 @@
 
    !                                  !!** 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 )
+   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 ::   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) 
+   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
 
    !! * Substitutions
@@ -80,51 +78,37 @@
       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( 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 +118,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 +196,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,10 +326,10 @@
       !! ** Purpose :   Initialization, namelist read
       !!----------------------------------------------------------------------
-      INTEGER ::  jn           ! local integers
-      INTEGER ::  inum, ierr   ! local integers
-      INTEGER ::  ios          ! Local integer output status for namelist read
+      INTEGER ::  inum, jn, ios, ierr           ! local integers
       !!
       NAMELIST/namptr/ ln_diaptr, ln_subbas
-      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: zmsk
+      !!----------------------------------------------------------------------
+
 
       REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
@@ -397,5 +339,5 @@
       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' )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist' )
       IF(lwm) WRITE ( numond, namptr )
 
@@ -406,51 +348,44 @@
          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
       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( 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           !
          !
       ENDIF 
@@ -471,57 +406,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 +471,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 +497,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 +509,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 +534,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 +546,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 +557,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 +601,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 +618,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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diatmb.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diatmb.F90	(revision 12705)
+++ 	(revision )
@@ -1,139 +1,0 @@
-MODULE diatmb 
-   !!======================================================================
-   !!                       ***  MODULE  diaharm  ***
-   !! Harmonic analysis of tidal constituents 
-   !!======================================================================
-   !! History :  3.6  !  08-2014  (E O'Dea)  Original code
-   !!            3.7  !  05-2016  (G. Madec)  use mbkt, mikt to account for ocean cavities
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain
-   !
-   USE in_out_manager  ! I/O units
-   USE iom             ! I/0 library
-   USE wet_dry
-
-   IMPLICIT NONE
-   PRIVATE
-
-   LOGICAL , PUBLIC ::   ln_diatmb     !: Top Middle and Bottom output 
-   PUBLIC   dia_tmb_init            ! routine called by nemogcm.F90
-   PUBLIC   dia_tmb                 ! routine called by diawri.F90
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dia_tmb_init 
-      !!---------------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_tmb_init  ***
-      !!     
-      !! ** Purpose :   Initialization of tmb namelist 
-      !!        
-      !! ** Method  :   Read namelist
-      !!---------------------------------------------------------------------------
-      INTEGER ::   ios                 ! Local integer output status for namelist read
-      !
-      NAMELIST/nam_diatmb/ ln_diatmb
-      !!----------------------------------------------------------------------
-      !
-      REWIND( numnam_ref )              ! Read Namelist nam_diatmb in reference namelist : TMB diagnostics
-      READ  ( numnam_ref, nam_diatmb, IOSTAT=ios, ERR= 901 )
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist' )
- 
-      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
-      READ  ( numnam_cfg, nam_diatmb, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist' )
-      IF(lwm) WRITE ( numond, nam_diatmb )
-
-      IF(lwp) THEN                   ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_tmb_init : Output Top, Middle, Bottom Diagnostics'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist nam_diatmb : set tmb outputs '
-         WRITE(numout,*) '      Switch for TMB diagnostics (T) or not (F)  ln_diatmb  = ', ln_diatmb
-      ENDIF
-      !
-   END SUBROUTINE dia_tmb_init
-
-
-   SUBROUTINE dia_calctmb( pfield, ptmb )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_tmb  ***
-      !!                   
-      !! ** Purpose :    Find the Top, Mid and Bottom fields of water Column
-      !!
-      !! ** Method  :    use mbkt, mikt to find surface, mid and bottom of 
-      !!              model levels due to potential existence of ocean cavities
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(in   ) ::   pfield   ! Input 3D field and mask
-      REAL(wp), DIMENSION(jpi, jpj,  3 ), INTENT(  out) ::   ptmb     ! top, middle, bottom extracted from pfield
-      !
-      INTEGER ::   ji, jj   ! Dummy loop indices
-      INTEGER ::   itop, imid, ibot   ! local integers
-      REAL(wp)::   zmdi = 1.e+20_wp   ! land value
-      !!---------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            itop = mikt(ji,jj)                        ! top    ocean 
-            ibot = mbkt(ji,jj)                        ! bottom ocean 
-            imid =  itop + ( ibot - itop + 1 ) / 2    ! middle ocean          
-            !                    
-            ptmb(ji,jj,1) = pfield(ji,jj,itop)*tmask(ji,jj,itop) + zmdi*( 1._wp-tmask(ji,jj,itop) )
-            ptmb(ji,jj,2) = pfield(ji,jj,imid)*tmask(ji,jj,imid) + zmdi*( 1._wp-tmask(ji,jj,imid) )
-            ptmb(ji,jj,3) = pfield(ji,jj,ibot)*tmask(ji,jj,ibot) + zmdi*( 1._wp-tmask(ji,jj,ibot) )
-         END DO
-      END DO
-      !
-   END SUBROUTINE dia_calctmb
-
-
-   SUBROUTINE dia_tmb
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE dia_tmb  ***
-      !! ** Purpose :   Write diagnostics for Top, Mid and Bottom of water Column
-      !!
-      !! ** Method  :  use mikt,mbkt to find surface, mid and bottom of model levels
-      !!      calls calctmb to retrieve TMB values before sending to iom_put
-      !!
-      !!--------------------------------------------------------------------
-      REAL(wp) ::   zmdi =1.e+20     ! land value
-      REAL(wp), DIMENSION(jpi,jpj,3) :: zwtmb    ! workspace 
-      !!--------------------------------------------------------------------
-      !
-      CALL dia_calctmb( tsn(:,:,:,jp_tem), zwtmb )
-      !ssh already output but here we output it masked
-      IF( ll_wd ) THEN
-         CALL iom_put( "sshnmasked", (sshn(:,:)+ssh_ref)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
-      ELSE
-         CALL iom_put( "sshnmasked", sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
-      ENDIF
-
-      CALL iom_put( "top_temp"  , zwtmb(:,:,1) )    ! tmb Temperature
-      CALL iom_put( "mid_temp"  , zwtmb(:,:,2) )    ! tmb Temperature
-      CALL iom_put( "bot_temp"  , zwtmb(:,:,3) )    ! tmb Temperature
-      !
-      CALL dia_calctmb( tsn(:,:,:,jp_sal), zwtmb )
-      CALL iom_put( "top_sal"   , zwtmb(:,:,1) )    ! tmb Salinity 
-      CALL iom_put( "mid_sal"   , zwtmb(:,:,2) )    ! tmb Salinity
-      CALL iom_put( "bot_sal"   , zwtmb(:,:,3) )    ! tmb Salinity
-      !
-      CALL dia_calctmb( un(:,:,:), zwtmb )
-      CALL iom_put( "top_u"     , zwtmb(:,:,1) )    ! tmb  U Velocity
-      CALL iom_put( "mid_u"     , zwtmb(:,:,2) )    ! tmb  U Velocity
-      CALL iom_put( "bot_u"     , zwtmb(:,:,3) )    ! tmb  U Velocity
-      !
-      CALL dia_calctmb( vn(:,:,:), zwtmb )
-      CALL iom_put( "top_v"     , zwtmb(:,:,1) )    ! tmb  V Velocity
-      CALL iom_put( "mid_v"     , zwtmb(:,:,2) )    ! tmb  V Velocity
-      CALL iom_put( "bot_v"     , zwtmb(:,:,3) )    ! tmb  V Velocity
-      !
-   END SUBROUTINE dia_tmb
-
-   !!======================================================================
-END MODULE diatmb
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90	(revision 12706)
@@ -48,5 +48,4 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager ! I/O manager
-   USE diatmb         ! Top,middle,bottom output
    USE dia25h         ! 25h Mean output
    USE iom            ! 
@@ -393,6 +392,4 @@
       CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
       !
-
-      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values 
           
       IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DOM/daymod.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DOM/daymod.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DOM/daymod.F90	(revision 12706)
@@ -93,6 +93,6 @@
       CALL ymds2ju( nyear, nmonth, nday, nhour*3600._wp+nminute*60._wp, fjulday )  
       IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < 0.1 / rday )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
-      IF( nn_time0*3600 - ndt05 .lt. 0 ) fjulday = fjulday + 1.                    ! move back to the day at nit000 (and not at nit000 - 1)
-
+      IF( nhour*3600 + nminute*60 - ndt05 .lt. 0 ) fjulday = fjulday + 1.                    ! move back to the day at nit000 (and not at nit000 - 1)
+      
       nsec1jan000 = 0
       CALL day_mth
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90	(revision 12706)
@@ -72,4 +72,5 @@
          IF(lwp) WRITE(numout,*) 'div_hor : horizontal velocity divergence '
          IF(lwp) WRITE(numout,*) '~~~~~~~   '
+         hdivn(:,:,:) = 0._wp    ! initialize hdivn for the halos at the first time step
       ENDIF
       !
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 12706)
@@ -226,28 +226,22 @@
             zcoef = atfp * rdt * r1_rau0
 
-            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            DO jk = 1, jpkm1
+               e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,jk) & 
+                             &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
+            END DO
 
             IF ( ln_rnf ) THEN
-               IF( ln_rnf_depth ) THEN
-                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
-                     DO jj = 1, jpj
-                        DO ji = 1, jpi
-                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
-                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
-                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
-                           ENDIF
-                        ENDDO
-                     ENDDO
-                  ENDDO
-               ELSE
-                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
-               ENDIF
-            END IF
-            !
-            ! ice shelf melting (deal separatly as it can be in depth)
-            ! PM: we could probably define a generic subroutine to do the in depth correction
-            !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...)
-            !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf)
-            IF ( ln_isf ) CALL isf_dynnxt( kt, atfp * rdt )
+               DO jk = 1, jpkm1
+                  e3t_b(:,:,jk) = e3t_b(:,:,jk) + zcoef * ( rnf_b(:,:) - rnf(:,:) ) * tmask(:,:,jk) & 
+                                &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
+               END DO
+            ENDIF
+
+            IF ( ln_isf ) THEN
+               DO jk = 1, jpkm1
+                  e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_b(:,:) - fwfisf(:,:) ) * tmask(:,:,jk) & 
+                                &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
+               END DO
+            ENDIF
             !
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90	(revision 12706)
@@ -47,5 +47,4 @@
    USE updtide         ! tide potential
    USE sbcwave         ! surface wave
-   USE diatmb          ! Top,middle,bottom output
 #if defined key_agrif
    USE agrif_oce_interp ! agrif
@@ -62,5 +61,4 @@
    USE iom             ! IOM library
    USE restart         ! only for lrst_oce
-   USE diatmb          ! Top,middle,bottom output
 
    USE iom   ! to remove
@@ -152,5 +150,5 @@
       REAL(wp) ::   r1_2dt_b, z1_hu, z1_hv          ! local scalars
       REAL(wp) ::   za0, za1, za2, za3              !   -      -
-      REAL(wp) ::   zmdi, zztmp, zldg               !   -      -
+      REAL(wp) ::   zztmp, zldg               !   -      -
       REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      -
       REAL(wp) ::   zun_save, zvn_save              !   -      -
@@ -175,6 +173,4 @@
       !                                         !* Allocate temporary arrays
       IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
-      !
-      zmdi=1.e+20                               !  missing data indicator for masking
       !
       zwdramp = r_rn_wdmin1               ! simplest ramp 
@@ -856,8 +852,6 @@
       IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
       !
-      IF( ln_diatmb ) THEN
-         CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
-         CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
-      ENDIF
+      CALL iom_put( "baro_u" , un_b )  ! Barotropic  U Velocity
+      CALL iom_put( "baro_v" , vn_b )  ! Barotropic  V Velocity
       !
    END SUBROUTINE dyn_spg_ts
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/ICB/icbthm.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/ICB/icbthm.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/ICB/icbthm.F90	(revision 12706)
@@ -20,4 +20,5 @@
    USE phycst         ! NEMO physical constants
    USE sbc_oce
+   USE lib_fortran, ONLY : DDPDD
 
    USE icb_oce        ! define iceberg arrays
@@ -55,5 +56,11 @@
       TYPE(iceberg), POINTER ::   this, next
       TYPE(point)  , POINTER ::   pt
+      !
+      COMPLEX(wp), DIMENSION(jpi,jpj) :: cicb_melt, cicb_hflx
       !!----------------------------------------------------------------------
+      !
+      !! initialiaze cicb_melt and cicb_heat
+      cicb_melt = CMPLX( 0.e0, 0.e0, wp ) 
+      cicb_hflx = CMPLX( 0.e0, 0.e0, wp ) 
       !
       z1_rday = 1._wp / rday
@@ -165,6 +172,12 @@
             z1_e1e2    = r1_e1e2t(ii,ij) * this%mass_scaling
             z1_dt_e1e2 = z1_dt * z1_e1e2
+            !
+            ! iceberg melt
+            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
             zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
-            berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt    * z1_e1e2    ! kg/m2/s
+            CALL DDPDD( CMPLX( zmelt * z1_e1e2, 0.e0, wp ), cicb_melt(ii,ij) )
+            !
+            ! iceberg heat flux
+            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
             !! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the
             !!     heat density of the icebergs is zero and the heat content flux to the ocean from iceberg
@@ -172,5 +185,7 @@
             zheat_hcflux = zmelt * pt%heat_density       ! heat content flux : kg/s x J/kg = J/s
             zheat_latent = - zmelt * rLfus               ! latent heat flux:  kg/s x J/kg = J/s
-            berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + ( zheat_hcflux + zheat_latent ) * z1_e1e2    ! W/m2
+            CALL DDPDD( CMPLX( ( zheat_hcflux + zheat_latent ) * z1_e1e2, 0.e0, wp ), cicb_hflx(ii,ij) )
+            !
+            ! diagnostics
             CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling,       &
                &                       zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
@@ -214,5 +229,8 @@
          !
       END DO
-
+      !
+      berg_grid%floating_melt = REAL(cicb_melt,wp)    ! kg/m2/s
+      berg_grid%calving_hflx  = REAL(cicb_hflx,wp)
+      !
       ! now use melt and associated heat flux in ocean (or not)
       !
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/IOM/iom.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/IOM/iom.F90	(revision 12706)
@@ -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
       !
@@ -795,4 +802,6 @@
       CHARACTER(LEN=100)    ::   clinfo    ! info character
       !---------------------------------------------------------------------
+      !
+      IF( iom_open_init == 0 )   RETURN   ! avoid to use iom_file(jf)%nfid that us not yet initialized
       !
       clinfo = '                    iom_close ~~~  '
@@ -1329,4 +1338,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 )  
@@ -1698,4 +1720,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
    !!----------------------------------------------------------------------
@@ -2040,12 +2073,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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfdyn.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfdyn.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfdyn.F90	(revision 12706)
@@ -115,9 +115,9 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
+      REWIND( numnam_ref )
       READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist' )
 
-      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
+      REWIND( numnam_cfg )
       READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist' )
@@ -417,5 +417,5 @@
             !
             zcmsmag   = (rn_csmc/rpi)**2                                            ! (C_smag/pi)^2
-            zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )         ! lower limit stability factor scaling
+            zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 12._wp * 12._wp * zcmsmag ) ! lower limit stability factor scaling
             zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )               ! upper limit stability factor scaling
             IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead 
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfslp.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfslp.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfslp.F90	(revision 12706)
@@ -208,4 +208,10 @@
                zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau )  )
                zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav )  )
+               !                                      ! Fred Dupont: add a correction for bottom partial steps:
+               !                                      !              max slope = 1/2 * e3 / e1
+               IF (ln_zps .AND. jk==mbku(ji,jj)) &
+                  zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , - 2._wp * e1u(ji,jj) / e3u_n(ji,jj,jk)* ABS( zau )  )
+               IF (ln_zps .AND. jk==mbkv(ji,jj)) &
+                  zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , - 2._wp * e2v(ji,jj) / e3v_n(ji,jj,jk)* ABS( zav )  )
                !                                      ! uslp and vslp output in zwz and zww, resp.
                zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
@@ -403,5 +409,4 @@
       REAL(wp) ::   zbeta0, ze3_e1, ze3_e2
       REAL(wp), DIMENSION(jpi,jpj)     ::   z1_mlbw
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zalbet
       REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) ::   zdxrho , zdyrho, zdzrho     ! Horizontal and vertical density gradients
       REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) ::   zti_mlb, ztj_mlb            ! for Griffies operator only
@@ -459,6 +464,6 @@
                      zdks = 0._wp
                   ENDIF
-                  zdzrho_raw = ( - rab_b(ji,jj,jk+kp,jp_tem) * zdkt & 
-                             &   + rab_b(ji,jj,jk+kp,jp_sal) * zdks &
+                  zdzrho_raw = ( - rab_b(ji,jj,jk   ,jp_tem) * zdkt & 
+                             &   + rab_b(ji,jj,jk   ,jp_sal) * zdks &
                              & ) / e3w_n(ji,jj,jk+kp)  
                   zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw )    ! force zdzrho >= repsln
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldftra.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldftra.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/LDF/ldftra.F90	(revision 12706)
@@ -152,8 +152,9 @@
       ! =================================
       !
-      REWIND( numnam_ref )              ! Namelist namtra_ldf in reference namelist : Lateral physics on tracers
+      REWIND( numnam_ref )
       READ  ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_ldf in reference namelist' )
-      REWIND( numnam_cfg )              ! Namelist namtra_ldf in configuration namelist : Lateral physics on tracers
+
+      REWIND( numnam_cfg )
       READ  ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist' )
@@ -510,9 +511,9 @@
       ENDIF
       !
-      REWIND( numnam_ref )              ! Namelist namtra_eiv in reference namelist : eddy induced velocity param.
+      REWIND( numnam_ref )
       READ  ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_eiv in reference namelist' )
       !
-      REWIND( numnam_cfg )              ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param.
+      REWIND( numnam_cfg )
       READ  ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist' )
@@ -770,7 +771,7 @@
             DO ji = 1, fs_jpim1   ! vector opt.
                zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   &
-                  &                                    * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * umask(ji,jj,jk)
+                  &                                    * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * wumask(ji,jj,jk)
                zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk) )   &
-                  &                                    * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * vmask(ji,jj,jk)
+                  &                                    * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * wvmask(ji,jj,jk)
             END DO
          END DO
@@ -851,4 +852,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 +886,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 +910,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( ln_diaptr )   CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
       !
       zztmp = 0.5_wp * 0.5
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/cpl_oasis3.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/cpl_oasis3.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/cpl_oasis3.F90	(revision 12706)
@@ -306,7 +306,13 @@
       ! End of definition phase
       !------------------------------------------------------------------
-      
+      !     
+#if defined key_agrif
+      IF( agrif_fixed() == Agrif_Nb_Fine_Grids() ) THEN
+#endif
       CALL oasis_enddef(nerror)
       IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef')
+#if defined key_agrif
+      ENDIF
+#endif
       !
       IF ( ltmp_wapatch ) THEN
@@ -357,7 +363,7 @@
                      WRITE(numout,*) 'oasis_put:  kstep ', kstep
                      WRITE(numout,*) 'oasis_put:   info ', kinfo
-                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(:,:,jc))
-                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(:,:,jc))
-                     WRITE(numout,*) '     -     Sum value is ', SUM(pdata(:,:,jc))
+                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(nldi:nlei,nldj:nlej,jc))
+                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(nldi:nlei,nldj:nlej,jc))
+                     WRITE(numout,*) '     -     Sum value is ', SUM(pdata(nldi:nlei,nldj:nlej,jc))
                      WRITE(numout,*) '****************'
                   ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/fldread.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/fldread.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/fldread.F90	(revision 12706)
@@ -833,6 +833,6 @@
 
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_z     ! depth of the data read in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_dz    ! thickness of the levels in bdy file
       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional)
       REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file
@@ -841,11 +841,11 @@
       INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number
       !!
-      INTEGER                                   ::   ipi             ! length of boundary data on local process
-      INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file
-      INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters
-      REAL(wp)                                  ::   zcoef
-      REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor
-      REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv)
-      REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth
+      INTEGER                  ::   ipi                 ! length of boundary data on local process
+      INTEGER                  ::   ipkb                ! number of vertical levels in boundary data file
+      INTEGER                  ::   ipkmax              ! number of vertical levels in boundary data file where no mask
+      INTEGER                  ::   jb, ji, jj, jk, jkb ! loop counters
+      REAL(wp)                 ::   zcoef, zi           ! 
+      REAL(wp)                 ::   ztrans, ztrans_new  ! transports
+      REAL(wp), DIMENSION(jpk) ::   zdepth, zdhalf      ! level and half-level depth
       !!---------------------------------------------------------------------
      
@@ -853,45 +853,20 @@
       ipkb = SIZE( pdta_read, 3 )
       
-      zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
-      !
-      WHERE( pdta_read == pfv )
-         pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data
-         pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct
-      ENDWHERE
-      
       DO jb = 1, ipi
          ji = idx_bdy(kbdy)%nbi(jb,kgrd)
          jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-         zh  = SUM(pdta_read_dz(jb,1,:) )
-         !
-         ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
+         !
+         ! --- max jk where input data /= FillValue --- !
+         ipkmax = 1
+         DO jkb = 2, ipkb
+            IF( pdta_read(jb,1,jkb) /= pfv )   ipkmax = MAX( ipkmax, jkb )
+         END DO
+         !
+         ! --- calculate depth at t,u,v points --- !
          SELECT CASE( kgrd )                         
-         CASE(1)
-            IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN
-               WRITE(ctmp1,"(I10.10)") jb 
-               CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
-               !   IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1),  ht_n(ji,jj), jb, jb, ji, jj
-            ENDIF
-         CASE(2)
-            IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN
-               WRITE(ctmp1,"(I10.10)") jb 
-               CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
-               !   IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1),  SUM(umask(ji,jj,:)), &
-               !      &                hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:)
-            ENDIF
-         CASE(3)
-            IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN
-               WRITE(ctmp1,"(I10.10)") jb
-               CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
-            ENDIF
-         END SELECT
-         !
-         SELECT CASE( kgrd )                         
-         CASE(1)
-            ! depth of T points:
+         CASE(1)            ! depth of T points:
             zdepth(:) = gdept_n(ji,jj,:)
-         CASE(2)
-            ! depth of U points: we must not use gdept_n as we don't want to do a communication
-            !   --> copy what is done for gdept_n in domvvl...
+         CASE(2)            ! depth of U points: we must not use gdept_n as we don't want to do a communication
+            !                 --> copy what is done for gdept_n in domvvl...
             zdhalf(1) = 0.0_wp
             zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1)
@@ -903,10 +878,9 @@
                zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
                zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1)
-               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  &
-                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk))
+               zdepth(jk) =       zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  &
+                  &         + (1.-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk))
             END DO
-         CASE(3)
-            ! depth of V points: we must not use gdept_n as we don't want to do a communication
-            !   --> copy what is done for gdept_n in domvvl...
+         CASE(3)            ! depth of V points: we must not use gdept_n as we don't want to do a communication
+            !                 --> copy what is done for gdept_n in domvvl...
             zdhalf(1) = 0.0_wp
             zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1)
@@ -918,76 +892,85 @@
                zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
                zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1)
-               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  &
-                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk))
+               zdepth(jk) =       zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  &
+                  &         + (1.-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk))
             END DO
          END SELECT
-         !
-         DO jk = 1, jpk                      
-            IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data
-               pdta(jb,1,jk) =  pdta_read(jb,1,1)
-            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data 
-               pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))
-            ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1
-               DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)
-                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &
-                     &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN   ! linear interpolation between 2 levels
+         !         
+         ! --- interpolate bdy data on the model grid --- !
+         DO jk = 1, jpk
+            IF(     zdepth(jk) <= pdta_read_z(jb,1,1)      ) THEN   ! above the first level of external data
+               pdta(jb,1,jk) = pdta_read(jb,1,1)
+            ELSEIF( zdepth(jk) >  pdta_read_z(jb,1,ipkmax) ) THEN   ! below the last level of external data /= FillValue
+               pdta(jb,1,jk) = pdta_read(jb,1,ipkmax)
+            ELSE                                                    ! inbetween: vertical interpolation between jkb & jkb+1
+               DO jkb = 1, ipkmax-1
+                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) ) <= 0._wp ) THEN ! linear interpolation between 2 levels
                      zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
-                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi
+                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) )
                   ENDIF
                END DO
             ENDIF
-         END DO   ! jpk
+         END DO
          !
       END DO   ! ipi
-      
-      IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term?
+
+      ! --- mask data and adjust transport --- !
+      SELECT CASE( kgrd )                         
+
+      CASE(1)                                 ! mask data (probably unecessary)
          DO jb = 1, ipi
             ji = idx_bdy(kbdy)%nbi(jb,kgrd)
             jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-            zh  = SUM(pdta_read_dz(jb,1,:) )
+            DO jk = 1, jpk                      
+               pdta(jb,1,jk) = pdta(jb,1,jk) * tmask(ji,jj,jk)
+            END DO
+         END DO
+         
+      CASE(2)                                 ! adjust the U-transport term
+         DO jb = 1, ipi
+            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
+            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
             ztrans = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
+            ENDDO
             ztrans_new = 0._wp
-            DO jkb = 1, ipkb                              ! calculate transport on input grid
-               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
-            ENDDO
             DO jk = 1, jpk                                ! calculate transport on model grid
-               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3u_n(ji,jj,jk ) * umask(ji,jj,jk)
+               ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)
             ENDDO
             DO jk = 1, jpk                                ! make transport correction
                IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
                   pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
-               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj)   * umask(ji,jj,jk)
+               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
+                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
                ENDIF
             ENDDO
          ENDDO
-      ENDIF
-      
-      IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term?
+
+      CASE(3)                                 ! adjust the V-transport term
          DO jb = 1, ipi
             ji = idx_bdy(kbdy)%nbi(jb,kgrd)
             jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-            zh  = SUM(pdta_read_dz(jb,1,:) )
             ztrans = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
+            ENDDO
             ztrans_new = 0._wp
-            DO jkb = 1, ipkb                              ! calculate transport on input grid
-               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
-            ENDDO
             DO jk = 1, jpk                                ! calculate transport on model grid
-               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v_n(ji,jj,jk ) * vmask(ji,jj,jk)
+               ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
             ENDDO
             DO jk = 1, jpk                                ! make transport correction
                IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
                   pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
-               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj)   * vmask(ji,jj,jk)
+               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
+                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
                ENDIF
             ENDDO
          ENDDO
-      ENDIF
-      
+      END SELECT
+
    END SUBROUTINE fld_bdy_interp
 
-
+   
    SUBROUTINE fld_rot( kt, sd )
       !!---------------------------------------------------------------------
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90	(revision 12706)
@@ -104,5 +104,5 @@
    !!              Ocean Surface Boundary Condition fields
    !!----------------------------------------------------------------------
-   INTEGER , PUBLIC ::  ncpl_qsr_freq            !: qsr coupling frequency per days from atmosphere
+   INTEGER , PUBLIC ::  ncpl_qsr_freq = 0        !: qsr coupling frequency per days from atmosphere (used by top)
    !
    LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress)
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcblk.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcblk.F90	(revision 12706)
@@ -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(:,:) ) 
+         CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average)
+         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(:,:) ) )
+         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 ) 
+          CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average)
+          CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
+          CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
+      ENDIF
       !
       IF(ln_ctl) THEN
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 12706)
@@ -574,4 +574,10 @@
       IF ( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl
 
+#if defined key_si3
+      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN 
+         IF( .NOT.srcv(jpr_ts_ice)%laction )  &
+            &   CALL ctl_stop( 'sbc_cpl_init: srcv(jpr_ts_ice)%laction should be set to true when ln_cndflx=T' )     
+      ENDIF
+#endif
       !                                                      ! ------------------------- !
       !                                                      !      Wave breaking        !    
@@ -863,5 +869,4 @@
       ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN  
          CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
-         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
       ENDIF
       ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
@@ -1041,9 +1046,4 @@
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
-      !
-      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
-      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
-         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
-      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq
       !
    END SUBROUTINE sbc_cpl_init
@@ -1111,4 +1111,12 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr
       !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN
+      !   cannot be done in the init phase when we use agrif as cpl_freq requires that oasis_enddef is done
+         ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
+         IF( ln_dm2dc .AND. ncpl_qsr_freq /= 86400 )   &
+            &   CALL ctl_stop( 'sbc_cpl_rcv: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
+         ncpl_qsr_freq = 86400 / ncpl_qsr_freq   ! used by top
+      ENDIF
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1244,16 +1252,4 @@
       IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
       !
-      !                                                      ! ================== !
-      !                                                      !   ice skin temp.   !
-      !                                                      ! ================== !
-#if defined key_si3
-      ! needed by Met Office
-      IF( srcv(jpr_ts_ice)%laction ) THEN 
-         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   tsfc_ice(:,:,:) = 0.0 
-         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   tsfc_ice(:,:,:) = -60.
-         ELSEWHERE                                        ;   tsfc_ice(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:)
-         END WHERE
-      ENDIF 
-#endif
       !                                                      ! ========================= ! 
       !                                                      ! Mean Sea Level Pressure   !   (taum) 
@@ -1635,11 +1631,11 @@
       !!                   sprecip           solid precipitation over the ocean  
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
-      !                                                !!           ! optional arguments, used only in 'mixed oce-ice' case
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
-      REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
+      REAL(wp), INTENT(in)   , DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
+      !                                                   !!           ! optional arguments, used only in 'mixed oce-ice' case or for Met-Office coupling
+      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
+      REAL(wp), INTENT(in)   , DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
+      REAL(wp), INTENT(inout), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] => inout for Met-Office
+      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
+      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
       !
       INTEGER  ::   ji, jj, jl   ! dummy loop index
@@ -1648,5 +1644,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice    !!gm , zfrqsr_tr_i
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu
       !!----------------------------------------------------------------------
       !
@@ -1774,11 +1770,12 @@
       IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
       IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
-      IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
-      IF( iom_use('precip') )       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
-      IF( iom_use('rain') )         CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
-      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , 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)  &
+      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
+      CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
+      IF ( iom_use('rain') ) CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
+      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)  &
          &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
       ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
@@ -1815,7 +1812,17 @@
 ! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
-         zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
-            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * ziceld(:,:)   &
-            &                                           + pist(:,:,1) * picefr(:,:) ) )
+         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
+            DO jl = 1, jpl
+               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    &
+                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
+                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+            END DO
+         ELSE
+            DO jl = 1, jpl
+               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    &
+                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
+                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+            END DO
+         ENDIF
       END SELECT
       !                                     
@@ -1902,14 +1909,14 @@
 #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) &
-           &                                                              * 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 ) &
-           &                                                              * ( 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 ) & 
-           &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice)
+      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_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 ) *  zsnw(:,:) )              ! heat flux from snow (over ice)
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
       !
@@ -1929,6 +1936,4 @@
             END DO
          ENDIF
-         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
-         zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
       CASE( 'oce and ice' )
          zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
@@ -1950,7 +1955,17 @@
 !       Create solar heat flux over ice using incoming solar heat flux and albedos
 !       ( see OASIS3 user guide, 5th edition, p39 )
-         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
-            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       &
-            &                     + palbi      (:,:,1) * picefr(:,:) ) )
+         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
+            DO jl = 1, jpl
+               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:,jl) * ( 1.- palbi(:,:,jl) )   &
+                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
+                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
+            END DO
+         ELSE
+            DO jl = 1, jpl
+               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:, 1) * ( 1.- palbi(:,:,jl) )   &
+                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
+                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
+            END DO
+         ENDIF
       CASE( 'none'      )       ! Not available as for now: needs additional coding  
       !                         ! since fields received, here zqsr_tot,  are not defined with none option
@@ -2012,8 +2027,14 @@
       !                                                      ! ========================= !
       CASE ('coupled')
-         qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
-         qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         IF( ln_mixcpl ) THEN
+            DO jl=1,jpl
+               qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:)
+               qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:)
+            ENDDO
+         ELSE
+            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
+            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         ENDIF
       END SELECT
-      !
       !                                                      ! ========================= !
       !                                                      !      Transmitted Qsr      !   [W/m2]
@@ -2022,9 +2043,13 @@
          !
          !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77)
+         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77)
          !
-         qtr_ice_top(:,:,:) = ztri * qsr_ice(:,:,:)
-         WHERE( phs(:,:,:) >= 0.0_wp )   qtr_ice_top(:,:,:) = 0._wp            ! snow fully opaque
-         WHERE( phi(:,:,:) <= 0.1_wp )   qtr_ice_top(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation
+         WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) )
+         ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm
+            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri
+         ELSEWHERE                                                         ! zero when hs>0
+            zqtr_ice_top(:,:,:) = 0._wp
+         END WHERE
          !     
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
@@ -2032,5 +2057,32 @@
          !                    ! ===> here we must receive the qtr_ice_top array from the coupler
          !                           for now just assume zero (fully opaque ice)
-         qtr_ice_top(:,:,:) = 0._wp
+         zqtr_ice_top(:,:,:) = 0._wp
+         !
+      ENDIF
+      !
+      IF( ln_mixcpl ) THEN
+         DO jl=1,jpl
+            qtr_ice_top(:,:,jl) = qtr_ice_top(:,:,jl) * xcplmask(:,:,0) + zqtr_ice_top(:,:,jl) * zmsk(:,:)
+         ENDDO
+      ELSE
+         qtr_ice_top(:,:,:) = zqtr_ice_top(:,:,:)
+      ENDIF
+      !                                                      ! ================== !
+      !                                                      !   ice skin temp.   !
+      !                                                      ! ================== !
+      ! needed by Met Office
+      IF( srcv(jpr_ts_ice)%laction ) THEN 
+         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   ztsu(:,:,:) =   0. + rt0 
+         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   ztsu(:,:,:) = -60. + rt0
+         ELSEWHERE                                        ;   ztsu(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) + rt0
+         END WHERE
+         !
+         IF( ln_mixcpl ) THEN
+            DO jl=1,jpl
+               pist(:,:,jl) = pist(:,:,jl) * xcplmask(:,:,0) + ztsu(:,:,jl) * zmsk(:,:)
+            ENDDO
+         ELSE
+            pist(:,:,:) = ztsu(:,:,:)
+         ENDIF
          !
       ENDIF
@@ -2195,5 +2247,5 @@
          CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
          END SELECT
-         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
+         CALL cpl_snd( jps_fice, isec, ztmp3, info )
       ENDIF
 
@@ -2255,5 +2307,5 @@
       !                                                      !      Ice melt ponds       ! 
       !                                                      ! ------------------------- !
-      ! needed by Met Office
+      ! needed by Met Office: 1) fraction of ponded ice 2) local/actual pond depth 
       IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN 
          SELECT CASE( sn_snd_mpnd%cldes)  
@@ -2261,12 +2313,12 @@
             SELECT CASE( sn_snd_mpnd%clcat )  
             CASE( 'yes' )  
-               ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl)
-               ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)  
+               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl)
+               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)  
             CASE( 'no' )  
                ztmp3(:,:,:) = 0.0  
                ztmp4(:,:,:) = 0.0  
                DO jl=1,jpl  
-                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip(:,:,jpl)  
-                 ztmp4(:,:,1) = ztmp4(:,:,1) + v_ip(:,:,jpl) 
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)  
+                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
                ENDDO  
             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )  
@@ -2306,5 +2358,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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90	(revision 12706)
@@ -236,4 +236,12 @@
 #endif
       !
+      !
+      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
@@ -536,4 +544,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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90	(revision 12706)
@@ -20,5 +20,5 @@
    USE sbc_oce        ! surface boundary condition variables
    USE eosbn2         ! Equation Of State
-   USE closea         ! closed seas
+   USE closea, ONLY: l_clo_rnf, clo_rnf ! closed seas
    !
    USE in_out_manager ! I/O manager
@@ -42,8 +42,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
@@ -64,4 +66,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)  
@@ -111,5 +114,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
@@ -117,5 +124,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
@@ -128,10 +142,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
@@ -238,6 +252,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
@@ -261,9 +275,9 @@
       !                                   ! ============
       !
-      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      REWIND( numnam_ref )
       READ  ( numnam_ref, namsbc_rnf, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_rnf in reference namelist' )
 
-      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      REWIND( numnam_cfg )
       READ  ( numnam_cfg, namsbc_rnf, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_rnf in configuration namelist' )
@@ -295,4 +309,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
       !
@@ -337,5 +366,5 @@
                IF( h_rnf(ji,jj) > 0._wp ) THEN
                   jk = 2
-                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
                   END DO
                   nk_rnf(ji,jj) = jk
@@ -394,5 +423,5 @@
                IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
                   jk = 2
-                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
                   END DO
                   nk_rnf(ji,jj) = jk
@@ -435,5 +464,5 @@
          !                                      !    - mixed upstream-centered (ln_traadv_cen2=T)
          !
-         IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
+         IF( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
             &                                              'be spread through depth by ln_rnf_depth'               )
          !
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcssr.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcssr.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcssr.F90	(revision 12706)
@@ -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_sssr_ice     ! 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_sssr_ice /= 1 ) THEN
+              ! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_sssr_ice .ne. 1
+              ! n.b. coefice is initialised and fixed to 1._wp if nn_sssr_ice = 1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     SELECT CASE ( nn_sssr_ice )
+                       CASE ( 0 )    ;  coefice(ji,jj) = 1._wp - fr_i(ji,jj)              ! no/reduced damping under ice
+                       CASE  DEFAULT ;  coefice(ji,jj) = 1._wp + ( nn_sssr_ice - 1 ) * fr_i(ji,jj) ! reinforced damping (x nn_sssr_ice) 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_sssr_ice
       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_sssr_ice    = ', nn_sssr_ice
+         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_sssr_ice=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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbc.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbc.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbc.F90	(revision 12706)
@@ -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' )
       !
@@ -133,9 +135,9 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnam_ref )              ! Namelist nambbc in reference namelist : Bottom momentum boundary condition
+      REWIND( numnam_ref )
       READ  ( numnam_ref, nambbc, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbc in reference namelist' )
       !
-      REWIND( numnam_cfg )              ! Namelist nambbc in configuration namelist : Bottom momentum boundary condition
+      REWIND( numnam_cfg )
       READ  ( numnam_cfg, nambbc, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbc in configuration namelist' )
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90	(revision 12706)
@@ -267,5 +267,5 @@
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
       REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
-      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
+      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d, zscale  !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
       !!----------------------------------------------------------------------
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90	(revision 12706)
@@ -73,5 +73,4 @@
    USE step_diu       ! diurnal bulk SST timestepping (called from here if run offline)
    USE crsini         ! initialise grid coarsening utility
-   USE diatmb         ! Top,middle,bottom output
    USE dia25h         ! 25h mean output
    USE sbc_oce , ONLY : lk_oasis
@@ -489,5 +488,4 @@
                            CALL     trd_init    ! Mixed-layer/Vorticity/Integral constraints trends
                            CALL dia_obs_init    ! Initialize observational data
-                           CALL dia_tmb_init    ! TMB outputs
                            CALL dia_25h_init    ! 25h mean  outputs
                            CALL dia_harm_init   ! tidal harmonics outputs
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/step.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 12706)
@@ -155,5 +155,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
@@ -189,10 +188,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
+      
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -206,7 +205,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
+      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics
       IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
                          CALL dia_wri ( kstp )        ! ocean model: outputs
@@ -245,7 +245,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
