Changeset 5836 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN

Timestamp:

2015-10-26T15:49:40+01:00 (8 years ago)

Author:

cetlod

Message:

merge the simplification branch onto the trunk, see ticket #1612

Location:

trunk/NEMOGCM/NEMO/OPA_SRC/DYN

Files:

• divcur.F90 (deleted)
• divhor.F90 (copied) (copied from branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90)
• dynadv.F90 (modified) (4 diffs)
• dynadv_cen2.F90 (modified) (3 diffs)
• dynadv_ubs.F90 (modified) (3 diffs)
• dynhpg.F90 (modified) (1 diff)
• dynldf.F90 (modified) (10 diffs)
• dynldf_bilap.F90 (deleted)
• dynldf_bilapg.F90 (deleted)
• dynldf_iso.F90 (modified) (21 diffs)
• dynldf_lap.F90 (deleted)
• dynldf_lap_blp.F90 (copied) (copied from branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap_blp.F90)
• dynnept.F90 (deleted)
• dynspg.F90 (modified) (7 diffs)
• dynspg_exp.F90 (modified) (1 diff)
• dynspg_flt.F90 (modified) (3 diffs)
• dynspg_oce.F90 (modified) (1 diff)
• dynspg_ts.F90 (modified) (28 diffs)
• dynvor.F90 (modified) (27 diffs)
• dynzad.F90 (modified) (8 diffs)
• dynzdf.F90 (modified) (7 diffs)
• dynzdf_imp.F90 (modified) (2 diffs)
• sshwzv.F90 (modified) (10 diffs)

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

@@ -76 +76 @@
       CASE ( 3 )   
                       CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme
-      !
-      CASE (-1 )                                            ! esopa: test all possibility with control print
-                      CALL dyn_keg     ( kt, nn_dynkeg )
-                      CALL dyn_zad     ( kt )
-                      CALL dyn_adv_cen2( kt )
-                      CALL dyn_adv_ubs ( kt )
       END SELECT
       !
@@ -126 +120 @@
       IF( ln_dynadv_cen2 )   ioptio = ioptio + 1
       IF( ln_dynadv_ubs  )   ioptio = ioptio + 1
-      IF( lk_esopa       )   ioptio =          1
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
@@ -139 +132 @@
       IF( ln_dynadv_cen2 )   nadv =  2
       IF( ln_dynadv_ubs  )   nadv =  3
-      IF( lk_esopa       )   nadv = -1
 
       IF(lwp) THEN                    ! Print the choice
@@ -151 +143 @@
          IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
          IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
-         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection formulation'
       ENDIF
       !

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90

@@ -90 +90 @@
          DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+               zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk)
+               zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk)
                !
                ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
@@ -114 +114 @@
       DO jk = 1, jpkm1                       ! ==================== !
          !                                         ! Vertical volume fluxesÊ
-         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
+         zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk)
          !
          IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
@@ -144 +144 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
-                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
+                  &  / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
                va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
-                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+                  &  / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
             END DO
          END DO

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

@@ -181 +181 @@
          DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+               zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk)
+               zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk)
                !
                ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
@@ -203 +203 @@
       DO jk = 1, jpkm1                       ! ==================== !
          !                                         ! Vertical volume fluxesÊ
-         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
+         zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk)
          !
          IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
@@ -233 +233 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
-                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
+                  &  / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
                va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
-                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+                  &  / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
             END DO
          END DO

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

@@ -1046 +1046 @@
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
-          zsshu_n(ji,jj) = (e12u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * &
-                         & r1_e12u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
-          zsshv_n(ji,jj) = (e12v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * &
-                         & r1_e12v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+          zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &
+                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
+          zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &
+                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
         END DO
       END DO

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

@@ -4 +4 @@
    !! Ocean physics:  lateral diffusivity trends 
    !!=====================================================================
-   !! History :  9.0  !  05-11  (G. Madec)  Original code (new step architecture)
+   !! History :  2.0  ! 2005-11  (G. Madec)  Original code (new step architecture)
+   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
+   !!                 !                                  add velocity dependent coefficient and optional read in file
    !!----------------------------------------------------------------------
 
@@ -14 +16 @@
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
-   USE ldfdyn_oce     ! ocean dynamics lateral physics
-   USE ldftra_oce     ! ocean tracers  lateral physics
-   USE ldfslp         ! lateral mixing: slopes of mixing orientation
-   USE dynldf_bilapg  ! lateral mixing            (dyn_ldf_bilapg routine)
-   USE dynldf_bilap   ! lateral mixing            (dyn_ldf_bilap  routine)
-   USE dynldf_iso     ! lateral mixing            (dyn_ldf_iso    routine)
-   USE dynldf_lap     ! lateral mixing            (dyn_ldf_lap    routine)
+   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
+   USE ldfslp         ! lateral diffusion: slopes of mixing orientation
+   USE dynldf_lap_blp ! lateral mixing   (dyn_ldf_lap & dyn_ldf_blp routines)
+   USE dynldf_iso     ! lateral mixing                 (dyn_ldf_iso routine )
    USE trd_oce        ! trends: ocean variables
-   USE trddyn         ! trend manager: dynamics   (trd_dyn        routine)
+   USE trddyn         ! trend manager: dynamics   (trd_dyn      routine)
    !
    USE prtctl         ! Print control
@@ -28 +27 @@
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -37 +36 @@
    PUBLIC   dyn_ldf_init  ! called by opa  module 
 
-   INTEGER ::   nldf = -2   ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals)
+   !                      ! Flag to control the type of lateral viscous operator
+   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in setting the operator
+   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral viscous trend)
+   !                          !!      laplacian     !    bilaplacian    !
+   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator
+   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       ! iso-neutral or geopotential operator
+
+   INTEGER ::   nldf   ! type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
 
    !! * Substitutions
@@ -43 +49 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -62 +68 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_ldf')
       !
-      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
+      IF( l_trddyn )   THEN                      ! temporary save of momentum trends
+         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
          ztrdu(:,:,:) = ua(:,:,:) 
          ztrdv(:,:,:) = va(:,:,:) 
@@ -70 +76 @@
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
       !
