Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/ISOMIP/EXP00/iodef.xml
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/ISOMIP/EXP00/iodef.xml	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/ISOMIP/EXP00/iodef.xml	(revision 5189)
@@ -32,5 +32,5 @@
       <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files -->
       <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->    
-      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->  
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/>  <!-- 5d files -->  
       <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files -->
 	<file id="file1" output_freq="1mo" name_suffix="_grid_T" description="ocean T grid variables" >
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 5189)
@@ -549,7 +549,6 @@
          &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
       IF( ln_zps .AND.        ln_isfcav)                            &
-         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-         &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-         &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
+         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+         &                                        rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
 
       rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 5189)
@@ -750,7 +750,6 @@
                &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
             IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
-               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF)
-               &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF)
+               &                                  rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level
 
 #if defined key_zdfkpp
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5189)
@@ -646,31 +646,4 @@
       CALL iom_close( inum )
       
-! need to be define for the extended grid south of -80S
-! some point are undefined but you need to have e1 and e2 .NE. 0
-      WHERE (e1t==0.0_wp)
-         e1t=1.0e2
-      END WHERE
-      WHERE (e1v==0.0_wp)
-         e1v=1.0e2
-      END WHERE
-      WHERE (e1u==0.0_wp)
-         e1u=1.0e2
-      END WHERE
-      WHERE (e1f==0.0_wp)
-         e1f=1.0e2
-      END WHERE
-      WHERE (e2t==0.0_wp)
-         e2t=1.0e2
-      END WHERE
-      WHERE (e2v==0.0_wp)
-         e2v=1.0e2
-      END WHERE
-      WHERE (e2u==0.0_wp)
-         e2u=1.0e2
-      END WHERE
-      WHERE (e2f==0.0_wp)
-         e2f=1.0e2
-      END WHERE
-       
     END SUBROUTINE hgr_read
     
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5189)
@@ -263,5 +263,5 @@
          END DO
       END DO
-      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
+      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector loop
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5189)
@@ -1538,5 +1538,5 @@
  
  ! remove single point "bay" on isf coast line in the ice shelf draft'
-         DO jk = 1, jpk
+         DO jk = 2, jpk
             WHERE (misfdep==0) misfdep=jpk
             zmask=0
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5189)
@@ -141,7 +141,6 @@
             &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level
          IF( ln_zps .AND.       ln_isfcav)                                 &
-            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-            &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF)
+            &                                            rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level
 #endif
          !   
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5189)
@@ -44,4 +44,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE iom
 
    IMPLICIT NONE
@@ -130,4 +131,10 @@
       INTEGER ::   ioptio = 0      ! temporary integer
       INTEGER ::   ios             ! Local integer output status for namelist read
+      !!
+      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
+      REAL(wp) ::   znad
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop, zrhd ! hypothesys on isf density
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_isf  ! density at bottom of ISF
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  ziceload     ! density at bottom of ISF
       !!
       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
@@ -191,6 +198,50 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
       ! 
-      ! initialisation of ice load
-      riceload(:,:)=0.0
+      ! initialisation of ice shelf load
+      IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0
+      IF (       ln_isfcav ) THEN
+         CALL wrk_alloc( jpi,jpj, 2,  ztstop) 
+         CALL wrk_alloc( jpi,jpj,jpk, zrhd  )
+         CALL wrk_alloc( jpi,jpj,     zrhdtop_isf, ziceload) 
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'  
+
+         ! To use density and not density anomaly
+         znad=1._wp
+         
+         ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
+         ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
+
+         ! compute density of the water displaced by the ice shelf 
+         DO jk = 1, jpk
+            CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk))
+         END DO
+      
+         ! compute rhd at the ice/oce interface (ice shelf side)
+         CALL eos(ztstop,risfdep,zrhdtop_isf)
+
+         ! Surface value + ice shelf gradient
+         ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
+         ziceload = 0._wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi   ! vector opt.
+               ikt=mikt(ji,jj)
+               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
+               DO jk=2,ikt-1
+                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) &
+                     &                              * (1._wp - tmask(ji,jj,jk))
+               END DO
+               IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
+                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
+            END DO
+         END DO
+         riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
+
+         CALL wrk_dealloc( jpi,jpj, 2,  ztstop) 
+         CALL wrk_dealloc( jpi,jpj,jpk, zrhd  )
+         CALL wrk_dealloc( jpi,jpj,     zrhdtop_isf, ziceload) 
+      END IF
       !
    END SUBROUTINE dyn_hpg_init
