Index: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynldf_lap_blp.F90	(revision 14037)
+++ NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynldf_lap_blp.F90	(revision 14062)
@@ -5,4 +5,5 @@
    !!======================================================================
    !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
+   !!           4.0  ! 2020-04  (A. Nasser, G. Madec)  Add symmetric mixing tensor 
    !!----------------------------------------------------------------------
 
@@ -19,5 +20,6 @@
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-
+   USE lib_mpp
+   
    IMPLICIT NONE
    PRIVATE
@@ -47,4 +49,6 @@
       !!
       !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
+      !!
+      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/ 
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
@@ -57,5 +61,6 @@
       REAL(wp) ::   zsign        ! local scalars
       REAL(wp) ::   zua, zva     ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zcur, zdiv
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms
       !!----------------------------------------------------------------------
       !
@@ -70,30 +75,78 @@
       ENDIF
       !
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         DO_2D( 0, 1, 0, 1 )
-            !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
-            zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask
-               &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
-               &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
-            !                                      ! ahm * div        (computed from 2 to jpi/jpj)
-            zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask
-               &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
-               &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
-         END_2D
+      SELECT CASE( nn_dynldf_typ )  
+      !              
+      CASE ( np_typ_rot )       !==  Vorticity-Divergence operator  ==!
          !
-         DO_2D( 0, 0, 0, 0 )                       ! - curl( curl) + grad( div )
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
-               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
-               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      )
+         ALLOCATE( zcur(jpi,jpj) , zdiv(jpi,jpj) )
+         !
+         DO jk = 1, jpkm1                                 ! Horizontal slab
+            !
+            DO_2D( 0, 1, 0, 1 )
+               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
+               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask
+                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
+               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
+               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask
+                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
+                  &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
+            END_2D
+            !
+            DO_2D( 0, 0, 0, 0 )                       ! - curl( curl) + grad( div )
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
+                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
+                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      )
                !
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use
-               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
-               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      )
-         END_2D
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use
+                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
+                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      )
+            END_2D
+            !
+         END DO                                           !   End of slab
+         !
+         DEALLOCATE( zcur , zdiv )
+         !
+      CASE ( np_typ_sym )       !==  Symmetric operator  ==!
+         !
+         ALLOCATE( zten(jpi,jpj) , zshe(jpi,jpj) )
+         !
+         DO jk = 1, jpkm1                                 ! Horizontal slab
+            !
+            DO_2D( 0, 1, 0, 1 )
+               !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
+               zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              &
+                  &     * (    e1f(ji-1,jj-1)    * r1_e2f(ji-1,jj-1)                                             &
+                  &         * ( pu(ji-1,jj  ,jk) * r1_e1u(ji-1,jj  )  - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) )   &
+                  &         +  e2f(ji-1,jj-1)    * r1_e1f(ji-1,jj-1)                                             &
+                  &         * ( pv(ji  ,jj-1,jk) * r1_e2v(ji  ,jj-1)  - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) )   ) 
+               !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask
+               zten(ji,jj)    = ahmt(ji,jj,jk)                                                       &
+                  &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         &
+                  &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   &
+                  &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         &
+                  &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )   
+            END_2D
+            !
+            DO_2D( 0, 0, 0, 0 )
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               &
+                  &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       &
+                  &            - zten(ji  ,jj  ) * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     &                                                    
+                  &        + (   zshe(ji  ,jj  ) * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
+                  &            - zshe(ji  ,jj-1) * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )   
+               !
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                               &
+                  &    * (   (   zshe(ji  ,jj  ) * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
+                  &            - zshe(ji-1,jj  ) * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj)     &
+                  &        - (   zten(ji  ,jj+1) * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                       &
+                  &            - zten(ji  ,jj  ) * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) )
+               !
+            END_2D
+            !
+         END DO
+         !
+         DEALLOCATE( zten , zshe )
+         !
+      END SELECT
       !
    END SUBROUTINE dyn_ldf_lap
