Changeset 13540 for NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA

Timestamp:

2020-09-29T12:41:06+02:00 (4 years ago)

Author:

andmirek

Message:

Ticket #2386: update to latest trunk

Location:

NEMO/branches/2020/r12377_ticket2386

Files:

• . (modified) (1 prop)
• src/OCE/TRA/eosbn2.F90 (modified) (17 diffs)
• src/OCE/TRA/traadv.F90 (modified) (2 diffs)
• src/OCE/TRA/traadv_cen.F90 (modified) (8 diffs)
• src/OCE/TRA/traadv_fct.F90 (modified) (28 diffs)
• src/OCE/TRA/traadv_mus.F90 (modified) (7 diffs)
• src/OCE/TRA/traadv_qck.F90 (modified) (15 diffs)
• src/OCE/TRA/traadv_ubs.F90 (modified) (15 diffs)
• src/OCE/TRA/traatf.F90 (modified) (6 diffs)
• src/OCE/TRA/traatf_qco.F90 (copied) (copied from NEMO/trunk/src/OCE/TRA/traatf_qco.F90)
• src/OCE/TRA/trabbc.F90 (modified) (2 diffs)
• src/OCE/TRA/trabbl.F90 (modified) (14 diffs)
• src/OCE/TRA/tradmp.F90 (modified) (4 diffs)
• src/OCE/TRA/traisf.F90 (modified) (4 diffs)
• src/OCE/TRA/traldf_iso.F90 (modified) (11 diffs)
• src/OCE/TRA/traldf_lap_blp.F90 (modified) (5 diffs)
• src/OCE/TRA/traldf_triad.F90 (modified) (13 diffs)
• src/OCE/TRA/tramle.F90 (modified) (12 diffs)
• src/OCE/TRA/tranpc.F90 (modified) (3 diffs)
• src/OCE/TRA/traqsr.F90 (modified) (8 diffs)
• src/OCE/TRA/trasbc.F90 (modified) (7 diffs)
• src/OCE/TRA/trazdf.F90 (modified) (7 diffs)
• src/OCE/TRA/zpshde.F90 (modified) (17 diffs)

NEMO/branches/2020/r12377_ticket2386

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev@HEAD      ext/AGRIF
+^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
 
 # SETTE
-^/utils/CI/sette@HEAD         sette
+^/utils/CI/sette@13507        sette

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev@HEAD      ext/AGRIF
+^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
 
 # SETTE
-^/utils/CI/sette@HEAD         sette
+^/utils/CI/sette@13507        sette

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/eosbn2.F90

@@ -180 +180 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -237 +238 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             !
             zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
@@ -273 +274 @@
       CASE( np_seos )                !==  simplified EOS  ==!
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
             zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
@@ -337 +338 @@
             END DO
             !
-            DO_3D_11_11( 1, jpkm1 )
+            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                !
                ! compute density (2*nn_sto_eos) times:
@@ -387 +388 @@
          ! Non-stochastic equation of state
          ELSE
-            DO_3D_11_11( 1, jpkm1 )
+            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
                !
                zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
@@ -425 +426 @@
       CASE( np_seos )                !==  simplified EOS  ==!
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
             zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
@@ -479 +480 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             !
             zh  = pdep(ji,jj) * r1_Z0                                  ! depth
@@ -514 +515 @@
       CASE( np_seos )                !==  simplified EOS  ==!
          !
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             !
             zt    = pts  (ji,jj,jp_tem)  - 10._wp
@@ -562 +563 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             !
             zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
@@ -615 +616 @@
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
             zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
@@ -669 +670 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             !
             zh  = pdep(ji,jj) * r1_Z0                                  ! depth
@@ -722 +723 @@
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             !
             zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
@@ -872 +873 @@
       IF( ln_timing )   CALL timing_start('bn2')
       !
-      DO_3D_11_11( 2, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90
          zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
             &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
@@ -920 +921 @@
       z1_T0   = 1._wp/40._wp
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          !
          zt  = ctmp   (ji,jj) * z1_T0
@@ -973 +974 @@
          !
          z1_S0 = 1._wp / 35.16504_wp
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity
             ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
@@ -1080 +1081 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             !
             zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
@@ -1139 +1140 @@
       CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==!
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
             zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv.F90

@@ -66 +66 @@
    INTEGER, PARAMETER ::   np_QCK     = 5   ! QUICK scheme
    
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -98 +99 @@
       IF( ln_wave .AND. ln_sdw )  THEN
          DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift
-            zuu(:,:,jk) = e2u  (:,:) * e3u(:,:,jk,Kmm) * ( uu(:,:,jk,Kmm) + usd(:,:,jk) )
-            zvv(:,:,jk) = e1v  (:,:) * e3v(:,:,jk,Kmm) * ( vv(:,:,jk,Kmm) + vsd(:,:,jk) )
-            zww(:,:,jk) = e1e2t(:,:)                 * ( ww(:,:,jk) + wsd(:,:,jk) )
+            zuu(:,:,jk) =   &
+               &  e2u  (:,:) * e3u(:,:,jk,Kmm) * ( uu(:,:,jk,Kmm) + usd(:,:,jk) )
+            zvv(:,:,jk) =   & 
+               &  e1v  (:,:) * e3v(:,:,jk,Kmm) * ( vv(:,:,jk,Kmm) + vsd(:,:,jk) )
+            zww(:,:,jk) =   & 
+               &  e1e2t(:,:)                 * ( ww(:,:,jk) + wsd(:,:,jk) )
          END DO
       ELSE

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_cen.F90

@@ -37 +37 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -103 +104 @@
          !
          CASE(  2  )                         !* 2nd order centered
-            DO_3D_10_10( 1, jpkm1 )
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
                zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
                zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
@@ -111 +112 @@
             ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero
             ztv(:,:,jpk) = 0._wp
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )          ! masked gradient
                ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
             END_3D
-            CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond.
+            CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond.
             !
-            DO_3D_00_10( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )           ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! C2 interpolation of T at u- & v-points (x2)
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
@@ -127 +128 @@
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v
             END_3D
+            CALL lbc_lnk_multi( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. )
             !
          CASE DEFAULT
-            CALL ctl_stop( 'traadv_fct: wrong value for nn_fct' )
+            CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
          END SELECT
          !
@@ -135 +137 @@
          !
          CASE(  2  )                         !* 2nd order centered
-            DO_3D_00_00( 2, jpk )
+            DO_3D( 0, 0, 0, 0, 2, jpk )
                zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) * wmask(ji,jj,jk)
             END_3D
@@ -141 +143 @@
          CASE(  4  )                         !* 4th order compact
             CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )      ! ztw = interpolated value of T at w-point
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zwz(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
             END_3D
@@ -149 +151 @@
          IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask)
             IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
-               DO_2D_11_11
+               DO_2D( 1, 1, 1, 1 )
                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
                END_2D
@@ -157 +159 @@
          ENDIF
          !               
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !--  Divergence of advective fluxes  --!
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    &
                &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
                &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    &
-               &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+               &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) &
+               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
-         !                             ! trend diagnostics
+         !                               ! trend diagnostics
          IF( l_trd ) THEN
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_fct.F90

@@ -46 +46 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -96 +97 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
+      !! -- init to 0
+      zwi(:,:,:) = 0._wp
+      zwx(:,:,:) = 0._wp
+      zwy(:,:,:) = 0._wp
+      zwz(:,:,:) = 0._wp
+      ztu(:,:,:) = 0._wp
+      ztv(:,:,:) = 0._wp
+      zltu(:,:,:) = 0._wp
+      zltv(:,:,:) = 0._wp
+      ztw(:,:,:) = 0._wp
       !
       l_trd = .FALSE.            ! set local switches
