Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynhpg.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynhpg.F90	(revision 14065)
@@ -44,5 +44,5 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE lbclnk          ! lateral boundary condition
+   USE lbclnk          ! lateral boundary condition 
    USE lib_mpp         ! MPP library
    USE eosbn2          ! compute density
@@ -154,5 +154,4 @@
       !!
       INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
-      REAL(wp) ::   znad
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zts_top, zrhd   ! hypothesys on isf density
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  zrhdtop_isf    ! density at bottom of ISF
@@ -206,5 +205,5 @@
       !
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
-      !
+      ! 
       IF(lwp) THEN
          WRITE(numout,*)
@@ -219,5 +218,5 @@
          WRITE(numout,*)
       ENDIF
-      !
+      !                          
       IF ( ln_hpg_djc ) THEN
          IF (ln_hpg_djc_vnh) THEN ! Von Neumann boundary condition
@@ -265,5 +264,5 @@
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
       REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
+      REAL(wp), DIMENSION(jpi,jpj) ::  zhpi, zhpj
       !!----------------------------------------------------------------------
       !
@@ -273,33 +272,28 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
       ENDIF
-
-      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
-
-      ! Surface value
-      DO_2D( 0, 0, 0, 0 )
+      !
+      zcoef0 = - grav * 0.5_wp            ! Local constant initialization
+      !
+      DO_2D( 0, 0, 0, 0 )                 ! Surface value
          zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
-         ! hydrostatic pressure gradient
-         zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
-         zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
-         ! add to the general momentum trend
-         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
-         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
+         !                                   ! hydrostatic pressure gradient
+         zhpi(ji,jj) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
+         zhpj(ji,jj) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
+         !                                   ! add to the general momentum trend
+         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj)
+         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj)
       END_2D
-
-      !
-      ! interior value (2=<jk=<jpkm1)
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      !
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! interior value (2=<jk=<jpkm1)
          zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
-         ! hydrostatic pressure gradient
-         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
-            &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    &
-            &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
-
-         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
-            &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    &
-            &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
-         ! add to the general momentum trend
-         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
-         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
+         !                                   ! hydrostatic pressure gradient
+         zhpi(ji,jj) = zhpi(ji,jj) + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )  &
+            &                                  - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
+
+         zhpj(ji,jj) = zhpj(ji,jj) + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )  &
+            &                                  - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
+         !                                   ! add to the general momentum trend
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj)
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj)
       END_3D
       !
@@ -411,7 +405,7 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !!
-      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
-      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
+      INTEGER  ::   ji, jj, jk, jii, jjj           ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, ztmp       ! local scalars
+      LOGICAL  ::   ll_tmp1, ll_tmp2               ! local logical variables
       REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zhpi, zhpj
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
@@ -427,7 +421,4 @@
       !
       zcoef0 = - grav * 0.5_wp
-      IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly
-      ELSE                    ;   znad = 1._wp         ! Variable volume: density
-      ENDIF
       !
       IF( ln_wd_il ) THEN
@@ -450,5 +441,5 @@
             zcpx(ji,jj) = 0._wp
           END IF
-
+   
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
@@ -471,52 +462,47 @@
         CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       END IF
-
-      ! Surface value
-      DO_2D( 0, 0, 0, 0 )
-         ! hydrostatic pressure gradient along s-surfaces
-         zhpi(ji,jj,1) =   &
-            &  zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    &
-            &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) &
-            &           * r1_e1u(ji,jj)
-         zhpj(ji,jj,1) =   &
-            &  zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    &
-            &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) &
-            &           * r1_e2v(ji,jj)
-         ! s-coordinate pressure gradient correction
-         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
+      !
+      DO_2D( 0, 0, 0, 0 )              ! Surface value
+         !                                   ! hydrostatic pressure gradient along s-surfaces
+         zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj)                      &
+            &          * (  e3w(ji+1,jj  ,1,Kmm) * rhd(ji+1,jj  ,1)  &
+            &             - e3w(ji  ,jj  ,1,Kmm) * rhd(ji  ,jj  ,1)  )
+         zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj)                      &
+            &          * (  e3w(ji  ,jj+1,1,Kmm) * rhd(ji  ,jj+1,1)  &
+            &             - e3w(ji  ,jj  ,1,Kmm) * rhd(ji  ,jj  ,1)  )
+         !                                   ! s-coordinate pressure gradient correction
+         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) )   &
             &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
-         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
+         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) )   &
             &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
          !
          IF( ln_wd_il ) THEN
             zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
-            zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
+            zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
             zuap = zuap * zcpx(ji,jj)
             zvap = zvap * zcpy(ji,jj)
          ENDIF
-         !
-         ! add to the general momentum trend
+         !                                   ! add to the general momentum trend
          puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap
          pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap
       END_2D
-
-      ! interior value (2=<jk=<jpkm1)
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         ! hydrostatic pressure gradient along s-surfaces
-         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
-            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
-            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
-         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
-            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
-            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
-         ! s-coordinate pressure gradient correction
-         zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
+      !
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior value (2=<jk=<jpkm1)
+         !                                   ! hydrostatic pressure gradient along s-surfaces
+         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)                         &
+            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )  &
+            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  )
+         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)                         &
+            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )  &
+            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  )
+         !                                   ! s-coordinate pressure gradient correction
+         zuap = -zcoef0 * ( rhd  (ji+1,jj  ,jk) + rhd  (ji,jj,jk) ) &
             &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj)
-         zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
+         zvap = -zcoef0 * ( rhd  (ji  ,jj+1,jk) + rhd  (ji,jj,jk) ) &
             &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj)
          !
          IF( ln_wd_il ) THEN
             zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
-            zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
+            zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
             zuap = zuap * zcpx(ji,jj)
             zvap = zvap * zcpy(ji,jj)
@@ -549,5 +535,5 @@
       !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
       !!      iceload is added
-      !!
+      !!      
       !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
       !!----------------------------------------------------------------------
@@ -556,6 +542,8 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !!
-      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      INTEGER  ::   ikt ,  ikti1,  iktj1   ! local integer
+      REAL(wp) ::   ze3w, ze3wi1, ze3wj1   ! local scalars
+      REAL(wp) ::   zcoef0, zuap, zvap     !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj
       REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top
@@ -565,7 +553,5 @@
       zcoef0 = - grav * 0.5_wp   ! Local constant initialization
       !
-      znad=1._wp                 ! To use density and not density anomaly
-      !
-      !                          ! iniitialised to 0. zhpi zhpi
+      !                          ! iniitialised to 0. zhpi zhpi 
       zhpi(:,:,:) = 0._wp   ;   zhpj(:,:,:) = 0._wp
 
@@ -581,54 +567,49 @@
       CALL eos( zts_top, risfdep, zrhdtop_oce )
 
-!==================================================================================
-!===== Compute surface value =====================================================
-!==================================================================================
+      !                     !===========================!
+      !                     !=====  surface value  =====!
+      !                     !===========================!
       DO_2D( 0, 0, 0, 0 )
-         ikt    = mikt(ji,jj)
-         iktp1i = mikt(ji+1,jj)
-         iktp1j = mikt(ji,jj+1)
-         ! 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 * e3w(ji+1,jj,iktp1i,Kmm)                                    &
-            &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
-            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &
-            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-            &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            )
-         zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    &
-            &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
-            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &
-            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-            &                                  + ( risfload(ji,jj+1) - risfload(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 )   &
+         ikt   = mikt(ji  ,jj  )   ;   ze3w   = e3w(ji  ,jj  ,ikt  ,Kmm)
+         ikti1 = mikt(ji+1,jj  )   ;   ze3wi1 = e3w(ji+1,jj  ,ikti1,Kmm)
+         iktj1 = mikt(ji  ,jj+1)   ;   ze3wj1 = e3w(ji  ,jj+1,iktj1,Kmm)
+         !                          ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
+         !                          ! we assume ISF is in isostatic equilibrium
+         zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * (   risfload(ji+1,jj) - risfload(ji,jj)  &
+            &                                       + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhdtop_oce(ji+1,jj) )     &
+            &                                                  - ze3w   * ( rhd(ji  ,jj,ikt  ) + zrhdtop_oce(ji  ,jj) ) )   )
+         zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * (   risfload(ji,jj+1) - risfload(ji,jj)  &
+            &                                       + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhdtop_oce(ji,jj+1) )      &
+            &                                                  - ze3w   * ( rhd(ji,jj  ,ikt  ) + zrhdtop_oce(ji,jj  ) ) )   ) 
+         !                          ! s-coordinate pressure gradient correction (=0 if z coordinate)
+         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) )   &
             &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
-         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
+         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) )   &
             &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
-         ! add to the general momentum trend
+         !                          ! add to the general momentum trend
          puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
          pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
       END_2D