-      CASE ( 0 )    ;   CALL dyn_ldf_lap    ( kt )      ! iso-level laplacian
-      CASE ( 1 )    ;   CALL dyn_ldf_iso    ( kt )      ! rotated laplacian (except dk[ dk[.] ] part)
-      CASE ( 2 )    ;   CALL dyn_ldf_bilap  ( kt )      ! iso-level bilaplacian
-      CASE ( 3 )    ;   CALL dyn_ldf_bilapg ( kt )      ! s-coord. horizontal bilaplacian
-      CASE ( 4 )                                        ! iso-level laplacian + bilaplacian
-         CALL dyn_ldf_lap    ( kt )
-         CALL dyn_ldf_bilap  ( kt )
-      CASE ( 5 )                                        ! rotated laplacian + bilaplacian (s-coord)
-         CALL dyn_ldf_iso    ( kt )
-         CALL dyn_ldf_bilapg ( kt )
+      CASE ( np_lap   )    ;   CALL dyn_ldf_lap  ( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian
+      CASE ( np_lap_i )    ;   CALL dyn_ldf_iso  ( kt )                         ! rotated      laplacian
+      CASE ( np_blp   )    ;   CALL dyn_ldf_blp  ( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian
       !
-      CASE ( -1 )                                       ! esopa: test all possibility with control print
-                        CALL dyn_ldf_lap    ( kt )
-                        CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf0 - Ua: ', mask1=umask,   &
-            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-                        CALL dyn_ldf_iso    ( kt )
-                        CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf1 - Ua: ', mask1=umask,   &
-            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-                        CALL dyn_ldf_bilap  ( kt )
-                        CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf2 - Ua: ', mask1=umask,   &
-            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-                        CALL dyn_ldf_bilapg ( kt )
-                        CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask,   &
-            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-      !
-      CASE ( -2 )                                       ! neither laplacian nor bilaplacian schemes used
-         IF( kt == nit000 ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'dyn_ldf : no lateral diffusion on momentum setup'
-            IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         ENDIF
       END SELECT
 
@@ -107 +86 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )
+         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
       ENDIF
       !                                          ! print sum trends (used for debugging)
@@ -126 +105 @@
       INTEGER ::   ioptio, ierr         ! temporary integers 
       !!----------------------------------------------------------------------
-    
+      !
       !                                   ! Namelist nam_dynldf: already read in ldfdyn module
-
+      !
       IF(lwp) THEN                        ! Namelist print
          WRITE(numout,*)
@@ -134 +113 @@
          WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '       Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
-         WRITE(numout,*) '          laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
-         WRITE(numout,*) '          bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
-         WRITE(numout,*) '          iso-level                   ln_dynldf_level = ', ln_dynldf_level
-         WRITE(numout,*) '          horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
-         WRITE(numout,*) '          iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
+         WRITE(numout,*) '          laplacian operator          ln_dynldf_lap = ', ln_dynldf_lap
+         WRITE(numout,*) '          bilaplacian operator        ln_dynldf_blp = ', ln_dynldf_blp
+         WRITE(numout,*) '          iso-level                   ln_dynldf_lev = ', ln_dynldf_lev
+         WRITE(numout,*) '          horizontal (geopotential)   ln_dynldf_hor = ', ln_dynldf_hor
+         WRITE(numout,*) '          iso-neutral                 ln_dynldf_iso = ', ln_dynldf_iso
       ENDIF
-
-      !                                   ! control the consistency
+      !                                   ! use of lateral operator or not
+      nldf = np_ERROR
       ioptio = 0
-      IF( ln_dynldf_lap   )   ioptio = ioptio + 1
-      IF( ln_dynldf_bilap )   ioptio = ioptio + 1
-      IF( ioptio <  1 ) CALL ctl_warn( '          neither laplacian nor bilaplacian operator set for dynamics' )
-      ioptio = 0
-      IF( ln_dynldf_level )   ioptio = ioptio + 1
-      IF( ln_dynldf_hor   )   ioptio = ioptio + 1
-      IF( ln_dynldf_iso   )   ioptio = ioptio + 1
-      IF( ioptio >  1 ) CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
-
-      IF( ln_dynldf_iso .AND. ln_traldf_hor ) CALL ctl_stop &
-      &   ( 'Not sensible to use geopotential diffusion for tracers with isoneutral diffusion for dynamics' )
-
-      !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals
-      ierr = 0
-      IF ( ln_dynldf_lap ) THEN      ! laplacian operator
-         IF ( ln_zco ) THEN                ! z-coordinate
-            IF ( ln_dynldf_level )   nldf = 0      ! iso-level  (no rotation)
-            IF ( ln_dynldf_hor   )   nldf = 0      ! horizontal (no rotation)
-            IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+      IF( ln_dynldf_lap )   ioptio = ioptio + 1
+      IF( ln_dynldf_blp )   ioptio = ioptio + 1
+      IF( ioptio >  1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on momentum' )
+      IF( ioptio == 0   )   nldf = np_no_ldf     ! No lateral mixing operator
+      !
+      IF( nldf /= np_no_ldf ) THEN        ! direction ==>> type of operator  
+         ioptio = 0
+         IF( ln_dynldf_lev )   ioptio = ioptio + 1
+         IF( ln_dynldf_hor )   ioptio = ioptio + 1
+         IF( ln_dynldf_iso )   ioptio = ioptio + 1
+         IF( ioptio >  1   )   CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
+         IF( ioptio == 0   )   CALL ctl_stop( '          use at least ONE direction (level/hor/iso)' )
+         !
+         !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals
+         ierr = 0
+         IF ( ln_dynldf_lap ) THEN      ! laplacian operator
+            IF ( ln_zco ) THEN                ! z-coordinate
+               IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation)
+            ENDIF
+            IF ( ln_zps ) THEN             ! z-coordinate with partial step
+               IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level              (no rotation)
+               IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level              (no rotation)
+               IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation)
+            ENDIF
+            IF ( ln_sco ) THEN             ! s-coordinate
+               IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_hor )   nldf = np_lap_i   ! horizontal             (   rotation)
+               IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation)
+            ENDIF
          ENDIF
-         IF ( ln_zps ) THEN             ! z-coordinate
-            IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed
-            IF ( ln_dynldf_hor   )   nldf = 0      ! horizontal (no rotation)
-            IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         !
+         IF( ln_dynldf_blp ) THEN          ! bilaplacian operator
+            IF ( ln_zco ) THEN                ! z-coordinate
+               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_hor )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
+            ENDIF
+            IF ( ln_zps ) THEN             ! z-coordinate with partial step
+               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation)
+               IF ( ln_dynldf_hor )   nldf = np_blp     ! iso-level              (no rotation)
+               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
+            ENDIF
+            IF ( ln_sco ) THEN             ! s-coordinate
+               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation)
+               IF ( ln_dynldf_hor )   ierr = 2          ! horizontal             (   rotation)
+               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
+            ENDIF
          ENDIF