@@ -128 +139 @@
       IF( ll_zAimp ) THEN
          ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
-         DO_3D_00_00( 1, jpkm1 )
-            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   &
+            &                               / e3t(ji,jj,jk,Krhs)
             zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
             zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
@@ -139 +151 @@
          !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
          !                    !* upstream tracer flux in the i and j direction 
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
             ! upstream scheme
             zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
@@ -148 +160 @@
             zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
          END_3D
-         !                    !* upstream tracer flux in the k direction *!
-         DO_3D_11_11( 2, jpkm1 )
+         !                               !* upstream tracer flux in the k direction *!
+         DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
             zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
             zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
             zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
          END_3D
-         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked)
-            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
-               DO_2D_11_11
+         IF( ln_linssh ) THEN               ! top ocean value (only in linear free surface as zwz has been w-masked)
+            IF( ln_isfcav ) THEN                        ! top of the ice-shelf cavities and at the ocean surface
+               DO_2D( 1, 1, 1, 1 )
                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                END_2D
-            ELSE                             ! no cavities: only at the ocean surface
-               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
+            ELSE                                        ! no cavities: only at the ocean surface
+               DO_2D( 1, 1, 1, 1 )
+                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
+               END_2D
             ENDIF
          ENDIF
          !               
-         DO_3D_00_00( 1, jpkm1 )
-            !                             ! total intermediate advective trends
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* trend and after field with monotonic scheme
+            !                               ! total intermediate advective trends
             ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
                &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
-            !                             ! update and guess with monotonic sheme
-            pt(ji,jj,jk,jn,Krhs) =                     pt(ji,jj,jk,jn,Krhs) +        ztra   / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
-            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
+            !                               ! update and guess with monotonic sheme
+            pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   &
+               &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
+            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) &
+               &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
          END_3D
          
@@ -178 +194 @@
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -184 +200 @@
                zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
             END_3D
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
                   &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
@@ -202 +218 @@
          !
          CASE(  2  )                   !- 2nd order centered
-            DO_3D_10_10( 1, jpkm1 )
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
                zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
                zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
@@ -211 +227 @@
             zltv(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1                 ! Laplacian
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )                 ! 1st derivative (gradient)
                   ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                   ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
                END_2D
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )                 ! 2nd derivative * 1/ 6
                   zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
                   zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
                END_2D
             END DO
-            CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
-            !
-            DO_3D_10_10( 1, jpkm1 )
+            CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )    ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
-               !                                                  ! C4 minus upstream advective fluxes 
+               !                                                        ! C4 minus upstream advective fluxes 
                zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
@@ -233 +249 @@
             ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
             ztv(:,:,jpk) = 0._wp
-            DO_3D_10_10( 1, jpkm1 )
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )    ! 1st derivative (gradient)
                ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
             END_3D
-            CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
-            !
-            DO_3D_00_00( 1, jpkm1 )
+            CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2)
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
@@ -255 +271 @@
          !
          CASE(  2  )                   !- 2nd order centered
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                   &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
@@ -262 +278 @@
          CASE(  4  )                   !- 4th order COMPACT
             CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
             END_3D
@@ -272 +288 @@
          !         
          IF ( ll_zAimp ) THEN
-            DO_3D_00_00( 1, jpkm1 )
-               !                             ! total intermediate advective trends
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !* trend and after field with monotonic scheme
+               !                                                ! total intermediate advective trends
                ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                   &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
@@ -282 +298 @@
             CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
             !
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -289 +305 @@
          END IF
          !
-         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. )
+         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp,  zwz, 'W',  1.0_wp )
          !
          !        !==  monotonicity algorithm  ==!
@@ -297 +313 @@
          !        !==  final trend with corrected fluxes  ==!
          !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
@@ -308 +324 @@
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -314 +330 @@
                zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
             END_3D
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
                   &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
@@ -374 +390 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikm1         ! local integer
-      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
-      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
-      !!----------------------------------------------------------------------
-      !
-      zbig  = 1.e+40_wp
-      zrtrn = 1.e-15_wp
-      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
+      REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
+      REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
+      REAL(dp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
+      !!----------------------------------------------------------------------
+      !
+      zbig  = 1.e+40_dp
+      zrtrn = 1.e-15_dp
+      zbetup(:,:,:) = 0._dp   ;   zbetdo(:,:,:) = 0._dp
 
       ! Search local extrema
@@ -393 +409 @@
       DO jk = 1, jpkm1
          ikm1 = MAX(jk-1,1)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
 
             ! search maximum in neighbourhood
@@ -423 +439 @@
          END_2D
       END DO
-      CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
+      CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp )   ! lateral boundary cond. (unchanged sign)
 
       ! 3. monotonic flux in the i & j direction (paa & pbb)
       ! ----------------------------------------
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
          zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
-         zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
+         zcu =       ( 0.5  + SIGN( 0.5_wp , paa(ji,jj,jk) ) )
          paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
 
          zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
          zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
-         zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
+         zcv =       ( 0.5  + SIGN( 0.5_wp , pbb(ji,jj,jk) ) )
          pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
 
-! monotonic flux in the k direction, i.e. pcc
-! -------------------------------------------
+      ! monotonic flux in the k direction, i.e. pcc
+      ! -------------------------------------------
          za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
          zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
-         zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
+         zc =       ( 0.5  + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) )
          pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
       END_3D
-      CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
+      CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1.0_wp , pbb, 'V', -1.0_wp )   ! lateral boundary condition (changed sign)
       !
    END SUBROUTINE nonosc
@@ -465 +481 @@
       !!----------------------------------------------------------------------
       
-      DO_3D_11_11( 3, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 3, jpkm1 )       !==  build the three diagonal matrix  ==!
          zwd (ji,jj,jk) = 4._wp
          zwi (ji,jj,jk) = 1._wp
@@ -479 +495 @@
       END_3D
       !
-      jk = 2                                          ! Switch to second order centered at top
-      DO_2D_11_11
+      jk = 2                                    ! Switch to second order centered at top
+      DO_2D( 1, 1, 1, 1 )
          zwd (ji,jj,jk) = 1._wp
          zwi (ji,jj,jk) = 0._wp
@@ -488 +504 @@
       !
       !                       !==  tridiagonal solve  ==!
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )           ! first recurrence
          zwt(ji,jj,2) = zwd(ji,jj,2)
       END_2D
-      DO_3D_11_11( 3, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 3, 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
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )           ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,2) = zwrm(ji,jj,2)
       END_2D
-      DO_3D_11_11( 3, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 3, jpkm1 )
          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
       END_3D
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )           ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
-      DO_3DS_11_11( jpk-2, 2, -1 )
+      DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 )
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D
@@ -530 +546 @@
       !                      !==  build the three diagonal matrix & the RHS  ==!
       !
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )    ! interior (from jk=3 to jpk-1)
          zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal
          zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal
@@ -549 +565 @@
       END IF
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )              ! 2nd order centered at top & bottom
          ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point
          ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point
