Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5203)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5204)
@@ -87,13 +87,11 @@
  
   SUBROUTINE sbc_isf(kt)
-    INTEGER, INTENT(in)          ::   kt         ! ocean time step
-    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
-    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
-    REAL(wp)                     ::   rmin
-    REAL(wp)                     ::   zhk
-    CHARACTER(len=256)           ::   cfisf, cvarzisf, cvarhisf   ! name for isf file
-    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
-    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
-    INTEGER           ::   ios           ! Local integer output status for namelist read
+    INTEGER, INTENT(in) :: kt         ! ocean time step
+    INTEGER             :: ji, jj, jk, inum, ierror
+    INTEGER             :: ikt, ikb   ! top and bottom level of the isf boundary layer
+    REAL(wp)            :: zhk
+    CHARACTER(len=256)  :: cvarzisf, cvarhisf   ! name for isf file
+    CHARACTER(LEN=32 )  :: cvarLeff                    ! variable name for efficient Length scale
+    INTEGER             :: ios           ! Local integer output status for namelist read
       !
       !!---------------------------------------------------------------------
@@ -125,4 +123,6 @@
          IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
          IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
+         IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
+         IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
          IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
          IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
@@ -150,7 +150,6 @@
             !: read effective lenght (BG03)
             IF (nn_isf == 2) THEN
-               ! Read Data and save some integral values
+               cvarLeff = 'soLeff' 
                CALL iom_open( sn_Leff_isf%clname, inum )
-               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
                CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
                CALL iom_close(inum)
@@ -237,6 +236,6 @@
             CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
             CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
-            CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
-            CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
+            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
+            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
             ! iom print
             CALL iom_put('ttbl',ttbl(:,:))
@@ -296,6 +295,6 @@
          ! 
          ! output
-         CALL iom_put('qisf'  , qisf)
-         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
+         IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf)
+         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl / soce )
       END IF
   
@@ -343,13 +342,10 @@
     INTEGER, INTENT ( in ) :: kt
 
-    INTEGER :: ji, jj, jk, jish  !temporary integer
-    INTEGER :: ijkmin
-    INTEGER :: ii, ij, ik 
-    INTEGER :: inum
-
-    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
-    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
-    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
-    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
+    INTEGER :: ji, jj, jk !temporary integer
+
+    REAL(wp) :: zt_sum    ! sum of the temperature between 200m and 600m
+    REAL(wp) :: zt_ave    ! averaged temperature between 200m and 600m
+    REAL(wp) :: zt_frz    ! freezing point temperature at depth z
+    REAL(wp) :: zpress    ! pressure to compute the freezing point in depth
     
     !!----------------------------------------------------------------------
@@ -360,23 +356,23 @@
     DO ji = 1, jpi
        DO jj = 1, jpj
-          ik = misfkt(ji,jj)
+          jk = misfkt(ji,jj)
           !! Initialize arrays to 0 (each step)
           zt_sum = 0.e0_wp
-          IF ( ik .GT. 1 ) THEN
-    ! 3. -----------the average temperature between 200m and 600m ---------------------
+          IF ( jk .GT. 1 ) THEN
+    ! 1. -----------the average temperature between 200m and 600m ---------------------
              DO jk = misfkt(ji,jj),misfkb(ji,jj)
              ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
              ! after verif with UNESCO, wrong sign in BG eq. 2
              ! Calculate freezing temperature
-                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
-                zt_frz = eos_fzp(tsb(ji,jj,ik,jp_sal), zpress) 
-                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
+                zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
+                zt_frz = eos_fzp(tsb(ji,jj,jk,jp_sal), zpress) 
+                zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
              ENDDO
              zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
     
-    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
+    ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
           ! For those corresponding to zonal boundary    
              qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
-                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
+                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk) 
              
              fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                  
@@ -407,161 +403,139 @@
       INTEGER, INTENT(in)          ::   kt         ! ocean time step
       !
-      LOGICAL :: ln_isomip = .true.
-      REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti
-      REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d 
-      !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz,zpress,zti
+      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
       REAL(wp) ::   zlamb1, zlamb2, zlamb3
       REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
       REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
-      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
-      REAL(wp) ::   zgammat, zgammas
-      REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==!
+      REAL(wp) ::   zeps = 1.e-20_wp        
+      REAL(wp) ::   zerr
       INTEGER  ::   ji, jj     ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
-      INTEGER  ::   ierror     ! return error code
-      LOGICAL  ::   lit=.TRUE.
       INTEGER  ::   nit
+      LOGICAL  ::   lit
       !!---------------------------------------------------------------------
       !
       ! coeficient for linearisation of tfreez
