New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 – NEMO

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r4624 r6225  
    1515   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme 
    1616   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     17   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity  
     18   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory 
    1719   !!---------------------------------------------------------------------- 
    1820 
    … …  
    2123   !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T) 
    2224   !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T) 
    23    !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 
    2425   !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T) 
    2526   !!   dyn_vor_init : set and control of the different vorticity option 
    … …  
    2930   USE dommsk         ! ocean mask 
    3031   USE dynadv         ! momentum advection (use ln_dynadv_vec value) 
    31    USE trdmod         ! ocean dynamics trends  
    32    USE trdmod_oce     ! ocean variables trends 
     32   USE trd_oce        ! trends: ocean variables 
     33   USE trddyn         ! trend manager: dynamics 
     34   ! 
    3335   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3436   USE prtctl         ! Print control 
    … …  
    4345 
    4446   PUBLIC   dyn_vor        ! routine called by step.F90 
    45    PUBLIC   dyn_vor_init   ! routine called by opa.F90 
     47   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90 
    4648 
    4749   !                                   !!* Namelist namdyn_vor: vorticity term 
    48    LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme 
    49    LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme 
    50    LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme 
    51    LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme 
    52  
    53    INTEGER ::   nvor = 0   ! type of vorticity trend used 
    54    INTEGER ::   ncor = 1   ! coriolis 
    55    INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term 
    56    INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term 
    57  
     50   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme    (ENE) 
     51   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme (ENS) 
     52   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                (MIX) 
     53   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme (EEN) 
     54   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
     55   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
     56 
     57   INTEGER ::   nvor_scheme        ! choice of the type of advection scheme 
     58   !                               ! associated indices: 
     59   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme 
     60   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 2   ! ENS scheme 
     61   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 3   ! MIX scheme 
     62   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme 
     63 
     64   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
     65   !                               ! associated indices: 
     66   INTEGER, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
     67   INTEGER, PARAMETER ::   np_RVO = 2         ! relative vorticity 
     68   INTEGER, PARAMETER ::   np_MET = 3         ! metric term 
     69   INTEGER, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity) 
     70   INTEGER, PARAMETER ::   np_CME = 5         ! Coriolis + metric term 
     71    
     72   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
     73   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8 
     74   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12 
     75    
    5876   !! * Substitutions 
    59 #  include "domzgr_substitute.h90" 
    6077#  include "vectopt_loop_substitute.h90" 
    6178   !!---------------------------------------------------------------------- 
    62    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     79   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    6380   !! $Id$ 
    6481   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    7390      !! ** Action : - Update (ua,va) with the now vorticity term trend 
    7491      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    75       !!               and planetary vorticity trends) ('key_trddyn') 
     92      !!               and planetary vorticity trends) and send them to trd_dyn  
     93      !!               for futher diagnostics (l_trddyn=T) 
    7694      !!---------------------------------------------------------------------- 
    7795      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    … …  
    84102      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    85103      ! 
    86       !                                          ! vorticity term  
    87       SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend 
    88       ! 
    89       CASE ( -1 )                                      ! esopa: test all possibility with control print 
    90          CALL vor_ene( kt, ntot, ua, va ) 
    91          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 
    92             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    93          CALL vor_ens( kt, ntot, ua, va ) 
    94          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 
    95             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    96          CALL vor_mix( kt ) 
    97          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 
    98             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    99          CALL vor_een( kt, ntot, ua, va ) 
    100          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 
    101             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    102          ! 
    103       CASE ( 0 )                                       ! energy conserving scheme 
    104          IF( l_trddyn )   THEN 
    105             ztrdu(:,:,:) = ua(:,:,:) 
    106             ztrdv(:,:,:) = va(:,:,:) 
    107             CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
    108             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    109             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    110             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
    111             ztrdu(:,:,:) = ua(:,:,:) 
    112             ztrdv(:,:,:) = va(:,:,:) 
    113             CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend 
    114             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    115             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    116             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    117             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
     104      SELECT CASE ( nvor_scheme )               !==  vorticity trend added to the general trend  ==! 
     105      ! 
     106      CASE ( np_ENE )                                 !* energy conserving scheme 
     107         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two 
     108            ztrdu(:,:,:) = ua(:,:,:) 
     109            ztrdv(:,:,:) = va(:,:,:) 
     110            CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     111            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     112            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     113            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
     114            ztrdu(:,:,:) = ua(:,:,:) 
     115            ztrdv(:,:,:) = va(:,:,:) 
     116            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend 
     117            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     118            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     119            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    118120         ELSE 
    119             CALL vor_ene( kt, ntot, ua, va )                ! total vorticity 
    120          ENDIF 
    121          ! 
    122       CASE ( 1 )                                       ! enstrophy conserving scheme 
    123          IF( l_trddyn )   THEN     
    124             ztrdu(:,:,:) = ua(:,:,:) 
    125             ztrdv(:,:,:) = va(:,:,:) 
    126             CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
    127             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    128             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    129             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
    130             ztrdu(:,:,:) = ua(:,:,:) 
    131             ztrdv(:,:,:) = va(:,:,:) 
    132             CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend 
    133             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    134             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    135             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    136             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
     121            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity trend 
     122         ENDIF 
     123         ! 
     124      CASE ( np_ENS )                                 !* enstrophy conserving scheme 
     125         IF( l_trddyn ) THEN                                ! trend diagnostics: splitthe trend in two     
     126            ztrdu(:,:,:) = ua(:,:,:) 
     127            ztrdv(:,:,:) = va(:,:,:) 
     128            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     129            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     130            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     131            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
     132            ztrdu(:,:,:) = ua(:,:,:) 
     133            ztrdv(:,:,:) = va(:,:,:) 
     134            CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend 
     135            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     136            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     137            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    137138         ELSE 
    138             CALL vor_ens( kt, ntot, ua, va )                ! total vorticity 
    139          ENDIF 
    140          ! 
    141       CASE ( 2 )                                       ! mixed ene-ens scheme 
    142          IF( l_trddyn )   THEN 
    143             ztrdu(:,:,:) = ua(:,:,:) 
    144             ztrdv(:,:,:) = va(:,:,:) 
     139            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity trend 
     140         ENDIF 
     141         ! 
     142      CASE ( np_MIX )                                 !* mixed ene-ens scheme 
     143         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two 
     144            ztrdu(:,:,:) = ua(:,:,:) 
     145            ztrdv(:,:,:) = va(:,:,:) 
     146            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens) 
     147            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     148            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     149            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
     150            ztrdu(:,:,:) = ua(:,:,:) 
     151            ztrdv(:,:,:) = va(:,:,:) 
     152            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene) 
     153            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     154            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     155            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
     156         ELSE 
    145157            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens) 
    146             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    147             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    148             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
    149             ztrdu(:,:,:) = ua(:,:,:) 
    150             ztrdv(:,:,:) = va(:,:,:) 
    151158            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene) 
    152             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    153             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    154             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    155             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
     159        ENDIF 
     160         ! 
     161      CASE ( np_EEN )                                 !* energy and enstrophy conserving scheme 
     162         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two 
     163            ztrdu(:,:,:) = ua(:,:,:) 
     164            ztrdv(:,:,:) = va(:,:,:) 
     165            CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     166            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     167            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     168            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
     169            ztrdu(:,:,:) = ua(:,:,:) 
     170            ztrdv(:,:,:) = va(:,:,:) 
     171            CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend 
     172            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     173            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     174            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    156175         ELSE 
    157             CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene) 
    158          ENDIF 
    159          ! 
    160       CASE ( 3 )                                       ! energy and enstrophy conserving scheme 
    161          IF( l_trddyn )   THEN 
    162             ztrdu(:,:,:) = ua(:,:,:) 
    163             ztrdv(:,:,:) = va(:,:,:) 
    164             CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
    165             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    166             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    167             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
    168             ztrdu(:,:,:) = ua(:,:,:) 
    169             ztrdv(:,:,:) = va(:,:,:) 
    170             CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend 
    171             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    172             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    173             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    174             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
    175          ELSE 
    176             CALL vor_een( kt, ntot, ua, va )                ! total vorticity 
     176            CALL vor_een( kt, ntot, ua, va )                ! total vorticity trend 
    177177         ENDIF 
    178178         ! 
    … …  
    198198      !! 
    199199      !! ** Method  :   Trend evaluated using now fields (centered in time)  
    200       !!      and the Sadourny (1975) flux form formulation : conserves the 
    201       !!      horizontal kinetic energy. 
    202       !!      The trend of the vorticity term is given by: 
    203       !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    204       !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ] 
    205       !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ] 
    206       !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    207       !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ] 
    208       !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ] 
    209       !!      Add this trend to the general momentum trend (ua,va): 
    210       !!          (ua,va) = (ua,va) + ( voru , vorv ) 
     200      !!       and the Sadourny (1975) flux form formulation : conserves the 
     201      !!       horizontal kinetic energy. 
     202      !!         The general trend of momentum is increased due to the vorticity  
     203      !!       term which is given by: 
     204      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ] 
     205      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ] 
     206      !!       where rvor is the relative vorticity 
    211207      !! 
    212208      !! ** Action : - Update (ua,va) with the now vorticity term trend 
    213       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    214       !!               and planetary vorticity trends) ('key_trddyn') 
    215209      !! 
    216210      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    217211      !!---------------------------------------------------------------------- 
    218       ! 
    219212      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    220213      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    … …  
    223216      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    224217      ! 
    225       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    226       REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars 
    227       REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz 
     218      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     219      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
     220      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz   ! 2D workspace 
    228221      !!---------------------------------------------------------------------- 
    229222      ! 
    … …  
    237230         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    238231      ENDIF 
    239  
    240       zfact2 = 0.5 * 0.5      ! Local constant initialization 
    241  
    242 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
     232      ! 
    243233      !                                                ! =============== 
    244234      DO jk = 1, jpkm1                                 ! Horizontal slab 
    245235         !                                             ! =============== 
    246236         ! 
    247          ! Potential vorticity and horizontal fluxes 
    248          ! ----------------------------------------- 
    249          SELECT CASE( kvor )      ! vorticity considered 
    250          CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis) 
    251          CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity 
    252          CASE ( 3 )                                                ! metric term 
     237         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     238         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     239            zwz(:,:) = ff(:,:)  
     240         CASE ( np_RVO )                           !* relative vorticity 
     241            DO jj = 1, jpjm1 
     242               DO ji = 1, fs_jpim1   ! vector opt. 
     243                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     244                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     245               END DO 
     246            END DO 
     247         CASE ( np_MET )                           !* metric term 
    253248            DO jj = 1, jpjm1 
    254249               DO ji = 1, fs_jpim1   ! vector opt. 
    255250                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    256251                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    257                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 
    258                END DO 
    259             END DO 
    260          CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity) 
    261          CASE ( 5 )                                                ! total (coriolis + metric) 
    262             DO jj = 1, jpjm1 
    263                DO ji = 1, fs_jpim1   ! vector opt. 
    264                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    265                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    266                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    267                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               & 
    268                        &       ) 
    269                END DO 
    270             END DO 
     252                       &     * 0.5 * r1_e1e2f(ji,jj) 
     253               END DO 
     254            END DO 
     255         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     256            DO jj = 1, jpjm1 
     257               DO ji = 1, fs_jpim1   ! vector opt. 
     258                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     259                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     260                     &                   * r1_e1e2f(ji,jj) 
     261               END DO 
     262            END DO 
     263         CASE ( np_CME )                           !* Coriolis + metric 
     264            DO jj = 1, jpjm1 
     265               DO ji = 1, fs_jpim1   ! vector opt. 
     266                  zwz(ji,jj) = ff(ji,jj)                                                                        & 
     267                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     268                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     269                       &     * 0.5 * r1_e1e2f(ji,jj) 
     270               END DO 
     271            END DO 
     272         CASE DEFAULT                                             ! error 
     273            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    271274         END SELECT 
     275         ! 
     276         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
     277            DO jj = 1, jpjm1 
     278               DO ji = 1, fs_jpim1   ! vector opt. 
     279                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     280               END DO 
     281            END DO 
     282         ENDIF 
    272283 
    273284         IF( ln_sco ) THEN 
    274             zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 
    275             zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    276             zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
     285            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
     286            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
     287            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
    277288         ELSE 
    278289            zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    279290            zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    280291         ENDIF 
    281  
    282          ! Compute and add the vorticity term trend 
    283          ! ---------------------------------------- 
     292         !                                   !==  compute and add the vorticity term trend  =! 
    284293         DO jj = 2, jpjm1 
    285294            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    288297               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    289298               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    290                pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    291                pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     299               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     300               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
    292301            END DO   
    293302         END DO   
    … …  
    302311 
    303312 
    304    SUBROUTINE vor_mix( kt ) 
    305       !!---------------------------------------------------------------------- 
    306       !!                 ***  ROUTINE vor_mix  *** 
    307       !! 
    308       !! ** Purpose :   Compute the now total vorticity trend and add it to 
    309       !!      the general trend of the momentum equation. 
    310       !! 
    311       !! ** Method  :   Trend evaluated using now fields (centered in time) 
    312       !!      Mixte formulation : conserves the potential enstrophy of a hori- 
    313       !!      zontally non-divergent flow for (rotzu x uh), the relative vor- 
    314       !!      ticity term and the horizontal kinetic energy for (f x uh), the 
    315       !!      coriolis term. the now trend of the vorticity term is given by: 
    316       !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    317       !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ] 
    318       !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ] 
    319       !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ] 
    320       !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ] 
    321       !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    322       !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ] 
    323       !!              +1/e1u  mj-1[ f          mi(e1v vn) ] 
    324       !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ] 
    325       !!              +1/e2v  mi-1[ f          mj(e2u un) ] 
    326       !!      Add this now trend to the general momentum trend (ua,va): 
    327       !!          (ua,va) = (ua,va) + ( voru , vorv ) 
    328       !! 
    329       !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 
    330       !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    331       !!               and planetary vorticity trends) ('key_trddyn') 
    332       !! 
    333       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    334       !!---------------------------------------------------------------------- 
    335       ! 
    336       INTEGER, INTENT(in) ::   kt   ! ocean timestep index 
    337       ! 
    338       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    339       REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars 
    340       REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      - 
    341       REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww 
    342       !!---------------------------------------------------------------------- 
    343       ! 
    344       IF( nn_timing == 1 )  CALL timing_start('vor_mix') 
    345       ! 
    346       CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )  
    347       ! 
    348       IF( kt == nit000 ) THEN 
    349          IF(lwp) WRITE(numout,*) 
    350          IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme' 
    351          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    352       ENDIF 
    353  
    354       zfact1 = 0.5 * 0.25      ! Local constant initialization 
    355       zfact2 = 0.5 * 0.5 
    356  
    357 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 
    358       !                                                ! =============== 
    359       DO jk = 1, jpkm1                                 ! Horizontal slab 
    360          !                                             ! =============== 
    361          ! 
    362          ! Relative and planetary potential vorticity and horizontal fluxes 
    363          ! ---------------------------------------------------------------- 
    364          IF( ln_sco ) THEN         
    365             IF( ln_dynadv_vec ) THEN 
    366                zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk) 
    367             ELSE                        
    368                DO jj = 1, jpjm1 
    369                   DO ji = 1, fs_jpim1   ! vector opt. 
    370                      zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    371                         &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    372                         &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) ) 
    373                   END DO 
    374                END DO 
    375             ENDIF 
    376             zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk) 
    377             zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    378             zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    379          ELSE 
    380             IF( ln_dynadv_vec ) THEN 
    381                zww(:,:) = rotn(:,:,jk) 
    382             ELSE                        
    383                DO jj = 1, jpjm1 
    384                   DO ji = 1, fs_jpim1   ! vector opt. 
    385                      zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    386                         &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    387                         &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) ) 
    388                   END DO 
    389                END DO 
    390             ENDIF 
    391             zwz(:,:) = ff (:,:) 
    392             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    393             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    394          ENDIF 
    395  
    396          ! Compute and add the vorticity term trend 
    397          ! ---------------------------------------- 
    398          DO jj = 2, jpjm1 
    399             DO ji = fs_2, fs_jpim1   ! vector opt. 
    400                zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 
    401                zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    402                zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 
    403                zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    404                ! enstrophy conserving formulation for relative vorticity term 
    405                zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 ) 
    406                zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 ) 
    407                ! energy conserving formulation for planetary vorticity term 
    408                zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    409                zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    410                ! mixed vorticity trend added to the momentum trends 
    411                ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua 
    412                va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva 
    413             END DO   
    414          END DO   
    415          !                                             ! =============== 
    416       END DO                                           !   End of slab 
    417       !                                                ! =============== 
    418       CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )  
    419       ! 
    420       IF( nn_timing == 1 )  CALL timing_stop('vor_mix') 
    421       ! 
    422    END SUBROUTINE vor_mix 
    423  
    424  
    425313   SUBROUTINE vor_ens( kt, kvor, pua, pva ) 
    426314      !!---------------------------------------------------------------------- 
    … …  
    434322      !!      potential enstrophy of a horizontally non-divergent flow. the 
    435323      !!      trend of the vorticity term is given by: 
    436       !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative: 
    437       !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ] 
    438       !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ] 
    439       !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    440       !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ] 
    441       !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ] 
     324      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ] 
     325      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ] 
    442326      !!      Add this trend to the general momentum trend (ua,va): 
    443327      !!          (ua,va) = (ua,va) + ( voru , vorv ) 
    444328      !! 
    445329      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 
    446       !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative  
    447       !!               and planetary vorticity trends) ('key_trddyn') 
    448330      !! 
    449331      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    450332      !!---------------------------------------------------------------------- 
    451       ! 
    452333      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    453334      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    … …  
    456337      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    457338      ! 
    458       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    459       REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars 
    460       REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww 
     339      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     340      REAL(wp) ::   zuav, zvau   ! local scalars 
     341      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace 
    461342      !!---------------------------------------------------------------------- 
    462343      ! 
    … …  
    470351         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    471352      ENDIF 
    472  
    473       zfact1 = 0.5 * 0.25      ! Local constant initialization 
    474  
    475 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
    476353      !                                                ! =============== 
    477354      DO jk = 1, jpkm1                                 ! Horizontal slab 
    478355         !                                             ! =============== 
    479356         ! 
    480          ! Potential vorticity and horizontal fluxes 
    481          ! ----------------------------------------- 
    482          SELECT CASE( kvor )      ! vorticity considered 
    483          CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis) 
    484          CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity 
    485          CASE ( 3 )                                                ! metric term 
     357         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     358         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     359            zwz(:,:) = ff(:,:)  
     360         CASE ( np_RVO )                           !* relative vorticity 
     361            DO jj = 1, jpjm1 
     362               DO ji = 1, fs_jpim1   ! vector opt. 
     363                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     364                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     365               END DO 
     366            END DO 
     367         CASE ( np_MET )                           !* metric term 
    486368            DO jj = 1, jpjm1 
    487369               DO ji = 1, fs_jpim1   ! vector opt. 
    488370                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    489371                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    490                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 
    491                END DO 
    492             END DO 
    493          CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity) 
    494          CASE ( 5 )                                                ! total (coriolis + metric) 
    495             DO jj = 1, jpjm1 
    496                DO ji = 1, fs_jpim1   ! vector opt. 
    497                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    498                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    499                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    500                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    501                        &       ) 
    502                END DO 
    503             END DO 
     372                       &     * 0.5 * r1_e1e2f(ji,jj) 
     373               END DO 
     374            END DO 
     375         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     376            DO jj = 1, jpjm1 
     377               DO ji = 1, fs_jpim1   ! vector opt. 
     378                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     379                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     380                     &                   * r1_e1e2f(ji,jj) 
     381               END DO 
     382            END DO 
     383         CASE ( np_CME )                           !* Coriolis + metric 
     384            DO jj = 1, jpjm1 
     385               DO ji = 1, fs_jpim1   ! vector opt. 
     386                  zwz(ji,jj) = ff(ji,jj)                                                                       & 
     387                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     388                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     389                       &     * 0.5 * r1_e1e2f(ji,jj) 
     390               END DO 
     391            END DO 
     392         CASE DEFAULT                                             ! error 
     393            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    504394         END SELECT 
    505395         ! 
    506          IF( ln_sco ) THEN 
    507             DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    508                DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    509                   zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 
    510                   zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 
    511                   zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 
    512                END DO 
    513             END DO 
     396         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==! 
     397            DO jj = 1, jpjm1 
     398               DO ji = 1, fs_jpim1   ! vector opt. 
     399                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     400               END DO 
     401            END DO 
     402         ENDIF 
     403         ! 
     404         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
     405            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
     406            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
     407            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
    514408         ELSE 
    515             DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    516                DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    517                   zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 
    518                   zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 
    519                END DO 
    520             END DO 
    521          ENDIF 
    522          ! 
    523          ! Compute and add the vorticity term trend 
    524          ! ---------------------------------------- 
     409            zwx(:,:) = e2u(:,:) * un(:,:,jk) 
     410            zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
     411         ENDIF 
     412         !                                   !==  compute and add the vorticity term trend  =! 
    525413         DO jj = 2, jpjm1 
    526414            DO ji = fs_2, fs_jpim1   ! vector opt. 
    527                zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   & 
    528                   &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) 
    529                zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
    530                   &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) 
     415               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  & 
     416                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) 
     417               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
     418                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) 
    531419               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    532420               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    … …  
    557445      !! 
    558446      !! ** Action : - Update (ua,va) with the now vorticity term trend 
    559       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    560       !!               and planetary vorticity trends) ('key_trddyn') 
    561447      !! 
    562448      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
    563449      !!---------------------------------------------------------------------- 
    564       ! 
    565450      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    566       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    567       !                                                           ! =nrvm (relative vorticity or metric) 
     451      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 
    568452      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    569453      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    570       !! 
    571       INTEGER  ::   ji, jj, jk                                    ! dummy loop indices 
    572       INTEGER  ::   ierr                                          ! local integer 
    573       REAL(wp) ::   zfac12, zua, zva                              ! local scalars 
    574       REAL(wp) ::   zmsk, ze3                                     ! local scalars 
    575       !                                                           !  3D workspace  
    576       REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz 
    577       REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse 
    578 #if defined key_vvl 
    579       REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T) 
    580 #else 
    581       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all 
    582 #endif 
     454      ! 
     455      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     456      INTEGER  ::   ierr         ! local integer 
     457      REAL(wp) ::   zua, zva     ! local scalars 
     458      REAL(wp) ::   zmsk, ze3    ! local scalars 
     459      ! 
     460      REAL(wp), POINTER, DIMENSION(:,:)   :: zwx, zwy, zwz, z1_e3f 
     461      REAL(wp), POINTER, DIMENSION(:,:)   :: ztnw, ztne, ztsw, ztse 
    583462      !!---------------------------------------------------------------------- 
    584463      ! 
    585464      IF( nn_timing == 1 )  CALL timing_start('vor_een') 
    586465      ! 
    587       CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        )  
    588       CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )  
    589 #if defined key_vvl 
    590       CALL wrk_alloc( jpi, jpj, jpk, ze3f                   ) 
    591 #endif 
     466      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f )  
     467      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   )  
    592468      ! 
    593469      IF( kt == nit000 ) THEN 
    … …  
    595471         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 
    596472         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    597 #if ! defined key_vvl 
    598          IF( .NOT.ALLOCATED(ze3f) ) THEN 
    599             ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr ) 
    600             IF( lk_mpp    )   CALL mpp_sum ( ierr ) 
    601             IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' ) 
    602          ENDIF 
    603          ze3f(:,:,:) = 0.d0 
    604 #endif 
    605473      ENDIF 
    606  
    607       IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 
    608          DO jk = 1, jpk 
    609             DO jj = 1, jpjm1 
    610                DO ji = 1, jpim1 
    611                   ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    612                      &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
    613                   zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    614                      &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
    615                   IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3 
    616                END DO 
    617             END DO 
    618          END DO 
    619          CALL lbc_lnk( ze3f, 'F', 1. ) 
    620       ENDIF 
    621  
    622       zfac12 = 1._wp / 12._wp    ! Local constant initialization 
    623  
    624        
    625 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 
     474      ! 
    626475      !                                                ! =============== 
    627476      DO jk = 1, jpkm1                                 ! Horizontal slab 
    628477         !                                             ! =============== 
    629           
    630          ! Potential vorticity and horizontal fluxes 
    631          ! ----------------------------------------- 
    632          SELECT CASE( kvor )      ! vorticity considered 
    633          CASE ( 1 )                                                ! planetary vorticity (Coriolis) 
    634             zwz(:,:) = ff(:,:)      * ze3f(:,:,jk) 
    635          CASE ( 2 )                                                ! relative  vorticity 
    636             zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 
    637          CASE ( 3 )                                                ! metric term 
     478         ! 
     479         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point 
     480         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     481            DO jj = 1, jpjm1 
     482               DO ji = 1, fs_jpim1   ! vector opt. 
     483                  ze3  = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     484                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  ) 
     485                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3 
     486                  ELSE                      ;   z1_e3f(ji,jj) = 0._wp 
     487                  ENDIF 
     488               END DO 
     489            END DO 
     490         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     491            DO jj = 1, jpjm1 
     492               DO ji = 1, fs_jpim1   ! vector opt. 
     493                  ze3  = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     494                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  ) 
     495                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
     496                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  ) 
     497                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3 
     498                  ELSE                      ;   z1_e3f(ji,jj) = 0._wp 
     499                  ENDIF 
     500               END DO 
     501            END DO 
     502         END SELECT 
     503         ! 
     504         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     505         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     506            DO jj = 1, jpjm1 
     507               DO ji = 1, fs_jpim1   ! vector opt. 
     508                  zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj) 
     509               END DO 
     510            END DO 
     511         CASE ( np_RVO )                           !* relative vorticity 
     512            DO jj = 1, jpjm1 
     513               DO ji = 1, fs_jpim1   ! vector opt. 
     514                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     515                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     516                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
     517               END DO 
     518            END DO 
     519         CASE ( np_MET )                           !* metric term 
    638520            DO jj = 1, jpjm1 
    639521               DO ji = 1, fs_jpim1   ! vector opt. 
    640522                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    641523                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    642                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) 
    643                END DO 
    644             END DO 
    645             CALL lbc_lnk( zwz, 'F', 1. ) 
    646         CASE ( 4 )                                                ! total (relative + planetary vorticity) 
    647             zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 
    648          CASE ( 5 )                                                ! total (coriolis + metric) 
    649             DO jj = 1, jpjm1 
    650                DO ji = 1, fs_jpim1   ! vector opt. 
    651                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    652                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    653                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    654                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    655                        &       ) * ze3f(ji,jj,jk) 
    656                END DO 
    657             END DO 
    658             CALL lbc_lnk( zwz, 'F', 1. ) 
     524                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
     525               END DO 
     526            END DO 
     527         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     528            DO jj = 1, jpjm1 
     529               DO ji = 1, fs_jpim1   ! vector opt. 
     530                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     531                     &                         - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     532                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj) 
     533               END DO 
     534            END DO 
     535         CASE ( np_CME )                           !* Coriolis + metric 
     536            DO jj = 1, jpjm1 
     537               DO ji = 1, fs_jpim1   ! vector opt. 
     538                  zwz(ji,jj) = (  ff(ji,jj)                                                                        & 
     539                       &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     540                       &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     541                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
     542               END DO 
     543            END DO 
     544         CASE DEFAULT                                             ! error 
     545            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    659546         END SELECT 
    660  
    661          zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    662          zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    663  
    664          ! Compute and add the vorticity term trend 
    665          ! ---------------------------------------- 
     547         ! 
     548         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
     549            DO jj = 1, jpjm1 
     550               DO ji = 1, fs_jpim1   ! vector opt. 
     551                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     552               END DO 
     553            END DO 
     554         ENDIF 
     555         ! 
     556         CALL lbc_lnk( zwz, 'F', 1. ) 
     557         ! 
     558         !                                   !==  horizontal fluxes  ==! 
     559         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
     560         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     561 
     562         !                                   !==  compute and add the vorticity term trend  =! 
    666563         jj = 2 
    667564         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
    668          DO ji = 2, jpi    
     565         DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
    669566               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
    670567               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
    … …  
    682579         DO jj = 2, jpjm1 
    683580            DO ji = fs_2, fs_jpim1   ! vector opt. 
    684                zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    685                   &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    686                zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    687                   &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     581               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     582                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     583               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
     584                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    688585               pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
    689586               pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
    … …  
    693590      END DO                                           !   End of slab 
    694591      !                                                ! =============== 
    695       CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        )  
    696       CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )  
    697 #if defined key_vvl 
    698       CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   ) 
    699 #endif 
     592      ! 
     593      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f )  
     594      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   )  
    700595      ! 
    701596      IF( nn_timing == 1 )  CALL timing_stop('vor_een') 
    … …  
    715610      INTEGER ::   ios             ! Local integer output status for namelist read 
    716611      !! 
    717       NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een 
     612      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk 
    718613      !!---------------------------------------------------------------------- 
    719614 
    … …  
    732627         WRITE(numout,*) '~~~~~~~~~~~~' 
    733628         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme' 
    734          WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene 
    735          WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens 
    736          WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix 
    737          WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een 
     629         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene 
     630         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens 
     631         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
     632         WRITE(numout,*) '           enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een 
     633         WRITE(numout,*) '              e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
     634         WRITE(numout,*) '           masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
    738635      ENDIF 
    739636 
     637!!gm  this should be removed when choosing a unique strategy for fmask at the coast 
    740638      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 
    741639      ! at angles with three ocean points and one land point 
     640      IF(lwp) WRITE(numout,*) 
     641      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat 
    742642      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 
    743643         DO jk = 1, jpk 
    … …  
    753653          ! 
    754654      ENDIF 
    755  
    756       ioptio = 0                     ! Control of vorticity scheme options 
    757       IF( ln_dynvor_ene )   ioptio = ioptio + 1 
    758       IF( ln_dynvor_ens )   ioptio = ioptio + 1 
    759       IF( ln_dynvor_mix )   ioptio = ioptio + 1 
    760       IF( ln_dynvor_een )   ioptio = ioptio + 1 
    761       IF( lk_esopa      )   ioptio =          1 
    762  
     655!!gm end 
     656 
     657      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme) 
     658      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENE   ;   ENDIF 
     659      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENS   ;   ENDIF 
     660      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_MIX   ;   ENDIF 
     661      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF 
     662      ! 
    763663      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 
    764  
    765       !                              ! Set nvor (type of scheme for vorticity) 
    766       IF( ln_dynvor_ene )   nvor =  0 
    767       IF( ln_dynvor_ens )   nvor =  1 
    768       IF( ln_dynvor_mix )   nvor =  2 
    769       IF( ln_dynvor_een )   nvor =  3 
    770       IF( lk_esopa      )   nvor = -1 
    771        
    772       !                              ! Set ncor, nrvm, ntot (type of vorticity) 
    773       IF(lwp) WRITE(numout,*) 
    774       ncor = 1 
     664      !                       
     665      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot) 
     666      ncor = np_COR 
    775667      IF( ln_dynadv_vec ) THEN      
    776668         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity' 
    777          nrvm = 2 
    778          ntot = 4 
     669         nrvm = np_RVO        ! relative vorticity 
     670         ntot = np_CRV        ! relative + planetary vorticity 
    779671      ELSE                         
    780672         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term' 
    781          nrvm = 3 
    782          ntot = 5 
     673         nrvm = np_MET        ! metric term 
     674         ntot = np_CME        ! Coriolis + metric term 
    783675      ENDIF 
    784676       
    785677      IF(lwp) THEN                   ! Print the choice 
    786678         WRITE(numout,*) 
    787          IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme' 
    788          IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme' 
    789          IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme' 
    790          IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme' 
    791          IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options' 
     679         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme' 
     680         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme' 
     681         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme' 
     682         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme' 
    792683      ENDIF 
    793684      ! 
Note: See TracChangeset for help on using the changeset viewer.