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

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90 @ 15574

Last change on this file since 15574 was 15574, checked in by techene, 3 years ago

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