@@ -566 +582 @@
       !                       !==  tridiagonal solver  ==!
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )           !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
          zwt(ji,jj,2) = zwd(ji,jj,2)
       END_2D
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, 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
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,2) = zwrm(ji,jj,2)
       END_2D
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
       END_3D
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
-      DO_3DS_00_00( jpk-2, 2, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D
@@ -622 +638 @@
       kstart =  1  + klev
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                         !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
          zwt(ji,jj,kstart) = pD(ji,jj,kstart)
       END_2D
-      DO_3D_00_00( kstart+1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )
          zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
       END_3D
       !
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                        !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
       END_2D
-      DO_3D_00_00( kstart+1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )
          pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
       END_3D
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                       !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
-      DO_3DS_00_00( jpk-2, kstart, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 )
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_mus.F90

@@ -47 +47 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -131 +132 @@
          zwx(:,:,jpk) = 0._wp                   ! bottom values
          zwy(:,:,jpk) = 0._wp  
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
             zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
             zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
          END_3D
          ! lateral boundary conditions   (changed sign)
-         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )
+         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )
          !                                !-- Slopes of tracer
          zslpx(:,:,jpk) = 0._wp                 ! bottom values
          zslpy(:,:,jpk) = 0._wp
-         DO_3D_01_01( 1, jpkm1 )
-            zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
-               &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
-            zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
-               &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
-         END_3D
-         !
-         DO_3D_01_01( 1, jpkm1 )
-            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
-               &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
-               &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
-            zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
-               &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
-               &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
-         END_3D
-         !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 1, 0, 1, 1, jpkm1 )
+            zslpx(ji,jj,jk) =                       ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
+               &            * ( 0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
+            zslpy(ji,jj,jk) =                       ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
+               &            * ( 0.25 + SIGN( 0.25_wp, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
+         END_3D
+         !
+         DO_3D( 0, 1, 0, 1, 1, jpkm1 )    !-- Slopes limitation
+            zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
+               &                                                     2.*ABS( zwx  (ji-1,jj,jk) ),   &
+               &                                                     2.*ABS( zwx  (ji  ,jj,jk) ) )
+            zslpy(ji,jj,jk) = SIGN( 1.0_wp, zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
+               &                                                     2.*ABS( zwy  (ji,jj-1,jk) ),   &
+               &                                                     2.*ABS( zwy  (ji,jj  ,jk) ) )
+         END_3D
+         !
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- MUSCL horizontal advective fluxes
             ! MUSCL fluxes
-            z0u = SIGN( 0.5, pU(ji,jj,jk) )
+            z0u = SIGN( 0.5_wp, pU(ji,jj,jk) )
             zalpha = 0.5 - z0u
             zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
@@ -165 +166 @@
             zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
             !
-            z0v = SIGN( 0.5, pV(ji,jj,jk) )
+            z0v = SIGN( 0.5_wp, pV(ji,jj,jk) )
             zalpha = 0.5 - z0v
             zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
@@ -172 +173 @@
             zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
          END_3D
-         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
-         !
-         DO_3D_00_00( 1, jpkm1 )
+         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! lateral boundary conditions   (changed sign)
+         !
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- Tracer advective trend
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
             &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
@@ -199 +200 @@
          !                                !-- Slopes of tracer
          zslpx(:,:,1) = 0._wp                   ! surface values
-         DO_3D_11_11( 2, jpkm1 )
-            zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
-               &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
-         END_3D
-         DO_3D_11_11( 2, jpkm1 )
-            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
-               &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
-               &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
-         END_3D
-         DO_3D_00_00( 1, jpk-2 )
-            z0w = SIGN( 0.5, pW(ji,jj,jk+1) )
+         DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+            zslpx(ji,jj,jk) =                        ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
+               &            * (  0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
+         END_3D
+         DO_3D( 1, 1, 1, 1, 2, jpkm1 )    !-- Slopes limitation
+            zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
+               &                                                     2.*ABS( zwx  (ji,jj,jk+1) ),   &
+               &                                                     2.*ABS( zwx  (ji,jj,jk  ) )  )
+         END_3D
+         DO_3D( 0, 0, 0, 0, 1, jpk-2 )    !-- vertical advective flux
+            z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) )
             zalpha = 0.5 + z0w
             zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
@@ -218 +219 @@
          IF( ln_linssh ) THEN                   ! top values, linear free surface only
             IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
-               DO_2D_11_11
+               DO_2D( 1, 1, 1, 1 )
                   zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
                END_2D
@@ -226 +227 @@
          ENDIF
          !
-         DO_3D_00_00( 1, jpkm1 )
-            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )     !-- vertical advective trend
+            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )   &
+               &                                      * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !                                ! send trends for diagnostic

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_qck.F90

@@ -41 +41 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -141 +142 @@
          !
 !!gm why not using a SHIFT instruction...
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
             zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
             zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
          END_3D
-         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions 
+         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions 
          
          !
          ! Horizontal advective fluxes
          ! ---------------------------
-         DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
             zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T 
          END_3D
          !
-         DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
             zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
             zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
@@ -163 +164 @@
          END_3D
          !--- Lateral boundary conditions 
-         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. )
+         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
 
          !--- QUICKEST scheme
@@ -169 +170 @@
          !
          ! Mask at the T-points in the x-direction (mask=0 or mask=1)
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
          END_3D
-         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions 
+         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions 
 
          !
@@ -178 +179 @@
          DO jk = 1, jpkm1  
             !
-            DO_2D_00_00
-               zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+            DO_2D( 0, 0, 0, 0 )
+               zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                !--- If the second ustream point is a land point
                !--- the flux is computed by the 1st order UPWIND scheme
@@ -188 +189 @@
          END DO
          !
-         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
+         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions
          !
          ! Computation of the trend
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             ! horizontal advective trends
@@ -232 +233 @@
             !                                             
             !--- Computation of the ustream and downstream value of the tracer and the mask
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                ! Upstream in the x-direction for the tracer
                zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
@@ -239 +240 @@
             END_2D
          END DO
-         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions 
+         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions 
 
          
@@ -246 +247 @@
          ! ---------------------------
          !
-         DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
             zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T 
          END_3D
          !
-         DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
             zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
             zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
@@ -260 +261 @@
 
          !--- Lateral boundary conditions 
-         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
+         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )
 
          !--- QUICKEST scheme
@@ -266 +267 @@
          !
          ! Mask at the T-points in the x-direction (mask=0 or mask=1)
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
          END_3D
-         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions 
+         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions 
          !
          ! Tracer flux on the x-direction
          DO jk = 1, jpkm1  
             !
-            DO_2D_00_00
-               zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+            DO_2D( 0, 0, 0, 0 )
+               zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                !--- If the second ustream point is a land point
                !--- the flux is computed by the 1st order UPWIND scheme
@@ -284 +285 @@
          END DO
          !
-         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
+         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions
          !
          ! Computation of the trend
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             ! horizontal advective trends
@@ -326 +327 @@
          !                                                       ! ===========
          !
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       !* Interior point   (w-masked 2nd order centered flux)
             zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk)
          END_3D
          IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
             IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
-               DO_2D_11_11
+               DO_2D( 1, 1, 1, 1 )
                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface 
                END_2D
@@ -339 +340 @@
          ENDIF
          !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !==  Tracer flux divergence added to the general trend  ==!
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
                &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
@@ -368 +369 @@
       !----------------------------------------------------------------------
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zc     = puc(ji,jj,jk)                         ! Courant number
          zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traadv_ubs.F90

@@ -39 +39 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -82 +83 @@
       !!
       !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 
-      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741. 
+      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731�1741.
       !!----------------------------------------------------------------------
       INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
@@ -123 +124 @@
          !                                                       ! ===========
          !                                              
-         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
-            DO_2D_10_10
+         DO jk = 1, jpkm1                !==  horizontal laplacian of before tracer ==!
+            DO_2D( 1, 0, 1, 0 )                   ! First derivative (masked gradient)
                zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
                zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
@@ -130 +131 @@
                ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
             END_2D
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )                   ! Second derivative (divergence)
                zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
                zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
@@ -137 +138 @@
             !                                    
          END DO         
-         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
+         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
          !    
-         DO_3D_10_10( 1, jpkm1 )
-            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2)
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )   !==  Horizontal advective fluxes  ==!     (UBS)
+            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )        ! upstream transport (x2)
             zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
             zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
@@ -155 +156 @@
          !
          DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
                   &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
-                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) &
+                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_2D
             !                                             
@@ -164 +166 @@
          !
          zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
-         !                                            ! and/or in trend diagnostic (l_trd=T) 
+         !                                                ! and/or in trend diagnostic (l_trd=T) 
          !                
          IF( l_trd ) THEN                  ! trend diagnostics
@@ -185 +187 @@
             IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
             !
-            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
-            DO_3D_11_11( 2, jpkm1 )
+            !                               !*  upstream advection with initial mass fluxes & intermediate update  ==!
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
                zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
                zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
                ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb)  ) * wmask(ji,jj,jk)
             END_3D
-            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
-               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
-                  DO_2D_11_11
+            IF( ln_linssh ) THEN                ! top ocean value (only in linear free surface as ztw has been w-masked)
+               IF( ln_isfcav ) THEN                   ! top of the ice-shelf cavities and at the ocean surface
+                  DO_2D( 1, 1, 1, 1 )
                      ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                   END_2D
-               ELSE                                ! no cavities: only at the ocean surface
+               ELSE                                   ! no cavities: only at the ocean surface
                   ztw(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
                ENDIF
             ENDIF
             !
-            DO_3D_00_00( 1, jpkm1 )
-               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* trend and after field with monotonic scheme
+               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
+                  &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak 
                zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
             END_3D
-            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
+            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1.0_wp )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
             !
             !                          !*  anti-diffusive flux : high order minus low order
-            DO_3D_11_11( 2, jpkm1 )
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
                ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                   &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
@@ -220 +223 @@
          CASE(  4  )                               ! 4th order COMPACT
             CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
             END_3D
@@ -227 +230 @@
          END SELECT
          !
-         DO_3D_00_00( 1, jpkm1 )
-            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !  final trend with corrected fluxes
+            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
+               &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !
-         IF( l_trd )  THEN       ! vertical advective trend diagnostics
-            DO_3D_00_00( 1, jpkm1 )
+         IF( l_trd )  THEN               ! vertical advective trend diagnostics
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )                 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
                zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
                   &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
@@ -270 +274 @@
       !!----------------------------------------------------------------------
       !
-      zbig  = 1.e+40_wp
+      zbig  = 1.e+38_wp
       zrtrn = 1.e-15_wp
       zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
@@ -282 +286 @@
       DO jk = 1, jpkm1     ! search maximum in neighbourhood
          ikm1 = MAX(jk-1,1)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
                &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
@@ -294 +298 @@
       DO jk = 1, jpkm1     ! search minimum in neighbourhood
          ikm1 = MAX(jk-1,1)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
                &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
@@ -306 +310 @@
       ! Positive and negative part of fluxes and beta terms
       ! ---------------------------------------------------
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          ! positive & negative part of the flux
          zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
@@ -318 +322 @@
       ! monotonic flux in the k direction, i.e. pcc
       ! -------------------------------------------
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
          zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
-         zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
+         zc = 0.5 * ( 1.e0 + SIGN( 1.0_wp, pcc(ji,jj,jk) ) )
          pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
       END_3D

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traatf.F90

@@ -58 +58 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -109 +110 @@
 #endif
       !                                              ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
+      CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp )
       !
       IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
@@ -155 +156 @@
          ENDIF
          !
-         CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., &
-                  &                    pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., &
-                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  )
+         CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kbb) , 'T', 1.0_wp, &
+                  &                    pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp, &
+                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp  )
          !
       ENDIF     
       !
       IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
-         zfact = 1._wp / rDt             
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact
-            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * zfact
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * r1_Dt
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * r1_Dt
          END DO
          CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
@@ -210 +210 @@
       DO jn = 1, kjpt
          !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             ztn = pt(ji,jj,jk,jn,Kmm)                                    
             ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers
@@ -229 +229 @@
       !! 
       !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
-      !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + rn_atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] )
-      !!                       /( e3t(Kmm)         + rn_atfp*[ e3t(Kbb)         - 2 e3t(Kmm)         + e3t(Kaa)    ] )
+      !!             pt(Kmm)  = ( e3t_Kmm*pt(Kmm) + rn_atfp*[ e3t_Kbb*pt(Kbb) - 2 e3t_Kmm*pt(Kmm) + e3t_Kaa*pt(Kaa) ] )
+      !!                       /( e3t_Kmm         + rn_atfp*[ e3t_Kbb         - 2 e3t_Kmm         + e3t_Kaa    ] )
       !!
       !! ** Action  : - pt(Kmm) ready for the next time step
@@ -275 +275 @@
       zfact2 = zfact1 * r1_rho0
       DO jn = 1, kjpt      
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             ze3t_b = e3t(ji,jj,jk,Kbb)
             ze3t_n = e3t(ji,jj,jk,Kmm)

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/trabbc.F90

@@ -46 +46 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -90 +91 @@
       ENDIF
       !                             !  Add the geothermal trend on temperature
-      DO_2D_00_00
-         pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm)
+      DO_2D( 0, 0, 0, 0 )
+         pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs)   &
+            &             + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm)
       END_2D
       !
-      CALL lbc_lnk( 'trabbc', pts(:,:,:,jp_tem,Krhs) , 'T', 1. )
+      CALL lbc_lnk( 'trabbc', pts(:,:,:,jp_tem,Krhs) , 'T', 1.0_wp )
       !
       IF( l_trdtra ) THEN        ! Send the trend for diagnostics

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/trabbl.F90

@@ -68 +68 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -125 +126 @@
             &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
-         CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1. , ahv_bbl, 'V', 1. )
+         CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1.0_wp , ahv_bbl, 'V', 1.0_wp )
          CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
          CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
@@ -138 +139 @@
             &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
-         CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1. , vtr_bbl, 'V', 1. )
+         CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp )
          CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
          CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
@@ -191 +192 @@
       DO jn = 1, kjpt                                     ! tracer loop
          !                                                ! ===========
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ik = mbkt(ji,jj)                             ! bottom T-level index
             zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
          END_2D
          !               
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )                               ! Compute the trend
             ik = mbkt(ji,jj)                            ! bottom T-level index
             pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
@@ -342 +343 @@
       ENDIF
       !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ik = mbkt(ji,jj)                             ! bottom T-level index
          zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
@@ -357 +358 @@
       IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
          !                                !-------------------!
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
             !                                                   ! i-direction
             za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
@@ -365 +366 @@
                &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
             !
-            zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
+            zsign  = SIGN(  0.5_wp, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
             ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
             !
@@ -375 +376 @@
                &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
             !
-            zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
+            zsign = SIGN(  0.5_wp, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
             ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
          END_2D
@@ -387 +388 @@
          !
          CASE( 1 )                                   != use of upper velocity
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
                !                                                  ! i-direction
                za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
@@ -395 +396 @@
                          - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
                !
-               zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
-               zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
+               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
+               zsigna= SIGN(  0.5_wp, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
                !
                !                                                          ! bbl velocity
@@ -407 +408 @@
                zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
                   &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
-               zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
-               zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
+               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
+               zsigna= SIGN(  0.5_wp, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
                !
                !                                                          ! bbl transport
@@ -416 +417 @@
          CASE( 2 )                                 != bbl velocity = F( delta rho )
             zgbbl = grav * rn_gambbl
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down
                !                                                  ! i-direction
                ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
@@ -504 +505 @@
       IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
       !
-      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
-      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
+      IF(lwp) THEN
+         IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
+         IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
+      ENDIF
       !
       !                             !* vertical index of  "deep" bottom u- and v-points
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )                 ! (the "shelf" bottom k-indices are mbku and mbkv)
          mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
          mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
@@ -514 +517 @@
       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
       zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp )  
-      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.) 
+      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp) 
       mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
       !
       !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
       mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
-            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            mgrhu(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
          ENDIF
          !
          IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
-            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            mgrhv(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
          ENDIF
       END_2D
       !
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )           !* bbl thickness at u- (v-) point; minimum of top & bottom e3u_0 (e3v_0)
          e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
          e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
       END_2D
-      CALL lbc_lnk_multi( 'trabbl', e3u_bbl_0, 'U', 1. , e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
+      CALL lbc_lnk_multi( 'trabbl', e3u_bbl_0, 'U', 1.0_wp , e3v_bbl_0, 'V', 1.0_wp )      ! lateral boundary conditions
       !
       !                             !* masked diffusive flux coefficients

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/tradmp.F90

@@ -112 +112 @@
       CASE( 0 )                        !*  newtonian damping throughout the water column  *!
          DO jn = 1, jpts
-            DO_3D_00_00( 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           &
                   &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
@@ -119 +119 @@
          !
       CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *!
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             IF( avt(ji,jj,jk) <= avt_c ) THEN
                pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
@@ -129 +129 @@
          !
       CASE ( 2 )                       !*  no damping in the mixed layer   *!
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
                pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
@@ -208 +208 @@
          !                          ! Read in mask from file
          CALL iom_open ( cn_resto, imask)
-         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto )
+         CALL iom_get  ( imask, jpdom_auto, 'resto', resto )
          CALL iom_close( imask )
       ENDIF

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traisf.F90

@@ -11 +11 @@
    !!----------------------------------------------------------------------
    USE isf_oce                                     ! Ice shelf variables
-   USE dom_oce , ONLY : e3t, r1_e1e2t            ! ocean space domain variables
+   USE dom_oce                                     ! ocean space domain variables
    USE isfutils, ONLY : debug                      ! debug option
    USE timing  , ONLY : timing_start, timing_stop  ! Timing
@@ -23 +23 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -107 +108 @@
       !
       ! update pts(:,:,:,:,Krhs)
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          !
          ikt = ktop(ji,jj)
@@ -140 +141 @@
       !
       DO jk = 1,jpk
-         ptsa(:,:,jk,jp_tem) = ptsa(:,:,jk,jp_tem) + ptsc(:,:,jk,jp_tem) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
-         ptsa(:,:,jk,jp_sal) = ptsa(:,:,jk,jp_sal) + ptsc(:,:,jk,jp_sal) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
+         ptsa(:,:,jk,jp_tem) =   &
+            &  ptsa(:,:,jk,jp_tem) + ptsc(:,:,jk,jp_tem) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
+         ptsa(:,:,jk,jp_sal) =   &
+            &  ptsa(:,:,jk,jp_sal) + ptsc(:,:,jk,jp_sal) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
       END DO
       !

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traldf_iso.F90

@@ -41 +41 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -140 +141 @@
       IF( kpass == 1 ) THEN                  !==  first pass only  ==!
          !
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             !
             zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
@@ -157 +158 @@
          !
          IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                akz(ji,jj,jk) = 0.25_wp * (                                                                     &
                   &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
@@ -166 +167 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               DO_3D_10_10( 2, jpkm1 )
-                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
+               DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+                  akz(ji,jj,jk) = 16._wp   &
+                     &   * ah_wslp2   (ji,jj,jk)   &
+                     &   * (  akz     (ji,jj,jk)   &
+                     &      + ah_wslp2(ji,jj,jk)   &
+                     &        / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               DO_3D_10_10( 2, jpkm1 )
+               DO_3D( 1, 0, 1, 0, 2, jpkm1 )
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                   zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
@@ -196 +200 @@
 
          ! Horizontal tracer gradient 
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )           ! bottom correction (partial bottom cell)
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)          
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)          
                   IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
@@ -225 +229 @@
             ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk)
             ENDIF
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )           !==  Horizontal fluxes
                zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
                zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
@@ -246 +250 @@
             END_2D
             !
-            DO_2D_00_00
-               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
-                  &                                                 + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
+            DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
+               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
+                  &       + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
                   &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_2D
@@ -262 +266 @@
          ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
          
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
             !
             zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
@@ -284 +288 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D_00_00( 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               &
@@ -293 +297 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D_00_00( 2, jpkm1 )
-                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk)                       &
-                     &           + ah_wslp2(ji,jj,jk)  * e1e2t(ji,jj)   &
+               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+                  ztfw(ji,jj,jk) =   &
+                     &  ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
                      &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
                END_3D
             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
-               DO_3D_00_00( 2, jpkm1 )
+               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                  &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
@@ -307 +311 @@
          ENDIF
          !         
-         DO_3D_00_00( 1, jpkm1 )
-            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
-               &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
+            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * r1_e1e2t(ji,jj)   &
+               &                                             / e3t(ji,jj,jk,Kmm)
          END_3D
          !

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traldf_lap_blp.F90

@@ -38 +38 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -98 +99 @@
       ELSE                    ;   zsign = -1._wp
       ENDIF
-      DO_3D_10_10( 1, jpkm1 )
+      DO_3D( 1, 0, 1, 0, 1, jpkm1 )
          zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)   !!gm   * umask(ji,jj,jk) pah masked!
          zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   !!gm   * vmask(ji,jj,jk)
@@ -107 +108 @@
          !                          ! =========== !    
          !                               
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )            !== First derivative (gradient)  ==!
             ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
             ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
          END_3D
-         IF( ln_zps ) THEN                ! set gradient at bottom/top ocean level
-            DO_2D_10_10
+         IF( ln_zps ) THEN                             ! set gradient at bottom/top ocean level
+            DO_2D( 1, 0, 1, 0 )                              ! bottom
                ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
                ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
             END_2D
-            IF( ln_isfcav ) THEN                ! top in ocean cavities only
-               DO_2D_10_10
+            IF( ln_isfcav ) THEN                             ! top in ocean cavities only
+               DO_2D( 1, 0, 1, 0 )
                   IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
                   IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
@@ -124 +125 @@
          ENDIF
          !
-         DO_3D_00_00( 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
                &                                      +    ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
@@ -199 +200 @@
       END SELECT
       !
-      CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
+      CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign)
       !                                               ! Partial top/bottom cell: GRADh( zlap )  
       IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traldf_triad.F90

@@ -41 +41 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -136 +137 @@
          DO ip = 0, 1                            ! i-k triads
             DO kp = 0, 1
-               DO_3D_10_10( 1, jpkm1 )
+               DO_3D( 1, 0, 1, 0, 1, jpkm1 )
                   ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
                   zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