-         IF ( ln_sco ) THEN             ! s-coordinate
-            IF ( ln_dynldf_level )   nldf = 0      ! iso-level  (no rotation)
-            IF ( ln_dynldf_hor   )   nldf = 1      ! horizontal (   rotation)
-            IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation)
-         ENDIF
-      ENDIF
-
-      IF( ln_dynldf_bilap ) THEN      ! bilaplacian operator
-         IF ( ln_zco ) THEN                ! z-coordinate
-            IF ( ln_dynldf_level )   nldf = 2      ! iso-level  (no rotation)
-            IF ( ln_dynldf_hor   )   nldf = 2      ! horizontal (no rotation)
-            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
-         ENDIF
-         IF ( ln_zps ) THEN             ! z-coordinate
-            IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed 
-            IF ( ln_dynldf_hor   )   nldf = 2      ! horizontal (no rotation)
-            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
-         ENDIF
-         IF ( ln_sco ) THEN             ! s-coordinate
-            IF ( ln_dynldf_level )   nldf = 2      ! iso-level  (no rotation)
-            IF ( ln_dynldf_hor   )   nldf = 3      ! horizontal (   rotation)
-            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
-         ENDIF
-      ENDIF
-      
-      IF( ln_dynldf_lap .AND. ln_dynldf_bilap ) THEN  ! mixed laplacian and bilaplacian operators
-         IF ( ln_zco ) THEN                ! z-coordinate
-            IF ( ln_dynldf_level )   nldf = 4      ! iso-level  (no rotation)
-            IF ( ln_dynldf_hor   )   nldf = 4      ! horizontal (no rotation)
-            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
-         ENDIF
-         IF ( ln_zps ) THEN             ! z-coordinate
-            IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed 
-            IF ( ln_dynldf_hor   )   nldf = 4      ! horizontal (no rotation)
-            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
-         ENDIF
-         IF ( ln_sco ) THEN             ! s-coordinate
-            IF ( ln_dynldf_level )   nldf = 4      ! iso-level  (no rotation)
-            IF ( ln_dynldf_hor   )   nldf = 5      ! horizontal (   rotation)
-            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
-         ENDIF
-      ENDIF
-
-      IF( lk_esopa )                 nldf = -1     ! esopa test
-
-      IF( ierr == 1 )   CALL ctl_stop( 'iso-level in z-coordinate - partial step, not allowed' )
-      IF( ierr == 2 )   CALL ctl_stop( 'isoneutral bilaplacian operator does not exist' )
-      IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
-         IF( .NOT.lk_ldfslp )   CALL ctl_stop( 'the rotation of the diffusive tensor require key_ldfslp' )
+         !
+         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' )
+         !
+         IF( nldf == np_lap_i )   l_ldfslp = .TRUE.      ! rotation require the computation of the slopes
+         !
       ENDIF
 
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nldf == -2 )   WRITE(numout,*) '              neither laplacian nor bilaplacian schemes used'
-         IF( nldf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
-         IF( nldf ==  0 )   WRITE(numout,*) '              laplacian operator'
-         IF( nldf ==  1 )   WRITE(numout,*) '              rotated laplacian operator'
-         IF( nldf ==  2 )   WRITE(numout,*) '              bilaplacian operator'
-         IF( nldf ==  3 )   WRITE(numout,*) '              rotated bilaplacian'
-         IF( nldf ==  4 )   WRITE(numout,*) '              laplacian and bilaplacian operators'
-         IF( nldf ==  5 )   WRITE(numout,*) '              rotated laplacian and bilaplacian operators'
+         IF( nldf == np_no_ldf )   WRITE(numout,*) '              NO lateral viscosity'
+         IF( nldf == np_lap    )   WRITE(numout,*) '              iso-level laplacian operator'
+         IF( nldf == np_lap_i  )   WRITE(numout,*) '              rotated laplacian operator with iso-level background'
+         IF( nldf == np_blp    )   WRITE(numout,*) '              iso-level bi-laplacian operator'
       ENDIF
       !

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

@@ -2 +2 @@
    !!======================================================================
    !!                     ***  MODULE  dynldf_iso  ***
-   !! Ocean dynamics:  lateral viscosity trend
+   !! Ocean dynamics:   lateral viscosity trend (rotated laplacian operator)
    !!======================================================================
    !! History :  OPA  !  97-07  (G. Madec)  Original code
@@ -8 +8 @@
    !!             -   !  2004-08  (C. Talandier) New trends organization
    !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion
+   !!            3.7  !  2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
+   !!                 !                                   add velocity dependent coefficient and optional read in file
    !!----------------------------------------------------------------------
-#if defined key_ldfslp   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_ldfslp'                      slopes of the direction of mixing
+
    !!----------------------------------------------------------------------
    !!   dyn_ldf_iso  : update the momentum trend with the horizontal part
@@ -19 +19 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE ldfdyn_oce      ! ocean dynamics lateral physics
-   USE ldftra_oce      ! ocean tracer   lateral physics
+   USE ldfdyn          ! lateral diffusion: eddy viscosity coef.
+   USE ldftra          ! lateral physics: eddy diffusivity
    USE zdf_oce         ! ocean vertical physics
    USE ldfslp          ! iso-neutral slopes 
    !
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
    USE wrk_nemo        ! Memory Allocation
@@ -42 +42 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "ldfdyn_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -83 +82 @@
       !!      horizontal fluxes associated with the rotated lateral mixing:
       !!      u-component:
-      !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ]
-      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
-      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
-      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
+      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ ub ]
+      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
+      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
+      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ]
       !!      v-component:
-      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
-      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
-      !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
-      !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
+      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ]
+      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ]
+      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
+      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
       !!      take the horizontal divergence of the fluxes:
       !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
@@ -106 +105 @@
       !!      of the rotated operator in dynzdf module
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars
-      REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      -
+      REAL(wp) ::   zmskt, zmskf                                     !   -      -
       REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
       REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
@@ -130 +128 @@
       ENDIF
 
-      ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
+!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
+      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
       IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
          !
@@ -136 +135 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
-                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
-                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
-                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
+                  uslp (ji,jj,jk) = - ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
+                  vslp (ji,jj,jk) = - ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
+                  wslpi(ji,jj,jk) = - ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
+                  wslpj(ji,jj,jk) = - ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
                END DO
             END DO
@@ -184 +183 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
-
-                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
-                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
-
-                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
+                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) * r1_e1t(ji,jj)
+
+                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
+                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
+
+                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
    
+                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    &
+                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
+                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         ELSE                   ! other coordinate system (zco or sco) : e3t
+            DO jj = 2, jpjm1
+               DO ji = fs_2, jpi   ! vector opt.
+                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * fse3t(ji,jj,jk) * r1_e1t(ji,jj)
+
+                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
+                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
+
+                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
+
                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
                      &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
@@ -196 +210 @@
                END DO
             END DO
-         ELSE                   ! other coordinate system (zco or sco) : e3t
-            DO jj = 2, jpjm1
-               DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
-
-                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
-                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
-
-                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
-
-                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
-                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
-                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
-               END DO
-            END DO
          ENDIF
 
@@ -216 +215 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
-
-               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
-                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
-
-               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
+               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * fse3f(ji,jj,jk) * r1_e2f(ji,jj)
+
+               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
+                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
+
+               zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
 
                zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
@@ -237 +236 @@
          DO jj = 2, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
-
-               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
-                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
-
-               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
-
-               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
-                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
-                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
+               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * fse3f(ji,jj,jk) * r1_e1f(ji,jj)
+
+               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
+                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
+
+               zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
+
+               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    &
+                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
+                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
             END DO
          END DO
@@ -254 +253 @@
             DO jj = 2, jpj
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
+                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) * r1_e2t(ji,jj)
+
+                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
+                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
+
+                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
+
+                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    &
+                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
+                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         ELSE                   ! other coordinate system (zco or sco) : e3t
+            DO jj = 2, jpj
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * fse3t(ji,jj,jk) * r1_e2t(ji,jj)
 
                   zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
                      &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
 
-                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
+                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
 
                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
@@ -266 +280 @@
                END DO
             END DO
-         ELSE                   ! other coordinate system (zco or sco) : e3t
-            DO jj = 2, jpj
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
-
-                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
-                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
-
-                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
-
-                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
-                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
-                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
-               END DO
-            END DO
          ENDIF
 
@@ -286 +285 @@
          ! Second derivative (divergence) and add to the general trend
          ! -----------------------------------------------------------
-
          DO jj = 2, jpjm1
-            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
-               ! volume elements
-               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-               ! horizontal component of isopycnal momentum diffusive trends
-               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
-                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
-               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
-                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
-               ! add the trends to the general trends
-               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
-               va (ji,jj,jk) = va (ji,jj,jk) + zvah
+            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj  )    &
+                  &                          + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj  ) - zivf(ji-1,jj)    &
+                  &                          + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
             END DO
          END DO
@@ -358 +349 @@
          DO jk = 2, jpkm1
             DO ji = 2, jpim1
-               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
-
+               zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk)
+               !
                zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
                zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
-
+               !
                zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
                              + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
-               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
-                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
+               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)   &
+                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ), 1. )
 
                zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
@@ -376 +367 @@
                                        +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
                ! update avmu (add isopycnal vertical coefficient to avmu)
