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/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/dynvor.F90

Last change on this file was 14834, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling
  • Property svn:keywords set to Id
File size: 64.4 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
17   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity
18   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
20   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
23   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity
24   !!            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(A2D(nn_hls))        ::   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( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
247         IF( kt == nit000 ) THEN
248            IF(lwp) WRITE(numout,*)
249            IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
250            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
251         ENDIF
252      ENDIF
253      !
254      !
255      SELECT CASE( kvor )                 !== relative vorticity considered  ==!
256      !
257      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used
258         ALLOCATE( zwz(A2D(nn_hls),jpk) )
259         DO jk = 1, jpkm1                                ! Horizontal slab
260            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
261               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
262                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
263            END_2D
264            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity
265               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
266                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
267               END_2D
268            ENDIF
269         END DO
270         IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
271         !
272      END SELECT
273
274      !                                                ! ===============
275      DO jk = 1, jpkm1                                 ! Horizontal slab
276         !                                             ! ===============
277         !
278         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
279         !
280         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
281            DO_2D( 0, 1, 0, 1 )
282               zwt(ji,jj) = ff_t(ji,jj) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
283            END_2D
284         CASE ( np_RVO )                           !* relative vorticity
285            DO_2D( 0, 1, 0, 1 )
286               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)       &
287                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk)   )   &
288                  &              * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
289            END_2D
290         CASE ( np_MET )                           !* metric term
291            DO_2D( 0, 1, 0, 1 )
292               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)       &
293                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   )   &
294                  &       * e3t(ji,jj,jk,Kmm)
295            END_2D
296         CASE ( np_CRV )                           !* Coriolis + relative vorticity
297            DO_2D( 0, 1, 0, 1 )
298               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)        &
299                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  )   &
300                  &       * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
301            END_2D
302         CASE ( np_CME )                           !* Coriolis + metric
303            DO_2D( 0, 1, 0, 1 )
304               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                               &
305                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)      &
306                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  )   &
307                    &     * e3t(ji,jj,jk,Kmm)
308            END_2D
309         CASE DEFAULT                                             ! error
310            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
311         END SELECT
312         !
313         !                                   !==  compute and add the vorticity term trend  =!
314         DO_2D( 0, 0, 0, 0 )
315            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
316               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
317               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
318               !
319            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
320               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &
321               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )
322         END_2D
323         !                                             ! ===============
324      END DO                                           !   End of slab
325      !                                                ! ===============
326      !
327      SELECT CASE( kvor )        ! deallocate zwz if necessary
328      CASE ( np_RVO , np_CRV )   ;   DEALLOCATE( zwz )
329      END SELECT
330      !
331   END SUBROUTINE vor_enT
332
333
334   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
335      !!----------------------------------------------------------------------
336      !!                  ***  ROUTINE vor_ene  ***
337      !!
338      !! ** Purpose :   Compute the now total vorticity trend and add it to
339      !!      the general trend of the momentum equation.
340      !!
341      !! ** Method  :   Trend evaluated using now fields (centered in time)
342      !!       and the Sadourny (1975) flux form formulation : conserves the
343      !!       horizontal kinetic energy.
344      !!         The general trend of momentum is increased due to the vorticity
345      !!       term which is given by:
346      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
347      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
348      !!       where rvor is the relative vorticity
349      !!
350      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
351      !!
352      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
353      !!----------------------------------------------------------------------
354      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
355      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
356      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
357      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
358      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
359      !
360      INTEGER  ::   ji, jj, jk           ! dummy loop indices
361      REAL(wp) ::   zx1, zy1, zx2, zy2, ze3f, zmsk   ! local scalars
362      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zwx, zwy, zwz   ! 2D workspace
363      !!----------------------------------------------------------------------
364      !
365      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
366         IF( kt == nit000 ) THEN
367            IF(lwp) WRITE(numout,*)
368            IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
369            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
370         ENDIF
371      ENDIF
372      !
373      !                                                ! ===============
374      DO jk = 1, jpkm1                                 ! Horizontal slab
375         !                                             ! ===============
376         !
377         SELECT CASE( kvor )                 !==  vorticity considered  ==!
378         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
379            DO_2D( 1, 0, 1, 0 )
380               zwz(ji,jj) = ff_f(ji,jj)
381            END_2D
382         CASE ( np_RVO )                           !* relative vorticity
383            DO_2D( 1, 0, 1, 0 )
384               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
385                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
386            END_2D
387            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
388               DO_2D( 1, 0, 1, 0 )
389                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
390               END_2D
391            ENDIF
392         CASE ( np_MET )                           !* metric term
393            DO_2D( 1, 0, 1, 0 )
394               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
395                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
396            END_2D
397         CASE ( np_CRV )                           !* Coriolis + relative vorticity
398            DO_2D( 1, 0, 1, 0 )
399               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
400                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
401            END_2D
402            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
403               DO_2D( 1, 0, 1, 0 )
404                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
405               END_2D
406            ENDIF
407         CASE ( np_CME )                           !* Coriolis + metric
408            DO_2D( 1, 0, 1, 0 )
409               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
410                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
411            END_2D
412         CASE DEFAULT                                             ! error
413            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
414         END SELECT
415         !
416#if defined key_qco   ||   defined key_linssh
417         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
418            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
419         END_2D
420#else
421         SELECT CASE( nn_e3f_typ  )           !==  potential vorticity  ==!
422         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
423            DO_2D( 1, 0, 1, 0 )
424               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
425                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
426                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
427                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
428               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
429               ELSE                       ;   zwz(ji,jj) = 0._wp
430               ENDIF
431            END_2D
432         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
433            DO_2D( 1, 0, 1, 0 )
434               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
435                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
436                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
437                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
438               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
439                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
440               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
441               ELSE                       ;   zwz(ji,jj) = 0._wp
442               ENDIF
443            END_2D
444         END SELECT
445#endif
446         !                                   !==  horizontal fluxes  ==!
447         DO_2D( 1, 1, 1, 1 )
448            zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
449            zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
450         END_2D
451         !
452         !                                   !==  compute and add the vorticity term trend  =!
453         DO_2D( 0, 0, 0, 0 )
454            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
455            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
456            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
457            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
458            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 )
459            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 )
460         END_2D
461         !                                             ! ===============
462      END DO                                           !   End of slab
463      !                                                ! ===============
464   END SUBROUTINE vor_ene
465
466
467   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
468      !!----------------------------------------------------------------------
469      !!                ***  ROUTINE vor_ens  ***
470      !!
471      !! ** Purpose :   Compute the now total vorticity trend and add it to
472      !!      the general trend of the momentum equation.
473      !!
474      !! ** Method  :   Trend evaluated using now fields (centered in time)
475      !!      and the Sadourny (1975) flux FORM formulation : conserves the
476      !!      potential enstrophy of a horizontally non-divergent flow. the
477      !!      trend of the vorticity term is given by:
478      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
479      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
480      !!      Add this trend to the general momentum trend:
481      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
482      !!
483      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
484      !!
485      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
486      !!----------------------------------------------------------------------
487      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
488      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
489      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
490      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
491      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
492      !
493      INTEGER  ::   ji, jj, jk   ! dummy loop indices
494      REAL(wp) ::   zuav, zvau, ze3f, zmsk   ! local scalars
495      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zwx, zwy, zwz   ! 2D workspace
496      !!----------------------------------------------------------------------
497      !
498      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
499         IF( kt == nit000 ) THEN
500            IF(lwp) WRITE(numout,*)
501            IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
502            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
503         ENDIF
504      ENDIF
505      !                                                ! ===============
506      DO jk = 1, jpkm1                                 ! Horizontal slab
507         !                                             ! ===============
508         !
509         SELECT CASE( kvor )                 !==  vorticity considered  ==!
510         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
511            DO_2D( 1, 0, 1, 0 )
512               zwz(ji,jj) = ff_f(ji,jj)
513            END_2D
514         CASE ( np_RVO )                           !* relative vorticity
515            DO_2D( 1, 0, 1, 0 )
516               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
517                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
518            END_2D
519            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
520               DO_2D( 1, 0, 1, 0 )
521                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
522               END_2D
523            ENDIF
524         CASE ( np_MET )                           !* metric term
525            DO_2D( 1, 0, 1, 0 )
526               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
527                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
528            END_2D
529         CASE ( np_CRV )                           !* Coriolis + relative vorticity
530            DO_2D( 1, 0, 1, 0 )
531               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
532                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
533            END_2D
534            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
535               DO_2D( 1, 0, 1, 0 )
536                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
537               END_2D
538            ENDIF
539         CASE ( np_CME )                           !* Coriolis + metric
540            DO_2D( 1, 0, 1, 0 )
541               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
542                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
543            END_2D
544         CASE DEFAULT                                             ! error
545            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
546         END SELECT
547         !
548         !
549#if defined key_qco   ||   defined key_linssh
550         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
551            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
552         END_2D
553#else
554         SELECT CASE( nn_e3f_typ )           !==  potential vorticity  ==!
555         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
556            DO_2D( 1, 0, 1, 0 )
557               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
558                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
559                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
560                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
561               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
562               ELSE                       ;   zwz(ji,jj) = 0._wp
563               ENDIF
564            END_2D
565         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
566            DO_2D( 1, 0, 1, 0 )
567               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
568                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
569                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
570                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
571               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
572                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
573               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
574               ELSE                       ;   zwz(ji,jj) = 0._wp
575               ENDIF
576            END_2D
577         END SELECT
578#endif
579         !                                   !==  horizontal fluxes  ==!
580         DO_2D( 1, 1, 1, 1 )
581            zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
582            zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
583         END_2D
584         !
585         !                                   !==  compute and add the vorticity term trend  =!
586         DO_2D( 0, 0, 0, 0 )
587            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
588               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
589            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
590               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
591            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
592            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
593         END_2D
594         !                                             ! ===============
595      END DO                                           !   End of slab
596      !                                                ! ===============
597   END SUBROUTINE vor_ens
598
599
600   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
601      !!----------------------------------------------------------------------
602      !!                ***  ROUTINE vor_een  ***
603      !!
604      !! ** Purpose :   Compute the now total vorticity trend and add it to
605      !!      the general trend of the momentum equation.
606      !!
607      !! ** Method  :   Trend evaluated using now fields (centered in time)
608      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
609      !!      both the horizontal kinetic energy and the potential enstrophy
610      !!      when horizontal divergence is zero (see the NEMO documentation)
611      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
612      !!
613      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
614      !!
615      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
616      !!----------------------------------------------------------------------
617      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
618      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
619      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
620      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
621      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
622      !
623      INTEGER  ::   ji, jj, jk   ! dummy loop indices
624      INTEGER  ::   ierr         ! local integer
625      REAL(wp) ::   zua, zva     ! local scalars
626      REAL(wp) ::   zmsk, ze3f   ! local scalars
627      REAL(wp), DIMENSION(A2D(nn_hls))       ::   z1_e3f
628#if defined key_loop_fusion
629      REAL(wp) ::   ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1
630      REAL(wp) ::   zwx, zwx_im1, zwx_jp1, zwx_im1_jp1
631      REAL(wp) ::   zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1
632#else
633      REAL(wp), DIMENSION(A2D(nn_hls))       ::   zwx , zwy
634      REAL(wp), DIMENSION(A2D(nn_hls))       ::   ztnw, ztne, ztsw, ztse
635#endif
636      REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) ::   zwz   ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
637      !!----------------------------------------------------------------------
638      !
639      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
640         IF( kt == nit000 ) THEN
641            IF(lwp) WRITE(numout,*)
642            IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
643            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
644         ENDIF
645      ENDIF
646      !
647      !                                                ! ===============
648      DO jk = 1, jpkm1                                 ! Horizontal slab
649         !                                             ! ===============
650         !
651#if defined key_qco   ||   defined key_linssh
652         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                 ! == reciprocal of e3 at F-point (key_qco)
653            z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
654         END_2D
655#else
656         SELECT CASE( nn_e3f_typ )           ! == reciprocal of e3 at F-point
657         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
658            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
659               ! round brackets added to fix the order of floating point operations
660               ! needed to ensure halo 1 - halo 2 compatibility
661               ze3f = (  (e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)    &
662                  &    +  e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk))   &
663                  &    + (e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)    &
664                  &    +  e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk))  )
665               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
666               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
667               ENDIF
668            END_2D
669         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
670            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
671               ! round brackets added to fix the order of floating point operations
672               ! needed to ensure halo 1 - halo 2 compatibility
673               ze3f = (  (e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)    &
674                  &    +  e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk))   &
675                  &    + (e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)    &
676                  &    +  e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk))  )
677               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
678                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
679               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
680               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
681               ENDIF
682            END_2D
683         END SELECT
684#endif
685         !
686         SELECT CASE( kvor )                 !==  vorticity considered  ==!
687         !
688         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
689            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
690               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
691            END_2D
692         CASE ( np_RVO )                           !* relative vorticity
693            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
694               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
695                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
696            END_2D
697            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
698               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
699                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
700               END_2D
701            ENDIF
702         CASE ( np_MET )                           !* metric term
703            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
704               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
705                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
706            END_2D
707         CASE ( np_CRV )                           !* Coriolis + relative vorticity
708            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
709            ! round brackets added to fix the order of floating point operations
710            ! needed to ensure halo 1 - halo 2 compatibility
711               zwz(ji,jj,jk) = (  ff_f(ji,jj) + ( ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
712                  &                               )                                                                  & ! bracket for halo 1 - halo 2 compatibility
713                  &                             - ( e1u(ji  ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)      &
714                  &                               )                                                                  & ! bracket for halo 1 - halo 2 compatibility
715                  &                             ) * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
716            END_2D
717            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
718               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
719                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
720               END_2D
721            ENDIF
722         CASE ( np_CME )                           !* Coriolis + metric
723            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
724               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
725                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
726            END_2D
727         CASE DEFAULT                                             ! error
728            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
729         END SELECT
730         !                                             ! ===============
731      END DO                                           !   End of slab
732      !                                                ! ===============
733      !
734      IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
735      !
736      !                                                ! ===============
737      !                                                ! Horizontal slab
738      !                                                ! ===============
739#if defined key_loop_fusion
740      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
741         !                                   !==  horizontal fluxes  ==!
742         zwx         = e2u(ji  ,jj  ) * e3u(ji  ,jj  ,jk,Kmm) * pu(ji  ,jj  ,jk)
743         zwx_im1     = e2u(ji-1,jj  ) * e3u(ji-1,jj  ,jk,Kmm) * pu(ji-1,jj  ,jk)
744         zwx_jp1     = e2u(ji  ,jj+1) * e3u(ji  ,jj+1,jk,Kmm) * pu(ji  ,jj+1,jk)
745         zwx_im1_jp1 = e2u(ji-1,jj+1) * e3u(ji-1,jj+1,jk,Kmm) * pu(ji-1,jj+1,jk)
746         zwy         = e1v(ji  ,jj  ) * e3v(ji  ,jj  ,jk,Kmm) * pv(ji  ,jj  ,jk)
747         zwy_ip1     = e1v(ji+1,jj  ) * e3v(ji+1,jj  ,jk,Kmm) * pv(ji+1,jj  ,jk)
748         zwy_jm1     = e1v(ji  ,jj-1) * e3v(ji  ,jj-1,jk,Kmm) * pv(ji  ,jj-1,jk)
749         zwy_ip1_jm1 = e1v(ji+1,jj-1) * e3v(ji+1,jj-1,jk,Kmm) * pv(ji+1,jj-1,jk)
750         !                                   !==  compute and add the vorticity term trend  =!
751         ztne     = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
752         ztnw     = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
753         ztnw_ip1 = zwz(ji  ,jj-1,jk) + zwz(ji  ,jj  ,jk) + zwz(ji+1,jj  ,jk)
754         ztse     = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
755         ztse_jp1 = zwz(ji  ,jj+1,jk) + zwz(ji  ,jj  ,jk) + zwz(ji-1,jj  ,jk)
756         ztsw_jp1 = zwz(ji  ,jj  ,jk) + zwz(ji-1,jj  ,jk) + zwz(ji-1,jj+1,jk)
757         ztsw_ip1 = zwz(ji+1,jj-1,jk) + zwz(ji  ,jj-1,jk) + zwz(ji  ,jj  ,jk)
758         !
759         zua = + r1_12 * r1_e1u(ji,jj) * (  ztne * zwy + ztnw_ip1 * zwy_ip1   &
760            &                             + ztse * zwy_jm1 + ztsw_ip1 * zwy_ip1_jm1 )
761         zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw_jp1 * zwx_im1_jp1 + ztse_jp1 * zwx_jp1   &
762            &                             + ztnw * zwx_im1 + ztne * zwx )
763         pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
764         pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
765      END_3D
766#else
767      DO jk = 1, jpkm1
768         !
769         !                                   !==  horizontal fluxes  ==!
770         DO_2D( 1, 1, 1, 1 )
771            zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
772            zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
773         END_2D
774         !
775         !                                   !==  compute and add the vorticity term trend  =!
776         DO_2D( 0, 1, 0, 1 )
777            ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
778            ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
779            ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
780            ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
781         END_2D
782         !
783         DO_2D( 0, 0, 0, 0 )
784            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
785               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
786            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
787               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
788            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
789            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
790         END_2D
791      END DO
792#endif
793         !                                             ! ===============
794         !                                             !   End of slab
795      !                                                ! ===============
796   END SUBROUTINE vor_een
797
798
799   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
800      !!----------------------------------------------------------------------
801      !!                ***  ROUTINE vor_eeT  ***
802      !!
803      !! ** Purpose :   Compute the now total vorticity trend and add it to
804      !!      the general trend of the momentum equation.
805      !!
806      !! ** Method  :   Trend evaluated using now fields (centered in time)
807      !!      and the Arakawa and Lamb (1980) vector form formulation using
808      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
809      !!      The change consists in
810      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
811      !!
812      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
813      !!
814      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
815      !!----------------------------------------------------------------------
816      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
817      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
818      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
819      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
820      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
821      !
822      INTEGER  ::   ji, jj, jk     ! dummy loop indices
823      INTEGER  ::   ierr           ! local integer
824      REAL(wp) ::   zua, zva       ! local scalars
825      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
826      REAL(wp), DIMENSION(A2D(nn_hls))       ::   zwx , zwy
827      REAL(wp), DIMENSION(A2D(nn_hls))       ::   ztnw, ztne, ztsw, ztse
828      REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined
829      !!----------------------------------------------------------------------
830      !
831      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
832         IF( kt == nit000 ) THEN
833            IF(lwp) WRITE(numout,*)
834            IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
835            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
836         ENDIF
837      ENDIF
838      !
839      !                                                ! ===============
840      DO jk = 1, jpkm1                                 ! Horizontal slab
841         !                                             ! ===============
842         !
843         !
844         SELECT CASE( kvor )                 !==  vorticity considered  ==!
845         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
846            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
847               zwz(ji,jj,jk) = ff_f(ji,jj)
848            END_2D
849         CASE ( np_RVO )                           !* relative vorticity
850            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
851               ! round brackets added to fix the order of floating point operations
852               ! needed to ensure halo 1 - halo 2 compatibility
853               zwz(ji,jj,jk) = (  (e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk))    &
854                  &             - (e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk))  ) &
855                  &          * r1_e1e2f(ji,jj)
856            END_2D
857            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
858               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
859                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
860               END_2D
861            ENDIF
862         CASE ( np_MET )                           !* metric term
863            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
864               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
865                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
866            END_2D
867         CASE ( np_CRV )                           !* Coriolis + relative vorticity
868            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
869               ! round brackets added to fix the order of floating point operations
870               ! needed to ensure halo 1 - halo 2 compatibility
871               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  (e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk))    &
872                  &                              - (e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk))  ) &
873                  &                           * r1_e1e2f(ji,jj)    )
874            END_2D
875            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
876               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
877                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
878               END_2D
879            ENDIF
880         CASE ( np_CME )                           !* Coriolis + metric
881            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
882               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
883                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
884            END_2D
885         CASE DEFAULT                                             ! error
886            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
887         END SELECT
888         !
889         !                                             ! ===============
890      END DO                                           !   End of slab
891      !                                                ! ===============
892      !
893      IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
894      !
895      !                                                ! ===============
896      DO jk = 1, jpkm1                                 ! Horizontal slab
897         !                                             ! ===============
898         !
899         !                                   !==  horizontal fluxes  ==!
900         DO_2D( 1, 1, 1, 1 )
901            zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
902            zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
903         END_2D
904         !
905         !                                   !==  compute and add the vorticity term trend  =!
906         DO_2D( 0, 1, 0, 1 )
907            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
908            ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
909            ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
910            ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
911            ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
912         END_2D
913         !
914         DO_2D( 0, 0, 0, 0 )
915            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
916               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
917            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
918               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
919            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
920            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
921         END_2D
922         !                                             ! ===============
923      END DO                                           !   End of slab
924      !                                                ! ===============
925   END SUBROUTINE vor_eeT
926
927
928   SUBROUTINE dyn_vor_init
929      !!---------------------------------------------------------------------
930      !!                  ***  ROUTINE dyn_vor_init  ***
931      !!
932      !! ** Purpose :   Control the consistency between cpp options for
933      !!              tracer advection schemes
934      !!----------------------------------------------------------------------
935      INTEGER ::   ji, jj, jk    ! dummy loop indices
936      INTEGER ::   ioptio, ios   ! local integer
937      REAL(wp) ::   zmsk    ! local scalars
938      !!
939      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
940         &                 ln_dynvor_een, nn_e3f_typ   , ln_dynvor_mix, ln_dynvor_msk
941      !!----------------------------------------------------------------------
942      !
943      IF(lwp) THEN
944         WRITE(numout,*)
945         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
946         WRITE(numout,*) '~~~~~~~~~~~~'
947      ENDIF
948      !
949      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
950901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
951      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
952902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
953      IF(lwm) WRITE ( numond, namdyn_vor )
954      !
955      IF(lwp) THEN                    ! Namelist print
956         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
957         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
958         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
959         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
960         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
961         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
962         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_e3f_typ = ', nn_e3f_typ
963         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
964         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
965      ENDIF
966
967!!gm  this should be removed when choosing a unique strategy for fmask at the coast
968      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
969      ! at angles with three ocean points and one land point
970      IF(lwp) WRITE(numout,*)
971      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
972      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
973         DO_3D( 1, 0, 1, 0, 1, jpk )
974            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
975               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
976         END_3D
977         !
978         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
979         !
980      ENDIF
981!!gm end
982
983      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
984      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
985      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
986      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
987      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
988      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
989      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
990      !
991      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
992      !
993      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
994      ncor = np_COR                       ! planetary vorticity
995      SELECT CASE( n_dynadv )
996      CASE( np_LIN_dyn )
997         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
998         nrvm = np_COR        ! planetary vorticity
999         ntot = np_COR        !     -         -
1000      CASE( np_VEC_c2  )
1001         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'
1002         nrvm = np_RVO        ! relative vorticity
1003         ntot = np_CRV        ! relative + planetary vorticity
1004      CASE( np_FLX_c2 , np_FLX_ubs  )
1005         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
1006         nrvm = np_MET        ! metric term
1007         ntot = np_CME        ! Coriolis + metric term
1008         !
1009         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
1010         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
1011            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
1012            DO_2D( 0, 0, 0, 0 )
1013               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
1014               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
1015            END_2D
1016            CALL lbc_lnk( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
1017            !
1018         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
1019            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
1020            DO_2D( 1, 0, 1, 0 )
1021               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
1022               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
1023            END_2D
1024            CALL lbc_lnk( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
1025         END SELECT
1026         !
1027      END SELECT
1028#if defined key_qco   ||   defined key_linssh
1029      SELECT CASE( nvor_scheme )    ! qco or linssh cases : pre-computed a specific e3f_0 for some vorticity schemes
1030      CASE( np_ENS , np_ENE , np_EEN , np_MIX )
1031         !
1032         ALLOCATE( e3f_0vor(jpi,jpj,jpk) )
1033         !
1034         SELECT CASE( nn_e3f_typ )
1035         CASE ( 0 )                        ! original formulation  (masked averaging of e3t divided by 4)
1036            DO_3D( 0, 0, 0, 0, 1, jpk )
1037               e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
1038                  &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
1039                  &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
1040                  &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) * 0.25_wp
1041            END_3D
1042         CASE ( 1 )                        ! new formulation  (masked averaging of e3t divided by the sum of mask)
1043            DO_3D( 0, 0, 0, 0, 1, jpk )
1044               zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk)   &
1045                  &  + tmask(ji,jj  ,jk) +tmask(ji+1,jj  ,jk)  )
1046               !
1047               IF( zmsk /= 0._wp ) THEN
1048                  e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
1049                     &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
1050                     &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
1051                     &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) / zmsk
1052               ELSE ; e3f_0vor(ji,jj,jk) = 0._wp
1053               ENDIF
1054            END_3D
1055         END SELECT
1056         !
1057         CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp )
1058         !                                 ! insure e3f_0vor /= 0
1059         WHERE( e3f_0vor(:,:,:) == 0._wp )   e3f_0vor(:,:,:) = e3f_0(:,:,:)
1060         !
1061      END SELECT
1062      !
1063#endif
1064      IF(lwp) THEN                   ! Print the choice
1065         WRITE(numout,*)
1066         SELECT CASE( nvor_scheme )
1067         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
1068         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
1069         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
1070                              IF( ln_dynadv_vec )   CALL ctl_warn('dyn_vor_init: ENT scheme may not work in vector form')
1071         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
1072         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
1073         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
1074         END SELECT
1075      ENDIF
1076      !
1077   END SUBROUTINE dyn_vor_init
1078
1079   !!==============================================================================
1080END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.