-!==================================================================================
-!===== Compute interior value =====================================================
-!==================================================================================
-      ! interior value (2=<jk=<jpkm1)
+      !   
+      !                     !=============================!
+      !                     !=====  interior values  =====!
+      !                     !=============================!
       DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         ! hydrostatic pressure gradient along s-surfaces
+         ze3w   = e3w(ji  ,jj  ,jk,Kmm)
+         ze3wi1 = e3w(ji+1,jj  ,jk,Kmm)
+         ze3wj1 = e3w(ji  ,jj+1,jk,Kmm)
+         !                          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
-            &           * (  e3w(ji+1,jj,jk,Kmm)                   &
-            &                  * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
-            &              - e3w(ji  ,jj,jk,Kmm)                   &
-            &                  * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
+            &           * (  ze3wi1 * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) * wmask(ji+1,jj,jk)   &
+            &              - ze3w   * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) ) * wmask(ji  ,jj,jk)   )
          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
-            &           * (  e3w(ji,jj+1,jk,Kmm)                   &
-            &                  * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
-            &              - e3w(ji,jj  ,jk,Kmm)                   &
-            &                  * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
-         ! s-coordinate pressure gradient correction
-         zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+            &           * (  ze3wj1 * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) * wmask(ji,jj+1,jk)   &
+            &              - ze3w   * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) ) * wmask(ji,jj  ,jk)   )
+         !                          ! s-coordinate pressure gradient correction
+         zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   &
             &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj)
-         zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+         zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   &
             &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj)
-         ! add to the general momentum trend
+         !                          ! add to the general momentum trend
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
@@ -651,5 +632,5 @@
       !!
       INTEGER  ::   ji, jj, jk          ! dummy loop indices
-      INTEGER  ::   iktb, iktt          ! jk indices at tracer points for top and bottom points
+      INTEGER  ::   iktb, iktt          ! jk indices at tracer points for top and bottom points 
       REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
       REAL(wp) ::   z_grav_10, z1_12
@@ -658,6 +639,5 @@
       LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rhd_opt
-
+ 
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzx, zdzy, zdzz                          ! Primitive grid differences ('delta_xyz')
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdz_i, zdz_j, zdz_k                       ! Harmonic average of primitive grid differences ('d_xyz')
@@ -665,5 +645,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrho_i, zdrho_j, zdrho_k                 ! Harmonic average of primitive rho differences ('d_rho')
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_rho_i, z_rho_j, z_rho_k                 ! Face intergrals
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays 
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
       !!----------------------------------------------------------------------
@@ -688,5 +668,5 @@
             zcpx(ji,jj) = 0._wp
           END IF
-
+   
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
@@ -721,20 +701,12 @@
       z1_12  = 1.0_wp / 12._wp
 
-
-!! To be removed once combined with KERNEL08.  KERNEL08 action removes znad from the hpg module, spg is calculated in spg moduel instead.  The final code will just have rhd(:,:,:), no rhd_opt(:,:,:)
-      IF( .NOT. ln_linssh ) THEN
-         rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option
-      ELSE
-         rhd_opt(:,:,:) = rhd(:,:,:)
-      END IF
-
       !----------------------------------------------------------------------------------------
-      !  1. compute and store elementary vertical differences in provisional arrays
+      !  1. compute and store elementary vertical differences in provisional arrays 
       !----------------------------------------------------------------------------------------
 
 !!bug gm   Not a true bug, but... zdzz=e3w  for zdzx, zdzy verify what it is really
 
-      DO_3D( 1, 1, 1, 1, 2, jpkm1 )
-         zdrhoz(ji,jj,jk) =   rhd_opt    (ji  ,jj  ,jk) - rhd_opt    (ji,jj,jk-1)
+      DO_3D( 1, 1, 1, 1, 2, jpkm1 )  
+         zdrhoz(ji,jj,jk) =   rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1)
          zdzz  (ji,jj,jk) = - gde3w(ji  ,jj  ,jk) + gde3w(ji,jj,jk-1)
       END_3D
@@ -745,9 +717,9 @@
       zep = 1.e-15
 
-!! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk)
+!! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 
       zdrho_k(:,:,:) = 0._wp
       zdz_k  (:,:,:) = 0._wp
 
-      DO_3D( 1, 1, 1, 1, 2, jpk-2 )
+      DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 
          cffw = 2._wp * zdrhoz(ji  ,jj  ,jk) * zdrhoz(ji,jj,jk+1)
          IF( cffw > zep) THEN
@@ -763,5 +735,5 @@
 
 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition
-      zdrho_k(:,:,1) = aco_bc_vrt * ( rhd_opt    (:,:,2) - rhd_opt    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2)
+      zdrho_k(:,:,1) = aco_bc_vrt * ( rhd    (:,:,2) - rhd    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2)
       zdz_k  (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k  (:,:,2)
 
@@ -769,6 +741,6 @@
          IF ( mbkt(ji,jj)>1 ) THEN
             iktb = mbkt(ji,jj)
-            zdrho_k(ji,jj,iktb) = aco_bc_vrt * (     rhd_opt(ji,jj,iktb) -     rhd_opt(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)
-            zdz_k  (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k  (ji,jj,iktb-1)
+            zdrho_k(ji,jj,iktb) = aco_bc_vrt * (     rhd(ji,jj,iktb) -     rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)
+            zdz_k  (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k  (ji,jj,iktb-1) 
          END IF
       END_2D
@@ -778,6 +750,6 @@
       !-------------------------------------------------------------
 
-!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights
-!! rho_k stores grav * FX / rho_0
+!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights  
+!! rho_k stores grav * FX / rho_0  
 
       !--------------------------------------------------------------
@@ -786,7 +758,7 @@
 ! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1)
       DO_2D( 0, 1, 0, 1)
-         z_rho_k(ji,jj,1) =  grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) )                        &
-            &                     * (  rhd_opt(ji,jj,1)                                        &
-            &                     + 0.5_wp * (   rhd_opt    (ji,jj,2) - rhd_opt    (ji,jj,1) ) &
+         z_rho_k(ji,jj,1) =  grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) )                        & 
+            &                     * (  rhd(ji,jj,1)                                        &
+            &                     + 0.5_wp * (   rhd    (ji,jj,2) - rhd    (ji,jj,1) ) &
             &                              * (   ssh   (ji,jj,Kmm) + gde3w(ji,jj,1) )          &
             &                              / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) )  )
@@ -798,5 +770,5 @@
 
       DO_3D( 0, 1, 0, 1, 2, jpkm1 )
-         z_rho_k(ji,jj,jk) = zcoef0 * (   rhd_opt    (ji,jj,jk) + rhd_opt    (ji,jj,jk-1) )                                   &
+         z_rho_k(ji,jj,jk) = zcoef0 * (   rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   &
             &                       * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) )                                               &
             &                       + z_grav_10 * (                                                                           &
@@ -804,20 +776,20 @@
             &   * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k  (ji,jj,jk) + zdz_k  (ji,jj,jk-1) ) )             &
             &   - ( zdz_k    (ji,jj,jk) - zdz_k    (ji,jj,jk-1) )                                                             &
-            &   * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) )   &
+            &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) )   &
             &                             )
       END_3D
 
       !----------------------------------------------------------------------------------------
-      !  5. compute and store elementary horizontal differences in provisional arrays
+      !  5. compute and store elementary horizontal differences in provisional arrays 
       !----------------------------------------------------------------------------------------
 
       DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-         zdrhox(ji,jj,jk) =   rhd_opt    (ji+1,jj  ,jk) - rhd_opt    (ji,jj,jk  )
+         zdrhox(ji,jj,jk) =   rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  )
          zdzx  (ji,jj,jk) = - gde3w(ji+1,jj  ,jk) + gde3w(ji,jj,jk  )
-         zdrhoy(ji,jj,jk) =   rhd_opt    (ji  ,jj+1,jk) - rhd_opt    (ji,jj,jk  )
+         zdrhoy(ji,jj,jk) =   rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  )
          zdzy  (ji,jj,jk) = - gde3w(ji  ,jj+1,jk) + gde3w(ji,jj,jk  )
       END_3D
 
-      CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. )
+      CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 
 
       !-------------------------------------------------------------------------
@@ -854,6 +826,6 @@
          ENDIF
       END_3D
-
-!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point
+      
+!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point      
 
       !----------------------------------------------------------------------------------
@@ -869,6 +841,6 @@
             ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
             IF (ji < jpi) THEN
-               IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp)  THEN
-                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd_opt    (ji+1,jj,jk) - rhd_opt    (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk)
+               IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp)  THEN  
+                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 
                   zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i  (ji+1,jj,jk)
                END IF
@@ -877,5 +849,5 @@
             IF (ji > 2) THEN
                IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN
-                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)
+                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)  
                   zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i  (ji-1,jj,jk)
                END IF
@@ -884,12 +856,12 @@
             IF (jj < jpj) THEN
                IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp)  THEN
-                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
+                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
                   zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j  (ji,jj+1,jk)
                END IF
-            END IF
+            END IF 
             ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
             IF (jj > 2) THEN