-               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
+               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0
             END DO
          END DO
@@ -384 +375 @@
          DO jk = 2, jpkm1
             DO ji = 2, jpim1
-               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
+               zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk)
 
                zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
@@ -398 +389 @@
                ! vertical flux on v field
                zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
-                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
-                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
-                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
+                  &                    +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
+                  &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
+                  &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
                ! update avmv (add isopycnal vertical coefficient to avmv)
-               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
+               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0
             END DO
          END DO
@@ -412 +403 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               ! volume elements
-               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-               ! part of the k-component of isopycnal momentum diffusive trends
-               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
-               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
-               ! add the trends to the general trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
-               va(ji,jj,jk) = va(ji,jj,jk) + zvav
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
             END DO
          END DO
@@ -432 +416 @@
    END SUBROUTINE dyn_ldf_iso
 
-# else
-   !!----------------------------------------------------------------------
-   !!   Dummy module                           NO rotation of mixing tensor
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
-   END SUBROUTINE dyn_ldf_iso
-#endif
-
    !!======================================================================
 END MODULE dynldf_iso

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

@@ -77 +77 @@
       !!             Note that as all external forcing a time averaging over a two rdt
       !!             period is used to prevent the divergence of odd and even time step.
-      !!
-      !! N.B. : When key_esopa is used all the scheme are tested, regardless 
-      !!        of the physical meaning of the results. 
       !!----------------------------------------------------------------------
       !
@@ -121 +118 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
-                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
+                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
                   spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
-                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
+                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
                END DO
             END DO
@@ -135 +132 @@
             DO jj = 2, jpjm1                         ! add tide potential forcing
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
+                  spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
+                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
                END DO 
             END DO
@@ -149 +146 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
-                  spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
+                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+                  spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
                END DO
             END DO
@@ -176 +173 @@
       CASE (  2 )   ;   CALL dyn_spg_flt( kt, kindic )      ! filtered
       !                                                    
-      CASE ( -1 )                                ! esopa: test all possibility with control print
-                        CALL dyn_spg_exp( kt )
-                        CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, &
-         &                            tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-                        CALL dyn_spg_ts ( kt )
-                        CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, &
-         &                           tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-                        CALL dyn_spg_flt( kt, kindic )
-                        CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, &
-         &                            tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       END SELECT
       !                    
@@ -248 +235 @@
       IF(lk_dynspg_flt)   ioptio = ioptio + 1
       !
-      IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   &
+      IF(  ioptio > 1  .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   &
            &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' )
       IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav )   &
            &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' )
       !
-      IF( lk_esopa     )   nspg = -1
       IF( lk_dynspg_exp)   nspg =  0
       IF( lk_dynspg_ts )   nspg =  1
       IF( lk_dynspg_flt)   nspg =  2
       !
-      IF( lk_esopa     )   nspg = -1
-      !
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nspg == -1 )   WRITE(numout,*) '     ESOPA test All scheme used'
          IF( nspg ==  0 )   WRITE(numout,*) '     explicit free surface'
          IF( nspg ==  1 )   WRITE(numout,*) '     free surface with time splitting scheme'
@@ -268 +251 @@
       ENDIF
 
-#if defined key_dynspg_flt || defined key_esopa
+#if defined key_dynspg_flt
       CALL solver_init( nit000 )   ! Elliptic solver initialisation
 #endif
-
-      !                        ! Control of timestep choice
-      IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN
-         IF( nn_cla == 1 )   CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )
-      ENDIF
-
       !               ! Control of hydrostatic pressure choice
-      IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN
-         CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )
-      ENDIF
+      IF( lk_dynspg_ts .AND. ln_dynhpg_imp )   CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init')

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

@@ -7 +7 @@
    !!            3.2  !  2009-06  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
    !!----------------------------------------------------------------------
-#if defined key_dynspg_exp   ||   defined key_esopa
+#if defined key_dynspg_exp
    !!----------------------------------------------------------------------
    !!   'key_dynspg_exp'                              explicit free surface

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

@@ -14 +14 @@
    !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
    !!            3.7  !  2014-04  (F. Roquet, G. Madec)  add some trends diag
-   !!----------------------------------------------------------------------
-#if defined key_dynspg_flt   ||   defined key_esopa  
+   !!             -   !  2014-12  (G. Madec) remove cross-land advection (cla)
+   !!----------------------------------------------------------------------
+#if defined key_dynspg_flt  
    !!----------------------------------------------------------------------
    !!   'key_dynspg_flt'                              filtered free surface
@@ -36 +37 @@
    USE bdydyn          ! ocean open boundary condition on dynamics
    USE bdyvol          ! ocean open boundary condition (bdy_vol routine)
-   USE cla             ! cross land advection
    USE trd_oce         ! trends: ocean variables
    USE trddyn          ! trend manager: dynamics
@@ -206 +206 @@
       CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces 
 #endif
-      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
 
       ! compute the next vertically averaged velocity (effect of the additional force not included)

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90

@@ -17 +17 @@
    
    !                                                       !!! Surface pressure gradient logicals
-#if   defined key_dynspg_exp  ||  defined key_esopa
+#if   defined key_dynspg_exp
    LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .TRUE.   !: Explicit free surface flag
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .FALSE.  !: Explicit free surface flag
 #endif
-#if   defined key_dynspg_ts   ||  defined key_esopa
+#if   defined key_dynspg_ts
    LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .TRUE.   !: Free surface with time splitting flag
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .FALSE.  !: Free surface with time splitting flag
 #endif
-#if   defined key_dynspg_flt  ||  defined key_esopa
+#if   defined key_dynspg_flt
    LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .TRUE.   !: Filtered free surface cst volume flag
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .FALSE.  !: Filtered free surface cst volume flag
 #endif
-
   !                                                                         !!! Time splitting scheme (key_dynspg_ts) 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshn_e, ssha_e   ! sea surface heigth (now, after)

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -12 +12 @@
    !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
    !!---------------------------------------------------------------------
-#if defined key_dynspg_ts   ||   defined key_esopa
+#if defined key_dynspg_ts
    !!----------------------------------------------------------------------
    !!   'key_dynspg_ts'         split explicit free surface
@@ -98 +98 @@
       ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
 
-      IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
-                                                    &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
+      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
+         &                          ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
 
       dyn_spg_ts_alloc = MAXVAL(ierr(:))
@@ -107 +107 @@
       !
    END FUNCTION dyn_spg_ts_alloc
+
 
    SUBROUTINE dyn_spg_ts( kt )
@@ -219 +220 @@
       !
       IF ( kt == nit000 .OR. lk_vvl ) THEN
-         IF ( ln_dynvor_een_old ) THEN
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
-                        &          ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp  
-                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
-               END DO
-            END DO
+         IF ( ln_dynvor_een ) THEN              !==  EEN scheme  ==!
+            SELECT CASE( nn_een_e3f )              !* ff/e3 at F-point
+            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
+                        &             ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp  
+                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj)
+                  END DO
+               END DO
+            CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     &
+                        &             ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   &
+                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
+                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
+                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj)
+                  END DO
+               END DO
+            END SELECT
             CALL lbc_lnk( zwz, 'F', 1._wp )
-            zwz(:,:) = ff(:,:) * zwz(:,:)
-
+            !
             ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
             DO jj = 2, jpj
-               DO ji = fs_2, jpi   ! vector opt.
+               DO ji = 2, jpi
                   ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
                   ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
@@ -239 +252 @@
                END DO
             END DO
-         ELSE IF ( ln_dynvor_een ) THEN
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     &
-                        &          ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   &
-                        &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
-                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
-                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
-               END DO
-            END DO
-            CALL lbc_lnk( zwz, 'F', 1._wp )
-            zwz(:,:) = ff(:,:) * zwz(:,:)
-
-            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-            DO jj = 2, jpj
-               DO ji = fs_2, jpi   ! vector opt.
-                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
-                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
-                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
-                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
-               END DO
-            END DO
-         ELSE
+            !
+         ELSE                                !== all other schemes (ENE, ENS, MIX)
             zwz(:,:) = 0._wp
-            zhf(:,:) = 0.
+            zhf(:,:) = 0._wp
             IF ( .not. ln_sco ) THEN
+
+!!gm  agree the JC comment  : this should be done in a much clear way
+
 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
 !     Set it to zero for the time being 
@@ -276 +271 @@
 
             DO jj = 1, jpjm1
-               zhf(:,jj) = zhf(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
             END DO
 
@@ -297 +292 @@
       ! If forward start at previous time step, and centered integration, 
       ! then update averaging weights:
-      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
+      IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
          ll_fw_start=.FALSE.
          CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
@@ -338 +333 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
-               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
-               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
+               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
+               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
+               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
                ! energy conserving formulation for planetary vorticity term
                zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
@@ -352 +347 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
-                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
+                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
                zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
-                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
                zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
                zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
@@ -360 +355 @@
          END DO
          !
-      ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN  ! enstrophy and energy conserving scheme
+      ELSEIF ( ln_dynvor_een ) THEN  ! enstrophy and energy conserving scheme
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
-                &                                      + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
-                &                                      + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
-                &                                      + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
-                &                                      + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
-                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
-                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
+               zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
+                &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
+                &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
+                &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
+               zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
+                &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
+                &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
+                &                                         + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
             END DO
          END DO
@@ -381 +376 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
-               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)
+               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
+               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)
             END DO
          END DO
@@ -431 +426 @@
             DO jj = 2, jpjm1              
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
-                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
+                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
+                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
                   zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
                   zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
@@ -441 +436 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
-                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
+                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
                   zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
-                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
+                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
                   zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
                   zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
@@ -549 +544 @@
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
                DO ji = 2, fs_jpim1   ! Vector opt.
-                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)     &
-                     &              * ( e12t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
-                     &              +   e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
-                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)     &
-                     &              * ( e12t(ji,jj  ) * zsshp2_e(ji,jj  )  &
-                     &              +   e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)     &
+                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
+                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
+                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)     &
+                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
+                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
                END DO
             END DO
@@ -602 +597 @@
          ! Sum over sub-time-steps to compute advective velocities
          za2 = wgtbtp2(jn)
-         zu_sum  (:,:) = zu_sum  (:,:) + za2 * zwx  (:,:) / e2u  (:,:)
-         zv_sum  (:,:) = zv_sum  (:,:) + za2 * zwy  (:,:) / e1v  (:,:)
+         zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
+         zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
          !
          ! Set next sea level:
@@ -609 +604 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
-                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e12t(ji,jj)
+                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
             END DO
          END DO
@@ -627 +622 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1      ! NO Vector Opt.
-                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)  &
-                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
-                     &              +   e12t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
-                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)  &
-                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
-                     &              +   e12t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
+                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)  &
+                     &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
+                     &              +   e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
+                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)  &
+                     &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
+                     &              +   e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
                END DO
             END DO
@@ -666 +661 @@
             DO jj = 2, jpjm1                            
                DO ji = 2, jpim1
-                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    &
-                     &      * ( e12t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
-                     &      +   e12t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
-                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e12v(ji  ,jj  )  &
-                     &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
-                     &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
+                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e1e2u(ji  ,jj)    &
+                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
+                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
+                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e1e2v(ji  ,jj  )  &
+                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
+                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
                   zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
                   zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
@@ -688 +683 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
-                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
-                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
+                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
+                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
+                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
                   zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
                   zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
@@ -701 +696 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
-                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
+                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
                   zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
-                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
                   zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
                   zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
@@ -709 +704 @@
             END DO
             !
-         ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==!
+         ELSEIF ( ln_dynvor_een ) THEN !==  energy and enstrophy conserving scheme  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
-                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
-                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
-                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
-                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
-                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
-                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
+                  zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
+                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
+                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
+                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
+                  zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
+                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
+                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
+                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
                END DO
             END DO
@@ -729 +724 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
+                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
+                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
                   zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
                   zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
@@ -745 +740 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ! Add surface pressure gradient
-               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
-               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
+               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
+               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
                zwx(ji,jj) = zu_spg
                zwy(ji,jj) = zv_spg
@@ -850 +845 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1      ! NO Vector Opt.
-               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj) &
-                  &              * ( e12t(ji  ,jj) * ssha(ji  ,jj)    &
-                  &              +   e12t(ji+1,jj) * ssha(ji+1,jj) )
-               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj) &
-                  &              * ( e12t(ji,jj  ) * ssha(ji,jj  )    &
-                  &              +   e12t(ji,jj+1) * ssha(ji,jj+1) )
+               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
+                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    &
+                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
+               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
+                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    &
+                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
             END DO
          END DO
@@ -1093 +1088 @@
          DO jj = 1, jpj
             DO ji =1, jpi
-               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
-               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
-               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
+               zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
+               zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
+               zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) )
             END DO
          END DO
@@ -1101 +1096 @@
          DO jj = 1, jpj
             DO ji =1, jpi
-               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
-               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
-               zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) )
-            END DO
-         END DO
-      ENDIF
-
-      zcmax = MAXVAL(zcu(:,:))
+               zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
+               zyr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
+               zcu(ji,jj) = SQRT( grav * ht(ji,jj) * (zxr2 + zyr2) )
+            END DO
+         END DO
+      ENDIF
+
+      zcmax = MAXVAL( zcu(:,:) )
       IF( lk_mpp )   CALL mpp_max( zcmax )
 
@@ -1114 +1109 @@
       IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
       
-      rdtbt = rdt / FLOAT(nn_baro)
+      rdtbt = rdt / REAL( nn_baro , wp )
       zcmax = zcmax * rdtbt
                      ! Print results
@@ -1194 +1189 @@
    !!======================================================================
 END MODULE dynspg_ts
-
-
-

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

@@ -16 +16 @@
    !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
    !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 
+   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory
    !!----------------------------------------------------------------------
 
@@ -22 +23 @@
    !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T)
    !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T)
-   !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
    !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T)
    !!   dyn_vor_init : set and control of the different vorticity option
@@ -32 +32 @@
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
+   !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
@@ -44 +45 @@
 
    PUBLIC   dyn_vor        ! routine called by step.F90
-   PUBLIC   dyn_vor_init   ! routine called by opa.F90
+   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
 
    !                                   !!* Namelist namdyn_vor: vorticity term
-   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation)
-
-   INTEGER ::   nvor = 0   ! type of vorticity trend used
-   INTEGER ::   ncor = 1   ! coriolis
-   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
-   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
-
+   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme    (ENE)
+   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme (ENS)
+   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                (MIX)
+   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme (EEN)
+   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
+   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
+
+   INTEGER ::   nvor_scheme        ! choice of the type of advection scheme
+   !                               ! associated indices:
+   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 2   ! ENS scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 3   ! MIX scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
+
+   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
+   !                               ! associated indices:
+   INTEGER, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
+   INTEGER, PARAMETER ::   np_RVO = 2         ! relative vorticity
+   INTEGER, PARAMETER ::   np_MET = 3         ! metric term
+   INTEGER, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
+   INTEGER, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
+   
+   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
+   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
+   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
+   
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -87 +103 @@
       IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
       !
-      !                                          ! vorticity term 
-      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
-      !
-      CASE ( -1 )                                      ! esopa: test all possibility with control print
-         CALL vor_ene( kt, ntot, ua, va )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_ens( kt, ntot, ua, va )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_mix( kt )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_een( kt, ntot, ua, va )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         !
-      CASE ( 0 )                                       ! energy conserving scheme
-         IF( l_trddyn )   THEN
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+      SELECT CASE ( nvor_scheme )               !==  vorticity trend added to the general trend  ==!
+      !
+      CASE ( np_ENE )                                 !* energy conserving scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -114 +115 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
-         ENDIF
-         !
-      CASE ( 1 )                                       ! enstrophy conserving scheme
-         IF( l_trddyn )   THEN    
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity trend
+         ENDIF
+         !
+      CASE ( np_ENS )                                 !* enstrophy conserving scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: splitthe trend in two    
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -132 +133 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
-         ENDIF
-         !
-      CASE ( 2 )                                       ! mixed ene-ens scheme
-         IF( l_trddyn )   THEN
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
+            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity trend
+         ENDIF
+         !
+      CASE ( np_MIX )                                 !* mixed ene-ens scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -150 +151 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
+            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
-         ENDIF
-         !
-      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
-         IF( l_trddyn )   THEN
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
+            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
+        ENDIF
+         !
+      CASE ( np_EEN )                                 !* energy and enstrophy conserving scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -168 +170 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
+            CALL vor_een( kt, ntot, ua, va )                ! total vorticity trend
          ENDIF
          !
@@ -197 +199 @@
       !!
       !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Sadourny (1975) flux form formulation : conserves the
-      !!      horizontal kinetic energy.
-      !!      The trend of the vorticity term is given by:
-      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
-      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
-      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
-      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
-      !!      Add this trend to the general momentum trend (ua,va):
-      !!          (ua,va) = (ua,va) + ( voru , vorv )
+      !!       and the Sadourny (1975) flux form formulation : conserves the
+      !!       horizontal kinetic energy.
+      !!         The general trend of momentum is increased due to the vorticity 
+      !!       term which is given by:
+      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ]
+      !!       where rvor is the relative vorticity
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
@@ -213 +211 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      !
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
@@ -220 +217 @@
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -234 +231 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
-
-      zfact2 = 0.5 * 0.5      ! Local constant initialization
-
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
+      !
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
          !
-         ! Potential vorticity and horizontal fluxes
-         ! -----------------------------------------
-         SELECT CASE( kvor )      ! vorticity considered
-         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
-         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
-         CASE ( 3 )                                                ! metric term
+         SELECT CASE( kvor )                 !==  vorticity considered  ==!
+         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
+            zwz(:,:) = ff(:,:) 
+         CASE ( np_RVO )                           !* relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_MET )                           !* metric term
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
-               END DO
-            END DO
-         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
-         CASE ( 5 )                                                ! total (coriolis + metric)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
-                       &       )
-               END DO
-            END DO
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CRV )                           !* Coriolis + relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &                   * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CME )                           !* Coriolis + metric
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj)                                                                        &
+                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE DEFAULT                                             ! error
+            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
          END SELECT
+         !
+         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
 
          IF( ln_sco ) THEN
@@ -276 +291 @@
             zwy(:,:) = e1v(:,:) * vn(:,:,jk)
          ENDIF
-
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
+         !                                   !==  compute and add the vorticity term trend  =!
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -285 +298 @@
                zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
                zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
+               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
             END DO  
          END DO  
@@ -299 +312 @@
 
 
-   SUBROUTINE vor_mix( kt )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE vor_mix  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time)
-      !!      Mixte formulation : conserves the potential enstrophy of a hori-
-      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
-      !!      ticity term and the horizontal kinetic energy for (f x uh), the
-      !!      coriolis term. the now trend of the vorticity term is given by:
-      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
-      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
-      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
-      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
-      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
-      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
-      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
-      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
-      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
-      !!      Add this now trend to the general momentum trend (ua,va):
-      !!          (ua,va) = (ua,va) + ( voru , vorv )
-      !!
-      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
-      !!
-      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
-      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('vor_mix')
-      !
-      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww ) 
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-
-      zfact1 = 0.5 * 0.25      ! Local constant initialization
-      zfact2 = 0.5 * 0.5
-
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         ! Relative and planetary potential vorticity and horizontal fluxes
-         ! ----------------------------------------------------------------
-         IF( ln_sco ) THEN        
-            IF( ln_dynadv_vec ) THEN
-               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
-            ELSE                       
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
-                  END DO
-               END DO
-            ENDIF
-            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         ELSE
-            IF( ln_dynadv_vec ) THEN
-               zww(:,:) = rotn(:,:,jk)
-            ELSE                       
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
-                  END DO
-               END DO
-            ENDIF
-            zwz(:,:) = ff (:,:)
-            zwx(:,:) = e2u(:,:) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
-         ENDIF
-
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
-               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
-               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-               ! enstrophy conserving formulation for relative vorticity term
-               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
-               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
-               ! energy conserving formulation for planetary vorticity term
-               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
-               ! mixed vorticity trend added to the momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
-            END DO  
-         END DO  
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('vor_mix')
-      !
-   END SUBROUTINE vor_mix
-
-
    SUBROUTINE vor_ens( kt, kvor, pua, pva )
       !!----------------------------------------------------------------------
@@ -429 +323 @@
       !!      potential enstrophy of a horizontally non-divergent flow. the
       !!      trend of the vorticity term is given by:
-      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
-      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
-      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
-      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
-      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
+      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
       !!      Add this trend to the general momentum trend (ua,va):
       !!          (ua,va) = (ua,va) + ( voru , vorv )
@@ -442 +332 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      !
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
@@ -449 +338 @@
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
       !
-      INTEGER  ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zuav, zvau   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -463 +352 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
-
-      zfact1 = 0.5 * 0.25      ! Local constant initialization
-
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
          !
-         ! Potential vorticity and horizontal fluxes
-         ! -----------------------------------------
-         SELECT CASE( kvor )      ! vorticity considered
-         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
-         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
-         CASE ( 3 )                                                ! metric term
+         SELECT CASE( kvor )                 !==  vorticity considered  ==!
+         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
+            zwz(:,:) = ff(:,:) 
+         CASE ( np_RVO )                           !* relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_MET )                           !* metric term
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
-               END DO
-            END DO
-         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
-         CASE ( 5 )                                                ! total (coriolis + metric)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
-                       &       )
-               END DO
-            END DO
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CRV )                           !* Coriolis + relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &                   * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CME )                           !* Coriolis + metric
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj)                                                                       &
+                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE DEFAULT                                             ! error
+            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
          END SELECT
          !
-         IF( ln_sco ) THEN
-            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
-               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
-                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
-                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
-                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
-               END DO
-            END DO
+         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
+         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
+            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
+            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
+            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
          ELSE
-            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
-               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
-                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
-                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
+            zwx(:,:) = e2u(:,:) * un(:,:,jk)
+            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
+         ENDIF
+         !                                   !==  compute and add the vorticity term trend  =!
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
-                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
-               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
-                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
+               zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
+                  &                       + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
+               zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
+                  &                       + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
                pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
                pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
