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/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/DYN/dynvor.F90 @ 14757

Last change on this file since 14757 was 14757, checked in by francesca, 3 years ago

Fortran 77 '.EQ.' operator replacement in conditional statements; [comm_cleanup] tags removal - ticket #2607

  • Property svn:keywords set to Id
File size: 60.6 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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
264                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
265               END_2D
266            ENDIF
267         END DO
268         IF (nn_hls==1) 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   ||   defined key_linssh
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   ||   defined key_linssh
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   ||   defined key_linssh
627         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                 ! == 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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
680            ! round brackets added to fix the order of floating point operations
681            ! needed to ensure halo 1 - halo 2 compatibility
682               zwz(ji,jj,jk) = (  ff_f(ji,jj) + ( ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
683                  &                               )                                                                  & ! bracket for halo 1 - halo 2 compatibility
684                  &                             - ( e1u(ji  ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)      &
685                  &                               )                                                                  & ! bracket for halo 1 - halo 2 compatibility
686                  &                             ) * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
687            END_2D
688            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
689               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
690                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
691               END_2D
692            ENDIF
693         CASE ( np_CME )                           !* Coriolis + metric
694            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
695               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
696                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
697            END_2D
698         CASE DEFAULT                                             ! error
699            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
700         END SELECT
701         !                                             ! ===============
702      END DO                                           !   End of slab
703      !                                                ! ===============
704      !
705      IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
706      !
707      !                                                ! ===============
708      DO jk = 1, jpkm1                                 ! Horizontal slab
709         !                                             ! ===============
710         !
711         !                                   !==  horizontal fluxes  ==!
712         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
713         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
714         !
715         !                                   !==  compute and add the vorticity term trend  =!
716         DO_2D( 0, 1, 0, 1 )
717            ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
718            ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
719            ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
720            ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
721         END_2D
722         !
723         DO_2D( 0, 0, 0, 0 )
724            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
725               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
726            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
727               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
728            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
729            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
730         END_2D
731         !                                             ! ===============
732      END DO                                           !   End of slab
733      !                                                ! ===============
734   END SUBROUTINE vor_een
735
736
737   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
738      !!----------------------------------------------------------------------
739      !!                ***  ROUTINE vor_eeT  ***
740      !!
741      !! ** Purpose :   Compute the now total vorticity trend and add it to
742      !!      the general trend of the momentum equation.
743      !!
744      !! ** Method  :   Trend evaluated using now fields (centered in time)
745      !!      and the Arakawa and Lamb (1980) vector form formulation using
746      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
747      !!      The change consists in
748      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
749      !!
750      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
751      !!
752      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
753      !!----------------------------------------------------------------------
754      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
755      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
756      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
757      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
758      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
759      !
760      INTEGER  ::   ji, jj, jk     ! dummy loop indices
761      INTEGER  ::   ierr           ! local integer
762      REAL(wp) ::   zua, zva       ! local scalars
763      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
764      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy
765      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
766      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined
767      !!----------------------------------------------------------------------
768      !
769      IF( kt == nit000 ) THEN
770         IF(lwp) WRITE(numout,*)
771         IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
772         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
773      ENDIF
774      !
775      !                                                ! ===============
776      DO jk = 1, jpkm1                                 ! Horizontal slab
777         !                                             ! ===============
778         !
779         !
780         SELECT CASE( kvor )                 !==  vorticity considered  ==!
781         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
782            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
783               zwz(ji,jj,jk) = ff_f(ji,jj)
784            END_2D
785         CASE ( np_RVO )                           !* relative vorticity
786            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
787               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
788                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
789                  &          * r1_e1e2f(ji,jj)
790            END_2D
791            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
792               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
793                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
794               END_2D
795            ENDIF
796         CASE ( np_MET )                           !* metric term
797            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
798               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
799                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
800            END_2D
801         CASE ( np_CRV )                           !* Coriolis + relative vorticity
802            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
803               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
804                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
805                  &                         * r1_e1e2f(ji,jj)    )
806            END_2D
807            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
808               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
809                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
810               END_2D
811            ENDIF
812         CASE ( np_CME )                           !* Coriolis + metric
813            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
814               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
815                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
816            END_2D
817         CASE DEFAULT                                             ! error
818            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
819         END SELECT
820         !
821         !                                             ! ===============
822      END DO                                           !   End of slab
823      !                                                ! ===============
824      !
825      IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
826      !
827      !                                                ! ===============
828      DO jk = 1, jpkm1                                 ! Horizontal slab
829         !                                             ! ===============
830         !
831         !                                   !==  horizontal fluxes  ==!
832         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
833         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
834         !
835         !                                   !==  compute and add the vorticity term trend  =!
836         DO_2D( 0, 1, 0, 1 )
837            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
838            ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
839            ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
840            ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
841            ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
842         END_2D
843         !
844         DO_2D( 0, 0, 0, 0 )
845            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
846               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
847            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
848               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
849            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
850            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
851         END_2D
852         !                                             ! ===============
853      END DO                                           !   End of slab
854      !                                                ! ===============
855   END SUBROUTINE vor_eeT
856
857
858   SUBROUTINE dyn_vor_init
859      !!---------------------------------------------------------------------
860      !!                  ***  ROUTINE dyn_vor_init  ***
861      !!
862      !! ** Purpose :   Control the consistency between cpp options for
863      !!              tracer advection schemes
864      !!----------------------------------------------------------------------
865      INTEGER ::   ji, jj, jk    ! dummy loop indices
866      INTEGER ::   ioptio, ios   ! local integer
867      REAL(wp) ::   zmsk    ! local scalars
868      !!
869      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
870         &                 ln_dynvor_een, nn_e3f_typ   , ln_dynvor_mix, ln_dynvor_msk
871      !!----------------------------------------------------------------------
872      !
873      IF(lwp) THEN
874         WRITE(numout,*)
875         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
876         WRITE(numout,*) '~~~~~~~~~~~~'
877      ENDIF
878      !
879      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
880901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
881      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
882902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
883      IF(lwm) WRITE ( numond, namdyn_vor )
884      !
885      IF(lwp) THEN                    ! Namelist print
886         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
887         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
888         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
889         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
890         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
891         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
892         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_e3f_typ = ', nn_e3f_typ
893         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
894         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
895      ENDIF
896
897!!gm  this should be removed when choosing a unique strategy for fmask at the coast
898      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
899      ! at angles with three ocean points and one land point
900      IF(lwp) WRITE(numout,*)
901      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
902      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
903         DO_3D( 1, 0, 1, 0, 1, jpk )
904            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
905               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
906         END_3D
907         !
908         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
909         !
910      ENDIF
911!!gm end
912
913      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
914      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
915      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
916      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
917      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
918      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
919      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
920      !
921      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
922      !
923      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
924      ncor = np_COR                       ! planetary vorticity
925      SELECT CASE( n_dynadv )
926      CASE( np_LIN_dyn )
927         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
928         nrvm = np_COR        ! planetary vorticity
929         ntot = np_COR        !     -         -
930      CASE( np_VEC_c2  )
931         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'
932         nrvm = np_RVO        ! relative vorticity
933         ntot = np_CRV        ! relative + planetary vorticity
934      CASE( np_FLX_c2 , np_FLX_ubs  )
935         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
936         nrvm = np_MET        ! metric term
937         ntot = np_CME        ! Coriolis + metric term
938         !
939         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
940         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
941            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
942            DO_2D( 0, 0, 0, 0 )
943               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
944               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
945            END_2D
946            CALL lbc_lnk( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
947            !
948         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
949            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
950            DO_2D( 1, 0, 1, 0 )
951               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
952               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
953            END_2D
954            CALL lbc_lnk( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
955         END SELECT
956         !
957      END SELECT
958#if defined key_qco   ||   defined key_linssh
959      SELECT CASE( nvor_scheme )    ! qco or linssh cases : pre-computed a specific e3f_0 for some vorticity schemes
960      CASE( np_ENS , np_ENE , np_EEN , np_MIX )
961         !
962         ALLOCATE( e3f_0vor(jpi,jpj,jpk) )
963         !
964         SELECT CASE( nn_e3f_typ )
965         CASE ( 0 )                        ! original formulation  (masked averaging of e3t divided by 4)
966            DO_3D( 0, 0, 0, 0, 1, jpk )
967               e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
968                  &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
969                  &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
970                  &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) * 0.25_wp
971            END_3D
972         CASE ( 1 )                        ! new formulation  (masked averaging of e3t divided by the sum of mask)
973            DO_3D( 0, 0, 0, 0, 1, jpk )
974               zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk)   &
975                  &  + tmask(ji,jj  ,jk) +tmask(ji+1,jj  ,jk)  )
976               !
977               IF( zmsk /= 0._wp ) THEN
978                  e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
979                     &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
980                     &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
981                     &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) / zmsk
982               ELSE ; e3f_0vor(ji,jj,jk) = 0._wp
983               ENDIF
984            END_3D
985         END SELECT
986         !
987         CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp )
988         !                                 ! insure e3f_0vor /= 0
989         WHERE( e3f_0vor(:,:,:) == 0._wp )   e3f_0vor(:,:,:) = e3f_0(:,:,:)
990         !
991      END SELECT
992      !
993#endif
994      IF(lwp) THEN                   ! Print the choice
995         WRITE(numout,*)
996         SELECT CASE( nvor_scheme )
997         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
998         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
999         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
1000                              IF( ln_dynadv_vec )   CALL ctl_warn('dyn_vor_init: ENT scheme may not work in vector form')
1001         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
1002         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
1003         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
1004         END SELECT
1005      ENDIF
1006      !
1007   END SUBROUTINE dyn_vor_init
1008
1009   !!==============================================================================
1010END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.