-      zlamb1=-0.0575
-      zlamb2=0.0901
+      zlamb1=-0.0575_wp
+      zlamb2= 0.0901_wp
       zlamb3=-7.61e-04
       IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
       !
-      CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
-
-      zcfac=0.0_wp 
-      IF (ln_conserve)  zcfac=1.0_wp
-      zpress(:,:)=0.0_wp
-      zgammat2d(:,:)=0.0_wp
-      zgammas2d(:,:)=0.0_wp
-      !
-      !
+      CALL wrk_alloc( jpi,jpj, zfrz  , zpress, zti     , zgammat, zgammas )
+      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx, zhtflx_b                   )
+
+      ! initialisation
+      IF (ln_conserve) THEN ;  zcfac=1.0_wp
+      ELSE                  ;  zcfac=0.0_wp
+      ENDIF
+      zpress (:,:)=0.0_wp
+      zgammat(:,:)=rn_gammat0 ; zgammas (:,:)=rn_gammas0;
+      zhtflx (:,:)=0.0_wp     ; zhtflx_b(:,:)=0.0_wp    ;
+      zfwflx (:,:)=0.0_wp
+      !
+      ! Crude approximation for pressure (but commonly used)
+      ! 1e-04 to convert from Pa to dBar
 !CDIR COLLAPSE
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ! Crude approximation for pressure (but commonly used)
-            ! 1e-04 to convert from Pa to dBar
             zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
-            !
          END DO
       END DO
 
-! Calculate in-situ temperature (ref to surface)
-      zti(:,:)=tinsitu( ttbl, stbl, zpress )
-! Calculate freezing temperature
-      zfrz(:,:)=eos_fzp( sss_m(:,:), zpress )
-
-      
-      zhtflx=0._wp ; zfwflx=0._wp
-      IF (nn_isfblk == 1) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (mikt(ji,jj) > 1 ) THEN
-                  nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp
-                  DO WHILE ( lit )
-! compute gamma
-                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
-! zhtflx is upward heat flux (out of ocean)
-                     zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
-! zwflx is upward water flux
-                     zfwflx = - zhtflx/lfusisf
-! test convergence and compute gammat
-                     IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
-
-                     nit = nit + 1
-                     IF (nit .GE. 100) THEN
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx
-                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-                     END IF
-! save gammat and compute zhtflx_b
-                     zgammat2d(ji,jj)=zgammat
-                     zhtflx_b = zhtflx
-                  END DO
-
-                  qisf(ji,jj) = - zhtflx
-! For genuine ISOMIP protocol this should probably be something like
-                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
-               ELSE
-                  fwfisf(ji,jj) = 0._wp
-                  qisf(ji,jj)   = 0._wp
-               END IF
-            !
+      ! Calculate in-situ temperature (ref to surface)
+      zti(:,:) = tinsitu( ttbl, stbl, zpress )
+
+      ! Calculate freezing temperature
+      zfrz(:,:) = eos_fzp( sss_m, zpress )
+
+      nit = 1 ; lit = .TRUE.
+      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
+         IF (nn_isfblk == 1) THEN ! ISOMIP case
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+            
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
+                  zfwflx(ji,jj) = - zhtflx(ji,jj)/lfusisf
+               END DO
             END DO
-         END DO
-
-      ELSE IF (nn_isfblk == 2 ) THEN
-
-! More complicated 3 equation thermodynamics as in MITgcm
-!CDIR COLLAPSE
-         DO jj = 2, jpj
-            DO ji = 2, jpi
-               IF (mikt(ji,jj) > 1 ) THEN
-                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
-                  DO WHILE ( lit )
-                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
-
-                     zeps1=rcp*rau0*zgammat
-                     zeps2=lfusisf*rau0*zgammas
-                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
-                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
-                     zeps6=zeps4-zti(ji,jj)
-                     zeps7=zeps4-tsurf
-                     zaqe=zlamb1 * (zeps1 + zeps3)
-                     zaqer=0.5/zaqe
-                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
-                     zcqe=zeps2*stbl(ji,jj)
-                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
-! Presumably zdis can never be negative because gammas is very small compared to gammat
-                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
-                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
-                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
+
+            ! Compute heat flux and upward fresh water flux
+            ! For genuine ISOMIP protocol this should probably be something like
+            qisf  (:,:) = - zhtflx(:,:)                                 * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * ( soce / MAX(stbl(:,:),zeps)) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         ELSE IF (nn_isfblk == 2 ) THEN
+            ! More complicated 3 equation thermodynamics as in MITgcm
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! compute coeficient to solve the 2nd order equation
+                  zeps1=rcp*rau0*zgammat(ji,jj)
+                  zeps2=lfusisf*rau0*zgammas(ji,jj)
+                  zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)
+                  zeps4=zlamb2+zlamb3*risfdep(ji,jj)
+                  zeps6=zeps4-zti(ji,jj)
+                  zeps7=zeps4-tsurf
+                  zaqe=zlamb1 * (zeps1 + zeps3)
+                  zaqer=0.5/MIN(zaqe,-zeps)
+                  zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
+                  zcqe=zeps2*stbl(ji,jj)
+                  zdis=zbqe*zbqe-4.0*zaqe*zcqe               
+
+                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
+                  ! compute s freeze
+                  IF (zsfrz .GE. 0.0_wp) THEN ; zsfrz=(-zbqe-SQRT(zdis))*zaqer
+                  ELSE                        ; zsfrz=(-zbqe+SQRT(zdis))*zaqer
+                  ENDIF
+
+                  ! compute t freeze
+                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
   
-! zfwflx is upward water flux
-                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
-! zhtflx is upward heat flux (out of ocean)
+                  ! zfwflx is upward water flux
+                  ! zhtflx is upward heat flux (out of ocean)
 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
-                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
-! zwflx is upward water flux
 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
-                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
-! test convergence and compute gammat
-                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
-
-                     nit = nit + 1
-                     IF (nit .GE. 51) THEN
-                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
-                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
-                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-                     END IF
-! save gammat and compute zhtflx_b
-                     zgammat2d(ji,jj)=zgammat
-                     zgammas2d(ji,jj)=zgammas
-                     zhtflx_b = zhtflx
-
-                  END DO
-! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
-                  qisf(ji,jj) = - zhtflx 
-! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
-                  fwfisf(ji,jj) = zfwflx 
-               ELSE
-                  fwfisf(ji,jj) = 0._wp
-                  qisf(ji,jj)   = 0._wp
-               ENDIF
-               !
+                  ! compute the upward water and heat flux
+                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
+                  zhtflx(ji,jj) = ( zgammat(ji,jj) * rau0 - zcfac * zfwflx(ji,jj) ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
+                  zfwflx(ji,jj) = ( zgammas(ji,jj) * rau0 - zcfac * zfwflx(ji,jj) ) * (zsfrz - stbl(ji,jj)) / MAX(stbl(ji,jj),zeps)
+               END DO
             END DO
-         END DO
-      ENDIF
-      ! lbclnk
-      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
-      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
+
+            ! compute heat and water flux
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         ENDIF
+
+         ! define if we need to iterate
+         IF ( nn_gammablk .LT. 2 ) THEN ; lit = .FALSE.
+         ELSE                           
+            ! check total number of iteration
+            IF (nit .GE. 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+            ELSE                   ; nit = nit + 1
+            ENDIF
+
+            ! compute error between 2 iterations
+            ! if needed save gammat and compute zhtflx_b for next iteration
+            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
+            IF ( zerr .LE. 0.01 ) THEN ; lit = .FALSE.
+            ELSE                       ; zhtflx_b(:,:) = zhtflx(:,:)
+            ENDIF
+         END IF
+      END DO
+      !
       ! output
-      CALL iom_put('isfgammat', zgammat2d)
-      CALL iom_put('isfgammas', zgammas2d)
-         !
-      !CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zqisf, zfwfisf  )
-      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
+      IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat)
+      IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas)
+      ! 
+      CALL wrk_dealloc( jpi,jpj, zfrz  , zpress, zti     , zgammat, zgammas )
+      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx, zhtflx_b                   )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
@@ -569,5 +543,5 @@
    END SUBROUTINE sbc_isf_cav
 
-   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
+   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf )
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange for heat flux
@@ -578,10 +552,10 @@
       !!                Jenkins et al., 2010, JPO, p2298-2312
       !!---------------------------------------------------------------------
-      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
-      INTEGER , INTENT(in)    :: ji,jj
-      LOGICAL , INTENT(inout) :: lit
+      REAL(wp), DIMENSION(:,:), INTENT(out) :: gt, gs
+      REAL(wp), DIMENSION(:,:), INTENT(in ) :: zqhisf, zqwisf
 
       INTEGER  :: ikt                 ! loop index
-      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
+      INTEGER  :: ji,jj
+      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity
       REAL(wp) :: zdku, zdkv                 ! U, V shear 
       REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
@@ -596,102 +570,96 @@
       REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
       REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
-      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
+      REAL(wp) ::   zeps     = 1.0e-20_wp        ! conversion: mm/s ==> m/s
       REAL(wp), DIMENSION(2) :: zts, zab
       !!---------------------------------------------------------------------
-      !
-      IF( nn_gammablk == 0 ) THEN
+      CALL wrk_alloc( jpi,jpj, zustar )
+      !
+      IF      ( nn_gammablk == 0 ) THEN
       !! gamma is constant (specified in namelist)
-         gt = rn_gammat0
-         gs = rn_gammas0
-         lit = .FALSE.
+         gt(:,:) = rn_gammat0
+         gs(:,:) = rn_gammas0
+
       ELSE IF ( nn_gammablk == 1 ) THEN
       !! gamma is assume to be proportional to u* 
-      !! WARNING in case of Losh 2008 tbl parametrization, 
-      !! you have to used the mean value of u in the boundary layer) 
-      !! not yet coded
       !! Jenkins et al., 2010, JPO, p2298-2312
-         ikt = mikt(ji,jj)
-      !! Compute U and V at T points
-   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
-   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
-          zut = utbl(ji,jj)
-          zvt = vtbl(ji,jj)
 
       !! compute ustar
-         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
-      !! Compute mean value over the TBL
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) )
 
       !! Compute gammats
-         gt = zustar * rn_gammat0
-         gs = zustar * rn_gammas0
-         lit = .FALSE.
+         gt(:,:) = zustar(:,:) * rn_gammat0
+         gs(:,:) = zustar(:,:) * rn_gammas0
+      
       ELSE IF ( nn_gammablk == 2 ) THEN
       !! gamma depends of stability of boundary layer
-      !! WARNING in case of Losh 2008 tbl parametrization, 
-      !! you have to used the mean value of u in the boundary layer) 
-      !! not yet coded
       !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
       !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
+      !! compute ustar
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) )
+
+      !! compute Pr and Sc number (can be improved)
+         zPr =   13.8
+         zSc = 2432.0
+
+      !! compute gamma mole
+         zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
+         zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
+
+      !! compute gamma
+         DO ji=2,jpi
+            DO jj=2,jpj
                ikt = mikt(ji,jj)
 
-      !! Compute U and V at T points
-               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
-               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
-
-      !! compute ustar
-               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
-               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
-                 gt = rn_gammat0
-                 gs = rn_gammas0
+               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think
+                  gt = rn_gammat0
+                  gs = rn_gammas0
                ELSE
       !! compute Rc number (as done in zdfric.F90)
-               zcoef = 0.5 / fse3w(ji,jj,ikt)
-               !                                            ! shear of horizontal velocity
-               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
-                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
-               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
-                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
-               !                                            ! richardson number (minimum value set to zero)
-               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
-
-      !! compute Pr and Sc number (can be improved)
-               zPr =   13.8
-               zSc = 2432.0
-
-      !! compute gamma mole
-               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
-               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
+                  zcoef = 0.5 / fse3w(ji,jj,ikt)
+                  !                                            ! shear of horizontal velocity
+                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
+                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
+                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
+                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
+                  !                                            ! richardson number (minimum value set to zero)
+                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
 
       !! compute bouyancy 
-               zts(jp_tem) = ttbl(ji,jj)
-               zts(jp_sal) = stbl(ji,jj)
-               zdep        = fsdepw(ji,jj,ikt)
-               !
-               CALL eos_rab( zts, zdep, zab )
+                  zts(jp_tem) = ttbl(ji,jj)
+                  zts(jp_sal) = stbl(ji,jj)
+                  zdep        = fsdepw(ji,jj,ikt)
+                  !
+                  CALL eos_rab( zts, zdep, zab )
                   !
       !! compute length scale 
-               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+                  zbuofdep = grav * ( zab(jp_tem) * zqhisf(ji,jj) - zab(jp_sal) * zqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       !! compute Monin Obukov Length
-               ! Maximum boundary layer depth
-               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
-               ! Compute Monin obukhov length scale at the surface and Ekman depth:
-               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
-               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+                  ! Maximum boundary layer depth
+                  zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp
+                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
+                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + epsln))
+                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
 
       !! compute eta* (stability parameter)
-               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
+                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp)))
 
       !! compute the sublayer thickness
-               zhnu = 5 * znu / zustar
+                  zhnu = 5 * znu / zustar(ji,jj)
+
       !! compute gamma turb
-               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
-               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
+                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
+                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
 
       !! compute gammats
-               gt = zustar / (zgturb + zgmolet)
-               gs = zustar / (zgturb + zgmoles)
+                  gt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+                  gs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
                END IF
+            END DO
+         END DO
+         CALL lbc_lnk(gt(:,:),'T',1.)
+         CALL lbc_lnk(gs(:,:),'T',1.)
       END IF
+      CALL wrk_dealloc( jpi,jpj, zustar )
 
    END SUBROUTINE
@@ -701,5 +669,5 @@
       !!                  ***  SUBROUTINE sbc_isf_tbl  ***
       !!
-      !! ** Purpose : compute mean T/S/U/V in the boundary layer 
+      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
       !!
       !!----------------------------------------------------------------------
@@ -726,29 +694,28 @@
       DO jj = 2,jpj
          DO ji = 2,jpi
-            IF (ssmask(ji,jj) == 1) THEN
-               ikt = mkt(ji,jj)
-               ikb = mkb(ji,jj)
-
-               ! level fully include in the ice shelf boundary layer
-               DO jk = ikt, ikb - 1
-                  ze3 = fse3t_n(ji,jj,jk)
-                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
-                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
-                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
-                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
-                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
-               END DO
-
-               ! level partially include in ice shelf boundary layer 
-               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
-               IF (cptin == 'T') &
-                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
-               IF (cptin == 'U') &
-                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
-               IF (cptin == 'V') &
-                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
-            END IF
+            ikt = mkt(ji,jj)
+            ikb = mkb(ji,jj)
+
+            ! level fully include in the ice shelf boundary layer
+            DO jk = ikt, ikb - 1
+               ze3 = fse3t_n(ji,jj,jk)
+               IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
+               IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
+                  &                                                       * r1_hisf_tbl(ji,jj) * ze3
+               IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
+                  &                                                       * r1_hisf_tbl(ji,jj) * ze3
+            END DO
+
+            ! level partially include in ice shelf boundary layer 
+            zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
+            IF (cptin == 'T') &
+                &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+            IF (cptin == 'U') &
+                &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
+            IF (cptin == 'V') &
+                &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
          END DO
       END DO
+      varout(:,:) = varout(:,:) * ssmask(:,:)
 
       CALL wrk_dealloc( jpi,jpj, mkt, mkb )      
@@ -777,7 +744,6 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt, ikb 
-      INTEGER  ::   nk_isf
-      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
-      REAL(wp)     ::   zfact     ! local scalar
+      REAL(wp) ::   zhk
+      REAL(wp) ::   zfact     ! local scalar
       !!----------------------------------------------------------------------
       !
@@ -795,5 +761,5 @@
                ! test on tmask useless ?????
                DO jk = ikt, mbkt(ji,jj)
-!                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -839,6 +805,4 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
       REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
-!      REAL(wp) :: fsatg
-!      REAL(wp) :: pfps, pfpt, pfphp 
       REAL(wp) :: zt, zs, zp, zh, zq, zxk
       INTEGER  :: ji, jj            ! dummy loop indices
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 5203)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5204)
@@ -268,6 +268,7 @@
          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 )    ! if level first is a p-step, ik.m1=1
+
+               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)
@@ -279,5 +280,5 @@
                   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) )
+                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                           ! case 2
                   zmaxu = -ze3wu / fse3w(ji,jj,iku)
@@ -285,5 +286,5 @@
                   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) )
+                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
                ENDIF
                !
@@ -294,5 +295,5 @@
                   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) )
+                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                           ! case 2
                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
@@ -300,6 +301,7 @@
                   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
+                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
+               ENDIF
+
             END DO
          END DO
@@ -314,4 +316,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
+
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
@@ -337,4 +340,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
+
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
@@ -342,12 +346,14 @@
                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
+               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
+               ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
+               ELSE                        ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( 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
          !
@@ -357,6 +363,6 @@
          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
+               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
@@ -366,36 +372,38 @@
                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+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) )
+                  zmaxu = ze3wu / fse3w(ji+1,jj,ikup1)
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
+                  ! gradient of tracers
+                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
+               ELSE                           ! case 2
+                  zmaxu = - ze3wu / fse3w(ji,jj,ikup1)
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,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) )
+                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( 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.
+                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikvp1)
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
+                  ! gradient of tracers
+                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
+               ELSE                           ! case 2
+                  zmaxv =  - ze3wv / fse3w(ji,jj,ikvp1)
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
+                  ! gradient of tracers
+                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( 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
@@ -403,7 +411,8 @@
       ! horizontal derivative of density anomalies (rd)
       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
-         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ;
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
+         pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp;
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
                iku = miku(ji,jj)
                ikv = mikv(ji,jj)
@@ -430,20 +439,21 @@
          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
+
+               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 ; 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
+               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
+               ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1
+               ELSE                      ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( 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
          !
       END IF  