@@ -553 +449 @@
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
       !!----------------------------------------------------------------------
-      !
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-      !                                                           ! =nrvm (relative vorticity or metric)
+      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric)
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
-      !!
-      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
-      INTEGER  ::   ierr                                          ! local integer
-      REAL(wp) ::   zfac12, zua, zva                              ! local scalars
-      REAL(wp) ::   zmsk, ze3                                     ! local scalars
-      !                                                           !  3D workspace 
-      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
-      REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse
-#if defined key_vvl
-      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
-#else
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
-#endif
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ierr         ! local integer
+      REAL(wp) ::   zua, zva     ! local scalars
+      REAL(wp) ::   zmsk, ze3    ! local scalars
+      !
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwx, zwy, zwz, z1_e3f
+      REAL(wp), POINTER, DIMENSION(:,:)   :: ztnw, ztne, ztsw, ztse
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('vor_een')
       !
-      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
-      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
-#if defined key_vvl
-      CALL wrk_alloc( jpi, jpj, jpk, ze3f                   )
-#endif
+      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
+      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
       !
       IF( kt == nit000 ) THEN
@@ -586 +472 @@
          IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-#if ! defined key_vvl
-         IF( .NOT.ALLOCATED(ze3f) ) THEN
-            ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
-            IF( lk_mpp    )   CALL mpp_sum ( ierr )
-            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
-         ENDIF
-         ze3f(:,:,:) = 0._wp
-#endif
       ENDIF
-
-      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
-
-         IF( ln_dynvor_een_old ) THEN ! original formulation
-            DO jk = 1, jpk
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
-                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3
-                  END DO
-               END DO
-            END DO
-         ELSE ! new formulation from NEMO 3.6
-            DO jk = 1, jpk
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
-                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
-                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
-                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
-                  END DO
-               END DO
-            END DO
-         ENDIF
-
-         CALL lbc_lnk( ze3f, 'F', 1. )
-      ENDIF
-
-      zfac12 = 1._wp / 12._wp    ! Local constant initialization
-
-      
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
+      !
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-         
-         ! Potential vorticity and horizontal fluxes
-         ! -----------------------------------------
-         SELECT CASE( kvor )      ! vorticity considered
-         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
-            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
-         CASE ( 2 )                                                ! relative  vorticity
-            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
-         CASE ( 3 )                                                ! metric term
+         !
+         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
+         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
+                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4.0_wp / ze3
+                  ELSE                      ;   z1_e3f(ji,jj) = 0.0_wp
+                  ENDIF
+               END DO
+            END DO
+         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
+                  zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
+                     &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
+                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3
+                  ELSE                      ;   z1_e3f(ji,jj) = 0.0_wp
+                  ENDIF
+               END DO
+            END DO
+         END SELECT
+         !
+         SELECT CASE( kvor )                 !==  vorticity considered  ==!
+         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_RVO )                           !* relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_MET )                           !* metric term
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
-               END DO
-            END DO
-            CALL lbc_lnk( zwz, 'F', 1. )
-        CASE ( 4 )                                                ! total (relative + planetary vorticity)
-            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
-         CASE ( 5 )                                                ! total (coriolis + metric)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
-                       &       ) * ze3f(ji,jj,jk)
-               END DO
-            END DO
-            CALL lbc_lnk( zwz, 'F', 1. )
+                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CRV )                           !* Coriolis + relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &                         - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CME )                           !* Coriolis + metric
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  ff(ji,jj)                                                                        &
+                       &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE DEFAULT                                             ! error
+            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
          END SELECT
-
+         !
+         CALL lbc_lnk( zwz, 'F', 1. )
+         !
+         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
+         !                                   !==  horizontal fluxes  ==!
          zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
          zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
 
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
+         !                                   !==  compute and add the vorticity term trend  =!
          jj = 2
          ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
-         DO ji = 2, jpi   
+         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
                ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
                ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
@@ -687 +580 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                  &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
+               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
+                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
+               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
+                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
                pua(ji,jj,jk) = pua(ji,jj,jk) + zua
                pva(ji,jj,jk) = pva(ji,jj,jk) + zva
@@ -698 +591 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
-      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
-#if defined key_vvl
-      CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   )
-#endif
+      !
+      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
+      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('vor_een')
@@ -720 +611 @@
       INTEGER ::   ios             ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
+      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk
       !!----------------------------------------------------------------------
 
@@ -737 +628 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
-         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
-         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
-         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
-         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
-         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
+         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
+         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
+         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
+         WRITE(numout,*) '           enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
+         WRITE(numout,*) '              e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
+         WRITE(numout,*) '           masked (=1) or unmasked(=0) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
       ENDIF
 
+!!gm  this should be removed when choosing a unique strategy for fmask at the coast
       ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
       ! at angles with three ocean points and one land point
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat
       IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
          DO jk = 1, jpk
@@ -759 +654 @@
           !
       ENDIF
-
-      ioptio = 0                     ! Control of vorticity scheme options
-      IF( ln_dynvor_ene )   ioptio = ioptio + 1
-      IF( ln_dynvor_ens )   ioptio = ioptio + 1
-      IF( ln_dynvor_mix )   ioptio = ioptio + 1
-      IF( ln_dynvor_een )   ioptio = ioptio + 1
-      IF( ln_dynvor_een_old )   ioptio = ioptio + 1
-      IF( lk_esopa      )   ioptio =          1
-
+!!gm end
+
+      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
+      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENE   ;   ENDIF
+      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENS   ;   ENDIF
+      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_MIX   ;   ENDIF
+      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF
+      !
       IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
-
-      !                              ! Set nvor (type of scheme for vorticity)
-      IF( ln_dynvor_ene )   nvor =  0
-      IF( ln_dynvor_ens )   nvor =  1
-      IF( ln_dynvor_mix )   nvor =  2
-      IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3
-      IF( lk_esopa      )   nvor = -1
-      
-      !                              ! Set ncor, nrvm, ntot (type of vorticity)
-      IF(lwp) WRITE(numout,*)
-      ncor = 1
+      !                      
+      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
+      ncor = np_COR
       IF( ln_dynadv_vec ) THEN     
          IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
-         nrvm = 2
-         ntot = 4
+         nrvm = np_RVO        ! relative vorticity
+         ntot = np_CRV        ! relative + planetary vorticity
       ELSE                        
          IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
-         nrvm = 3
-         ntot = 5
+         nrvm = np_MET        ! metric term
+         ntot = np_CME        ! Coriolis + metric term
       ENDIF
       
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
-         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
-         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
-         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
-         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
-         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
+         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme'
+         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme'
+         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme'
+         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme'
       ENDIF
       !

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

@@ -49 +49 @@
       !!
       !! ** Method  :   The now vertical advection of momentum is given by:
-      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
-      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
+      !!         w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ]
+      !!         w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ]
       !!      Add this trend to the general trend (ua,va):
       !!         (ua,va) = (ua,va) + w dz(u,v)
@@ -85 +85 @@
          DO jj = 2, jpj                   ! vertical fluxes 
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
@@ -121 +121 @@
             DO ji = fs_2, fs_jpim1       ! vector opt.
                !                         ! vertical momentum advective trends
-               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
-               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
+               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
                !                         ! add the trends to the general momentum trends
                ua(ji,jj,jk) = ua(ji,jj,jk) + zua
@@ -146 +146 @@
       !
    END SUBROUTINE dyn_zad