-               IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN
-                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)
+               IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 
+                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 
                   zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j  (ji,jj-1,jk)
                END IF
@@ -903,11 +875,11 @@
 
       !--------------------------------------------------------------
-      ! 7. Calculate integrals on side faces
+      ! 7. Calculate integrals on side faces  
       !-------------------------------------------------------------
 
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
 ! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height)
-         z_rho_i(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji+1,jj,jk) + rhd_opt    (ji,jj,jk) )                                       &
-             &                    * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )
+         z_rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                       &
+             &                    * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                    
          IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN
             z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * (                                                               &
@@ -915,10 +887,10 @@
              &   * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i  (ji+1,jj,jk) + zdz_i  (ji,jj,jk) ) )              &
              &   - (   zdz_i    (ji+1,jj,jk) - zdz_i    (ji,jj,jk) )                                                            &
-             &   * (   rhd_opt    (ji+1,jj,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) )  &
+             &   * (   rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) )  &
              &                                               )
          END IF
-
-         z_rho_j(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji,jj+1,jk) + rhd_opt    (ji,jj,jk) )                                       &
-             &                    * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )
+  
+         z_rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                       &
+             &                    * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                  
          IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN
             z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * (                                                               &
@@ -926,5 +898,5 @@
              &   * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j  (ji,jj+1,jk) + zdz_j  (ji,jj,jk) ) )              &
              &   - (   zdz_j    (ji,jj+1,jk) - zdz_j    (ji,jj,jk) )                                                            &
-             &   * (   rhd_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) )  &
+             &   * (   rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) )  &
              &                                                 )
          END IF
@@ -932,7 +904,7 @@
 
       !--------------------------------------------------------------
-      ! 8. Integrate in the vertical
+      ! 8. Integrate in the vertical   
       !-------------------------------------------------------------
-
+      !
       ! ---------------
       !  Surface value
@@ -943,5 +915,5 @@
          IF( ln_wd_il ) THEN
            zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
-           zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
+           zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
          ENDIF
          ! add to the general momentum trend
@@ -963,5 +935,5 @@
          IF( ln_wd_il ) THEN
            zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
-           zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
+           zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
          ENDIF
          ! add to the general momentum trend
@@ -1000,7 +972,8 @@
       REAL(wp) :: zrhdt1
       REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zpgu, zpgv   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zsshu_n, zsshv_n
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdept, zrhh
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
-      REAL(wp), DIMENSION(jpi,jpj)   ::   zsshu_n, zsshv_n
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
       !!----------------------------------------------------------------------
@@ -1015,6 +988,14 @@
       zcoef0 = - grav
       znad = 1._wp
-      IF( ln_linssh )   znad = 0._wp
-
+      IF( ln_linssh )   znad = 1._wp
+      !
+      ! ---------------
+      !  Surface pressure gradient to be removed
+      ! ---------------
+      DO_2D( 0, 0, 0, 0 )
+         zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
+         zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
+      END_2D
+      !
       IF( ln_wd_il ) THEN
          ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
@@ -1034,10 +1015,10 @@
             zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
                         &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) )
-
+           
              zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
           ELSE
             zcpx(ji,jj) = 0._wp
           END IF
-
+   
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
@@ -1082,5 +1063,5 @@
       ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
       DO_2D( 1, 1, 1, 1 )
-         zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad
+         zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm)
       END_2D
 
@@ -1120,12 +1101,12 @@
 !!gm BUG ?    if it is ssh at u- & v-point then it should be:
 !          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * &
-!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
+!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
 !          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * &
-!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
+!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
 !!gm not this:
        zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * &
-                      & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
+                      & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
        zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * &
-                      & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
+                      & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
       END_2D
 
@@ -1133,6 +1114,6 @@
 
       DO_2D( 0, 0, 0, 0 )
-       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)
-       zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad)
+       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 
+       zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) )
       END_2D
 
@@ -1216,5 +1197,5 @@
             zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
          ENDIF
-         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 
       ENDIF
 
@@ -1272,9 +1253,9 @@
          ENDIF
          IF( ln_wd_il ) THEN
-            zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)
-            zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)
-         ENDIF
-
-         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
+            zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 
+            zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 
+         ENDIF
+
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk)
       ENDIF
          !
@@ -1307,5 +1288,5 @@
       !!----------------------------------------------------------------------
       !
-!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!
+!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!  
       jpi   = size(fsp,1)
       jpj   = size(fsp,2)
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg.F90	(revision 14065)
@@ -82,6 +82,7 @@
       INTEGER  ::   ji, jj, jk                   ! dummy loop indices
       REAL(wp) ::   z2dt, zg_2, zintp, zgrho0r, zld   ! local scalars
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
+      REAL(wp)             , DIMENSION(jpi,jpj) ::   zpgu, zpgv   ! 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zpice
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   ztrdu, ztrdv
       !!----------------------------------------------------------------------
       !
@@ -99,6 +100,6 @@
          !
          DO_2D( 0, 0, 0, 0 )
-            spgu(ji,jj) = 0._wp
-            spgv(ji,jj) = 0._wp
+            zpgu(ji,jj) = 0._wp
+            zpgv(ji,jj) = 0._wp
          END_2D
          !
@@ -106,7 +107,7 @@
             zg_2 = grav * 0.5
             DO_2D( 0, 0, 0, 0 )                       ! gradient of Patm using inverse barometer ssh
-               spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
+               zpgu(ji,jj) = zpgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
                   &                                + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
-               spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
+               zpgv(ji,jj) = zpgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
                   &                                + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
             END_2D
@@ -121,6 +122,6 @@
             !
             DO_2D( 0, 0, 0, 0 )                      ! add tide potential forcing
-               spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
-               spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
+               zpgu(ji,jj) = zpgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
+               zpgv(ji,jj) = zpgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
             END_2D
             !
@@ -128,6 +129,6 @@
                zld = rn_scal_load * grav
                DO_2D( 0, 0, 0, 0 )                   ! add scalar approximation for load potential
-                  spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
-                  spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
+                  zpgu(ji,jj) = zpgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
+                  zpgv(ji,jj) = zpgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
                END_2D
             ENDIF
@@ -140,6 +141,6 @@
             zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrho0r
             DO_2D( 0, 0, 0, 0 )
-               spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
-               spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
+               zpgu(ji,jj) = zpgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+               zpgv(ji,jj) = zpgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
             END_2D
             DEALLOCATE( zpice )
@@ -148,12 +149,12 @@
          IF( ln_wave .and. ln_bern_srfc ) THEN          !== Add J terms: depth-independent Bernoulli head
             DO_2D( 0, 0, 0, 0 )
-               spgu(ji,jj) = spgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
-               spgv(ji,jj) = spgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
+               zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
+               zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
             END_2D
          ENDIF
          !
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )       !== Add all terms to the general trend
-            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
-            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
+            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zpgu(ji,jj)
+            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zpgv(ji,jj)
          END_3D
          !
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_exp.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_exp.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_exp.F90	(revision 14065)
@@ -60,4 +60,5 @@
       !!
       INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   zpgu, zpgv   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -67,22 +68,18 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (explicit free surface)'
          !
-         spgu(:,:) = 0._wp   ;   spgv(:,:) = 0._wp
+         zpgu(:,:) = 0._wp   ;   zpgv(:,:) = 0._wp
          !
          IF( .NOT.ln_linssh .AND. lwp ) WRITE(numout,*) '      non linear free surface: spg is included in dynhpg'
       ENDIF
-
-      IF( ln_linssh ) THEN          !* linear free surface : add the surface pressure gradient trend
-         !
-         DO_2D( 0, 0, 0, 0 )                 ! now surface pressure gradient
-            spgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
-            spgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
-         END_2D
-         !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! Add it to the general trend
-            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
-            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
-         END_3D
-         !
-      ENDIF
+      !
+      DO_2D( 0, 0, 0, 0 )
+         zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
+         zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
+      END_2D
+      !
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zpgu(ji,jj)
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zpgv(ji,jj)
+      END_3D
       !
    END SUBROUTINE dyn_spg_exp
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90	(revision 14065)
@@ -1,5 +1,5 @@
 MODULE dynspg_ts
 
-   !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !
+   !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 
 
    !!======================================================================
@@ -23,5 +23,5 @@
 
    !!----------------------------------------------------------------------
-   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme
+   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme 
    !!   dyn_spg_ts_init: initialisation of the time-splitting scheme
    !!   ts_wgt         : set time-splitting weights for temporal averaging (or not)
@@ -50,5 +50,5 @@
    USE agrif_oce
 #endif
-#if defined key_asminc
+#if defined key_asminc   
    USE asminc          ! Assimilation increment
 #endif
@@ -66,5 +66,5 @@
    PRIVATE
 
-   PUBLIC dyn_spg_ts        ! called by dyn_spg
+   PUBLIC dyn_spg_ts        ! called by dyn_spg 
    PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init
 
@@ -122,10 +122,10 @@
       !! ** Purpose : - Compute the now trend due to the explicit time stepping
       !!              of the quasi-linear barotropic system, and add it to the
-      !!              general momentum trend.
+      !!              general momentum trend. 
       !!
       !! ** Method  : - split-explicit schem (time splitting) :
       !!      Barotropic variables are advanced from internal time steps
       !!      "n"   to "n+1" if ln_bt_fw=T
-      !!      or from
+      !!      or from 
       !!      "n-1" to "n+1" if ln_bt_fw=F
       !!      thanks to a generalized forward-backward time stepping (see ref. below).
@@ -136,9 +136,9 @@
       !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
       !!      These are used to advect tracers and are compliant with discrete
-      !!      continuity equation taken at the baroclinic time steps. This
+      !!      continuity equation taken at the baroclinic time steps. This 
       !!      ensures tracers conservation.
       !!      - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component.
       !!
-      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
+      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 
       !!---------------------------------------------------------------------
       INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
@@ -148,5 +148,5 @@
       !
       INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
-      LOGICAL  ::   ll_fw_start           ! =T : forward integration
+      LOGICAL  ::   ll_fw_start           ! =T : forward integration 
       LOGICAL  ::   ll_init               ! =T : special startup of 2d equations
       INTEGER  ::   noffset               ! local integers  : time offset for bdy update
@@ -162,7 +162,9 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
       REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes
+#if defined key_qco 
       REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
-      !
-      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.
+#endif
+      !
+      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True. 
 
       INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point
@@ -179,14 +181,14 @@
       IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
       !
-      zwdramp = r_rn_wdmin1               ! simplest ramp
+      zwdramp = r_rn_wdmin1               ! simplest ramp 
 !     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
-      !                                         ! inverse of baroclinic time step
-      r1_Dt_b = 1._wp / rDt
-      !
-      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart
+      !                                         ! inverse of baroclinic time step 
+      r1_Dt_b = 1._wp / rDt 
+      !
+      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart 
       ll_fw_start = .FALSE.
       !                                         ! time offset in steps for bdy data update
       IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_e
-      ELSE                       ;   noffset =   0
+      ELSE                       ;   noffset =   0 
       ENDIF
       !
@@ -215,5 +217,5 @@
             ! and the averaging weights. We don't have an easy way of telling whether we did
             ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false.
-            ! at the end of the first timestep) so just do this in all cases.
+            ! at the end of the first timestep) so just do this in all cases. 
             ll_fw_start = .FALSE.
             CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
@@ -225,8 +227,9 @@
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
       ! -----------------------------------------------------------------------------
-      !
+      !      
       !
       !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends)
       !                                   !  ---------------------------  !
+#if defined key_qco 
       DO jk = 1 , jpk
          ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
@@ -236,4 +239,8 @@
       zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)
       zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm)
+#else
+      zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)
+      zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm)
+#endif 
       !
       !
@@ -243,39 +250,20 @@
          vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk)
       END DO
-
+      
 !!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
 !!gm  Is it correct to do so ?   I think so...
-
-      !                                   !=  remove 2D Coriolis and pressure gradient trends  =!
-      !                                   !  -------------------------------------------------  !
+      
+      !                                   !=  remove 2D Coriolis trend  =!
+      !                                   !  --------------------------  !
       !
       IF( kt == nit000 .OR. .NOT. ln_linssh )   CALL dyn_cor_2D_init( Kmm )   ! Set zwz, the barotropic Coriolis force coefficient
-      !       ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes
+      !                      ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes
       !
       !                                         !* 2D Coriolis trends
-      zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes
+      zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes 
       zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
       !
       CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
-         &                                                                     zu_trd, zv_trd   )   ! ==>> out
-      !
-      IF( .NOT.ln_linssh ) THEN                 !* surface pressure gradient   (variable volume only)
-         !
-         IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg
-            CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy
-            DO_2D( 0, 0, 0, 0 )                                ! SPG with the application of W/D gravity filters
-               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj ,Kmm) )   &
-                  &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
-               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj ,Kmm) )   &
-                  &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth
-            END_2D
-         ELSE                                      ! now suface pressure gradient
-            DO_2D( 0, 0, 0, 0 )
-               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e1u(ji,jj)
-               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e2v(ji,jj)
-            END_2D
-         ENDIF
-         !
-      ENDIF
+         &                                                                          zu_trd, zv_trd   )   ! ==>> out
       !
       DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend
@@ -287,4 +275,5 @@
       !                                   !  -----------------------------------------------------------  !
       CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc,  zCdU_u, zCdU_v )      ! also provide the barotropic drag coefficients
+      !
       !                                   !=  Add atmospheric pressure forcing  =!
       !                                   !  ----------------------------------  !
@@ -319,5 +308,5 @@
             zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
          END_2D
-      ENDIF
+      ENDIF  
       !
       !              !----------------!
@@ -330,8 +319,8 @@
       ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
          zztmp = r1_rho0 * r1_2
-         zssh_frc(:,:) = zztmp * (  emp(:,:)        + emp_b(:,:)                    &
-                &                 - rnf(:,:)        - rnf_b(:,:)                    &
-                &                 + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)             &
-                &                 + fwfisf_par(:,:) + fwfisf_par_b(:,:)             )
+         zssh_frc(:,:) = zztmp * (   emp(:,:)        + emp_b(:,:)          &
+            &                      - rnf(:,:)        - rnf_b(:,:)          &
+            &                      + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)   &
+            &                      + fwfisf_par(:,:) + fwfisf_par_b(:,:)   )
       ENDIF
       !                                   !=  Add Stokes drift divergence  =!   (if exist)
@@ -369,11 +358,11 @@
       !
       ! -----------------------------------------------------------------------
-      !  Phase 2 : Integration of the barotropic equations
+      !  Phase 2 : Integration of the barotropic equations 
       ! -----------------------------------------------------------------------
       !
       !                                             ! ==================== !
       !                                             !    Initialisations   !
-      !                                             ! ==================== !
-      ! Initialize barotropic variables:
+      !                                             ! ==================== !  
+      ! Initialize barotropic variables:      
       IF( ll_init )THEN
          sshbb_e(:,:) = 0._wp
@@ -391,26 +380,26 @@
       ENDIF
       !
-      IF( ln_bt_fw ) THEN                 ! FORWARD integration: start from NOW fields
-         sshn_e(:,:) =    pssh (:,:,Kmm)
-         un_e  (:,:) =    puu_b(:,:,Kmm)
+      IF( ln_bt_fw ) THEN                 ! FORWARD integration: start from NOW fields                    
+         sshn_e(:,:) =    pssh (:,:,Kmm)            
+         un_e  (:,:) =    puu_b(:,:,Kmm)            
          vn_e  (:,:) =    pvv_b(:,:,Kmm)
          !
-         hu_e  (:,:) =    hu(:,:,Kmm)
-         hv_e  (:,:) =    hv(:,:,Kmm)
-         hur_e (:,:) = r1_hu(:,:,Kmm)
+         hu_e  (:,:) =    hu(:,:,Kmm)       
+         hv_e  (:,:) =    hv(:,:,Kmm) 
+         hur_e (:,:) = r1_hu(:,:,Kmm)    
          hvr_e (:,:) = r1_hv(:,:,Kmm)
       ELSE                                ! CENTRED integration: start from BEFORE fields
          sshn_e(:,:) =    pssh (:,:,Kbb)
-         un_e  (:,:) =    puu_b(:,:,Kbb)
+         un_e  (:,:) =    puu_b(:,:,Kbb)         
          vn_e  (:,:) =    pvv_b(:,:,Kbb)
          !
-         hu_e  (:,:) =    hu(:,:,Kbb)
-         hv_e  (:,:) =    hv(:,:,Kbb)
-         hur_e (:,:) = r1_hu(:,:,Kbb)
+         hu_e  (:,:) =    hu(:,:,Kbb)       
+         hv_e  (:,:) =    hv(:,:,Kbb) 
+         hur_e (:,:) = r1_hu(:,:,Kbb)    
          hvr_e (:,:) = r1_hv(:,:,Kbb)
       ENDIF
       !
       ! Initialize sums:
-      puu_b (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)
+      puu_b (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)          
       pvv_b (:,:,Kaa) = 0._wp
       pssh  (:,:,Kaa) = 0._wp       ! Sum for after averaged sea level
@@ -419,9 +408,9 @@
       !
       IF( ln_wd_dl ) THEN
-         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)
-         zvwdmask(:,:) = 0._wp  !
-         zuwdav2 (:,:) = 0._wp
-         zvwdav2 (:,:) = 0._wp
-      END IF
+         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary) 
+         zvwdmask(:,:) = 0._wp  ! 
+         zuwdav2 (:,:) = 0._wp 
+         zvwdav2 (:,:) = 0._wp   
+      END IF 
 
       !                                             ! ==================== !
@@ -443,9 +432,9 @@
          !
          !                       !* Set extrapolation coefficients for predictor step:
-         IF ((jn<3).AND.ll_init) THEN      ! Forward
-           za1 = 1._wp
-           za2 = 0._wp
-           za3 = 0._wp
-         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
+         IF ((jn<3).AND.ll_init) THEN      ! Forward           
+           za1 = 1._wp                                          
+           za2 = 0._wp                        
+           za3 = 0._wp                        
+         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105 
            za1 =  1.781105_wp              ! za1 =   3/2 +   bet
            za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
@@ -467,5 +456,5 @@
             !--------------------------------------------------------------------------------!
             zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
-
+            
             ! set wetting & drying mask at tracer points for this barotropic mid-step
             IF( ln_wd_dl )   CALL wad_tmsk( zsshp2_e, ztwdmask )
@@ -495,5 +484,5 @@
                     &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
             END_2D
-#endif
+#endif               
             !
          ENDIF
@@ -504,5 +493,5 @@
          !                             ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
          IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
-         !
+         !      
          !                             ! resulting flux at mid-step (not over the full domain)
          zhU(1:jpim1,1:jpj  ) = e2u(1:jpim1,1:jpj  ) * ua_e(1:jpim1,1:jpj  ) * zhup2_e(1:jpim1,1:jpj  )   ! not jpi-column
@@ -515,9 +504,9 @@
          IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e)    !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
 
-         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where
+         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where 
             !                          ! the direction of the flow is from dry cells
             CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask )   ! not jpi colomn for U, not jpj row for V
             !
-         ENDIF
+         ENDIF    
          !
          !
@@ -543,5 +532,5 @@
          un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
          vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
-         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
+         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 
          IF ( ln_wd_dl_bc ) THEN
             zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
@@ -549,5 +538,5 @@
          END IF
          !
-         !
+         !  
          ! Sea Surface Height at u-,v-points (vvl case only)
          IF( .NOT.ln_linssh ) THEN
@@ -569,5 +558,5 @@
 #endif
          ENDIF
-         !
+         !         
          ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
          !--            m+1/2           m+1              m               m-1              m-2     --!
@@ -626,8 +615,8 @@
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
             DO_2D( 0, 0, 0, 0 )
-               ua_e(ji,jj) = (                                 un_e(ji,jj)   &
+               ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
                          &     + rDt_e * (                   zu_spg(ji,jj)   &
                          &                                 + zu_trd(ji,jj)   &
-                         &                                 + zu_frc(ji,jj) ) &
+                         &                                 + zu_frc(ji,jj) ) & 
                          &   ) * ssumask(ji,jj)
 
@@ -657,5 +646,5 @@
                z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
                !
-               ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      &
+               ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
                     &            + rDt_e * (  zhu_bck        * zu_spg (ji,jj)  &   !
                     &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   !
@@ -675,5 +664,5 @@
             END_2D
          ENDIF
-
+       
          IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
             hu_e (2:jpim1,2:jpjm1) =    hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
@@ -692,5 +681,5 @@
          !                                                 ! open boundaries
          IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
-#if defined key_agrif
+#if defined key_agrif                                                           
          IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
 #endif
@@ -711,16 +700,16 @@
          !                                             !* Sum over whole bt loop
          !                                             !  ----------------------
-         za1 = wgtbtp1(jn)
+         za1 = wgtbtp1(jn)                                    
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
-            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:)
-            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:)
+            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) 
+            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) 
          ELSE                                       ! Sum transports
-            IF ( .NOT.ln_wd_dl ) THEN
+            IF ( .NOT.ln_wd_dl ) THEN  
                puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:)
                pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:)
-            ELSE
+            ELSE 
                puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
                pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
-            END IF
+            END IF 
          ENDIF
          !                                          ! Sum sea level
@@ -740,5 +729,5 @@
                zun_save = un_adv(ji,jj)
                zvn_save = vn_adv(ji,jj)
-               !                          ! apply the previously computed correction
+               !                          ! apply the previously computed correction 
                un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
                vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
@@ -781,5 +770,5 @@
                &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
          END_2D
-#endif
+#endif   
          CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
          !
@@ -796,5 +785,5 @@
 
 
-      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
+      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)  
       DO jk = 1, jpkm1
          puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk)
@@ -802,14 +791,14 @@
       END DO
 
-      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
+      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 
          DO jk = 1, jpkm1
             puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
-                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
-            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
-                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
+                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 
+            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
+                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)  
          END DO
-      END IF
-
-
+      END IF 
+
+      
       CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
       CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
@@ -829,5 +818,5 @@
          vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
       ENDIF
-#endif
+#endif      
       !                                   !* write time-spliting arrays in the restart
       IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
@@ -850,8 +839,8 @@
       LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
       LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
-      INTEGER, INTENT(inout) :: jpit      ! cycle length
+      INTEGER, INTENT(inout) :: jpit      ! cycle length    
       REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, & ! Primary weights
                                                          zwgt2    ! Secondary weights
-
+      
       INTEGER ::  jic, jn, ji                      ! temporary integers
       REAL(wp) :: za1, za2
@@ -862,5 +851,5 @@
 
       ! Set time index when averaged value is requested
-      IF (ll_fw) THEN
+      IF (ll_fw) THEN 
          jic = nn_e
       ELSE
@@ -870,6 +859,6 @@
       ! Set primary weights:
       IF (ll_av) THEN
-           ! Define simple boxcar window for primary weights
-           ! (width = nn_e, centered around jic)
+           ! Define simple boxcar window for primary weights 
+           ! (width = nn_e, centered around jic)     
          SELECT CASE ( nn_bt_flt )
               CASE( 0 )  ! No averaging
@@ -879,5 +868,5 @@
               CASE( 1 )  ! Boxcar, width = nn_e
                  DO jn = 1, 3*nn_e
-                    za1 = ABS(float(jn-jic))/float(nn_e)
+                    za1 = ABS(float(jn-jic))/float(nn_e) 
                     IF (za1 < 0.5_wp) THEN
                       zwgt1(jn) = 1._wp
@@ -888,5 +877,5 @@
               CASE( 2 )  ! Boxcar, width = 2 * nn_e
                  DO jn = 1, 3*nn_e
-                    za1 = ABS(float(jn-jic))/float(nn_e)
+                    za1 = ABS(float(jn-jic))/float(nn_e) 
                     IF (za1 < 1._wp) THEN
                       zwgt1(jn) = 1._wp
@@ -901,5 +890,5 @@
          jpit = jic
       ENDIF
-
+    
       ! Set secondary weights
       DO jn = 1, jpit
@@ -930,17 +919,17 @@
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
          !                                   ! ---------------
          IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN    !* Read the restart file
-            CALL iom_get( numror, jpdom_auto, 'ub2_b'  , ub2_b  (:,:), cd_type = 'U', psgn = -1._wp )
-            CALL iom_get( numror, jpdom_auto, 'vb2_b'  , vb2_b  (:,:), cd_type = 'V', psgn = -1._wp )
-            CALL iom_get( numror, jpdom_auto, 'un_bf'  , un_bf  (:,:), cd_type = 'U', psgn = -1._wp )
-            CALL iom_get( numror, jpdom_auto, 'vn_bf'  , vn_bf  (:,:), cd_type = 'V', psgn = -1._wp )
+            CALL iom_get( numror, jpdom_auto, 'ub2_b'  , ub2_b  (:,:), cd_type = 'U', psgn = -1._wp )   
+            CALL iom_get( numror, jpdom_auto, 'vb2_b'  , vb2_b  (:,:), cd_type = 'V', psgn = -1._wp ) 
+            CALL iom_get( numror, jpdom_auto, 'un_bf'  , un_bf  (:,:), cd_type = 'U', psgn = -1._wp )   
+            CALL iom_get( numror, jpdom_auto, 'vn_bf'  , vn_bf  (:,:), cd_type = 'V', psgn = -1._wp ) 
             IF( .NOT.ln_bt_av ) THEN
-               CALL iom_get( numror, jpdom_auto, 'sshbb_e'  , sshbb_e(:,:), cd_type = 'T', psgn =  1._wp )
-               CALL iom_get( numror, jpdom_auto, 'ubb_e'    ,   ubb_e(:,:), cd_type = 'U', psgn = -1._wp )
+               CALL iom_get( numror, jpdom_auto, 'sshbb_e'  , sshbb_e(:,:), cd_type = 'T', psgn =  1._wp )   
+               CALL iom_get( numror, jpdom_auto, 'ubb_e'    ,   ubb_e(:,:), cd_type = 'U', psgn = -1._wp )   
                CALL iom_get( numror, jpdom_auto, 'vbb_e'    ,   vbb_e(:,:), cd_type = 'V', psgn = -1._wp )
-               CALL iom_get( numror, jpdom_auto, 'sshb_e'   ,  sshb_e(:,:), cd_type = 'T', psgn =  1._wp )
-               CALL iom_get( numror, jpdom_auto, 'ub_e'     ,    ub_e(:,:), cd_type = 'U', psgn = -1._wp )
+               CALL iom_get( numror, jpdom_auto, 'sshb_e'   ,  sshb_e(:,:), cd_type = 'T', psgn =  1._wp ) 
+               CALL iom_get( numror, jpdom_auto, 'ub_e'     ,    ub_e(:,:), cd_type = 'U', psgn = -1._wp )   
                CALL iom_get( numror, jpdom_auto, 'vb_e'     ,    vb_e(:,:), cd_type = 'V', psgn = -1._wp )
             ENDIF
@@ -948,5 +937,5 @@
             ! Read time integrated fluxes
             IF ( .NOT.Agrif_Root() ) THEN
-               CALL iom_get( numror, jpdom_auto, 'ub2_i_b'  , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp )
+               CALL iom_get( numror, jpdom_auto, 'ub2_i_b'  , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp )   
                CALL iom_get( numror, jpdom_auto, 'vb2_i_b'  , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp )
             ELSE
@@ -974,5 +963,5 @@
          !
          IF (.NOT.ln_bt_av) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
             CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
             CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
@@ -1017,5 +1006,5 @@
       ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
       IF( ln_bt_auto )   nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
-
+      
       rDt_e = rn_Dt / REAL( nn_e , wp )
       zcmax = zcmax * rDt_e
@@ -1053,5 +1042,5 @@
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_e'
-         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e'
+         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e' 
          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
       END SELECT
@@ -1071,5 +1060,5 @@
       ENDIF
       IF( zcmax>0.9_wp ) THEN
-         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
+         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )          
       ENDIF
       !
@@ -1082,5 +1071,5 @@
    END SUBROUTINE dyn_spg_ts_init
 
-
+   
    SUBROUTINE dyn_cor_2D_init( Kmm )
       !!---------------------------------------------------------------------
@@ -1108,5 +1097,5 @@
             DO_2D( 0, 0, 0, 0 )
                zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1)   &
-                    &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp
+                    &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp  
                IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
             END_2D
@@ -1145,5 +1134,5 @@
          !
       END SELECT
-
+      
    END SUBROUTINE dyn_cor_2d_init
 
@@ -1170,8 +1159,8 @@
                !
             zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
-               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   &
-               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   )
-         END_2D
-         !
+               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
+               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
+         END_2D
+         !         
       CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
          DO_2D( 0, 0, 0, 0 )
@@ -1195,5 +1184,5 @@
          END_2D
          !
-      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)
+      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
          DO_2D( 0, 0, 0, 0 )
             zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
@@ -1215,9 +1204,9 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
-      !!
-      !! ** Purpose :   set wetting & drying mask at tracer points
-      !!              for the current barotropic sub-step
-      !!
-      !! ** Method  :   ???
+      !!                    
+      !! ** Purpose :   set wetting & drying mask at tracer points 
+      !!              for the current barotropic sub-step 
+      !!
+      !! ** Method  :   ??? 
       !!
       !! ** Action  :  ptmsk : wetting & drying t-mask
@@ -1229,16 +1218,16 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_wd_dl_rmp ) THEN
+      IF( ln_wd_dl_rmp ) THEN     
          DO_2D( 1, 1, 1, 1 )
-            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN
-               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN
+            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
+               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN 
                ptmsk(ji,jj) = 1._wp
             ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
                ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
-            ELSE
+            ELSE 
                ptmsk(ji,jj) = 0._wp
             ENDIF
          END_2D
-      ELSE
+      ELSE  
          DO_2D( 1, 1, 1, 1 )
             IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
@@ -1254,9 +1243,9 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
-      !!
-      !! ** Purpose :   set wetting & drying mask at tracer points
-      !!              for the current barotropic sub-step
-      !!
-      !! ** Method  :   ???
+      !!                    
+      !! ** Purpose :   set wetting & drying mask at tracer points 
+      !!              for the current barotropic sub-step 
+      !!
+      !! ** Method  :   ??? 
       !!
       !! ** Action  :  ptmsk : wetting & drying t-mask
@@ -1270,6 +1259,6 @@
       !
       DO_2D( 1, 1, 1, 0 )   ! not jpi-column
-         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj)
-         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)
+         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
+         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
          ENDIF
          phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
@@ -1279,7 +1268,7 @@
       DO_2D( 1, 0, 1, 1 )   ! not jpj-row
          IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
-         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)
-         ENDIF
-         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
+         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
+         ENDIF
+         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
          pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
       END_2D
@@ -1292,5 +1281,5 @@
       !!                   ***  ROUTINE  wad_sp  ***
       !!
-      !! ** Purpose :
+      !! ** Purpose : 
       !!----------------------------------------------------------------------
       INTEGER  ::   ji ,jj               ! dummy loop indices
@@ -1325,5 +1314,5 @@
               &      MAX(   pshn(ji,jj)              ,  pshn(ji,jj+1) ) >                &
               &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
-
+         
          IF(ll_tmp1) THEN
             zcpy(ji,jj) = 1.0_wp
@@ -1337,17 +1326,17 @@
          ENDIF
       END_2D
-
+            
    END SUBROUTINE wad_spg
-
+     
 
    SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_drg_init  ***
-      !!
-      !! ** Purpose : - add the baroclinic top/bottom drag contribution to
+      !!                    
+      !! ** Purpose : - add the baroclinic top/bottom drag contribution to 
       !!              the baroclinic part of the barotropic RHS
       !!              - compute the barotropic drag coefficients
       !!
-      !! ** Method  :   computation done over the INNER domain only
+      !! ** Method  :   computation done over the INNER domain only 
       !!----------------------------------------------------------------------
       INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
@@ -1366,5 +1355,5 @@
       !
       IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
-
+         
          DO_2D( 0, 0, 0, 0 )
             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
@@ -1381,16 +1370,16 @@
       !
       IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
-
+         
          DO_2D( 0, 0, 0, 0 )
-            ikbu = mbku(ji,jj)
-            ikbv = mbkv(ji,jj)
+            ikbu = mbku(ji,jj)       
+            ikbv = mbkv(ji,jj)    
             zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
             zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
          END_2D
       ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
-
+         
          DO_2D( 0, 0, 0, 0 )
-            ikbu = mbku(ji,jj)
-            ikbv = mbkv(ji,jj)
+            ikbu = mbku(ji,jj)       
+            ikbv = mbkv(ji,jj)    
             zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
             zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
@@ -1401,11 +1390,11 @@
          zztmp = -1._wp / rDt_e
          DO_2D( 0, 0, 0, 0 )
-            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 &
+            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
                  &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
-            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 &
+            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
                  &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
          END_2D
       ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
-
+         
          DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
@@ -1419,5 +1408,5 @@
          !
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
-
+            
             DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
@@ -1427,5 +1416,5 @@
             END_2D
          ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
-
+            
             DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
@@ -1437,5 +1426,5 @@
          !
          !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
-
+         
          DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
@@ -1458,6 +1447,6 @@
       !                             ! set Half-step back interpolation coefficient
       IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
-         za0 = 1._wp
-         za1 = 0._wp
+         za0 = 1._wp                        
+         za1 = 0._wp                           
          za2 = 0._wp
          za3 = 0._wp
@@ -1466,7 +1455,7 @@
          za1 =-0.1666666666666_wp                 ! za1 = gam
          za2 = 0.0833333333333_wp                 ! za2 = eps
-         za3 = 0._wp
-      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
-         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion
+         za3 = 0._wp              
+      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
+         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion  
             za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
             za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
@@ -1480,5 +1469,5 @@
             za2 = zgamma
             za3 = zepsilon
-         ENDIF
+         ENDIF 
       ENDIF
    END SUBROUTINE ts_bck_interp
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/IOM/iom.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/IOM/iom.F90	(revision 14065)
@@ -2280,5 +2280,5 @@
       !
       CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig0(Nis0)-1,jbegin=mjg0(Njs0)-1,ni=Ni_0,nj=Nj_0)
-      CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
+      CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni=jpi, data_jbegin = -nn_hls, data_nj=jpj)
 !don't define lon and lat for restart reading context.
       IF ( .NOT.ldrxios ) &
@@ -2289,7 +2289,9 @@
          ! mask land points, keep values on coast line -> specific mask for U, V and W points
          SELECT CASE ( cdgrd )
-         CASE('T')   ;   zmask(:,:,:)       = tmask(:,:,:)
-         CASE('U')   ;   zmask(2:jpim1,:,:) = tmask(2:jpim1,:,:) + tmask(3:jpi,:,:)
-         CASE('V')   ;   zmask(:,2:jpjm1,:) = tmask(:,2:jpjm1,:) + tmask(:,3:jpj,:)
+         CASE('T')   ;   zmask( :     , :     ,:) = tmask(:,:,:)
+         CASE('U')   ;   zmask(2:jpim1, :     ,:) = tmask(2:jpim1, :     ,:) + tmask(3:jpi  , :   ,:)
+         CASE('V')   ;   zmask( :     ,2:jpjm1,:) = tmask( :     ,2:jpjm1,:) + tmask( :     ,3:jpj,:)
+         CASE('F')   ;   zmask(2:jpim1,2:jpjm1,:) = tmask(2:jpim1,2:jpjm1,:) + tmask(2:jpim1,3:jpj,:)    &
+            &                                     + tmask(3:jpi  ,2:jpjm1,:) + tmask(3:jpi  ,3:jpj,:)
          CASE('W')   ;   zmask(:,:,2:jpk  ) = tmask(:,:,1:jpkm1) + tmask(:,:,2:jpk)   ;   zmask(:,:,1) = tmask(:,:,1)
          END SELECT
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isf_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isf_oce.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isf_oce.F90	(revision 14065)
@@ -1,6 +1,6 @@
 MODULE isf_oce
    !!======================================================================
-   !!                       ***  MODULE  sbcisf  ***
-   !! Surface module :  compute iceshelf melt and heat flux
+   !!                       ***  MODULE  isf_oce  ***
+   !! Ice shelves :  ice shelves variables defined in memory
    !!======================================================================
    !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
@@ -146,4 +146,5 @@
    END SUBROUTINE isf_alloc_par
 
+   
    SUBROUTINE isf_alloc_cav()
       !!---------------------------------------------------------------------
@@ -173,4 +174,5 @@
    END SUBROUTINE isf_alloc_cav
 
+   
    SUBROUTINE isf_alloc_cpl()
       !!---------------------------------------------------------------------
@@ -184,14 +186,14 @@
       ierr = 0
       !
-      ALLOCATE( risfcpl_ssh(jpi,jpj), risfcpl_tsc(jpi,jpj,jpk,jpts), risfcpl_vol(jpi,jpj,jpk), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      risfcpl_tsc(:,:,:,:) = 0.0 ; risfcpl_vol(:,:,:) = 0.0 ; risfcpl_ssh(:,:) = 0.0
-
-      IF ( ln_isfcpl_cons) THEN
-         ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) ,risfcpl_cons_ssh(jpi,jpj), STAT=ialloc )
+      ALLOCATE( risfcpl_ssh(jpi,jpj) , risfcpl_tsc(jpi,jpj,jpk,jpts) , risfcpl_vol(jpi,jpj,jpk) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      risfcpl_tsc(:,:,:,:) = 0._wp ; risfcpl_vol(:,:,:) = 0._wp ; risfcpl_ssh(:,:) = 0._wp
+
+      IF ( ln_isfcpl_cons ) THEN
+         ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) , risfcpl_cons_ssh(jpi,jpj) , STAT=ialloc )
          ierr = ierr + ialloc
          !
-         risfcpl_cons_tsc(:,:,:,:) = 0.0 ; risfcpl_cons_vol(:,:,:) = 0.0 ; risfcpl_cons_ssh(:,:) = 0.0
+         risfcpl_cons_tsc(:,:,:,:) = 0._wp ; risfcpl_cons_vol(:,:,:) = 0._wp ; risfcpl_cons_ssh(:,:) = 0._wp
          !
       END IF
@@ -202,4 +204,5 @@
    END SUBROUTINE isf_alloc_cpl
 
+   
    SUBROUTINE isf_dealloc_cpl()
       !!---------------------------------------------------------------------
@@ -213,5 +216,5 @@
       ierr = 0
       !
-      DEALLOCATE( risfcpl_ssh, risfcpl_tsc, risfcpl_vol, STAT=ialloc )
+      DEALLOCATE( risfcpl_ssh , risfcpl_tsc , risfcpl_vol , STAT=ialloc )
       ierr = ierr + ialloc
       !
@@ -221,4 +224,5 @@
    END SUBROUTINE isf_dealloc_cpl
 
+   
    SUBROUTINE isf_alloc()
       !!---------------------------------------------------------------------
@@ -233,19 +237,19 @@
       ierr = 0       ! set to zero if no array to be allocated
       !
-      ALLOCATE(fwfisf_par(jpi,jpj)  , fwfisf_par_b(jpi,jpj), &
-         &     fwfisf_cav(jpi,jpj)  , fwfisf_cav_b(jpi,jpj), &
-         &     fwfisf_oasis(jpi,jpj),            STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risf_par_tsc(jpi,jpj,jpts), risf_par_tsc_b(jpi,jpj,jpts), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risf_cav_tsc(jpi,jpj,jpts), risf_cav_tsc_b(jpi,jpj,jpts), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risfload(jpi,jpj), STAT=ialloc)
-      ierr = ierr + ialloc
-      !
-      ALLOCATE( mskisf_cav(jpi,jpj), STAT=ialloc)
+      ALLOCATE( fwfisf_par  (jpi,jpj) , fwfisf_par_b(jpi,jpj) ,     &
+         &      fwfisf_cav  (jpi,jpj) , fwfisf_cav_b(jpi,jpj) ,     &
+         &      fwfisf_oasis(jpi,jpj)                         , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( risfload(jpi,jpj) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( mskisf_cav(jpi,jpj) , STAT=ialloc )
       ierr = ierr + ialloc
       !
@@ -254,13 +258,13 @@
       !
       ! initalisation of fwf and tsc array to 0
-      risfload(:,:)       = 0.0_wp
-      fwfisf_oasis(:,:)   = 0.0_wp
-      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp
-      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp
-      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp
-      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp
-      !
-
+      risfload    (:,:)   = 0._wp
+      fwfisf_oasis(:,:)   = 0._wp
+      fwfisf_par  (:,:)   = 0._wp   ;   fwfisf_par_b  (:,:)   = 0._wp
+      fwfisf_cav  (:,:)   = 0._wp   ;   fwfisf_cav_b  (:,:)   = 0._wp
+      risf_cav_tsc(:,:,:) = 0._wp   ;   risf_cav_tsc_b(:,:,:) = 0._wp
+      risf_par_tsc(:,:,:) = 0._wp   ;   risf_par_tsc_b(:,:,:) = 0._wp
+      !
    END SUBROUTINE isf_alloc
-
+   
+   !!======================================================================
 END MODULE isf_oce
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfload.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfload.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfload.F90	(revision 14065)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE  isfload  ***
-   !! isfload module :  compute ice shelf load (needed for the hpg)
+   !! Ice Shelves :   compute ice shelf load (needed for the hpg)
    !!======================================================================
    !! History :  4.1  !  2019-09  (P. Mathiot) original code
@@ -8,5 +8,5 @@
 
    !!----------------------------------------------------------------------
-   !!   isfload      : compute ice shelf load
+   !!   isf_load      : compute ice shelf load
    !!----------------------------------------------------------------------
 
@@ -23,9 +23,14 @@
    PRIVATE
 
-   PUBLIC isf_load
+   PUBLIC   isf_load   ! called by isfstp.F90
+   !
    !! * Substitutions
 #  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
-
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -37,8 +42,6 @@
       !!
       !!--------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pisfload
-      !!-------------------------- IN  -------------------------------------
-      INTEGER,                      INTENT(in)    :: Kmm           ! ocean time level index
+      INTEGER,                      INTENT(in   ) ::   Kmm        ! ocean time level index      
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pisfload   ! ice shelf load
       !!----------------------------------------------------------------------
       !
@@ -46,5 +49,5 @@
       !               the smaller the residual flow is, the better it is.
       !
-      ! ice shelf cavity
+      ! type of ice shelf cavity
       SELECT CASE ( cn_isfload )
       CASE ( 'uniform' )
@@ -56,5 +59,6 @@
    END SUBROUTINE isf_load
 
-   SUBROUTINE isf_load_uniform( Kmm, pisfload )
+   
+   SUBROUTINE isf_load_uniform( Kmm, pload )
       !!--------------------------------------------------------------------
       !!                  ***  SUBROUTINE isf_load  ***
@@ -67,18 +71,13 @@
       !!
       !!--------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pisfload
-      !!-------------------------- IN  -------------------------------------
-      INTEGER,                      INTENT(in)    :: Kmm           ! ocean time level index
-      !!--------------------------------------------------------------------
+      INTEGER,                      INTENT(in   ) ::   Kmm     ! ocean time level index      
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pload   ! ice shelf load
+      !
       INTEGER  :: ji, jj, jk
       INTEGER  :: ikt
-      REAL(wp)                          :: znad        ! 
       REAL(wp), DIMENSION(jpi,jpj)      :: zrhdtop_isf ! water density    displaced by the ice shelf (at the interface)
       REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts_top     ! water properties displaced by the ice shelf   
       REAL(wp), DIMENSION(jpi,jpj,jpk)  :: zrhd        ! water density    displaced by the ice shelf
       !!----------------------------------------------------------------------
-      !
-      znad = 1._wp                     !- To use density and not density anomaly
       !
       !                                !- assume water displaced by the ice shelf is at T=rn_isfload_T and S=rn_isfload_S (rude)
@@ -87,4 +86,5 @@
       DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
          CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) )
+!!st ==>> CALL eos( zts_top(:,:,:), gdept_0(:,:,jk), zrhd(:,:,jk) )
       END DO
       !
@@ -93,23 +93,21 @@
       !
       !                                !- Surface value + ice shelf gradient
-      pisfload(:,:) = 0._wp                       ! compute pressure due to ice shelf load 
+      pload(:,:) = 0._wp                      ! compute pressure due to ice shelf load 
       DO_2D( 1, 1, 1, 1 )
          ikt = mikt(ji,jj)
          !
          IF ( ikt > 1 ) THEN
+            !                                 ! top layer of the ice shelf
+            pload(ji,jj) = pload(ji,jj)   &
+               &         + zrhd (ji,jj,1) * e3w(ji,jj,1,Kmm)
             !
-            ! top layer of the ice shelf
-            pisfload(ji,jj) = pisfload(ji,jj) + (znad + zrhd(ji,jj,1) )   &
-               &                                * e3w(ji,jj,1,Kmm)
-            !
-            ! core layers of the ice shelf
-            DO jk = 2, ikt-1
-               pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk))   &
-                  &                                * e3w(ji,jj,jk,Kmm)
+            DO jk = 2, ikt-1                  ! core layers of the ice shelf
+               pload(ji,jj) = pload(ji,jj) + (zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk))   &
+                  &                        *   e3w(ji,jj,jk,Kmm)
             END DO
-            !
-            ! deepest part of the ice shelf (between deepest T point and ice/ocean interface
-            pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
-               &                                              * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) )
+            !                                 ! deepest part of the ice shelf (between deepest T point and ice/ocean interface
+            pload(ji,jj) = pload(ji,jj) + ( zrhdtop_isf(ji,jj) +  zrhd(ji,jj,ikt-1)     )   &
+               &                        * (     risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) )
+!!st ==>>     &                        * (     risfdep(ji,jj) - gdept_0(ji,jj,ikt-1) )
             !
          END IF
@@ -117,4 +115,5 @@
       !
    END SUBROUTINE isf_load_uniform
-
+   
+   !!======================================================================
 END MODULE isfload
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfstp.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfstp.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ISF/isfstp.F90	(revision 14065)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE  isfstp  ***
-   !! Surface module :  compute iceshelf load, melt and heat flux
+   !! Ice Shelves :  compute iceshelf load, melt and heat flux
    !!======================================================================
    !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
@@ -42,8 +42,7 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
-
 CONTAINS
  
-  SUBROUTINE isf_stp( kt, Kmm )
+   SUBROUTINE isf_stp( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isf_stp  ***
@@ -58,11 +57,12 @@
       !!              - compute fluxes
       !!              - write restart variables
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
-      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
-      !!----------------------------------------------------------------------
-      INTEGER :: jk                               ! loop index
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t    ! e3t 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt    ! ocean time step
+      INTEGER, INTENT(in) ::   Kmm   ! ocean time level index
+      !
+      INTEGER :: jk                              ! loop index
+#if defined key_qco
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t   ! 3D workspace
+#endif
       !!---------------------------------------------------------------------
       !
@@ -83,11 +83,15 @@
          ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
          rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
+#if defined key_qco
          DO jk = 1, jpk
             ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
          END DO 
-         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         CALL isf_tbl_lvl( ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
+#else
+         CALL isf_tbl_lvl( ht(:,:),  e3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
+#endif
          !
          ! 1.3: compute ice shelf melt
-         CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav)
+         CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav )
          !
       END IF
@@ -108,11 +112,15 @@
          ! by simplicity, we assume the top level where param applied do not change with time (done in init part)
          rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
+#if defined key_qco
          DO jk = 1, jpk
             ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
          END DO
-         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         CALL isf_tbl_lvl( ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+#else
+         CALL isf_tbl_lvl( ht(:,:),  e3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+#endif
          !
          ! 2.3: compute ice shelf melt
-         CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par)
+         CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par )
          !
       END IF
@@ -128,5 +136,6 @@
    END SUBROUTINE isf_stp
 
-   SUBROUTINE isf_init(Kbb, Kmm, Kaa)
+   
+   SUBROUTINE isf_init( Kbb, Kmm, Kaa )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isfstp_init  ***
@@ -142,34 +151,25 @@
       !!              - call cav/param/isfcpl init routine
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa      ! ocean time level indices
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time level indices
+      !!----------------------------------------------------------------------
       !
       ! constrain: l_isfoasis need to be known
       !
-      ! Read namelist
-      CALL isf_nam()
-      !
-      ! Allocate public array
-      CALL isf_alloc()
-      !
-      ! check option compatibility
-      CALL isf_ctl()
-      !
-      ! compute ice shelf load
-      IF ( ln_isfcav ) CALL isf_load( Kmm, risfload )
+      CALL isf_nam()                                              ! Read namelist
+      !
+      CALL isf_alloc()                                            ! Allocate public array
+      !
+      CALL isf_ctl()                                              ! check option compatibility
+      !
+      IF( ln_isfcav ) CALL isf_load( Kmm, risfload )              ! compute ice shelf load
       !
       ! terminate routine now if no ice shelf melt formulation specify
-      IF ( ln_isf ) THEN
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation melt in the cavity
-         IF ( ln_isfcav_mlt ) CALL isf_cav_init()
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation parametrised melt
-         IF ( ln_isfpar_mlt ) CALL isf_par_init()
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation ice sheet coupling
-         IF( ln_isfcpl ) CALL isfcpl_init(Kbb, Kmm, Kaa)
+      IF( ln_isf ) THEN
+         !
+         IF( ln_isfcav_mlt )   CALL isf_cav_init()                ! initialisation melt in the cavity
+         !
+         IF( ln_isfpar_mlt )   CALL isf_par_init()                ! initialisation parametrised melt
+         !
+         IF( ln_isfcpl     )   CALL isfcpl_init( Kbb, Kmm, Kaa )  ! initialisation ice sheet coupling
          !
       END IF
@@ -177,4 +177,5 @@
   END SUBROUTINE isf_init
 
+  
   SUBROUTINE isf_ctl()
       !!---------------------------------------------------------------------
@@ -283,5 +284,6 @@
       END IF
    END SUBROUTINE isf_ctl
-   !
+
+   
    SUBROUTINE isf_nam
       !!---------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/oce.F90	(revision 14063)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/oce.F90	(revision 14065)
@@ -50,6 +50,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_i_b, vb2_i_b         !: Half step time integrated fluxes 
 #endif
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   spgu, spgv               !: horizontal surface pressure gradient
 
    !! interpolated gradient (only used in zps case)
@@ -91,18 +89,17 @@
       !
       ierr(:) = 0 
-      ALLOCATE( uu   (jpi,jpj,jpk,jpt)  , vv   (jpi,jpj,jpk,jpt) ,                              &          
-         &      ww   (jpi,jpj,jpk)      , hdiv(jpi,jpj,jpk)      ,                              &
-         &      ts   (jpi,jpj,jpk,jpts,jpt),                                                    &
-         &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) ,                             &
-         &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)      ,                             &
-         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) )
+      ALLOCATE( uu   (jpi,jpj,jpk,jpt)  , vv   (jpi,jpj,jpk,jpt)           ,     &          
+         &      ww   (jpi,jpj,jpk)      , hdiv(jpi,jpj,jpk)                ,     &
+         &      ts   (jpi,jpj,jpk,jpts,jpt)                                ,     &
+         &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts)          ,     &
+         &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)               ,     &
+         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)               , STAT=ierr(1) )
          !
-      ALLOCATE( ssh(jpi,jpj,jpt)  , uu_b(jpi,jpj,jpt) , vv_b(jpi,jpj,jpt) , &
-         &      spgu  (jpi,jpj)   , spgv(jpi,jpj)                     ,     &
-         &      gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts)                ,     &
-         &      gru(jpi,jpj)      , grv(jpi,jpj)                      ,     &
-         &      gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts)                ,     &
-         &      grui(jpi,jpj)     , grvi(jpi,jpj)                     ,     &
-         &      riceload(jpi,jpj)                                     , STAT=ierr(2) )
+      ALLOCATE( ssh (jpi,jpj,jpt)  , uu_b(jpi,jpj,jpt) , vv_b(jpi,jpj,jpt) ,     &
+         &      gtsu(jpi,jpj,jpts) , gtsv(jpi,jpj,jpts)                    ,     &
+         &      gru (jpi,jpj)      , grv (jpi,jpj)                         ,     &
+         &      gtui(jpi,jpj,jpts) , gtvi(jpi,jpj,jpts)                    ,     &
+         &      grui(jpi,jpj)      , grvi(jpi,jpj)                         ,     &
+         &      riceload(jpi,jpj)                                          , STAT=ierr(2) )
          !
       ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(3) )