@@ -465,78 +516,35 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj, zrhd
+      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j   ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, znad                    ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop
-      REAL(wp), POINTER, DIMENSION(:,:)     ::  ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj, 2, ztstop) 
-      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
-      CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
-      !
-     IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
-      ENDIF
-
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce, zpshpi, zpshpj
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj,  2, ztstop) 
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj)
+      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce, zpshpi, zpshpj) 
+      !
       ! Local constant initialization
       zcoef0 = - grav * 0.5_wp
+   
       ! To use density and not density anomaly
-!      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
-!      ELSE                 ;     znad = 0._wp         ! Fixed volume
-!      ENDIF
       znad=1._wp
+
       ! iniitialised to 0. zhpi zhpi 
       zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
 
-!==================================================================================     
-!=====Compute iceload and contribution of the half first wet layer ================= 
-!===================================================================================
-
-      ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
-      ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
-
-      ! compute density of the water displaced by the ice shelf 
-      zrhd = rhd ! save rhd
-      DO jk = 1, jpk
-           zdept(:,:)=gdept_1d(jk)
-           CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))
-      END DO
-      WHERE ( tmask(:,:,:) == 1._wp)
-        rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd
-      END WHERE
-      
-      ! compute rhd at the ice/oce interface (ice shelf side)
-      CALL eos(ztstop,risfdep,zrhdtop_isf)
-
-      ! compute rhd at the ice/oce interface (ocean side)
-      DO ji=1,jpi
-        DO jj=1,jpj
-          ikt=mikt(ji,jj)
-          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1)
-          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2)
-        END DO
-      END DO
-      CALL eos(ztstop,risfdep,zrhdtop_oce)
-      !
-      ! Surface value + ice shelf gradient
-      ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
-      ziceload = 0._wp
-      DO jj = 1, jpj
-         DO ji = 1, jpi   ! vector opt.
-            ikt=mikt(ji,jj)
-            ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
-            DO jk=2,ikt-1
-               ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) &
-                  &                              * (1._wp - tmask(ji,jj,jk))
-            END DO
-            IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) &
-                               &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
-         END DO
-      END DO
-      riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
-      ! compute zp from z=0 to first T wet point (correction due to zps not yet applied)
+      ! compute rhd at the ice/oce interface (ocean side) usefull to reduce unrealistic water current if homogene water
+!  TO BE DISCUSS do we compute density at w point for the ocean/ice atmosphere
+!      DO ji=1,jpi
+!        DO jj=1,jpj
+!          ikt=mikt(ji,jj)
+!          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1)
+!          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2)
+!        END DO
+!      END DO
+!      CALL eos(ztstop,risfdep,zrhdtop_oce)
+
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -544,14 +552,20 @@
             ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
             ! we assume ISF is in isostatic equilibrium
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    &
-               &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   &
-               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
-               &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              ) 
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    &
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   &
-               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
-               &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             ) 
+            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,iktp1i) * ( znad + rhd(ji+1,jj  ,iktp1i) )   &
+               &                                  - fse3w(ji  ,jj  ,ikt   ) * ( znad + rhd(ji  ,jj  ,ikt   ) )   &
+               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                       ) 
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,iktp1j) * ( znad + rhd(ji  ,jj+1,iktp1j) )   &
+               &                                  - fse3w(ji  ,jj  ,ikt   ) * ( znad + rhd(ji  ,jj  ,ikt   ) )   &
+               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj) )                      ) 
+!            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    &
+!               &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   &
+!               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &
+!               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
+!               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                              ) 
+!            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    &
+!               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   &
+!               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
+!               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
+!               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj) )                             ) 
             ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
@@ -564,141 +578,34 @@
          END DO
       END DO
-!==================================================================================     
-!===== Compute partial cell contribution for the top cell ========================= 
-!==================================================================================
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            iku = miku(ji,jj) ; 
-            zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp
-            ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
-            ! u direction
-            IF ( iku .GT. 1 ) THEN
-               ! case iku
-               zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            &
-                      &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       &
-                      &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )
-               ! corrective term ( = 0 if z coordinate )
-               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)
-               ! zhpi will be added in interior loop
-               ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap
-               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
-               IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku)
-
-               ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)
-               zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &    
-                  &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          &
-                  &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   &
-                  &              - fse3w(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          &
-                  &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   )
-               zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
-            END IF
-               
-            ! v direction
-            ikv = mikv(ji,jj)
-            ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-            IF ( ikv .GT. 1 ) THEN
-               ! case ikv
-               zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            &
-                     &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           &
-                     &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )
-               ! corrective term ( = 0 if z coordinate )
-               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)
-               ! zhpi will be added in interior loop
-               va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap
-               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
-               IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv) 
-               
-               ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)
-               zhpjint        =  zcoef0 / e2v(ji,jj)                                                              &
-                  &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         &
-                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    &
-                  &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         &
-                  &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  )
-               zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
-            END IF
-         END DO
-      END DO
 
 !==================================================================================     
 !===== Compute interior value ===================================================== 
 !==================================================================================
