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 @ 14801

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

add loop fusion to DYN and TRA modules - ticket #2607

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