@@ -156 +157 @@
          DO jp = 0, 1                            ! j-k triads 
             DO kp = 0, 1
-               DO_3D_10_10( 1, jpkm1 )
+               DO_3D( 1, 0, 1, 0, 1, jpkm1 )
                   ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,Kmm)
                   zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
@@ -178 +179 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
-               DO_3D_10_10( 2, jpkm1 )
-                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
+               DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+                  akz(ji,jj,jk) = 16._wp           &
+                     &   * ah_wslp2   (ji,jj,jk)   &
+                     &   * (  akz     (ji,jj,jk)   &
+                     &      + ah_wslp2(ji,jj,jk)   &
+                     &        / ( e3w (ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
-               DO_3D_10_10( 2, jpkm1 )
+               DO_3D( 1, 0, 1, 0, 2, jpkm1 )
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                   zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
@@ -207 +211 @@
          zftv(:,:,:) = 0._wp
          !
-         DO_3D_10_10( 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
-            DO_2D_10_10
+            DO_2D( 1, 0, 1, 0 )                    ! bottom level
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
-               DO_2D_10_10
+               DO_2D( 1, 0, 1, 0 )
                   IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
                   IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
@@ -242 +246 @@
                DO ip = 0, 1              !==  Horizontal & vertical fluxes
                   DO kp = 0, 1
-                     DO_2D_10_10
+                     DO_2D( 1, 0, 1, 0 )
                         ze1ur = r1_e1u(ji,jj)
                         zdxt  = zdit(ji,jj,jk) * ze1ur
@@ -263 +267 @@
                DO jp = 0, 1
                   DO kp = 0, 1
-                     DO_2D_10_10
+                     DO_2D( 1, 0, 1, 0 )
                         ze2vr = r1_e2v(ji,jj)
                         zdyt  = zdjt(ji,jj,jk) * ze2vr
@@ -285 +289 @@
                DO ip = 0, 1               !==  Horizontal & vertical fluxes
                   DO kp = 0, 1
-                     DO_2D_10_10
+                     DO_2D( 1, 0, 1, 0 )
                         ze1ur = r1_e1u(ji,jj)
                         zdxt  = zdit(ji,jj,jk) * ze1ur
@@ -306 +310 @@
                DO jp = 0, 1
                   DO kp = 0, 1
-                     DO_2D_10_10
+                     DO_2D( 1, 0, 1, 0 )
                         ze2vr = r1_e2v(ji,jj)
                         zdyt  = zdjt(ji,jj,jk) * ze2vr
@@ -325 +329 @@
             ENDIF
             !                             !==  horizontal divergence and add to the general trend  ==!
-            DO_2D_00_00
-               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
+            DO_2D( 0, 0, 0, 0 )
+               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
+                  &                       + zsign * (  zftu(ji-1,jj  ,jk) - zftu(ji,jj,jk)       &
                   &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
                   &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
@@ -335 +340 @@
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
-            DO_3D_10_00( 2, jpkm1 )
+            DO_3D( 1, 0, 0, 0, 2, jpkm1 )
                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
@@ -343 +348 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
-               DO_3D_10_00( 2, jpkm1 )
+               DO_3D( 1, 0, 0, 0, 2, jpkm1 )
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
                      &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                END_3D
             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
-               DO_3D_10_00( 2, jpkm1 )
+               DO_3D( 1, 0, 0, 0, 2, jpkm1 )
                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
@@ -356 +361 @@
          ENDIF
          !
-         DO_3D_00_00( 1, jpkm1 )
-            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
+            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
+            &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
                &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
          END_3D

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/tramle.F90

@@ -49 +49 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -99 +100 @@
       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
-         DO_3DS_11_11( jpkm1, nlb10, -1 )
+         DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
             IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
          END_3D
@@ -109 +110 @@
       zbm (:,:) = 0._wp
       zn2 (:,:) = 0._wp
-      DO_3D_11_11( 1, ikmax )
+      DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
          zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
          zmld(ji,jj) = zmld(ji,jj) + zc
@@ -118 +119 @@
       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
       CASE ( 0 )                                               != min of the 2 neighbour MLDs
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
             zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
          END_2D
       CASE ( 1 )                                               != average of the 2 neighbour MLDs
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
             zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
          END_2D
       CASE ( 2 )                                               != max of the 2 neighbour MLDs
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
             zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
@@ -145 +146 @@
       !
       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
                &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
@@ -156 +157 @@
          !
       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
                &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
@@ -166 +167 @@
       !
       IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
             IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
@@ -173 +174 @@
       !
       !                                      !==  structure function value at uw- and vw-points  ==!
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
          zhv(ji,jj) = 1._wp / zhv(ji,jj)
@@ -181 +182 @@
       zpsi_vw(:,:,:) = 0._wp
       !
-      DO_3D_10_10( 2, ikmax )
+      DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0
          zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
          zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj)
@@ -195 +196 @@
       !                                      !==  transport increased by the MLE induced transport ==!
       DO jk = 1, ikmax
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )                      ! CAUTION pu,pv must be defined at row/column i=1 / j=1
             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
          END_2D
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
                &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
@@ -282 +283 @@
             IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
             z1_t2 = 1._wp / ( rn_time * rn_time )
-            DO_2D_01_01
+            DO_2D( 0, 1, 0, 1 )                      ! "coriolis+ time^-1" at u- & v-points
                zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp
                zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp
@@ -288 +289 @@
                rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 )
             END_2D
-            CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1. , rfv, 'V', 1. )
+            CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1.0_wp , rfv, 'V', 1.0_wp )
             !
          ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/tranpc.F90

@@ -35 +35 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -102 +103 @@
          inpcc = 0
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )                                ! interior column only
             !
             IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
@@ -309 +310 @@
          ENDIF
          !
-         CALL lbc_lnk_multi( 'tranpc', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
+         CALL lbc_lnk_multi( 'tranpc', pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp )
          !
          IF( lwp .AND. l_LB_debug ) THEN

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traqsr.F90

@@ -63 +63 @@
    REAL(wp) ::   xsi1r   ! inverse of rn_si1
    !
-   REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
+   REAL(wp) , PUBLIC, DIMENSION(3,61)   ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
 
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -110 +111 @@
       REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
       REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
-      REAL(wp) ::   zz0 , zz1                !    -         -
-      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
-      REAL(wp) ::   zlogc, zlogc2, zlogc3 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d
+      REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !    -         -
+      REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
+      REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: ze0, ze1, ze2, ze3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d
       !!----------------------------------------------------------------------
       !
@@ -138 +138 @@
             IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
             z1_2 = 0.5_wp
-            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux
+            CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux
          ELSE                                           ! No restart or restart not found: Euler forward time stepping
             z1_2 = 1._wp
@@ -160 +160 @@
       CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
          !
-         ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , &
-            &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , &
-            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
+         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,   &
+            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,   &
+            &      ztmp3d(jpi,jpj,nksr + 1)                     )
          !
          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
             CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
