Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90	(revision 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90	(revision 5969)
@@ -20,5 +20,5 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean
+   USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean
    USE sbcrnf          ! river runoff 
    USE sbcisf          ! ice shelf
@@ -89,9 +89,9 @@
       END DO
       !
-      IF( ln_rnf                     )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field)
+      IF( ln_rnf )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field)
       !
-      IF( ln_divisf .AND. nn_isf > 0 )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field)
+      IF( ln_isf )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field)
       !
-      CALL lbc_lnk( hdivn, 'T', 1. )                       !==  lateral boundary cond.  ==!   (no sign change)
+      CALL lbc_lnk( hdivn, 'T', 1. )                !==  lateral boundary cond.  ==!   (no sign change)
       !
       IF( nn_timing == 1 )  CALL timing_stop('div_hor')
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 5969)
@@ -231,5 +231,5 @@
       IF(  ioptio > 1  .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   &
            &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme' )
-      IF( lk_dynspg_exp .AND. ln_isfcav )   &
+      IF( ln_dynspg_exp .AND. ln_isfcav )   &
            &   CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' )
       !
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5969)
@@ -231,5 +231,4 @@
                   END DO
                END DO
-            END DO
             CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
                DO jj = 1, jpjm1
@@ -244,6 +243,5 @@
             END SELECT
             CALL lbc_lnk( zwz, 'F', 1._wp )
-            zwz(:,:) = ff(:,:) * zwz(:,:)
-
+            !
             ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
             DO jj = 2, jpj
@@ -572,10 +570,10 @@
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
                DO ji = 2, fs_jpim1   ! Vector opt.
-                  zwx(ji,jj) = z1_2 * ssumask(ji,jj)  * r1_e12u(ji,jj)     &
-                     &              * ( e12t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
-                     &              +   e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
-                  zwy(ji,jj) = z1_2 * ssvmask(ji,jj)  * r1_e12v(ji,jj)     &
-                     &              * ( e12t(ji,jj  ) * zsshp2_e(ji,jj  )  &
-                     &              +   e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+                  zwx(ji,jj) = z1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
+                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
+                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
+                  zwy(ji,jj) = z1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
+                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
+                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
                END DO
             END DO
@@ -650,10 +648,10 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1      ! NO Vector Opt.
-                  zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e12u(ji,jj)    &
-                     &              * ( e12t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
-                     &              +   e12t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
-                  zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e12v(ji,jj)    &
-                     &              * ( e12t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
-                     &              +   e12t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
+                  zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
+                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
+                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
+                  zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
+                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
+                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
                END DO
             END DO
@@ -689,8 +687,8 @@
             DO jj = 2, jpjm1                            
                DO ji = 2, jpim1
-                  zx1 = z1_2 * ssumask(ji  ,jj) *  r1_e12u(ji  ,jj)    &
+                  zx1 = z1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
                      &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
                      &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
-                  zy1 = z1_2 * ssvmask(ji  ,jj) *  r1_e12v(ji  ,jj  )  &
+                  zy1 = z1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
                      &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
                      &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
@@ -765,6 +763,6 @@
          !
          ! Add top stresses:
-         zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * zun_e(:,:) * hur_e(:,:)
-         zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
+         zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:)
+         zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
          !
          ! Surface pressure trend:
@@ -780,5 +778,5 @@
          !
          ! Set next velocities:
-         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form
+         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN 	! Vector form
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
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 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5969)
@@ -82,6 +82,5 @@
       ENDIF
       
-      !                             ! Vertical momentum advection at uw and vw-points
-      DO jk = 2, jpk                
+      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
          DO jj = 2, jpj                   ! vertical fluxes 
             DO ji = fs_2, jpi             ! vector opt.
@@ -89,6 +88,6 @@
             END DO
          END DO
-         DO jj = 2, jpjm1                 ! interior vertical momentum advection at w-point
-            DO ji = fs_2, fs_jpim1
+         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) )
@@ -96,12 +95,25 @@
          END DO   
       END DO
-      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
+      !
+      ! 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 jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
@@ -194,9 +206,5 @@
          DO jj = 2, jpj                   
             DO ji = fs_2, jpi             ! vector opt.
-<<<<<<< .working
-               zww(ji,jj,jk) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk)
-=======
                zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
->>>>>>> .merge-right.r5936
             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 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5969)
@@ -118,4 +118,5 @@
       REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
       REAL(wp) ::   zck, zfk,      zbw             !   -      -
+      REAL(wp) ::   zdepu, zdepv                   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zslpml_hmlpu, zslpml_hmlpv
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww
@@ -177,4 +178,5 @@
       CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml
 
+
       ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
       ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
@@ -183,8 +185,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)  &
-                  &                                  - 0.5 * ( risfdep(ji,jj) + risfdep(ji+1,jj  ) )
-               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)  &
-                  &                                  - 0.5 * ( risfdep(ji,jj) + risfdep(ji  ,jj+1) )
+               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)       &
+                  &                                  - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji+1,jj  ) ) )
+               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)       &
+                  &                                  - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji  ,jj+1) ) )
             END DO
          END DO
@@ -193,5 +195,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
-               zslpml_hmlpu(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
+               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
             END DO
          END DO
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5969)
@@ -91,5 +91,5 @@
          IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' )
          !
-         area = glob_sum( e12t(:,:) * tmask(:,:,1))           ! interior global domain surface
+         area = glob_sum( e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface
          ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
          ! and in case of no melt, it can generate HSSW.
@@ -133,5 +133,5 @@
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
-            a_fwb   = glob_sum( e12t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
+            a_fwb   = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
@@ -159,6 +159,6 @@
             ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
             !
-            zsurf_neg = glob_sum( e12t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
-            zsurf_pos = glob_sum( e12t(:,:)*ztmsk_pos(:,:) )
+            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
+            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
             !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges) 
             z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
@@ -172,9 +172,9 @@
             ENDIF
             !
-            zsum_fwf   = glob_sum( e12t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
+            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
 !!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
             z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
             !                                                  ! weight to respect erp field 2D structure 
-            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e12t(:,:) )
+            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
             z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
             !                                                  ! final correction term to apply
@@ -191,10 +191,10 @@
                IF( z_fwf < 0._wp ) THEN
                   WRITE(numout,*)'   z_fwf < 0'
-                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv'
+                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
                ELSE
                   WRITE(numout,*)'   z_fwf >= 0'
-                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv'
+                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
                ENDIF
-               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e12t(:,:) )*1.e-9,' Sv'
+               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
                WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
                WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
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 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5969)
@@ -388,5 +388,5 @@
                ! For those corresponding to zonal boundary    
                qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
-                           & * r1_e12t(ji,jj) * tmask(ji,jj,jk)
+                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
              
                fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                  
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 5968)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5969)
@@ -234,5 +234,4 @@
             END DO
             IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
-            ! PM: IF needed to avoid erasing surface value (jk=1)
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
@@ -333,5 +332,5 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
@@ -348,5 +347,5 @@
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
                            &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
-                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk)
+                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * wmask(ji,jj,jk) / fse3w(ji,jj,jk)
                      END DO
                   END DO
@@ -356,5 +355,5 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      &
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk)                      &
                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