-
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            iku=miku(ji,jj); ikv=mikv(ji,jj)
-            DO jk = 2, jpkm1
+      ! interior value (2=<jk=<jpkm1)
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                ! hydrostatic pressure gradient along s-surfaces
-               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
-               zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              &
-                  &                            + zcoef0 / e1u(ji,jj)                                                           &
-                  &                                     * ( fse3w(ji+1,jj  ,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   &
-                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   ) 
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
+                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
+                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
+                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
+                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
                ! s-coordinate pressure gradient correction
-               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
-               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    &
-                  &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1)
-               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
-
-               ! hydrostatic pressure gradient along s-surfaces
-               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
-               zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              &
-                  &                            + zcoef0 / e2v(ji,jj)                                                           &
-                  &                                     * ( fse3w(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   &
-                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   )
-               ! s-coordinate pressure gradient correction
-               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
-               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     &
-                  &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1)
+               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
+               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
                ! add to the general momentum trend
-               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk)
+               ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-
-!==================================================================================     
-!===== Compute bottom cell contribution (partial cell) ============================ 
-!==================================================================================
-
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-            iku = mbku(ji,jj)
-            ikv = mbkv(ji,jj)
-
-            IF (iku .GT. 1) THEN
-               ! remove old value (interior case)
-               zuap            = -zcoef0 * ( rhd   (ji+1,jj  ,iku) + rhd   (ji,jj,iku) + 2._wp * znad )   &
-                     &                   * ( fsde3w(ji+1,jj  ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj)
-               ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap
-               ! put new value
-               ! -zpshpi to avoid double contribution of the partial step in the top layer 
-               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj)
-               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 
-               ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap
-            END IF
-            ! v direction
-            IF (ikv .GT. 1) THEN
-               ! remove old value (interior case)
-               zvap            = -zcoef0 * ( rhd   (ji  ,jj+1,ikv) + rhd   (ji,jj,ikv) + 2._wp * znad )   &
-                     &                   * ( fsde3w(ji  ,jj+1,ikv) - fsde3w(ji,jj,ikv) )   / e2v(ji,jj)
-               va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap
-               ! put new value
-               ! -zpshpj to avoid double contribution of the partial step in the top layer 
-               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj)
-               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)   
-               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap
-            END IF
-         END DO
-      END DO
-     
-      ! set back to original density value into the ice shelf cell (maybe useless because it is masked)
-      rhd = zrhd
-      !
-      CALL wrk_dealloc( jpi,jpj,2, ztstop)
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
-      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)
+      !
+      CALL wrk_dealloc( jpi,jpj,2  , ztstop)
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj)
+      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce, zpshpi, zpshpj)
       !
    END SUBROUTINE hpg_isf
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5189)
@@ -82,38 +82,26 @@
       ENDIF
       
-      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
+      !                             ! Vertical momentum advection at uw and vw-points
+      DO jk = 2, jpk                
          DO jj = 2, jpj                   ! vertical fluxes 
-            DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+            DO ji = fs_2, jpi
+               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
-         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
-            DO ji = fs_2, fs_jpim1        ! vector opt.
-               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
-               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
+         DO jj = 2, jpjm1                 ! interior vertical momentum advection at w-point
+            DO ji = fs_2, fs_jpim1
+               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )
+               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )
             END DO  
          END DO   
       END DO
