Index: NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13469)
+++ NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13470)
@@ -253,5 +253,5 @@
          IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg
             CALL wad_spg( sshn, zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
                   &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
@@ -260,5 +260,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 * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
@@ -268,5 +268,5 @@
       ENDIF
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
           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)
@@ -281,5 +281,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)
@@ -287,5 +287,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)
@@ -299,5 +299,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_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
@@ -305,5 +305,5 @@
       ELSE
          zztmp = r1_rau0 * 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_n(ji,jj)
             zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
@@ -443,10 +443,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 )
                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 )
                zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
                     &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
@@ -508,5 +508,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) - rdtbt * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
@@ -533,5 +533,5 @@
          ! 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 +553,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)
@@ -571,5 +571,5 @@
          ! 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 +580,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 +598,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)   & 
                          &     + rdtbt * (                   zu_spg(ji,jj)   &
@@ -613,5 +613,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 +637,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 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
                   va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
@@ -703,5 +703,5 @@
       IF (ln_bt_fw) THEN
          IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) 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 +734,5 @@
       ELSE
          ! At this stage, ssha 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) * ssha(ji  ,jj)      &
@@ -969,5 +969,5 @@
       ! Max courant number for ext. grav. waves
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
          zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
@@ -1093,5 +1093,5 @@
          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_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
                     &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp  
@@ -1099,5 +1099,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_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
                     &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
@@ -1110,5 +1110,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  )
@@ -1119,5 +1119,5 @@
       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_n(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
@@ -1152,5 +1152,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)   )      &
@@ -1171,5 +1171,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
@@ -1194,5 +1194,5 @@
       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) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) )
             z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) )
@@ -1207,5 +1207,5 @@
          !         
       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)
@@ -1218,5 +1218,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)
@@ -1228,5 +1228,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  ) &
@@ -1262,5 +1262,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 
@@ -1273,5 +1273,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
@@ -1301,5 +1301,5 @@
       !!----------------------------------------------------------------------
       !
-      DO_2D_11_10
+      DO_2D( 1, 1, 1, 0 )
          IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
          ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
@@ -1309,5 +1309,5 @@
       END_2D
       !
-      DO_2D_10_11
+      DO_2D( 1, 0, 1, 1 )
          IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
          ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
@@ -1331,5 +1331,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
       !!----------------------------------------------------------------------
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
               &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
@@ -1397,10 +1397,10 @@
       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) )
@@ -1412,5 +1412,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)    
@@ -1420,5 +1420,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)    
@@ -1430,5 +1430,5 @@
       IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
          zztmp = -1._wp / rdtbt
-         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_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
@@ -1438,5 +1438,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_n(ji,jj) * 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_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
@@ -1450,5 +1450,5 @@
          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)
@@ -1458,5 +1458,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)
@@ -1468,5 +1468,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_n(ji,jj) * 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_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
