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.
dynvor.F90 in NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90 @ 13696

Last change on this file since 13696 was 13696, checked in by techene, 4 years ago

#2555 use same e3f definition as in EEN in ENS and ENE instead of previous ugly fix and change time splitting accordingly change namelist vorticity scheme param nn_een_e3f into nn_e3f_typ
  • Property svn:keywords set to Id
File size: 58.4 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            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
19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
20   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
23   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity
24   !!----------------------------------------------------------------------
25
26   !!----------------------------------------------------------------------
27   !!   dyn_vor       : Update the momentum trend with the vorticity trend
28   !!       vor_enT   : energy conserving scheme at T-pt  (ln_dynvor_enT=T)
29   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
30   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
31   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
32   !!       vor_eeT   : energy conserving at T-pt         (ln_dynvor_eeT=T)
33   !!   dyn_vor_init  : set and control of the different vorticity option
34   !!----------------------------------------------------------------------
35   USE oce            ! ocean dynamics and tracers
36   USE dom_oce        ! ocean space and time domain
37   USE dommsk         ! ocean mask
38   USE dynadv         ! momentum advection
39   USE trd_oce        ! trends: ocean variables
40   USE trddyn         ! trend manager: dynamics
41   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
42   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
43   !
44   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
45   USE prtctl         ! Print control
46   USE in_out_manager ! I/O manager
47   USE lib_mpp        ! MPP library
48   USE timing         ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   dyn_vor        ! routine called by step.F90
54   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
55
56   !                                   !!* Namelist namdyn_vor: vorticity term
57   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
58   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
59   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
60   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
61   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
62   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
63   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
64   INTEGER, PUBLIC ::   nn_e3f_typ      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
65
66   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
67   !                                       ! associated indices:
68   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
69   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
70   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
71   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
72   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
73   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
74
75   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
76   !                               ! associated indices:
77   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
78   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
79   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
80   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
81   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
82
83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
84   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
85   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
86   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
87   !
88   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   e3f_0vor   ! e3f used in EEN, ENE and ENS cases (key_qco only)
89   
90   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
91   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
92   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
93   
94   !! * Substitutions
95#  include "do_loop_substitute.h90"
96#  include "domzgr_substitute.h90"
97
98   !!----------------------------------------------------------------------
99   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
100   !! $Id$
101   !! Software governed by the CeCILL license (see ./LICENSE)
102   !!----------------------------------------------------------------------
103CONTAINS
104
105   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
106      !!----------------------------------------------------------------------
107      !!
108      !! ** Purpose :   compute the lateral ocean tracer physics.
109      !!
110      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
111      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
112      !!               and planetary vorticity trends) and send them to trd_dyn
113      !!               for futher diagnostics (l_trddyn=T)
114      !!----------------------------------------------------------------------
115      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index
116      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices
117      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
118      !
119      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
120      !!----------------------------------------------------------------------
121      !
122      IF( ln_timing )   CALL timing_start('dyn_vor')
123      !
124      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
125         !
126         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
127         !
128         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
129         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
130         SELECT CASE( nvor_scheme )
131         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
132            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
133         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
134            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
135         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
136            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
137         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
138            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
139         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
140            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
141         END SELECT
142         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
143         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
144         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
145         !
146         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
147            ztrdu(:,:,:) = puu(:,:,:,Krhs)
148            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
149            SELECT CASE( nvor_scheme )
150            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
151            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
152            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
153            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
154            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
155            END SELECT
156            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
157            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
158            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
159         ENDIF
160         !
161         DEALLOCATE( ztrdu, ztrdv )
162         !
163      ELSE              !==  total vorticity trend added to the general trend  ==!
164         !
165         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
166         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
167                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
168            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
169         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
170                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
171            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
172         CASE( np_ENE )                        !* energy conserving scheme
173                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
174            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
175         CASE( np_ENS )                        !* enstrophy conserving scheme
176                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
177            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
178         CASE( np_MIX )                        !* mixed ene-ens scheme
179                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
180                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
181            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
182         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
183                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
184            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
185         END SELECT
186         !
187      ENDIF
188      !
189      !                       ! print sum trends (used for debugging)
190      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
191         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
192      !
193      IF( ln_timing )   CALL timing_stop('dyn_vor')
194      !
195   END SUBROUTINE dyn_vor
196
197
198   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
199      !!----------------------------------------------------------------------
200      !!                  ***  ROUTINE vor_enT  ***
201      !!
202      !! ** Purpose :   Compute the now total vorticity trend and add it to
203      !!      the general trend of the momentum equation.
204      !!
205      !! ** Method  :   Trend evaluated using now fields (centered in time)
206      !!       and t-point evaluation of vorticity (planetary and relative).
207      !!       conserves the horizontal kinetic energy.
208      !!         The general trend of momentum is increased due to the vorticity
209      !!       term which is given by:
210      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
211      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
212      !!       where rvor is the relative vorticity at f-point
213      !!
214      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
215      !!----------------------------------------------------------------------
216      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
217      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
218      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
219      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
220      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
221      !
222      INTEGER  ::   ji, jj, jk           ! dummy loop indices
223      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
224      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace
225      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwz      ! 3D workspace
226      !!----------------------------------------------------------------------
227      !
228      IF( kt == nit000 ) THEN
229         IF(lwp) WRITE(numout,*)
230         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
231         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
232      ENDIF
233      !
234      !
235      SELECT CASE( kvor )                 !== relative vorticity considered  ==!
236      !
237      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used
238         ALLOCATE( zwz(jpi,jpj,jpk) )
239         DO jk = 1, jpkm1                                ! Horizontal slab
240            DO_2D( 1, 0, 1, 0 )
241               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
242                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
243            END_2D
244            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity
245               DO_2D( 1, 0, 1, 0 )
246                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
247               END_2D
248            ENDIF
249         END DO
250         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
251         !
252      END SELECT
253
254      !                                                ! ===============
255      DO jk = 1, jpkm1                                 ! Horizontal slab
256         !                                             ! ===============
257         !
258         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
259         !
260         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
261            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
262         CASE ( np_RVO )                           !* relative vorticity
263            DO_2D( 0, 1, 0, 1 )
264               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)       &
265                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk)   )   &
266                  &              * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
267            END_2D
268         CASE ( np_MET )                           !* metric term
269            DO_2D( 0, 1, 0, 1 )
270               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)       &
271                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   )   &
272                  &       * e3t(ji,jj,jk,Kmm)
273            END_2D
274         CASE ( np_CRV )                           !* Coriolis + relative vorticity
275            DO_2D( 0, 1, 0, 1 )
276               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)        &
277                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  )   &
278                  &       * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
279            END_2D
280         CASE ( np_CME )                           !* Coriolis + metric
281            DO_2D( 0, 1, 0, 1 )
282               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                               &
283                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)      &
284                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  )   &
285                    &     * e3t(ji,jj,jk,Kmm)
286            END_2D
287         CASE DEFAULT                                             ! error
288            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
289         END SELECT
290         !
291         !                                   !==  compute and add the vorticity term trend  =!
292         DO_2D( 0, 0, 0, 0 )
293            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
294               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
295               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
296               !
297            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
298               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
299               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
300         END_2D
301         !                                             ! ===============
302      END DO                                           !   End of slab
303      !                                                ! ===============
304      !
305      SELECT CASE( kvor )        ! deallocate zwz if necessary
306      CASE ( np_RVO , np_CRV )   ;   DEALLOCATE( zwz )
307      END SELECT
308      !
309   END SUBROUTINE vor_enT
310
311
312   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
313      !!----------------------------------------------------------------------
314      !!                  ***  ROUTINE vor_ene  ***
315      !!
316      !! ** Purpose :   Compute the now total vorticity trend and add it to
317      !!      the general trend of the momentum equation.
318      !!
319      !! ** Method  :   Trend evaluated using now fields (centered in time)
320      !!       and the Sadourny (1975) flux form formulation : conserves the
321      !!       horizontal kinetic energy.
322      !!         The general trend of momentum is increased due to the vorticity
323      !!       term which is given by:
324      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
325      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
326      !!       where rvor is the relative vorticity
327      !!
328      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
329      !!
330      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
331      !!----------------------------------------------------------------------
332      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
333      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
334      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
335      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
336      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
337      !
338      INTEGER  ::   ji, jj, jk           ! dummy loop indices
339      REAL(wp) ::   zx1, zy1, zx2, zy2, ze3f, zmsk   ! local scalars
340      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
341      !!----------------------------------------------------------------------
342      !
343      IF( kt == nit000 ) THEN
344         IF(lwp) WRITE(numout,*)
345         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
346         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
347      ENDIF
348      !
349      !                                                ! ===============
350      DO jk = 1, jpkm1                                 ! Horizontal slab
351         !                                             ! ===============
352         !
353         SELECT CASE( kvor )                 !==  vorticity considered  ==!
354         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
355            zwz(:,:) = ff_f(:,:) 
356         CASE ( np_RVO )                           !* relative vorticity
357            DO_2D( 1, 0, 1, 0 )
358               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
359                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
360            END_2D
361            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
362               DO_2D( 1, 0, 1, 0 )
363                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
364               END_2D
365            ENDIF
366         CASE ( np_MET )                           !* metric term
367            DO_2D( 1, 0, 1, 0 )
368               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
369                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
370            END_2D
371         CASE ( np_CRV )                           !* Coriolis + relative vorticity
372            DO_2D( 1, 0, 1, 0 )
373               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
374                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
375            END_2D
376            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
377               DO_2D( 1, 0, 1, 0 )
378                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
379               END_2D
380            ENDIF
381         CASE ( np_CME )                           !* Coriolis + metric
382            DO_2D( 1, 0, 1, 0 )
383               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
384                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
385            END_2D
386         CASE DEFAULT                                             ! error
387            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
388         END SELECT
389         !
390#if defined key_qco
391         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
392            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
393         END_2D
394#else
395         SELECT CASE( nn_e3f_typ  )           !==  potential vorticity  ==!
396         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
397            DO_2D( 1, 0, 1, 0 )
398               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
399                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
400                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
401                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
402               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
403               ELSE                       ;   zwz(ji,jj) = 0._wp
404               ENDIF
405            END_2D
406         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
407            DO_2D( 1, 0, 1, 0 )
408               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
409                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
410                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
411                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
412               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
413                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
414               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
415               ELSE                       ;   zwz(ji,jj) = 0._wp
416               ENDIF
417            END_2D
418         END SELECT
419#endif
420         !                                   !==  horizontal fluxes  ==!
421         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
422         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
423         !
424         !                                   !==  compute and add the vorticity term trend  =!
425         DO_2D( 0, 0, 0, 0 )
426            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
427            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
428            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
429            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
430            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
431            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
432         END_2D
433         !                                             ! ===============
434      END DO                                           !   End of slab
435      !                                                ! ===============
436   END SUBROUTINE vor_ene
437
438
439   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
440      !!----------------------------------------------------------------------
441      !!                ***  ROUTINE vor_ens  ***
442      !!
443      !! ** Purpose :   Compute the now total vorticity trend and add it to
444      !!      the general trend of the momentum equation.
445      !!
446      !! ** Method  :   Trend evaluated using now fields (centered in time)
447      !!      and the Sadourny (1975) flux FORM formulation : conserves the
448      !!      potential enstrophy of a horizontally non-divergent flow. the
449      !!      trend of the vorticity term is given by:
450      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
451      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
452      !!      Add this trend to the general momentum trend:
453      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
454      !!
455      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
456      !!
457      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
458      !!----------------------------------------------------------------------
459      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
460      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
461      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
462      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
463      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
464      !
465      INTEGER  ::   ji, jj, jk   ! dummy loop indices
466      REAL(wp) ::   zuav, zvau, ze3f, zmsk   ! local scalars
467      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
468      !!----------------------------------------------------------------------
469      !
470      IF( kt == nit000 ) THEN
471         IF(lwp) WRITE(numout,*)
472         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
473         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
474      ENDIF
475      !                                                ! ===============
476      DO jk = 1, jpkm1                                 ! Horizontal slab
477         !                                             ! ===============
478         !
479         SELECT CASE( kvor )                 !==  vorticity considered  ==!
480         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
481            zwz(:,:) = ff_f(:,:) 
482         CASE ( np_RVO )                           !* relative vorticity
483            DO_2D( 1, 0, 1, 0 )
484               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
485                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
486            END_2D
487            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
488               DO_2D( 1, 0, 1, 0 )
489                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
490               END_2D
491            ENDIF
492         CASE ( np_MET )                           !* metric term
493            DO_2D( 1, 0, 1, 0 )
494               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
495                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
496            END_2D
497         CASE ( np_CRV )                           !* Coriolis + relative vorticity
498            DO_2D( 1, 0, 1, 0 )
499               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
500                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
501            END_2D
502            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
503               DO_2D( 1, 0, 1, 0 )
504                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
505               END_2D
506            ENDIF
507         CASE ( np_CME )                           !* Coriolis + metric
508            DO_2D( 1, 0, 1, 0 )
509               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
510                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
511            END_2D
512         CASE DEFAULT                                             ! error
513            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
514         END SELECT
515         !
516         !
517#if defined key_qco
518         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
519            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
520         END_2D
521#else
522         SELECT CASE( nn_e3f_typ )           !==  potential vorticity  ==!
523         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
524            DO_2D( 1, 0, 1, 0 )
525               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
526                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
527                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
528                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
529               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
530               ELSE                       ;   zwz(ji,jj) = 0._wp
531               ENDIF
532            END_2D
533         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
534            DO_2D( 1, 0, 1, 0 )
535               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
536                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
537                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
538                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
539               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
540                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
541               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
542               ELSE                       ;   zwz(ji,jj) = 0._wp
543               ENDIF
544            END_2D
545         END SELECT
546#endif
547         !                                   !==  horizontal fluxes  ==!
548         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
549         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
550         !
551         !                                   !==  compute and add the vorticity term trend  =!
552         DO_2D( 0, 0, 0, 0 )
553            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
554               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
555            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
556               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
557            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
558            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
559         END_2D
560         !                                             ! ===============
561      END DO                                           !   End of slab
562      !                                                ! ===============
563   END SUBROUTINE vor_ens
564
565
566   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
567      !!----------------------------------------------------------------------
568      !!                ***  ROUTINE vor_een  ***
569      !!
570      !! ** Purpose :   Compute the now total vorticity trend and add it to
571      !!      the general trend of the momentum equation.
572      !!
573      !! ** Method  :   Trend evaluated using now fields (centered in time)
574      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
575      !!      both the horizontal kinetic energy and the potential enstrophy
576      !!      when horizontal divergence is zero (see the NEMO documentation)
577      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
578      !!
579      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
580      !!
581      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
582      !!----------------------------------------------------------------------
583      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
584      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
585      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
586      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
587      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
588      !
589      INTEGER  ::   ji, jj, jk   ! dummy loop indices
590      INTEGER  ::   ierr         ! local integer
591      REAL(wp) ::   zua, zva     ! local scalars
592      REAL(wp) ::   zmsk, ze3f   ! local scalars
593      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy , z1_e3f
594      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
595      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
596      !!----------------------------------------------------------------------
597      !
598      IF( kt == nit000 ) THEN
599         IF(lwp) WRITE(numout,*)
600         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
601         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
602      ENDIF
603      !
604      !                                                ! ===============
605      DO jk = 1, jpkm1                                 ! Horizontal slab
606         !                                             ! ===============
607         !
608#if defined key_qco
609         DO_2D( 1, 0, 1, 0 )                 ! == reciprocal of e3 at F-point (key_qco)
610            z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
611         END_2D
612#else
613         SELECT CASE( nn_e3f_typ )           ! == reciprocal of e3 at F-point
614         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
615            DO_2D( 1, 0, 1, 0 )
616               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
617                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
618                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
619                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
620               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
621               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
622               ENDIF
623            END_2D
624         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
625            DO_2D( 1, 0, 1, 0 )
626               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
627                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
628                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
629                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
630               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
631                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
632               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
633               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
634               ENDIF
635            END_2D
636         END SELECT
637#endif
638         !
639         SELECT CASE( kvor )                 !==  vorticity considered  ==!
640         !
641         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
642            DO_2D( 1, 0, 1, 0 )
643               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
644            END_2D
645         CASE ( np_RVO )                           !* relative vorticity
646            DO_2D( 1, 0, 1, 0 )
647               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
648                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
649            END_2D
650            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
651               DO_2D( 1, 0, 1, 0 )
652                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
653               END_2D
654            ENDIF
655         CASE ( np_MET )                           !* metric term
656            DO_2D( 1, 0, 1, 0 )
657               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
658                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
659            END_2D
660         CASE ( np_CRV )                           !* Coriolis + relative vorticity
661            DO_2D( 1, 0, 1, 0 )
662               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
663                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
664                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
665            END_2D
666            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
667               DO_2D( 1, 0, 1, 0 )
668                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
669               END_2D
670            ENDIF
671         CASE ( np_CME )                           !* Coriolis + metric
672            DO_2D( 1, 0, 1, 0 )
673               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
674                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
675            END_2D
676         CASE DEFAULT                                             ! error
677            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
678         END SELECT
679         !                                             ! ===============
680      END DO                                           !   End of slab
681      !                                                ! ===============
682      !
683      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
684      !
685      !                                                ! ===============
686      DO jk = 1, jpkm1                                 ! Horizontal slab
687         !                                             ! ===============
688         !
689         !                                   !==  horizontal fluxes  ==!
690         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
691         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
692         !
693         !                                   !==  compute and add the vorticity term trend  =!
694         DO_2D( 0, 1, 0, 1 )
695            ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
696            ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
697            ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
698            ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
699         END_2D
700         !
701         DO_2D( 0, 0, 0, 0 )
702            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
703               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
704            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
705               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
706            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
707            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
708         END_2D
709         !                                             ! ===============
710      END DO                                           !   End of slab
711      !                                                ! ===============
712   END SUBROUTINE vor_een
713
714
715   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
716      !!----------------------------------------------------------------------
717      !!                ***  ROUTINE vor_eeT  ***
718      !!
719      !! ** Purpose :   Compute the now total vorticity trend and add it to
720      !!      the general trend of the momentum equation.
721      !!
722      !! ** Method  :   Trend evaluated using now fields (centered in time)
723      !!      and the Arakawa and Lamb (1980) vector form formulation using
724      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
725      !!      The change consists in
726      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
727      !!
728      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
729      !!
730      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
731      !!----------------------------------------------------------------------
732      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
733      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
734      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
735      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
736      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
737      !
738      INTEGER  ::   ji, jj, jk     ! dummy loop indices
739      INTEGER  ::   ierr           ! local integer
740      REAL(wp) ::   zua, zva       ! local scalars
741      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
742      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy 
743      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
744      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
745      !!----------------------------------------------------------------------
746      !
747      IF( kt == nit000 ) THEN
748         IF(lwp) WRITE(numout,*)
749         IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
750         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
751      ENDIF
752      !
753      !                                                ! ===============
754      DO jk = 1, jpkm1                                 ! Horizontal slab
755         !                                             ! ===============
756         !
757         !
758         SELECT CASE( kvor )                 !==  vorticity considered  ==!
759         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
760            DO_2D( 1, 0, 1, 0 )
761               zwz(ji,jj,jk) = ff_f(ji,jj)
762            END_2D
763         CASE ( np_RVO )                           !* relative vorticity
764            DO_2D( 1, 0, 1, 0 )
765               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
766                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
767                  &          * r1_e1e2f(ji,jj)
768            END_2D
769            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
770               DO_2D( 1, 0, 1, 0 )
771                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
772               END_2D
773            ENDIF
774         CASE ( np_MET )                           !* metric term
775            DO_2D( 1, 0, 1, 0 )
776               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
777                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
778            END_2D
779         CASE ( np_CRV )                           !* Coriolis + relative vorticity
780            DO_2D( 1, 0, 1, 0 )
781               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
782                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
783                  &                         * r1_e1e2f(ji,jj)    )
784            END_2D
785            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
786               DO_2D( 1, 0, 1, 0 )
787                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
788               END_2D
789            ENDIF
790         CASE ( np_CME )                           !* Coriolis + metric
791            DO_2D( 1, 0, 1, 0 )
792               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
793                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
794            END_2D
795         CASE DEFAULT                                             ! error
796            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
797         END SELECT
798         !
799         !                                             ! ===============
800      END DO                                           !   End of slab
801      !                                                ! ===============
802      !
803      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
804      !
805      !                                                ! ===============
806      DO jk = 1, jpkm1                                 ! Horizontal slab
807         !                                             ! ===============
808         !
809         !                                   !==  horizontal fluxes  ==!
810         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
811         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
812         !
813         !                                   !==  compute and add the vorticity term trend  =!
814         DO_2D( 0, 1, 0, 1 )
815            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
816            ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
817            ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
818            ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
819            ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
820         END_2D
821         !
822         DO_2D( 0, 0, 0, 0 )
823            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
824               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
825            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
826               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
827            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
828            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
829         END_2D
830         !                                             ! ===============
831      END DO                                           !   End of slab
832      !                                                ! ===============
833   END SUBROUTINE vor_eeT
834
835
836   SUBROUTINE dyn_vor_init
837      !!---------------------------------------------------------------------
838      !!                  ***  ROUTINE dyn_vor_init  ***
839      !!
840      !! ** Purpose :   Control the consistency between cpp options for
841      !!              tracer advection schemes
842      !!----------------------------------------------------------------------
843      INTEGER ::   ji, jj, jk    ! dummy loop indices
844      INTEGER ::   ioptio, ios   ! local integer
845      REAL(wp) ::   zmsk    ! local scalars
846      !!
847      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
848         &                 ln_dynvor_een, nn_e3f_typ   , ln_dynvor_mix, ln_dynvor_msk
849      !!----------------------------------------------------------------------
850      !
851      IF(lwp) THEN
852         WRITE(numout,*)
853         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
854         WRITE(numout,*) '~~~~~~~~~~~~'
855      ENDIF
856      !
857      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
858901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
859      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
860902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
861      IF(lwm) WRITE ( numond, namdyn_vor )
862      !
863      IF(lwp) THEN                    ! Namelist print
864         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
865         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
866         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
867         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
868         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
869         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
870         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_e3f_typ = ', nn_e3f_typ
871         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
872         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
873      ENDIF
874
875!!gm  this should be removed when choosing a unique strategy for fmask at the coast
876      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
877      ! at angles with three ocean points and one land point
878      IF(lwp) WRITE(numout,*)
879      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
880      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
881         DO_3D( 1, 0, 1, 0, 1, jpk )
882            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
883               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
884         END_3D
885         !
886         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
887         !
888      ENDIF
889!!gm end
890
891      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
892      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
893      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
894      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
895      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
896      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
897      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
898      !
899      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
900      !                     
901      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
902      ncor = np_COR                       ! planetary vorticity
903      SELECT CASE( n_dynadv )
904      CASE( np_LIN_dyn )
905         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
906         nrvm = np_COR        ! planetary vorticity
907         ntot = np_COR        !     -         -
908      CASE( np_VEC_c2  )
909         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
910         nrvm = np_RVO        ! relative vorticity
911         ntot = np_CRV        ! relative + planetary vorticity         
912      CASE( np_FLX_c2 , np_FLX_ubs  )
913         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
914         nrvm = np_MET        ! metric term
915         ntot = np_CME        ! Coriolis + metric term
916         !
917         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
918         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
919            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
920            DO_2D( 0, 0, 0, 0 )
921               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
922               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
923            END_2D
924            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
925            !
926         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
927            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
928            DO_2D( 1, 0, 1, 0 )
929               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
930               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
931            END_2D
932            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
933         END SELECT
934         !
935      END SELECT
936!!st ADD !! pense a recalculer le hf_0
937#if defined key_qco
938      SELECT CASE( nvor_scheme )    ! qco case: pre-computed a specific e3f_0 for some vorticity schemes
939      CASE( np_ENS , np_ENE , np_EEN , np_MIX )
940         !
941         ALLOCATE( e3f_0vor(jpi,jpj,jpk) )
942         !
943         SELECT CASE( nn_e3f_typ )
944         CASE ( 0 )                        ! original formulation  (masked averaging of e3t divided by 4)
945            DO_3D( 0, 0, 0, 0, 1, jpk )
946               e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
947                  &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
948                  &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
949                  &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) * 0.25_wp
950            END_3D
951         CASE ( 1 )                        ! new formulation  (masked averaging of e3t divided by the sum of mask)
952            DO_3D( 0, 0, 0, 0, 1, jpk )
953               zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk)   &
954                  &  + tmask(ji,jj  ,jk) +tmask(ji+1,jj  ,jk)  )
955               !
956               IF( zmsk /= 0._wp ) THEN
957                  e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
958                     &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
959                     &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
960                     &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) / zmsk
961               ENDIF
962            END_3D
963         END SELECT
964         !
965         CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp )
966         !                                 ! insure e3f_0vor /= 0
967         WHERE( e3f_0vor(:,:,:) == 0._wp )   e3f_0vor(:,:,:) = e3f_0(:,:,:)
968         !
969      END SELECT
970      !
971#endif
972!!st END     
973      IF(lwp) THEN                   ! Print the choice
974         WRITE(numout,*)
975         SELECT CASE( nvor_scheme )
976         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
977         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
978         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
979                              IF( ln_dynadv_vec )   CALL ctl_warn('dyn_vor_init: ENT scheme may not work in vector form')
980         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
981         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
982         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
983         END SELECT         
984      ENDIF
985      !
986   END SUBROUTINE dyn_vor_init
987
988   !!==============================================================================
989END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.