-      !
-      ! Surface and bottom advective fluxes set to zero
-      IF ( ln_isfcav ) THEN
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1           ! vector opt.
-               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
-               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
-               zwuw(ji,jj,jpk) = 0._wp
-               zwvw(ji,jj,jpk) = 0._wp
-            END DO
-         END DO
-      ELSE
-         DO jj = 2, jpjm1        
-            DO ji = fs_2, fs_jpim1           ! vector opt.
-               zwuw(ji,jj, 1 ) = 0._wp
-               zwvw(ji,jj, 1 ) = 0._wp
-               zwuw(ji,jj,jpk) = 0._wp
-               zwvw(ji,jj,jpk) = 0._wp
-            END DO  
-         END DO
-      END IF
+      DO jj = 2, jpjm1                    ! First/last level advective fluxes
+         DO ji = fs_2, fs_jpim1           ! w is zero at k=1 and k=jpk
+            zwuw(ji,jj, 1 ) = 0._wp
+            zwvw(ji,jj, 1 ) = 0._wp
+            zwuw(ji,jj,jpk) = 0._wp
+            zwvw(ji,jj,jpk) = 0._wp
+         END DO
+      END DO
 
       DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
@@ -205,5 +193,5 @@
          DO jj = 2, jpj                   
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5189)
@@ -162,5 +162,5 @@
             !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
             !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1
-            !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
+            !                                !          tmask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
             !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster
             zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              &
@@ -188,12 +188,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
-                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp)
-                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp)
-                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
-                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp)
-                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp)
-               ENDDO
-            ENDDO
+                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) - 0.5 * ( risfdep(ji,jj) + risfdep(ji+1,jj  ) )
+                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) - 0.5 * ( risfdep(ji,jj) + risfdep(ji  ,jj+1) )
+               END DO
+            END DO
          ELSE
             DO jj = 2, jpjm1
@@ -201,7 +197,9 @@
                   zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
                   zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
-               ENDDO
-            ENDDO
+               END DO
+            END DO
          END IF
+! PM + GM can be optimized by using : zuslp_hmlpu(ji,jj)= uslpml(ji,jj) / zhmlpu (ji,jj)
+
          DO jk = 2, jpkm1                            !* Slopes at u and v points
             DO jj = 2, jpjm1
@@ -220,8 +218,8 @@
                   zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) )
                   ! thickness of water column between surface and level k at u/v point
-                  zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              &
-                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) )
-                  zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              &
-                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) )
+                  zdepu = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji+1,jj  ,jk) )                            &
+                                   - ( risfdep(ji,jj)    + risfdep(ji+1,jj)    ) - fse3u(ji,jj,miku(ji,jj)) )
+                  zdepv = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji,jj+1,jk) )                              &
+                                   - ( risfdep(ji,jj)    + risfdep(ji,jj+1)    ) - fse3v(ji,jj,mikv(ji,jj)) )
                   !
                   zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          &
@@ -315,5 +313,5 @@
                   !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
                   zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
-                  zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp )
+                  zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - fsdepw(ji,jj,mikt(ji,jj)), 10._wp )
                   zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) &
                      &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
@@ -426,6 +424,6 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt. 
-               uslp(ji,jj,1) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) )  * umask(ji,jj,1) 
-               vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) )  * vmask(ji,jj,1) 
+               uslp(ji,jj,1)  = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) )  * umask(ji,jj,1) 
+               vslp(ji,jj,1)  = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) )  * vmask(ji,jj,1) 
                wslpi(ji,jj,1) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5 
                wslpj(ji,jj,1) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5 
@@ -436,6 +434,6 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt. 
-                  uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
-                  vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
+                  uslp(ji,jj,jk)  = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
+                  vslp(ji,jj,jk)  = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) &
                     &                              * wmask(ji,jj,jk) * 0.5 
@@ -752,8 +750,6 @@
             DO ji = 1, jpi
                ik = nmln(ji,jj) - 1
-               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN
-                  omlmask(ji,jj,jk) = 1._wp
-               ELSE
-                  omlmask(ji,jj,jk) = 0._wp
+               IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp
+               ELSE                ; omlmask(ji,jj,jk) = 0._wp
                ENDIF
             END DO
@@ -777,6 +773,6 @@
             !
             !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point)
-            iku = MIN(  MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
-            ikv = MIN(  MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
+            iku = MIN(  MAX( nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
+            ikv = MIN(  MAX( nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
             zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) )
             zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) )
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 5189)
@@ -851,5 +851,5 @@
                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+                  &            / fse3w(ji,jj,jk) * wmask(ji,jj,jk)
             END DO
          END DO
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5189)
@@ -108,12 +108,12 @@
       REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
       REAL(wp) ::  zcoef0, zbtr, ztra            !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw 
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zdkt, zdk1t, z2d
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdit, zdjt , ztfw 
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso')
       !
