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/cfgs/penzAM98/MY_SRC – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC/dynvor.F90 @ 13562

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

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

Download in other formats: