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_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynvor.F90 @ 14062

Last change on this file since 14062 was 14062, checked in by ayoung, 3 years ago

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