+
 
    SUBROUTINE dyn_zad_zts ( kt )
@@ -182 +183 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
       !
-      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
-      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 
+      CALL wrk_alloc( jpi,jpj,jpk,     zwuw, zwvw, zww ) 
+      CALL wrk_alloc( jpi,jpj,jpk,3,   zus , zvs ) 
       !
       IF( kt == nit000 ) THEN
@@ -205 +206 @@
          DO jj = 2, jpj                   
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
       END DO
-      !
-      ! Surface and bottom advective fluxes set to zero
-      DO jj = 2, jpjm1        
+
+      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
          DO ji = fs_2, fs_jpim1           ! vector opt.
+ !!gm missing ISF boundary condition
             zwuw(ji,jj, 1 ) = 0._wp
             zwvw(ji,jj, 1 ) = 0._wp
@@ -251 +252 @@
                DO ji = fs_2, fs_jpim1       ! vector opt.
                   !                         ! vertical momentum advective trends
-                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
-                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
+                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
                   zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
                   zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
@@ -283 +284 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
-      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,     zwuw, zwvw, zww ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,3,   zus , zvs ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

@@ -9 +9 @@
 
    !!----------------------------------------------------------------------
-   !!   dyn_zdf      : Update the momentum trend with the vertical diffusion
-   !!   dyn_zdf_init : initializations of the vertical diffusion scheme
+   !!   dyn_zdf       : Update the momentum trend with the vertical diffusion
+   !!   dyn_zdf_init  : initializations of the vertical diffusion scheme
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE zdf_oce         ! ocean vertical physics variables
-
-   USE dynzdf_exp      ! vertical diffusion: explicit (dyn_zdf_exp     routine)
-   USE dynzdf_imp      ! vertical diffusion: implicit (dyn_zdf_imp     routine)
-
-   USE ldfdyn_oce      ! ocean dynamics: lateral physics
-   USE trd_oce         ! trends: ocean variables
-   USE trddyn          ! trend manager: dynamics
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! MPP library
-   USE prtctl          ! Print control
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables 
+   USE zdf_oce        ! ocean vertical physics variables
+   USE dynzdf_exp     ! vertical diffusion: explicit (dyn_zdf_exp     routine)
+   USE dynzdf_imp     ! vertical diffusion: implicit (dyn_zdf_imp     routine)
+   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE prtctl         ! Print control
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -61 +60 @@
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf')
+      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf')
       !
       !                                          ! set time step
@@ -79 +78 @@
       CASE ( 1 )   ;   CALL dyn_zdf_imp( kt, r2dt )      ! implicit scheme
       !
-      CASE ( -1 )                                        ! esopa: test all possibility with control print
-                       CALL dyn_zdf_exp( kt, r2dt )
-                       CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask,               &
-                          &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-                       CALL dyn_zdf_imp( kt, r2dt )
-                       CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask,               &
-                          &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       END SELECT
 
@@ -96 +88 @@
       !                                          ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
-            &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf')
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         !
+      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf')
       !
    END SUBROUTINE dyn_zdf
@@ -114 +106 @@
       USE zdftke
       USE zdfgls
-      USE zdfkpp
       !!----------------------------------------------------------------------
       !
@@ -123 +114 @@
       !
       ! Force implicit schemes
-      IF( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp )   nzdf = 1   ! TKE, GLS or KPP physics
-      IF( ln_dynldf_iso                           )   nzdf = 1   ! iso-neutral lateral physics
-      IF( ln_dynldf_hor .AND. ln_sco              )   nzdf = 1   ! horizontal lateral physics in s-coordinate
-      !
-      IF( lk_esopa )    nzdf = -1                   ! Esopa key: All schemes used
+      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE or GLS physics
+      IF( ln_dynldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
+      IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
       !
       IF(lwp) THEN                                  ! Print the choice
@@ -133 +122 @@
          WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
          IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
          IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

@@ -209 +209 @@
       !-----------------------------------------------------------------------
       !
-      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-      DO jk = 2, jpkm1
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -309 +308 @@
       !-----------------------------------------------------------------------
       !
-      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-      DO jk = 2, jpkm1        
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -20 +20 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE domvvl          ! Variable volume
-   USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE restart         ! only for lrst_oce
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE phycst
-   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-   USE lib_mpp         ! MPP library
+   USE divhor          ! horizontal divergence
+   USE phycst          ! physical constants
    USE bdy_oce
    USE bdy_par         
@@ -36 +31 @@
    USE asminc          ! Assimilation increment
 #endif
+   USE in_out_manager  ! I/O manager
+   USE restart         ! only for lrst_oce
+   USE prtctl          ! Print control
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+   USE lib_mpp         ! MPP library
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
@@ -66 +66 @@
       !!      by the time step.
       !!
-      !! ** action  :   ssha    : after sea surface height
+      !! ** action  :   ssha, after sea surface height
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
-      INTEGER, INTENT(in) ::   kt                      ! time step
+      INTEGER, INTENT(in) ::   kt   ! time step
       ! 
-      INTEGER             ::   jk                      ! dummy loop indice
-      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
-      !
-      CALL wrk_alloc( jpi, jpj, zhdiv ) 
+      INTEGER  ::   jk            ! dummy loop indice
+      REAL(wp) ::   z2dt, zcoef   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt')
+      !
+      CALL wrk_alloc( jpi,jpj,   zhdiv ) 
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         !
-      ENDIF
-      !
-      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
+      ENDIF
+      !
+      CALL div_hor( kt )                              ! Horizontal divergence
       !
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
@@ -106 +103 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      z1_rau0 = 0.5_wp * r1_rau0
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+      zcoef = 0.5_wp * r1_rau0
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
 
 #if ! defined key_dynspg_ts
       ! These lines are not necessary with time splitting since
       ! boundary condition on sea level is set during ts loop
-#if defined key_agrif
+# if defined key_agrif
       CALL agrif_ssh( kt )
-#endif
-#if defined key_bdy
-      IF (lk_bdy) THEN
-         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
-         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
-      ENDIF
-#endif
+# endif
+# if defined key_bdy
+      IF( lk_bdy ) THEN
+         CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
+         CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
+      ENDIF
+# endif
 #endif
 
 #if defined key_asminc
-      !                                                ! Include the IAU weighted SSH increment
-      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN     ! Include the IAU weighted SSH increment
          CALL ssh_asm_inc( kt )
          ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
       ENDIF
 #endif
-
       !                                           !------------------------------!
       !                                           !           outputs            !
@@ -159 +154 @@
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt           ! time step
+      INTEGER, INTENT(in) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   z1_2dt       ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      !
-      INTEGER             ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp)            ::   z1_2dt       ! local scalars
-      !!----------------------------------------------------------------------
-      
-      IF( nn_timing == 1 )  CALL timing_start('wzv')
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('wzv')
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
@@ -177 +170 @@
          !
          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
       ENDIF
       !                                           !------------------------------!
@@ -193 +185 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+                  zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
                END DO
             END DO
@@ -216 +208 @@
 
 #if defined key_bdy
-      IF (lk_bdy) THEN
+      IF( lk_bdy ) THEN
          DO jk = 1, jpkm1
             wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
@@ -224 +216 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
-
-
+      !
    END SUBROUTINE wzv
+
 
    SUBROUTINE ssh_swp( kt )
@@ -265 +257 @@
          sshb(:,:) = sshn(:,:)                           ! before <-- now
          sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
+         !
       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered