Index: /NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzad.F90
===================================================================
--- /NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzad.F90	(revision 14427)
+++ /NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzad.F90	(revision 14428)
@@ -4,6 +4,7 @@
    !! Ocean dynamics : vertical advection trend
    !!======================================================================
-   !! History :  OPA  ! 1991-01  (G. Madec) Original code
-   !!   NEMO     0.5  ! 2002-07  (G. Madec) Free form, F90
+   !! History :  OPA  !  1991-01  (G. Madec) Original code
+   !!   NEMO     0.5  !  2002-07  (G. Madec) Free form, F90
+   !!            4.5  !  2021-01  (S. Techene, G. Madec) memory optimization
    !!----------------------------------------------------------------------
 
@@ -54,13 +55,12 @@
       !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
-      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step inedx
-      INTEGER                             , INTENT( in )  ::  Kmm, Krhs        ! ocean time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
+      INTEGER                             , INTENT(in   ) ::    kt, Kmm, Krhs   ! ocean time-step & time-level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::    puu, pvv        ! ocean velocities and RHS of momentum equation
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zua, zva     ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zww
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwuw, zwvw
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
+      REAL(wp) ::   zWf, zWfi, zzWfu, zzWdzU   ! local scalars
+      REAL(wp) ::        zWfj, zzWfv, zzWdzV   !   -      -
+      REAL(wp), DIMENSION(A2D(0))             ::   zWdzU, zWdzV   ! 2D inner workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   ! 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -71,5 +71,5 @@
          IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme'
       ENDIF
-
+      !
       IF( l_trddyn )   THEN           ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
@@ -77,34 +77,40 @@
          ztrdv(:,:,:) = pvv(:,:,:,Krhs)
       ENDIF
-
-      DO jk = 2, jpkm1                ! Vertical momentum advection at level w and u- and v- vertical
-         DO_2D( 0, 1, 0, 1 )              ! vertical fluxes
-          IF( ln_vortex_force ) THEN
-            zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
-          ELSE
-            zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
-          ENDIF
-         END_2D
-         DO_2D( 0, 0, 0, 0 )              ! vertical momentum advection at w-point
-            zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) )
-            zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) )
-         END_2D
-      END DO
       !
-      ! Surface and bottom advective fluxes set to zero
-      DO_2D( 0, 0, 0, 0 )
-         zwuw(ji,jj, 1 ) = 0._wp
-         zwvw(ji,jj, 1 ) = 0._wp
-         zwuw(ji,jj,jpk) = 0._wp
-         zwvw(ji,jj,jpk) = 0._wp
-      END_2D
+      !                                !==  vertical momentum advection ==!   at u- and v-points
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Vertical momentum advection at u- and v-points
-         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
-            &                                      / e3u(ji,jj,jk,Kmm)
-         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj)   &
-            &                                      / e3v(ji,jj,jk,Kmm)
+      zWdzU(A2D(0)) = 0._wp                  ! set surface (jk=1) vertical advection to zero
+      zWdzV(A2D(0)) = 0._wp
+      !
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
+         !                                ! vertical transport at jk+1 uw/vw-level (x2): 2*mi/j[e1e2t*(We)]
+         IF( ln_vortex_force ) THEN             ! We = ww+wsd
+            zWf  = e1e2t(ji  ,jj  ) * ( ww(ji  ,jj  ,jk+1) + wsd(ji  ,jj  ,jk+1) )
+            zWfi = e1e2t(ji+1,jj  ) * ( ww(ji+1,jj  ,jk+1) + wsd(ji+1,jj  ,jk+1) )
+            zWfj = e1e2t(ji  ,jj+1) * ( ww(ji  ,jj+1,jk+1) + wsd(ji  ,jj+1,jk+1) )
+         ELSE                                   ! We = ww
+            zWf  = e1e2t(ji  ,jj  ) *   ww(ji  ,jj  ,jk+1)
+            zWfi = e1e2t(ji+1,jj  ) *   ww(ji+1,jj  ,jk+1)
+            zWfj = e1e2t(ji  ,jj+1) *   ww(ji  ,jj+1,jk+1)
+         ENDIF
+         zzWfu = zWfi + zWf                     ! averaging at uw- and vw-points (x2)
+         zzWfv = zWfj + zWf 
+         !                                ! vertical advection at jk+1 uw-level (x4): zzWfu/v*dk+1[u/v] 
+         zzWdzU = zzWfu * ( puu(ji,jj,jk,Kmm) - puu(ji,jj,jk+1,Kmm) )
+         zzWdzV = zzWfv * ( pvv(ji,jj,jk,Kmm) - pvv(ji,jj,jk+1,Kmm) )
+         !
+         !                                ! vertical advection  at jk u/v-level 
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
+            &                                    * ( zWdzU(ji,jj) + zzWdzU )
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
+            &                                    * ( zWdzV(ji,jj) + zzWdzV )
+         !
+         IF( jk /= jpkm1 ) THEN             ! save for next level computation 
+            zWdzU(ji,jj) = zzWdzU 
+            zWdzV(ji,jj) = zzWdzV
+         ENDIF
+         !
       END_3D
-
+      !
       IF( l_trddyn ) THEN             ! save the vertical advection trends for diagnostic
          ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