+            !
+            ! Separation in R-G-B depending on the surface Chl
+            ! perform and store as many of the 2D calculations as possible
+            ! before the 3D loop (use the temporary 2D arrays to replace the
+            ! most expensive calculations)
+            !
+            DO_2D( 0, 0, 0, 0 )
+                       ! zlogc = log(zchl)
+               zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )     
+                       ! zc1 : log(zCze)  = log (1.12  * zchl**0.803)
+               zc1   = 0.113328685307 + 0.803 * zlogc                         
+                       ! zc2 : log(zCtot) = log(40.6  * zchl**0.459)
+               zc2   = 3.703768066608 + 0.459 * zlogc                           
+                       ! zc3 : log(zze)   = log(568.2 * zCtot**(-0.746))
+               zc3   = 6.34247346942  - 0.746 * zc2                           
+                       ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293))
+               IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2        
+               !   
+               ze0(ji,jj) = zlogc                                                 ! ze0 = log(zchl)
+               ze1(ji,jj) = EXP( zc1 )                                            ! ze1 = zCze
+               ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi
+               ze3(ji,jj) = EXP( - zc3 )                                          ! ze3 = 1/zze
+            END_2D
+            
+!
+            DO_3D( 0, 0, 0, 0, 1, nksr + 1 )
+               ! zchl    = ALOG( ze0(ji,jj) )
+               zlogc = ze0(ji,jj)
+               !
+               zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
+               zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
+               zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
+               ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
+               !
+               zCze   = ze1(ji,jj)
+               zrdpsi = ze2(ji,jj)                                                 ! 1/zdelpsi
+               zpsi   = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm)                           ! gdepw/zze
+               !
+               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
+               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) )
+               ! Convert chlorophyll value to attenuation coefficient look-up table index
+               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
+            END_3D
+         ELSE                                !* constant chlorophyll
+            zchl = 0.05
+            ! NB. make sure constant value is such that: 
+            zchl = MIN( 10. , MAX( 0.03, zchl ) )
+            ! Convert chlorophyll value to attenuation coefficient look-up table index
+            zlui = 41 + 20.*LOG10(zchl) + 1.e-15
             DO jk = 1, nksr + 1
-               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
-                  DO ji = 2, jpim1
-                     zchl    = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
-                     zCtot   = 40.6  * zchl**0.459
-                     zze     = 568.2 * zCtot**(-0.746)
-                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
-                     zpsi    = gdepw(ji,jj,jk,Kmm) / zze
-                     !
-                     zlogc   = LOG( zchl )
-                     zlogc2  = zlogc * zlogc
-                     zlogc3  = zlogc * zlogc * zlogc
-                     zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
-                     zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
-                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
-                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
-                     zCze    = 1.12  * (zchl)**0.803 
-                     !
-                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
-                  END DO
-                  !
-               END DO
+               ztmp3d(:,:,jk) = zlui 
             END DO
-         ELSE                                !* constant chrlorophyll
-           DO jk = 1, nksr + 1
-              zchl3d(:,:,jk) = 0.05 
-            ENDDO
          ENDIF
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
-         DO_2D_00_00
-            ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
-            ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
-            ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
-            ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
-            zea(ji,jj,1) =          qsr(ji,jj)
+         DO_2D( 0, 0, 0, 0 )
+            ze0(ji,jj) = rn_abs * qsr(ji,jj)
+            ze1(ji,jj) = zcoef  * qsr(ji,jj)
+            ze2(ji,jj) = zcoef  * qsr(ji,jj)
+            ze3(ji,jj) = zcoef  * qsr(ji,jj)
+            ! store the surface SW radiation; re-use the surface ztmp3d array 
+            ! since the surface attenuation coefficient is not used
+            ztmp3d(ji,jj,1) =       qsr(ji,jj)
          END_2D
          !
-         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
-            DO_2D_00_00
-               zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
-               irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-               zekb(ji,jj) = rkrgb(1,irgb)
-               zekg(ji,jj) = rkrgb(2,irgb)
-               zekr(ji,jj) = rkrgb(3,irgb)
-            END_2D
-
-            DO_2D_00_00
-               zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       )
-               zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) )
-               zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) )
-               zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) )
-               ze0(ji,jj,jk) = zc0
-               ze1(ji,jj,jk) = zc1
-               ze2(ji,jj,jk) = zc2
-               ze3(ji,jj,jk) = zc3
-               zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
-            END_2D
-         END DO
-         !
-         DO_3D_00_00( 1, nksr )
-            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
+         !                                    !* interior equi-partition in R-G-B depending on vertical profile of Chl
+         DO_3D( 0, 0, 0, 0, 2, nksr + 1 )
+            ze3t = e3t(ji,jj,jk-1,Kmm)
+            irgb = NINT( ztmp3d(ji,jj,jk) )
+            zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r )
+            zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
+            zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
+            zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
+            ze0(ji,jj) = zc0
+            ze1(ji,jj) = zc1
+            ze2(ji,jj) = zc2
+            ze3(ji,jj) = zc3
+            ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
          END_3D
          !
-         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
+         DO_3D( 0, 0, 0, 0, 1, nksr )          !* now qsr induced heat content
+            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
+         END_3D
+         !
+         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 
          !
       CASE( np_2BD  )            !==  2-bands fluxes  ==!
@@ -236 +256 @@
          zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
          zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
-         DO_3D_00_00( 1, nksr )
+         DO_3D( 0, 0, 0, 0, 1, nksr )             ! solar heat absorbed at T-point in the top 400m 
             zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
             zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
@@ -245 +265 @@
       !
       !                          !-----------------------------!
-      DO_3D_00_00( 1, nksr )
+      !                          !  update to the temp. trend  !
+      !                          !-----------------------------!
+      DO_3D( 0, 0, 0, 0, 1, nksr )
          pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
-            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
+            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) )   &
+            &                             / e3t(ji,jj,jk,Kmm)
       END_3D
       !
       ! sea-ice: store the 1st ocean level attenuation coefficient
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
          ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
@@ -396 +419 @@
          IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
          !
+         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
+         !                                   
+         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
+         !
+         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
+         !
       END SELECT
       !
@@ -402 +431 @@
       ! 1st ocean level attenuation coefficient (used in sbcssm)
       IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
+         CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
       ELSE
          fraqsr_1lev(:,:) = 1._wp   ! default : no penetration

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/trasbc.F90

@@ -43 +43 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -111 +112 @@
             zfact = 0.5_wp
             sbc_tsc(:,:,:) = 0._wp
-            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before heat content sbc trend
-            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content sbc trend
+            CALL iom_get( numror, jpdom_auto, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before heat content sbc trend
+            CALL iom_get( numror, jpdom_auto, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content sbc trend
          ELSE                                   ! No restart or restart not found: Euler forward time stepping
             zfact = 1._wp
@@ -123 +124 @@
       ENDIF
       !                             !==  Now sbc tracer content fields  ==!
-      DO_2D_01_00
+      DO_2D( 0, 1, 0, 0 )
          sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
          sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
       END_2D
       IF( ln_linssh ) THEN                !* linear free surface  
-         DO_2D_01_00
+         DO_2D( 0, 1, 0, 0 )                    !==>> add concentration/dilution effect due to constant volume cell
             sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
             sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
-         END_2D
+         END_2D                                 !==>> output c./d. term
          IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
          IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
@@ -137 +138 @@
       !
       DO jn = 1, jpts               !==  update tracer trend  ==!
-         DO_2D_01_00
-            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,Kmm)
+         DO_2D( 0, 1, 0, 0 )
+            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) )    &
+               &                                                / e3t(ji,jj,1,Kmm)
          END_2D
       END DO
@@ -155 +157 @@
       IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff 
          zfact = 0.5_wp
-         DO_2D_01_00
+         DO_2D( 0, 1, 0, 0 )
             IF( rnf(ji,jj) /= 0._wp ) THEN
                zdep = zfact / h_rnf(ji,jj)