-      CALL wrk_alloc( jpi, jpj,      z2d ) 
-      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
+      CALL wrk_alloc( jpi, jpj,      zdkt, zdk1t, z2d  ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt , ztfw ) 
       !
 
@@ -168,34 +168,19 @@
          !!   II - horizontal trend  (full)
          !!----------------------------------------------------------------------
-!!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t ) 
+!CDIR PARALLEL DO PRIVATE( zdk1t ) 
+         !                                                ! ===============
+         DO jk = 1, jpkm1                                 ! Horizontal slab
+            !                                             ! ===============
             ! 1. Vertical tracer gradient at level jk and jk+1
             ! ------------------------------------------------
-         ! 
-         ! interior value 
-         DO jk = 2, jpkm1               
-            DO jj = 1, jpj
-               DO ji = 1, jpi   ! vector opt.
-                  zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)
-                  !
-                  zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
-         zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2)
-         zdkt (:,:,1) = zdk1t(:,:,1)
-         IF ( ln_isfcav ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi   ! vector opt.
-                  ikt = mikt(ji,jj) ! surface level
-                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)
-                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)
-               END DO
-            END DO
-         END IF
-
-         ! 2. Horizontal fluxes
-         ! --------------------   
-         DO jk = 1, jpkm1
+            ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
+            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)
+            !
+            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)
+            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
+            ENDIF
+
+            ! 2. Horizontal fluxes
+            ! --------------------   
             DO jj = 1 , jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -203,9 +188,9 @@
                   zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
                   !
-                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
-                     &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. )
-                  !
-                  zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
-                     &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
+                  zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
+                     &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. )
+                  !
+                  zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   &
+                     &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. )
                   !
                   zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
@@ -213,9 +198,9 @@
                   !
                   zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
-                     &              + zcof1 * (  zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk)      &
-                     &                         + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk)
+                     &              + zcof1 * (  zdkt (ji+1,jj  ) + zdk1t(ji,jj)      &
+                     &                         + zdk1t(ji+1,jj  ) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
                   zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
-                     &              + zcof2 * (  zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk)      &
-                     &                         + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                  
+                     &              + zcof2 * (  zdkt (ji  ,jj+1) + zdk1t(ji,jj)      &
+                     &                         + zdk1t(ji  ,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                  
                END DO
             END DO
@@ -322,6 +307,6 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, z2d ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
+      CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d  ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , ztfw ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso')
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 5189)
@@ -222,5 +222,5 @@
       !
       IF( nn_isf > 0 ) THEN
-         zfact = 0.5e0
+         zfact = 0.5_wp
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1
@@ -235,12 +235,13 @@
                ! compute tfreez for the temperature correction (we add water at freezing temperature)
 !                  zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
-                  zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
+                  zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
                ! compute trend
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
-                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
-                     &               - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                &
+                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)               &
+                     &                     - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &
                      &           * r1_hisf_tbl(ji,jj)
-                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
-                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                                &
+                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal))               &
+                     &           * r1_hisf_tbl(ji,jj)
                END DO
    
@@ -248,12 +249,14 @@
                ! compute tfreez for the temperature correction (we add water at freezing temperature)
 !               zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04
-               zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress )
+               zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress )
                ! compute trend
-               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
-                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
-                  &                  - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
+               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 &
+                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)               &
+                  &                        - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
                   &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
-               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
-                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
+               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                                 &
+                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal))               &
+                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
+
             END DO
          END DO
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5189)
@@ -196,7 +196,6 @@
    END SUBROUTINE zps_hde
    !
-   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   &
-      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  &
-      &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi )
+   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
+      &                          prd, pgru, pgrv, pgrui, pgrvi )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zps_hde  ***
@@ -252,11 +251,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom)
       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top)
       !
       INTEGER  ::   ji, jj, jn      ! Dummy loop indices
@@ -269,8 +262,8 @@
       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf')
       !
-      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
+      pgtu (:,:,:)=0.0_wp ; pgtv(:,:,:) =0.0_wp ;
       pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ;
-      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ;
-      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ;
+      zti  (:,:,:)=0.0_wp ; ztj  (:,:,:)=0.0_wp ;
+      zhi  (:,:  )=0.0_wp ; zhj  (:,:  )=0.0_wp ;
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
@@ -280,43 +273,132 @@
                iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
                ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
