New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynvor.F90 in NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90 @ 13005

Last change on this file since 13005 was 13005, checked in by gm, 4 years ago
ADE and more options to AM98 config
Property svn:keywords set to `Id`
File size: 71.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: