Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/DYN/dynspg_ts.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/DYN/dynspg_ts.F90	(revision 13540)
@@ -87,4 +87,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -161,4 +162,5 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
       REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes
+      REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
       !
       REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True. 
@@ -227,6 +229,11 @@
       !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends)
       !                                   !  ---------------------------  !
-      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)
+      DO jk = 1 , jpk
+         ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
+         ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+      END DO
+      !
+      zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)
+      zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm)
       !
       !
@@ -250,5 +257,5 @@
       zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
       !
-      CALL dyn_cor_2d( hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
+      CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
          &                                                                     zu_trd, zv_trd   )   ! ==>> out
       !
@@ -257,5 +264,5 @@
          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_00_00
+            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
@@ -264,5 +271,5 @@
             END_2D
          ELSE                                      ! now suface pressure gradient
-            DO_2D_00_00
+            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) 
@@ -272,5 +279,5 @@
       ENDIF
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend
           zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
           zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
@@ -284,5 +291,5 @@
       IF( ln_apr_dyn ) THEN
          IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_frc(ji,jj) = zu_frc(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
                zv_frc(ji,jj) = zv_frc(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
@@ -290,5 +297,5 @@
          ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
             zztmp = grav * r1_2
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  &
                     &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
@@ -302,5 +309,5 @@
       !                                   !  ----------------------------------  !
       IF( ln_bt_fw ) THEN                        ! Add wind forcing
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
@@ -308,5 +315,5 @@
       ELSE
          zztmp = r1_rho0 * r1_2
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
@@ -468,10 +475,10 @@
             !
             !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk)
-            DO_2D_11_10
+            DO_2D( 1, 1, 1, 0 )   ! not jpi-column
                zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
                     &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
                     &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
             END_2D
-            DO_2D_10_11
+            DO_2D( 1, 0, 1, 1 )   ! not jpj-row
                zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
                     &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
@@ -508,5 +515,5 @@
          !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
          !-------------------------------------------------------------------------!
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
             ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
@@ -514,4 +521,10 @@
          !
          CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp )
+         !
+         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
+         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
+#if defined key_agrif
+         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
+#endif
          !
          !                             ! Sum over sub-time-steps to compute advective velocities
@@ -525,13 +538,8 @@
          END IF
          !
-         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
-         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
-#if defined key_agrif
-         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
-#endif
          !  
          ! Sea Surface Height at u-,v-points (vvl case only)
          IF( .NOT.ln_linssh ) THEN                                
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
                   &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
@@ -553,5 +561,5 @@
          !                             ! Surface pressure gradient
          zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
             zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
@@ -567,9 +575,9 @@
          ! at each time step. We however keep them constant here for optimization.
          ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
-         CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
+         CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
          !
          ! Add tidal astronomical forcing if defined
          IF ( ln_tide .AND. ln_tide_pot ) THEN
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
@@ -580,5 +588,5 @@
 !jth do implicitly instead
          IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
@@ -598,5 +606,5 @@
          !------------------------------------------------------------------------------------------------------------------------!
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
                          &     + rDt_e * (                   zu_spg(ji,jj)   &
@@ -613,5 +621,5 @@
             !
          ELSE                           !* Flux form
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
                !                    ! backward interpolated depth used in spg terms at jn+1/2
@@ -637,5 +645,5 @@
 !jth implicit bottom friction:
          IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                   ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
                   va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
@@ -643,9 +651,12 @@
          ENDIF
        
-         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
+         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)
             hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
             hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
             hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
+         ENDIF
+         !
+         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
             CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
                  &                         , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
@@ -654,6 +665,4 @@
             CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  )
          ENDIF
-         !
-         !
          !                                                 ! open boundaries
          IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
@@ -703,5 +712,5 @@
       IF (ln_bt_fw) THEN
          IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zun_save = un_adv(ji,jj)
                zvn_save = vn_adv(ji,jj)
@@ -734,5 +743,5 @@
       ELSE
          ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
                &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      &
@@ -891,21 +900,21 @@
          !                                   ! ---------------
          IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN    !* Read the restart file
-            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
-            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
-            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
-            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
+            CALL iom_get( numror, jpdom_auto, 'ub2_b'  , ub2_b  (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
+            CALL iom_get( numror, jpdom_auto, 'vb2_b'  , vb2_b  (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 
+            CALL iom_get( numror, jpdom_auto, 'un_bf'  , un_bf  (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
+            CALL iom_get( numror, jpdom_auto, 'vn_bf'  , vn_bf  (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 
             IF( .NOT.ln_bt_av ) THEN
-               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
-               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'sshbb_e'  , sshbb_e(:,:), cd_type = 'T', psgn =  1._wp, ldxios = lrxios )   
+               CALL iom_get( numror, jpdom_auto, 'ubb_e'    ,   ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
+               CALL iom_get( numror, jpdom_auto, 'vbb_e'    ,   vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'sshb_e'   ,  sshb_e(:,:), cd_type = 'T', psgn =  1._wp, ldxios = lrxios ) 
+               CALL iom_get( numror, jpdom_auto, 'ub_e'     ,    ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
+               CALL iom_get( numror, jpdom_auto, 'vb_e'     ,    vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
             ENDIF
 #if defined key_agrif
             ! Read time integrated fluxes
             IF ( .NOT.Agrif_Root() ) THEN
-               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'ub2_i_b'  , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
+               CALL iom_get( numror, jpdom_auto, 'vb2_i_b'  , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
             ENDIF
 #endif
@@ -966,5 +975,5 @@
       ! Max courant number for ext. grav. waves
       !
-      DO_2D_11_11
+      DO_2D( 0, 0, 0, 0 )
          zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
          zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
@@ -972,5 +981,5 @@
       END_2D
       !
-      zcmax = MAXVAL( zcu(:,:) )
+      zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
       CALL mpp_max( 'dynspg_ts', zcmax )
 
@@ -1088,8 +1097,8 @@
       !
       SELECT CASE( nvor_scheme )
-      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
+      CASE( np_EEN )                != EEN scheme using e3f energy & enstrophy scheme
          SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
                     &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp  
@@ -1097,5 +1106,5 @@
             END_2D
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      &
                     &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   &
@@ -1108,5 +1117,5 @@
          !
          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
             ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
@@ -1115,7 +1124,7 @@
          END_2D
          !
-      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
+      CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme
          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-         DO_2D_01_01
+         DO_2D( 0, 1, 0, 1 )
             z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
             ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
@@ -1150,5 +1159,5 @@
             !
             !zhf(:,:) = hbatf(:,:)
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )
                zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
                     &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
@@ -1169,5 +1178,5 @@
          CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
          ! JC: TBC. hf should be greater than 0 
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
          END_2D
@@ -1179,5 +1188,5 @@
 
 
-   SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
+   SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE dyn_cor_2d  ***
@@ -1187,23 +1196,23 @@
       INTEGER  ::   ji ,jj                             ! dummy loop indices
       REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: phu, phv, punb, pvnb, zhU, zhV
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pht, phu, phv, punb, pvnb, zhU, zhV
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
       !!----------------------------------------------------------------------
       SELECT CASE( nvor_scheme )
       CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
             z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
             zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
-               &               * (  e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
-               &                  + e1e2t(ji  ,jj)*ht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
+               &               * (  e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
+               &                  + e1e2t(ji  ,jj)*pht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
                !
             zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
-               &               * (  e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
-               &                  + e1e2t(ji,jj  )*ht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
+               &               * (  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_00_00
+         DO_2D( 0, 0, 0, 0 )
             zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
             zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
@@ -1216,5 +1225,5 @@
          !
       CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
               &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
@@ -1226,5 +1235,5 @@
          !
       CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
              &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
@@ -1260,5 +1269,5 @@
       !
       IF( ln_wd_dl_rmp ) THEN     
-         DO_2D_11_11
+         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 
@@ -1271,5 +1280,5 @@
          END_2D
       ELSE  
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
             ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
@@ -1299,5 +1308,5 @@
       !!----------------------------------------------------------------------
       !
-      DO_2D_11_10
+      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)  
@@ -1307,5 +1316,5 @@
       END_2D
       !
-      DO_2D_10_11
+      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)  
@@ -1329,5 +1338,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
       !!----------------------------------------------------------------------
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
               &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
@@ -1396,12 +1405,12 @@
       !                    !==  Set the barotropic drag coef.  ==!
       !
-      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
          
-         DO_2D_00_00
+         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) )
             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
          END_2D
       ELSE                          ! bottom friction only
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
@@ -1413,5 +1422,5 @@
       IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
             ikbv = mbkv(ji,jj)    
@@ -1421,5 +1430,5 @@
       ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
          
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
             ikbv = mbkv(ji,jj)    
@@ -1431,5 +1440,5 @@
       IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
          zztmp = -1._wp / rDt_e
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             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  )
@@ -1439,5 +1448,5 @@
       ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
          
-         DO_2D_00_00
+         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)
             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
@@ -1447,9 +1456,9 @@
       !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
       !
-      IF( ln_isfcav ) THEN
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN
          !
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
             
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
                iktv = mikv(ji,jj)
@@ -1459,5 +1468,5 @@
          ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
             
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
                iktv = mikv(ji,jj)
@@ -1469,5 +1478,5 @@
          !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
          
-         DO_2D_00_00
+         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)
             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