+               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
+               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
+               !
+               ! i- direction
+               IF( ze3wu >= 0._wp ) THEN      ! case 1
+                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
+                  ! interpolated values of tracers
+                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
+                  ! gradient of  tracers
+                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
+               ELSE                           ! case 2
+                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
+                  ! interpolated values of tracers
+                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
+                  ! gradient of tracers
+                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
+               ENDIF
+               !
+               ! j- direction
+               IF( ze3wv >= 0._wp ) THEN      ! case 1
+                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
+                  ! interpolated values of tracers
+                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
+               ELSE                           ! case 2
+                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
+                  ! interpolated values of tracers
+                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
+         !
+      END DO
+
+      ! horizontal derivative of density anomalies (rd)
+      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
+         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
+         !
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               iku = mbku(ji,jj)
+               ikv = mbkv(ji,jj)
+               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
+               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
+
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2
+               ENDIF
+
+            END DO
+         END DO
+
+         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
+         ! step and store it in  zri, zrj for each  case
+         CALL eos( zti, zhi, zri )
+         CALL eos( ztj, zhj, zrj )
+
+         ! Gradient of density at the last level 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               iku = mbku(ji,jj)
+               ikv = mbkv(ji,jj)
+               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
+               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
+
+               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
+               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
+               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions
+         !
+      END IF
+         ! (ISH)  compute grui and gruvi
+      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1
+               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1
+               !
                ! (ISF) case partial step top and bottom in adjacent cell in vertical
                ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
-               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
-               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-               !
+               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
+               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv) 
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
-                  ! interpolated values of tracers
-                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
+                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1)
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) )
+                  ! gradient of tracers
+                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
+               ELSE                           ! case 2
+                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1)
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) )
                   ! gradient of  tracers
-                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
-               ELSE                           ! case 2
-                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
-                  ! interpolated values of tracers
-                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
-                  ! gradient of tracers
-                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
+                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
                ENDIF
                !
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
-                  ! interpolated values of tracers
-                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
-                  ! gradient of tracers
-                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
-               ELSE                           ! case 2
-                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
-                  ! interpolated values of tracers
-                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
-                  ! gradient of tracers
-                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
+                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1)
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) )
+                  ! gradient of tracers
+                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
+               ELSE                           ! case 2
+                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1)
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) )
+                  ! gradient of tracers
+                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
+               ENDIF
+            END DO!!
+         END DO!!
+         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
          !
       END DO
@@ -324,24 +406,23 @@
       ! horizontal derivative of density anomalies (rd)
       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
-         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
-         pgzu(:,:)=0.0_wp   ; pgzv(:,:)=0.0_wp ;
-         pmru(:,:)=0.0_wp   ; pmru(:,:)=0.0_wp ;
-         pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ;
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               iku = mbku(ji,jj)
-               ikv = mbkv(ji,jj)
-               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
-               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2
-               ENDIF
-            END DO
-         END DO
-         
+         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ;
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               iku = miku(ji,jj)
+               ikv = mikv(ji,jj)
+               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
+               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv) 
+
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2
+               ENDIF
+
+            END DO
+         END DO
+
          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
          ! step and store it in  zri, zrj for each  case
@@ -352,166 +433,20 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
-               ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points
-               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
-               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-               IF( ze3wu >= 0._wp ) THEN 
-                  pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku)
-                  pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1
-                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1 
-                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
-                                * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) &
-                                   - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
-               ELSE  
-                  pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu)
-                  pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2
-                  pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2
-                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
-                                * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) &
-                                   -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN
-                  pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 
-                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1
-                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1
-                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
-                                * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) &
-                                   - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
-               ELSE 
-                  pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv)
-                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2
-                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2
-                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
-                                * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) &
-                                   -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions
-         !
-      END IF
-         ! (ISH)  compute grui and gruvi
-      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1
-               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1
-               !
-               ! (ISF) case partial step top and bottom in adjacent cell in vertical
-               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
-               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
-               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
-               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
-               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-               ! i- direction
-               IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1)
-                  ! interpolated values of tracers
-                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) )
-                  ! gradient of tracers
-                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
-               ELSE                           ! case 2
-                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1)
-                  ! interpolated values of tracers
-                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) )
-                  ! gradient of  tracers
-                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
-               ENDIF
-               !
-               ! j- direction
-               IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1)
-                  ! interpolated values of tracers
-                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) )
-                  ! gradient of tracers
-                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
-               ELSE                           ! case 2
-                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1)
-                  ! interpolated values of tracers
-                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) )
-                  ! gradient of tracers
-                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
-               ENDIF
-            END DO!!
-         END DO!!
-         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
-         !
-      END DO
-
-      ! horizontal derivative of density anomalies (rd)
-      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
-         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ;
-         pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ;
-         pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ;
-         pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ;
-
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               iku = miku(ji,jj)
-               ikv = mikv(ji,jj)
-               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
-               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2
-               ENDIF
-            END DO
-         END DO
-
-         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
-         ! step and store it in  zri, zrj for each  case
-         CALL eos( zti, zhi, zri )  
-         CALL eos( ztj, zhj, zrj )
-
-         ! Gradient of density at the last level 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
                iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1
                ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1
-               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
-               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-               IF( ze3wu >= 0._wp ) THEN
-                 pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)
-                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1
-                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1 
-                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  &
-                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
-                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
-               ELSE
-                 pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
-                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2
-                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2
-                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   &
-                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
-                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN
-                 pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 
-                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1
-                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1
-                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  & 
-                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
-                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
-                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco
-               ELSE
-                 pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)
-                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2
-                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2
-                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   &
-                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
-                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
-               ENDIF
+               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
+               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv) 
+
+               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) )          ! i: 1
+               ELSE                      ; pgrui(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) )      ! i: 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1
+               ELSE                      ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2
+               ENDIF
+
             END DO
          END DO
          CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions
          !
       END IF  
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5189)
@@ -171,4 +171,5 @@
             END DO
          END DO
+         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
          IF ( ln_isfcav ) THEN
             DO jj = 2, jpjm1
@@ -202,7 +203,7 @@
                END DO
             END DO
+            CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition
          END IF
          !
-         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
          !
          IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
@@ -295,4 +296,19 @@
          bfrua(:,:) = - bfrcoef2d(:,:)
          bfrva(:,:) = - bfrcoef2d(:,:)
+         !
+         IF ( ln_isfcav ) THEN
+            IF(ln_tfr2d) THEN
+               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+               CALL iom_open('tfr_coef.nc',inum)
+               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
+               CALL iom_close(inum)
+               tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+            ELSE
+               tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable
+            ENDIF
+            !
+            tfrua(:,:) = - tfrcoef2d(:,:)
+            tfrva(:,:) = - tfrcoef2d(:,:)
+         END IF
          !
       CASE( 2 )
@@ -336,4 +352,17 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
+         
+         IF ( ln_isfcav ) THEN
+            IF(ln_tfr2d) THEN
+               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+               CALL iom_open('tfr_coef.nc',inum)
+               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! bfrcoef2d is used as tmp array
+               CALL iom_close(inum)
+               !
+               tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+            ELSE
+               tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable
+            ENDIF
+         END IF
          !
          IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all
@@ -346,4 +375,14 @@
                END DO
             END DO
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ikbt = mikt(ji,jj)
+                     ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp
+                     tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
+                     tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max)
+                  END DO
+               END DO
+            END IF
          ENDIF
          !
@@ -398,4 +437,28 @@
              zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  )
              zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  )
+! (ISF)
+             IF ( ln_isfcav ) THEN
+                ikbu = miku(ji,jj)       ! deepest ocean level at u- and v-points
+                ikbv = mikv(ji,jj)
+                zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt
+                zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt
+                IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN
+                   IF( ln_ctl ) THEN
+                      WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu
+                      WRITE(numout,*) 'BFR ', ABS( tfrcoef2d(ji,jj) ), zfru
+                   ENDIF
+                   ictu = ictu + 1
+                ENDIF
+                IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN
+                   IF( ln_ctl ) THEN
+                      WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv
+                      WRITE(numout,*) 'BFR ', tfrcoef2d(ji,jj), zfrv
+                   ENDIF
+                   ictv = ictv + 1
+                ENDIF
+                zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  )
+                zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  )
+             END IF
+! END ISF
          END DO
       END DO
@@ -405,4 +468,6 @@
          CALL mpp_min( zminbfr )
          CALL mpp_max( zmaxbfr )
+         IF ( ln_isfcav) CALL mpp_min( zmintfr )
+         IF ( ln_isfcav) CALL mpp_max( zmaxtfr )
       ENDIF
       IF( .NOT.ln_bfrimp) THEN
@@ -411,6 +476,6 @@
          WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
          WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
-         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
          WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
+         IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
       ENDIF
       ENDIF
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 5189)
@@ -31,5 +31,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer
 
    REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
@@ -111,11 +111,10 @@
       END DO
       !
-      ! w-level of the turbocline
+      ! w-level of the turbocline and mixing layer (iom_use)
       imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
       DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
          DO jj = 1, jpj
             DO ji = 1, jpi
-               imkt = mikt(ji,jj)
-               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline 
+               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline 
             END DO
          END DO
@@ -127,12 +126,18 @@
             iikn = nmln(ji,jj)
             imkt = mikt(ji,jj)
-            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! Turbocline depth 
-            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth
-            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
+            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
+            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth
+            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
          END DO
       END DO
       IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode
-         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
-         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth
+         IF ( iom_use("mldr10_1") ) THEN
+            IF( .NOT. ln_isfcav ) CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
+            IF(       ln_isfcav ) CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
+         END IF
+         IF ( iom_use("mldkz5") ) THEN
+            IF( .NOT. ln_isfcav ) CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
+            IF(       ln_isfcav ) CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
+         END IF
       ENDIF
       
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 5189)
@@ -45,7 +45,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtsu, gtsv   !: horizontal gradient of T, S bottom u-point 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gru , grv    !: horizontal gradient of rd at bottom u-point
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aru , arv    
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzu , gzv    
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3ru, ge3rv   !: horizontal gradient of T, S and rd at top v-point  
 
    !! (ISF) interpolated gradient (only used for ice shelf case) 
@@ -53,7 +50,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtui, gtvi   !: horizontal gradient of T, S and rd at top u-point 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   grui, grvi   !: horizontal gradient of T, S and rd at top v-point  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   arui, arvi   !: horizontal average  of rd          at top v-point  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzui, gzvi   !: horizontal gradient of z           at top v-point  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3rui, ge3rvi   !: horizontal gradient of T, S and rd at top v-point  
    !! (ISF) ice load
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   riceload
@@ -104,12 +98,6 @@
          &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       &
          &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     &
-         &     aru(jpi,jpj)      , arv(jpi,jpj)      ,                     &
-         &     gzu(jpi,jpj)      , gzv(jpi,jpj)      ,                     &
          &     gru(jpi,jpj)      , grv(jpi,jpj)      ,                     &
-         &     ge3ru(jpi,jpj)    , ge3rv(jpi,jpj)    ,                     &
          &     gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts),                     &
-         &     arui(jpi,jpj)     , arvi(jpi,jpj)     ,                     &
-         &     gzui(jpi,jpj)     , gzvi(jpi,jpj)     ,                     &
-         &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     &
          &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     &
          &     riceload(jpi,jpj),                             STAT=ierr(2) )
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5189)
@@ -150,7 +150,6 @@
             &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
          IF( ln_zps .AND.       ln_isfcav)                               &
-            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-            &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-            &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
+            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
          IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
                          CALL ldf_slp_grif( kstp )
@@ -181,11 +180,10 @@
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
-            IF( ln_zps .AND. .NOT. ln_isfcav)                               &
-               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
-               &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
-            IF( ln_zps .AND.       ln_isfcav)                               &
-               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND. .NOT. ln_isfcav )                               &
+               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient
+               &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND.       ln_isfcav )                                          &
+               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+               &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
 
                                   ua(:,:,:) = 0.e0             ! set dynamics trends to zero
@@ -266,7 +264,6 @@
                &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
             IF( ln_zps .AND.       ln_isfcav)                                &
-               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+               &                                           rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
       ELSE                                                  ! centered hpg  (eos then time stepping)
          IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
@@ -276,7 +273,6 @@
                &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
          IF( ln_zps .AND.       ln_isfcav)                                   & 
-               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+               &                                           rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
          ENDIF
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 5189)
@@ -86,5 +86,5 @@
             &            CALL zps_hde    ( kstp, jptra, trn, gtru, gtrv )   ! Partial steps: now horizontal gradient of passive
          IF( ln_zps .AND.        ln_isfcav)        &
-            &            CALL zps_hde_isf( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! Partial steps: now horizontal gradient of passive
+            &            CALL zps_hde_isf( kstp, jptra, trn, gtru, gtrv, gtrui, gtrvi )  ! Partial steps: isf case
                                                                 ! tracers at the bottom ocean level
          !
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 5188)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 5189)
@@ -146,5 +146,5 @@
         &    CALL zps_hde    ( nit000, jptra, trn, gtru, gtrv  )  ! Partial steps: before horizontal gradient
       IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav )   &
-        &    CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )       ! tracers at the bottom ocean level
+        &    CALL zps_hde_isf( nit000, jptra, trn, gtru, gtrv, gtrui, gtrvi )       ! tracers at the bottom ocean level
 
 
