Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/trazdf.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/trazdf.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/trazdf.F90	(revision 13710)
@@ -35,4 +35,6 @@
    PUBLIC   tra_zdf_imp   ! called by trczdf.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -77,5 +79,5 @@
       ! JMM avoid negative salinities near river outlet ! Ugly fix
       ! JMM : restore negative salinities to small salinities:
-!!$   WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp
+!!$      WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp
 !!gm
 
@@ -95,6 +97,6 @@
       ENDIF
       !                                          ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
-         &                       tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
+         &                                  tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_zdf')
@@ -154,19 +156,11 @@
             IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
                IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)  
-                        END DO
-                     END DO
-                  END DO
+                  DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)  
+                  END_3D
                ELSE                          ! standard or triad iso-neutral operator
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
-                        END DO
-                     END DO
-                  END DO
+                  DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
+                  END_3D
                ENDIF
             ENDIF
@@ -174,26 +168,18 @@
             ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
             IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection
-               DO jk = 1, jpkm1
-                  DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
-                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm)
-                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
-                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   &
-                           &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
-                        zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
-                        zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp )
-                    END DO
-                  END DO
-               END DO
+               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+                  zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm)
+                  zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   &
+                     &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
+                  zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
+                  zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp )
+               END_3D
             ELSE
-               DO jk = 1, jpkm1
-                  DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt.
-                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm)
-                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
-                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
-                    END DO
-                  END DO
-               END DO
+               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+                  zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm)
+                  zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
+               END_3D
             ENDIF
             !
@@ -217,46 +203,28 @@
             !   used as a work space array: its value is modified.
             !
-            DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
-               DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction)
-                  zwt(ji,jj,1) = zwd(ji,jj,1)
-               END DO
-            END DO
-            DO jk = 2, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
-                  END DO
-               END DO
-            END DO
+            DO_2D( 0, 0, 0, 0 )
+               zwt(ji,jj,1) = zwd(ji,jj,1)
+            END_2D
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
+            END_3D
             !
          ENDIF 
          !         
-         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
-            DO ji = fs_2, fs_jpim1
-               pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
-            END DO
-         END DO
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side
-                  pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
-               END DO
-            END DO
-         END DO
+         DO_2D( 0, 0, 0, 0 )
+            pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
+         END_2D
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side
+            pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
+         END_3D
          !
-         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
-            DO ji = fs_2, fs_jpim1
-               pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
-            END DO
-         END DO
-         DO jk = jpk-2, 1, -1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   &
-                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_2D( 0, 0, 0, 0 )
+            pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
+         END_2D
+         DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
+            pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   &
+               &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
+         END_3D
          !                                            ! ================= !
       END DO                                          !  end tracer loop  !