@@ -180 +182 @@
           !
          IF( ln_linssh ) THEN 
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
                pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim
@@ -186 +188 @@
             END_2D
          ELSE
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) )
                pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/trazdf.F90

@@ -37 +37 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -84 +85 @@
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) &
-               &          / (e3t(:,:,jk,Kmm)*rDt) ) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) &
-              &           / (e3t(:,:,jk,Kmm)*rDt) ) - ztrds(:,:,jk)
+            ztrdt(:,:,jk) = (   (  pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa)     &
+               &                 - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb)  )  &
+               &              / (  e3t(:,:,jk,Kmm)*rDt  )   )                 &
+               &          - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = (   (  pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa)     &
+               &                 - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb)  )  &
+               &             / (   e3t(:,:,jk,Kmm)*rDt  )   )                 &
+               &          - ztrds(:,:,jk)
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
-         CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1. )
+         CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1.0_wp , ztrds, 'T', 1.0_wp )
 !!gm
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt )
@@ -156 +161 @@
             IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
                IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
-                  DO_3D_00_00( 2, jpkm1 )
+                  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_3D_00_00( 2, jpkm1 )
+                  DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                      zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
                   END_3D
@@ -168 +173 @@
             ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
             IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection
-               DO_3D_00_00( 1, jpkm1 )
+               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)
@@ -177 +182 @@
                END_3D
             ELSE
-               DO_3D_00_00( 1, jpkm1 )
+               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)
@@ -203 +208 @@
             !   used as a work space array: its value is modified.
             !
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )      !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) ! done one for all passive tracers (so included in the IF instruction)
                zwt(ji,jj,1) = zwd(ji,jj,1)
             END_2D
-            DO_3D_00_00( 2, jpkm1 )
+            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
@@ -212 +217 @@
          ENDIF 
          !         
-         DO_2D_00_00
-            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)
+         DO_2D( 0, 0, 0, 0 )         !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+            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_00_00( 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
+         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_2D_00_00
+         DO_2D( 0, 0, 0, 0 )         !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
             pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
          END_2D
-         DO_3DS_00_00( jpk-2, 1, -1 )
+         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)

NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/zpshde.F90

@@ -32 +32 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -65 +66 @@
       !!              ___ |   |   |           ___  |   |   |
       !!                  
-      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
-      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
-      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
+      !!      case 1->   e3w(i+1,:,:,Kmm) >= e3w(i,:,:,Kmm) ( and e3w(:,j+1,:,Kmm) >= e3w(:,j,:,Kmm) ) then
+      !!          t~ = t(i+1,j  ,k) + (e3w(i+1,j,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Ti+1)/e3w(i+1,j,k,Kmm)
+      !!        ( t~ = t(i  ,j+1,k) + (e3w(i,j+1,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Tj+1)/e3w(i,j+1,k,Kmm)  )
       !!          or
-      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
-      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
-      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
+      !!      case 2->   e3w(i+1,:,:,Kmm) <= e3w(i,:,:,Kmm) ( and e3w(:,j+1,:,Kmm) <= e3w(:,j,:,Kmm) ) then
+      !!          t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i+1,j,k,Kmm)) * dk(Ti)/e3w(i,j,k,Kmm)
+      !!        ( t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i,j+1,k,Kmm)) * dk(Tj)/e3w(i,j,k,Kmm) )
       !!          Idem for di(s) and dj(s)          
       !!
@@ -106 +107 @@
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
          !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
             ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
-!!gm BUG ? when applied to before fields, e3w(:,:,:,Kbb) should be used....
+!!gm BUG ? when applied to before fields, e3w(:,:,k,Kbb) should be used....
             ze3wu = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
             ze3wv = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
@@ -145 +146 @@
       END DO
       !
-      CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1. , pgtv(:,:,:), 'V', -1. )   ! Lateral boundary cond.
+      CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
       !                
       IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
          pgru(:,:) = 0._wp
          pgrv(:,:) = 0._wp                ! depth of the partial step level
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -166 +167 @@
          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
          !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )              ! Gradient of density at the last level
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -178 +179 @@
             ENDIF
          END_2D
-         CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1. , pgrv , 'V', -1. )   ! Lateral boundary conditions
+         CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions
          !
       END IF
@@ -214 +215 @@
       !!              ___ |   |   |           ___  |   |   |
       !!                  
-      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
-      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
-      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
+      !!      case 1->   e3w(i+1,j,k,Kmm) >= e3w(i,j,k,Kmm) ( and e3w(i,j+1,k,Kmm) >= e3w(i,j,k,Kmm) ) then
+      !!          t~ = t(i+1,j  ,k) + (e3w(i+1,j  ,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Ti+1)/e3w(i+1,j  ,k,Kmm)
+      !!        ( t~ = t(i  ,j+1,k) + (e3w(i  ,j+1,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Tj+1)/e3w(i  ,j+1,k,Kmm)  )
       !!          or
-      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
-      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
-      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
+      !!      case 2->   e3w(i+1,j,k,Kmm) <= e3w(i,j,k,Kmm) ( and e3w(i,j+1,k,Kmm) <= e3w(i,j,k,Kmm) ) then
+      !!          t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i+1,j  ,k,Kmm)) * dk(Ti)/e3w(i,j,k,Kmm)
+      !!        ( t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i  ,j+1,k,Kmm)) * dk(Tj)/e3w(i,j,k,Kmm) )
       !!          Idem for di(s) and dj(s)          
       !!
@@ -261 +262 @@
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
          !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
 
             iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
@@ -301 +302 @@
       END DO
       !
-      CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1. , pgtv(:,:,:), 'V', -1. )   ! Lateral boundary cond.
+      CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
 
       ! horizontal derivative of density anomalies (rd)
@@ -307 +308 @@
          pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
          !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
 
             iku = mbku(ji,jj)
@@ -328 +329 @@
          CALL eos( ztj, zhj, zrj )
 
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )            ! Gradient of density at the last level
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -343 +344 @@
          END_2D
 
-         CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1. , pgrv , 'V', -1. )   ! Lateral boundary conditions
+         CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions
          !
       END IF
@@ -350 +351 @@
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
             iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1
             ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1
@@ -356 +357 @@
             ! (ISF) case partial step top and bottom in adjacent cell in vertical
             ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
-            ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
+            ! in this case e3w(i,j,k,Kmm) - e3w(i,j+1,k,Kmm) is not the distance between Tj~ and Tj
             ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
             ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
@@ -394 +395 @@
          !
       END DO
-      CALL lbc_lnk_multi( 'zpshde', pgtui(:,:,:), 'U', -1. , pgtvi(:,:,:), 'V', -1. )   ! Lateral boundary cond.
+      CALL lbc_lnk_multi( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond.
 
       IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
          !
          pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp;
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )
 
             iku = miku(ji,jj)
@@ -419 +420 @@
          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
          !
-         DO_2D_10_10
+         DO_2D( 1, 0, 1, 0 )              ! Gradient of density at the last level
             iku = miku(ji,jj) 
             ikv = mikv(ji,jj) 
@@ -433 +434 @@
 
          END_2D
-         CALL lbc_lnk_multi( 'zpshde', pgrui, 'U', -1. , pgrvi, 'V', -1. )   ! Lateral boundary conditions
+         CALL lbc_lnk_multi( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp )   ! Lateral boundary conditions
          !
       END IF  

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev@HEAD      ext/AGRIF
+^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
 
 # SETTE
-^/utils/CI/sette@HEAD         sette
+^/utils/CI/sette@13507        sette