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_tam.F90 in branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

Context Navigation

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/DYN/dynvor_tam.F90 @ 4580

Last change on this file since 4580 was 4580, checked in by pabouttier, 10 years ago
Correct a wrong variable name in dynvor_tam (ze3f instead of ze3f_ad)
Property svn:executable set to ``*
File size: 83.3 KB

Line
1	MODULE dynvor_tam
2	#if defined key_tam
3	!!======================================================================
4	!! * MODULE dynvor_tam *
5	!! Ocean dynamics: Update the momentum trend with the relative and
6	!! planetary vorticity trends
7	!! Tangent and Adjoint Module
8	!!======================================================================
9	!! History of the drect module:
10	!! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code
11	!! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code
12	!! 6.0 ! 1996-01 (G. Madec) s-coord, suppress work arrays
13	!! 8.5 ! 2002-08 (G. Madec) F90: Free form and module
14	!! NEMO 1.0 ! 2004-02 (G. Madec) vor_een: Original code
15	!! - ! 2003-08 (G. Madec) add vor_ctl
16	!! - ! 2005-11 (G. Madec) add dyn_vor (new step architecture)
17	!! 2.0 ! 2006-11 (G. Madec) flux form advection: add metric term
18	!! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme
19	!! History of the TAM module:
20	!! 9.0 ! 2008-06 (A. Vidard) Skeleton
21	!! 9.0 ! 2009-01 (A. Vidard) TAM of the 06-11 version
22	!! 9.0 ! 2010-01 (F. Vigilant) Add een TAM option
23	!! NEMO 3.2 ! 2010-04 (F. Vigilant) 3.2 version
24	!! NEMO 3.4 ! 2012-07 (P.-A. Bouttier) 3.4 version
25	!!----------------------------------------------------------------------
26
27	!!----------------------------------------------------------------------
28	!! dyn_vor : Update the momentum trend with the vorticity trend
29	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
30	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
31	!! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
32	!! vor_een : energy and enstrophy conserving (ln_dynvor_een=T)
33	!! vor_ctl : set and control of the different vorticity option
34	!!----------------------------------------------------------------------
35	USE par_oce
36	USE oce
37	USE oce_tam
38	USE divcur
39	USE divcur_tam
40	USE dom_oce
41	USE dynadv
42	USE dynvor
43	USE lbclnk
44	USE lbclnk_tam
45	USE in_out_manager
46	USE gridrandom
47	USE dotprodfld
48	USE tstool_tam
49	USE dommsk
50	USE lib_mpp
51	USE wrk_nemo
52	USE timing
53
54	IMPLICIT NONE
55	PRIVATE
56
57	PUBLIC dyn_vor_init_tam
58	PUBLIC dyn_vor_tan ! routine called by step_tam.F90
59	PUBLIC dyn_vor_adj ! routine called by step_tam.F90
60	PUBLIC dyn_vor_adj_tst ! routine called by the tst.F90
61
62	INTEGER :: nvor = 0 ! type of vorticity trend used
63	INTEGER :: ncor = 1 ! coriolis
64	INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term
65	INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
66
67	!! * Substitutions
68	# include "domzgr_substitute.h90"
69	# include "vectopt_loop_substitute.h90"
70
71	CONTAINS
72
73	SUBROUTINE dyn_vor_tan( kt )
74	!!----------------------------------------------------------------------
75	!!
76	!! ** Purpose of the direct routine:
77	!! compute the lateral ocean tracer physics.
78	!!
79	!! ** Action : - Update (ua,va) with the now vorticity term trend
80	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
81	!! and planetary vorticity trends) ('key_trddyn')
82	!!----------------------------------------------------------------------
83	!!
84	INTEGER, INTENT( in ) :: kt ! ocean time-step index
85	!!----------------------------------------------------------------------
86	!
87	IF( nn_timing == 1 ) CALL timing_start('dyn_vor_tan')
88	!
89	!! ! vorticity term
90	SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend
91	!
92	CASE ( -1 ) ! esopa: test all possibility with control print
93	CALL vor_ene_tan( kt, ntot, ua_tl, va_tl )
94	CALL vor_ens_tan( kt, ntot, ua_tl, va_tl )
95	CALL vor_mix_tan( kt )
96	CALL vor_een_tan( kt, ntot, ua_tl, va_tl )
97	!
98	CASE ( 0 ) ! energy conserving scheme
99	CALL vor_ene_tan( kt, ntot, ua_tl, va_tl ) ! total vorticity
100	!
101	CASE ( 1 ) ! enstrophy conserving scheme
102	CALL vor_ens_tan( kt, ntot, ua_tl, va_tl ) ! total vorticity
103	!
104	CASE ( 2 ) ! mixed ene-ens scheme
105	CALL vor_mix_tan( kt ) ! total vorticity (mix=ens-ene)
106	!
107	CASE ( 3 ) ! energy and enstrophy conserving scheme
108	CALL vor_een_tan( kt, ntot, ua_tl, va_tl ) ! total vorticity
109	!
110	END SELECT
111	!
112	IF( nn_timing == 1 ) CALL timing_stop('dyn_vor_tan')
113	!
114	END SUBROUTINE dyn_vor_tan
115	SUBROUTINE vor_ene_tan( kt, kvor, pua_tl, pva_tl )
116	!!----------------------------------------------------------------------
117	!! * ROUTINE vor_ene *
118	!!
119	!! ** Purpose : Compute the now total vorticity trend and add it to
120	!! the general trend of the momentum equation.
121	!!
122	!! ** Method : Trend evaluated using now fields (centered in time)
123	!! and the Sadourny (1975) flux form formulation : conserves the
124	!! horizontal kinetic energy.
125	!! The trend of the vorticity term is given by:
126	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
127	!! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ]
128	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ]
129	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
130	!! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ]
131	!! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ]
132	!! Add this trend to the general momentum trend (ua,va):
133	!! (ua,va) = (ua,va) + ( voru , vorv )
134	!!
135	!! ** Action : - Update (ua,va) with the now vorticity term trend
136	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
137	!! and planetary vorticity trends) ('key_trddyn')
138	!!
139	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
140	!!----------------------------------------------------------------------
141	INTEGER , INTENT(in ) :: kt ! ocean time-step index
142	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
143	! ! =nrvm (relative vorticity or metric)
144	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_tl ! total u-trend
145	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_tl ! total v-trend
146	!!
147	INTEGER :: ji, jj, jk ! dummy loop indices
148	REAL(wp) :: zx1, zy1, zfact2 ! temporary scalars
149	REAL(wp) :: zx1tl, zy1tl ! temporary scalars
150	REAL(wp) :: zx2, zy2 ! " "
151	REAL(wp) :: zx2tl, zy2tl ! " "
152	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! temporary 2D workspace
153	REAL(wp), POINTER, DIMENSION(:,:) :: zwxtl, zwytl, zwztl ! temporary 2D workspace
154	!!----------------------------------------------------------------------
155	!
156	IF( nn_timing == 1 ) CALL timing_start('vor_ene_tan')
157	!
158	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
159	CALL wrk_alloc( jpi, jpj, zwxtl, zwytl, zwztl )
160	!
161	IF( kt == nit000 ) THEN
162	IF(lwp) WRITE(numout,*)
163	IF(lwp) WRITE(numout,*) 'dyn:vor_ene_tan : vorticity term: energy conserving scheme'
164	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
165	ENDIF
166
167	zfact2 = 0.5 * 0.5 ! Local constant initialization
168
169	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
170	! ! ===============
171	DO jk = 1, jpkm1 ! Horizontal slab
172	! ! ===============
173	!
174	! Potential vorticity and horizontal fluxes
175	! -----------------------------------------
176	SELECT CASE( kvor ) ! vorticity considered
177	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ; zwztl(:,:) = 0.
178	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ; zwztl(:,:) = rotn_tl(:,:,jk) ! relative vorticity
179	CASE ( 3 ) ! metric term
180	DO jj = 1, jpjm1
181	DO ji = 1, fs_jpim1 ! vector opt.
182	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
183	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
184	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
185	zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
186	& - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
187	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
188	END DO
189	END DO
190	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ; zwztl(:,:) = rotn_tl(:,:,jk) ! total (relative + planetary vorticity)
191	CASE ( 5 ) ! total (coriolis + metric)
192	DO jj = 1, jpjm1
193	DO ji = 1, fs_jpim1 ! vector opt.
194	zwz(ji,jj) = ( ff (ji,jj) &
195	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
196	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
197	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
198	& )
199	zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
200	& - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
201	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
202	END DO
203	END DO
204	END SELECT
205
206	IF( ln_sco ) THEN
207	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
208	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
209	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
210	zwztl(:,:) = zwztl(:,:) / fse3f(:,:,jk)
211	zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk)
212	zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk)
213	ELSE
214	zwx(:,:) = e2u(:,:) * un(:,:,jk)
215	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
216	zwxtl(:,:) = e2u(:,:) * un_tl(:,:,jk)
217	zwytl(:,:) = e1v(:,:) * vn_tl(:,:,jk)
218	ENDIF
219
220	! Compute and add the vorticity term trend
221	! ----------------------------------------
222	DO jj = 2, jpjm1
223	DO ji = fs_2, fs_jpim1 ! vector opt.
224	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
225	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
226	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
227	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
228	zy1tl = zwytl(ji,jj-1) + zwytl(ji+1,jj-1)
229	zy2tl = zwytl(ji,jj ) + zwytl(ji+1,jj )
230	zx1tl = zwxtl(ji-1,jj) + zwxtl(ji-1,jj+1)
231	zx2tl = zwxtl(ji ,jj) + zwxtl(ji ,jj+1)
232	pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwztl(ji ,jj-1) * zy1 + zwztl(ji,jj) * zy2 ) &
233	& + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1tl + zwz(ji,jj) * zy2tl )
234	pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwztl(ji-1,jj ) * zx1 + zwztl(ji,jj) * zx2 ) &
235	& - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1tl + zwz(ji,jj) * zx2tl )
236	END DO
237	END DO
238	! ! ===============
239	END DO ! End of slab
240	! ! ===============
241	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
242	CALL wrk_dealloc( jpi, jpj, zwxtl, zwytl, zwztl )
243	!
244	IF( nn_timing == 1 ) CALL timing_stop('vor_ene_tan')
245	!
246	END SUBROUTINE vor_ene_tan
247
248
249	SUBROUTINE vor_mix_tan( kt )
250	!!----------------------------------------------------------------------
251	!! * ROUTINE vor_mix *
252	!!
253	!! ** Purpose : Compute the now total vorticity trend and add it to
254	!! the general trend of the momentum equation.
255	!!
256	!! ** Method : Trend evaluated using now fields (centered in time)
257	!! Mixte formulation : conserves the potential enstrophy of a hori-
258	!! zontally non-divergent flow for (rotzu x uh), the relative vor-
259	!! ticity term and the horizontal kinetic energy for (f x uh), the
260	!! coriolis term. the now trend of the vorticity term is given by:
261	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
262	!! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
263	!! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]
264	!! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
265	!! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]
266	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
267	!! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]
268	!! +1/e1u mj-1[ f mi(e1v vn) ]
269	!! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]
270	!! +1/e2v mi-1[ f mj(e2u un) ]
271	!! Add this now trend to the general momentum trend (ua,va):
272	!! (ua,va) = (ua,va) + ( voru , vorv )
273	!!
274	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
275	!! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
276	!! and planetary vorticity trends) ('key_trddyn')
277	!!
278	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
279	!!----------------------------------------------------------------------
280	INTEGER, INTENT(in) :: kt ! ocean timestep index
281	!!
282	INTEGER :: ji, jj, jk ! dummy loop indices
283	REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! temporary scalars
284	REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! " "
285	REAL(wp) :: zuatl, zcuatl, zx1tl, zy1tl ! temporary scalars
286	REAL(wp) :: zvatl, zcvatl, zx2tl, zy2tl ! " "
287	REAL(wp), POINTER, DIMENSION(:,:) :: zwxtl, zwytl, zwztl, zwwtl ! temporary 3D workspace
288	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww ! temporary 3D workspace
289	!!----------------------------------------------------------------------
290	!
291	IF( nn_timing == 1 ) CALL timing_start('vor_mix_tan')
292	!
293	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )
294	CALL wrk_alloc( jpi, jpj, zwxtl, zwytl, zwztl, zwwtl )
295	!
296	IF( kt == nit000 ) THEN
297	IF(lwp) WRITE(numout,*)
298	IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
299	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
300	ENDIF
301
302	zfact1 = 0.5 * 0.25 ! Local constant initialization
303	zfact2 = 0.5 * 0.5
304
305	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
306	! ! ===============
307	DO jk = 1, jpkm1 ! Horizontal slab
308	! ! ===============
309	!
310	! Relative and planetary potential vorticity and horizontal fluxes
311	! ----------------------------------------------------------------
312	IF( ln_sco ) THEN
313	IF( ln_dynadv_vec ) THEN
314	zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
315	zwwtl(:,:) = rotn_tl(:,:,jk) / fse3f(:,:,jk)
316	ELSE
317	DO jj = 1, jpjm1
318	DO ji = 1, fs_jpim1 ! vector opt.
319	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
320	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
321	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
322	zwwtl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
323	& - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
324	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
325	END DO
326	END DO
327	ENDIF
328	zwz(:,:) = ff (:,:) / fse3f(:,:,jk)
329	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
330	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
331	zwztl(:,:) = 0._wp
332	zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk)
333	zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk)
334	ELSE
335	IF( ln_dynadv_vec ) THEN
336	zww(:,:) = rotn(:,:,jk)
337	zwwtl(:,:) = rotn_tl(:,:,jk)
338	ELSE
339	DO jj = 1, jpjm1
340	DO ji = 1, fs_jpim1 ! vector opt.
341	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
342	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
343	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
344	zwwtl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
345	& - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
346	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
347	END DO
348	END DO
349	ENDIF
350	zwz(:,:) = ff (:,:)
351	zwx(:,:) = e2u(:,:) * un(:,:,jk)
352	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
353	zwztl(:,:) = 0.0_wp
354	zwxtl(:,:) = e2u(:,:) * un_tl(:,:,jk)
355	zwytl(:,:) = e1v(:,:) * vn_tl(:,:,jk)
356	ENDIF
357
358	! Compute and add the vorticity term trend
359	! ----------------------------------------
360	DO jj = 2, jpjm1
361	DO ji = fs_2, fs_jpim1 ! vector opt.
362	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
363	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
364	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
365	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
366	zy1tl = ( zwytl(ji,jj-1) + zwytl(ji+1,jj-1) ) / e1u(ji,jj)
367	zy2tl = ( zwytl(ji,jj ) + zwytl(ji+1,jj ) ) / e1u(ji,jj)
368	zx1tl = ( zwxtl(ji-1,jj) + zwxtl(ji-1,jj+1) ) / e2v(ji,jj)
369	zx2tl = ( zwxtl(ji ,jj) + zwxtl(ji ,jj+1) ) / e2v(ji,jj)
370	! enstrophy conserving formulation for relative vorticity term
371	zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
372	zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )
373	zuatl = zfact1 * ( zwwtl(ji ,jj-1) + zwwtl(ji,jj) ) * ( zy1 + zy2 ) &
374	& + zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1tl + zy2tl )
375	zvatl =-zfact1 * ( zwwtl(ji-1,jj ) + zwwtl(ji,jj) ) * ( zx1 + zx2 ) &
376	& - zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1tl + zx2tl )
377	! energy conserving formulation for planetary vorticity term
378	zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
379	zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
380	zcuatl = zfact2 * ( zwztl(ji ,jj-1) * zy1 + zwztl(ji,jj) * zy2 ) &
381	& + zfact2 * ( zwz(ji ,jj-1) * zy1tl + zwz(ji,jj) * zy2tl )
382	zcvatl =-zfact2 * ( zwztl(ji-1,jj ) * zx1 + zwztl(ji,jj) * zx2 ) &
383	& -zfact2 * ( zwz(ji-1,jj ) * zx1tl + zwz(ji,jj) * zx2tl )
384	! mixed vorticity trend added to the momentum trends
385	ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + zcuatl + zuatl
386	va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + zcvatl + zvatl
387	END DO
388	END DO
389	! ! ===============
390	END DO ! End of slab
391	! ! ===============
392	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )
393	CALL wrk_dealloc( jpi, jpj, zwxtl, zwytl, zwztl, zwwtl )
394	!
395	IF( nn_timing == 1 ) CALL timing_stop('vor_mix_tan')
396	!
397	END SUBROUTINE vor_mix_tan
398	SUBROUTINE vor_ens_tan( kt, kvor, pua_tl, pva_tl )
399	!!----------------------------------------------------------------------
400	!! * ROUTINE vor_ens_tan *
401	!!
402	!! ** Purpose of the direct routine:
403	!! Compute the now total vorticity trend and add it to
404	!! the general trend of the momentum equation.
405	!!
406	!! ** Method of the direct routine:
407	!! Trend evaluated using now fields (centered in time)
408	!! and the Sadourny (1975) flux FORM formulation : conserves the
409	!! potential enstrophy of a horizontally non-divergent flow. the
410	!! trend of the vorticity term is given by:
411	!! * s-coordinate (ln_sco=T), the e3. are inside the derivative:
412	!! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
413	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
414	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
415	!! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ]
416	!! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ]
417	!! Add this trend to the general momentum trend (ua,va):
418	!! (ua,va) = (ua,va) + ( voru , vorv )
419	!!
420	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
421	!! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
422	!! and planetary vorticity trends) ('key_trddyn')
423	!!
424	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
425	!!----------------------------------------------------------------------
426	INTEGER , INTENT(in ) :: kt ! ocean time-step index
427	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
428	! ! =nrvm (relative vorticity or metric)
429	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_tl ! total u-trend
430	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_tl ! total v-trend
431	!!
432	INTEGER :: ji, jj, jk ! dummy loop indices
433	REAL(wp) :: zfact1, zuav, zvau ! temporary scalars
434	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! temporary 3D workspace
435	REAL(wp) :: zuavtl, zvautl ! temporary scalars
436	REAL(wp), POINTER, DIMENSION(:,:) :: zwxtl, zwytl, zwztl ! temporary 3D workspace
437	!!----------------------------------------------------------------------
438	!
439	IF( nn_timing == 1 ) CALL timing_start('vor_ens_tan')
440	!
441	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
442	CALL wrk_alloc( jpi, jpj, zwxtl, zwytl, zwztl )
443	!
444	IF( kt == nit000 ) THEN
445	IF(lwp) WRITE(numout,*)
446	IF(lwp) WRITE(numout,*) 'dyn_vor_ens_tan : vorticity term: enstrophy conserving scheme'
447	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
448	ENDIF
449	! Local constant initialization
450	zfact1 = 0.5 * 0.25
451
452	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz zwxtl, zwytl, zwztl)
453	! ! ===============
454	DO jk = 1, jpkm1 ! Horizontal slab
455	! ! ===============
456	! Potential vorticity and horizontal fluxes
457	! -----------------------------------------
458	SELECT CASE( kvor ) ! vorticity considered
459	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
460	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
461	CASE ( 3 ) ! metric term
462	DO jj = 1, jpjm1
463	DO ji = 1, fs_jpim1 ! vector opt.
464	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
465	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
466	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
467	END DO
468	END DO
469	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
470	CASE ( 5 ) ! total (coriolis + metric)
471	DO jj = 1, jpjm1
472	DO ji = 1, fs_jpim1 ! vector opt.
473	zwz(ji,jj) = ( ff (ji,jj) &
474	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
475	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
476	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
477	& )
478	END DO
479	END DO
480	END SELECT
481
482	IF( ln_sco ) THEN
483	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
484	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
485	zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
486	zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
487	zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
488	END DO
489	END DO
490	ELSE
491	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
492	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
493	zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
494	zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
495	END DO
496	END DO
497	ENDIF
498
499	! ===================
500	! Tangent counterpart
501	! ===================
502	! Potential vorticity and horizontal fluxes
503	! -----------------------------------------
504	SELECT CASE( kvor ) ! vorticity considered
505	CASE ( 1 ) ; zwztl(:,:) = 0.0_wp ! planetary vorticity (Coriolis)
506	CASE ( 2 ,4) ; zwztl(:,:) = rotn_tl(:,:,jk) ! relative vorticity
507	CASE ( 3 ,5 ) ! metric term
508	DO jj = 1, jpjm1
509	DO ji = 1, fs_jpim1 ! vector opt.
510	zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) &
511	& * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
512	& - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) &
513	& * ( e1u(ji ,jj+1) - e1u(ji,jj) ) &
514	& ) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
515	END DO
516	END DO
517	END SELECT
518
519	IF( ln_sco ) THEN
520	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
521	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
522	zwztl(ji,jj) = zwztl(ji,jj) / fse3f(ji,jj,jk)
523	zwxtl(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un_tl(ji,jj,jk)
524	zwytl(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn_tl(ji,jj,jk)
525	END DO
526	END DO
527	ELSE
528	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
529	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
530	zwxtl(ji,jj) = e2u(ji,jj) * un_tl(ji,jj,jk)
531	zwytl(ji,jj) = e1v(ji,jj) * vn_tl(ji,jj,jk)
532	END DO
533	END DO
534	ENDIF
535
536	! Compute and add the vorticity term trend
537	! ----------------------------------------
538	DO jj = 2, jpjm1
539	DO ji = fs_2, fs_jpim1 ! vector opt.
540	zuav = zfact1 / e1u(ji,jj) * ( zwy( ji ,jj-1) + zwy( ji+1,jj-1) &
541	& + zwy( ji ,jj ) + zwy( ji+1,jj ) )
542	zvau =-zfact1 / e2v(ji,jj) * ( zwx( ji-1,jj ) + zwx( ji-1,jj+1) &
543	& + zwx( ji ,jj ) + zwx( ji ,jj+1) )
544	zuavtl = zfact1 / e1u(ji,jj) * ( zwytl(ji ,jj-1) + zwytl(ji+1,jj-1) &
545	& + zwytl(ji ,jj ) + zwytl(ji+1,jj ) )
546	zvautl =-zfact1 / e2v(ji,jj) * ( zwxtl(ji-1,jj ) + zwxtl(ji-1,jj+1) &
547	& + zwxtl(ji ,jj ) + zwxtl(ji ,jj+1) )
548	pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) &
549	& + zuavtl * ( zwz( ji,jj-1) + zwz( ji,jj) ) &
550	& + zuav * ( zwztl(ji,jj-1) + zwztl(ji,jj) )
551	pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) &
552	& + zvautl * ( zwz( ji-1,jj) + zwz( ji,jj) ) &
553	& + zvau * ( zwztl(ji-1,jj) + zwztl(ji,jj) )
554	END DO
555	END DO
556	! ! ===============
557	END DO ! End of slab
558	! ! ===============
559	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
560	CALL wrk_dealloc( jpi, jpj, zwxtl, zwytl, zwztl )
561	!
562	IF( nn_timing == 1 ) CALL timing_stop('vor_ens_tan')
563	!
564	END SUBROUTINE vor_ens_tan
565
566	SUBROUTINE vor_een_tan( kt, kvor, pua_tl, pva_tl )
567	!!----------------------------------------------------------------------
568	!! * ROUTINE vor_een_tan *
569	!!
570	!! ** Purpose : Compute the now total vorticity trend and add it to
571	!! the general trend of the momentum equation.
572	!!
573	!! ** Method : Trend evaluated using now fields (centered in time)
574	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
575	!! both the horizontal kinetic energy and the potential enstrophy
576	!! when horizontal divergence is zero (see the NEMO documentation)
577	!! Add this trend to the general momentum trend (ua,va).
578	!!
579	!! ** Action : - Update (ua,va) with the now vorticity term trend
580	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
581	!! and planetary vorticity trends) ('key_trddyn')
582	!!
583	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
584	!!----------------------------------------------------------------------
585	INTEGER , INTENT(in ) :: kt ! ocean time-step index
586	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
587	! ! =nrvm (relative vorticity or metric)
588	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_tl ! total u-trend
589	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_tl ! total v-trend
590	!!
591	INTEGER :: ji, jj, jk ! dummy loop indices
592	INTEGER :: ierr
593	REAL(wp) :: zfac12, zua, zva ! temporary scalars
594	REAL(wp) :: zuatl, zvatl ! temporary scalars
595	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! temporary 2D workspace
596	REAL(wp), POINTER, DIMENSION(:,:) :: ztnw, ztne, ztsw, ztse ! temporary 3D workspace
597	REAL(wp), POINTER, DIMENSION(:,:) :: zwxtl, zwytl, zwztl ! temporary 2D workspace
598	REAL(wp), POINTER, DIMENSION(:,:) :: ztnwtl, ztnetl, ztswtl, ztsetl ! temporary 3D workspace
599	#if defined key_vvl
600	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3f_tl
601	#else
602	REAL(wp), POINTER, DIMENSION(:,:,:), SAVE :: ze3f_tl
603	#endif
604	!!----------------------------------------------------------------------
605	IF( nn_timing == 1 ) CALL timing_start('vor_een_tan')
606	!
607	CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz )
608	CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
609	CALL wrk_alloc( jpi, jpj, zwxtl , zwytl , zwztl )
610	CALL wrk_alloc( jpi, jpj, ztnwtl, ztnetl, ztswtl, ztsetl )
611
612	zuatl = 0.0_wp ; zvatl = 0.0_wp
613	zwx (:,:) = 0.0_wp ; zwy (:,:) = 0.0_wp ; zwz (:,:) = 0.0_wp
614	ztnw(:,:) = 0.0_wp ; ztne(:,:) = 0.0_wp ; ztsw(:,:) = 0.0_wp ; ztse(:,:) = 0.0_wp
615	zwxtl (:,:) = 0.0_wp ; zwytl (:,:) = 0.0_wp ; zwztl (:,:) = 0.0_wp
616	ztnwtl(:,:) = 0.0_wp ; ztnetl(:,:) = 0.0_wp ; ztswtl(:,:) = 0.0_wp ; ztsetl(:,:) = 0.0_wp
617	#if defined key_vvl
618	CALL ctl_stop( 'STOP', 'key_vvl not available for tam' )
619	#endif
620	IF( kt == nit000 ) THEN
621	IF(lwp) WRITE(numout,*)
622	IF(lwp) WRITE(numout,*) 'dyn:vor_een_tam : vorticity term: energy and enstrophy conserving scheme'
623	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
624	IF( .NOT.lk_vvl ) THEN
625	ALLOCATE( ze3f_tl(jpi,jpj,jpk) , STAT=ierr )
626	IF( lk_mpp ) CALL mpp_sum ( ierr )
627	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
628	ze3f_tl = 0._wp
629	ENDIF
630	ENDIF
631
632	IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t)
633	DO jk = 1, jpk
634	DO jj = 1, jpjm1
635	DO ji = 1, jpim1
636	ze3f_tl(ji,jj,jk) = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
637	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) ) * 0.25_wp
638	IF( ze3f_tl(ji,jj,jk) /= 0.0_wp ) ze3f_tl(ji,jj,jk) = 1.0_wp / ze3f_tl(ji,jj,jk)
639	END DO
640	END DO
641	END DO
642	CALL lbc_lnk( ze3f_tl, 'F', 1._wp )
643	ENDIF
644
645	zfac12 = 1.0_wp / 12.0_wp ! Local constant initialization
646
647
648	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
649	! ! ===============
650	DO jk = 1, jpkm1 ! Horizontal slab
651	! ! ===============
652
653	! Potential vorticity and horizontal fluxes
654	! -----------------------------------------
655	SELECT CASE( kvor ) ! vorticity considered
656	CASE ( 1 )
657	zwz(:,:) = ff(:,:) * ze3f_tl(:,:,jk) ! planetary vorticity (Coriolis)
658	zwztl(:,:) = 0.0_wp
659	CASE ( 2 )
660	zwz(:,:) = rotn(:,:,jk) * ze3f_tl(:,:,jk) ! relative vorticity
661	zwztl(:,:) = rotn_tl(:,:,jk) * ze3f_tl(:,:,jk)
662	CASE ( 3 ) ! metric term
663	DO jj = 1, jpjm1
664	DO ji = 1, fs_jpim1 ! vector opt.
665	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
666	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
667	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_tl(ji,jj,jk)
668	END DO
669	END DO
670	CALL lbc_lnk( zwz, 'F', 1._wp )
671	DO jj = 1, jpjm1
672	DO ji = 1, fs_jpim1 ! vector opt.
673	zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
674	& - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
675	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_tl(ji,jj,jk)
676	END DO
677	END DO
678	CALL lbc_lnk( zwztl, 'F', 1._wp )
679	CASE ( 4 )
680	zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f_tl(:,:,jk) ! total (relative + planetary vorticity)
681	zwztl(:,:) = ( rotn_tl(:,:,jk) ) * ze3f_tl(:,:,jk)
682	CASE ( 5 ) ! total (coriolis + metric)
683	DO jj = 1, jpjm1
684	DO ji = 1, fs_jpim1 ! vector opt.
685	zwz(ji,jj) = ( ff (ji,jj) &
686	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
687	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
688	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
689	& ) * ze3f_tl(ji,jj,jk)
690	END DO
691	END DO
692	CALL lbc_lnk( zwz, 'F', 1._wp )
693	DO jj = 1, jpjm1
694	DO ji = 1, fs_jpim1 ! vector opt.
695	zwztl(ji,jj) = ( ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
696	& - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
697	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
698	& ) * ze3f_tl(ji,jj,jk)
699	END DO
700	END DO
701	CALL lbc_lnk( zwztl, 'F', 1._wp )
702	END SELECT
703
704	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
705	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
706
707	zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk)
708	zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk)
709
710	! Compute and add the vorticity term trend
711	! ----------------------------------------
712	jj=2
713	ztne(1,:) = 0.0_wp ; ztnw(1,:) = 0.0_wp ; ztse(1,:) = 0.0_wp ; ztsw(1,:) = 0.0_wp
714	ztnetl(1,:) = 0.0_wp ; ztnwtl(1,:) = 0.0_wp ; ztsetl(1,:) = 0.0_wp ; ztswtl(1,:) = 0.0_wp
715	DO ji = 2, jpi
716	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
717	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
718	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
719	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
720
721	ztnetl(ji,jj) = zwztl(ji-1,jj ) + zwztl(ji ,jj ) + zwztl(ji ,jj-1)
722	ztnwtl(ji,jj) = zwztl(ji-1,jj-1) + zwztl(ji-1,jj ) + zwztl(ji ,jj )
723	ztsetl(ji,jj) = zwztl(ji ,jj ) + zwztl(ji ,jj-1) + zwztl(ji-1,jj-1)
724	ztswtl(ji,jj) = zwztl(ji ,jj-1) + zwztl(ji-1,jj-1) + zwztl(ji-1,jj )
725	END DO
726	DO jj = 3, jpj
727	DO ji = fs_2, jpi ! vector opt.
728	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
729	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
730	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
731	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
732
733	ztnetl(ji,jj) = zwztl(ji-1,jj ) + zwztl(ji ,jj ) + zwztl(ji ,jj-1)
734	ztnwtl(ji,jj) = zwztl(ji-1,jj-1) + zwztl(ji-1,jj ) + zwztl(ji ,jj )
735	ztsetl(ji,jj) = zwztl(ji ,jj ) + zwztl(ji ,jj-1) + zwztl(ji-1,jj-1)
736	ztswtl(ji,jj) = zwztl(ji ,jj-1) + zwztl(ji-1,jj-1) + zwztl(ji-1,jj )
737	END DO
738	END DO
739	DO jj = 2, jpjm1
740	DO ji = fs_2, fs_jpim1 ! vector opt.
741	zuatl = + zfac12 / e1u(ji,jj) * ( ztnetl(ji,jj ) * zwy(ji ,jj ) + ztne(ji,jj ) * zwytl(ji ,jj ) &
742	& + ztnwtl(ji+1,jj) * zwy(ji+1,jj ) + ztnw(ji+1,jj) * zwytl(ji+1,jj ) &
743	& + ztsetl(ji,jj ) * zwy(ji ,jj-1) + ztse(ji,jj ) * zwytl(ji ,jj-1) &
744	& + ztswtl(ji+1,jj) * zwy(ji+1,jj-1) + ztsw(ji+1,jj) * zwytl(ji+1,jj-1))
745
746	zvatl = - zfac12 / e2v(ji,jj) * ( ztswtl(ji,jj+1) * zwx(ji-1,jj+1) + ztsw(ji,jj+1) * zwxtl(ji-1,jj+1) &
747	& + ztsetl(ji,jj+1) * zwx(ji ,jj+1) + ztse(ji,jj+1) * zwxtl(ji ,jj+1) &
748	& + ztnwtl(ji,jj ) * zwx(ji-1,jj ) + ztnw(ji,jj ) * zwxtl(ji-1,jj ) &
749	& + ztnetl(ji,jj ) * zwx(ji ,jj ) + ztne(ji,jj ) * zwxtl(ji ,jj ) )
750	pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) + zuatl
751	pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) + zvatl
752	END DO
753	END DO
754	! ! ===============
755	END DO ! End of slab
756	! ! ===============
757	CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz )
758	CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
759	CALL wrk_dealloc( jpi, jpj, zwxtl , zwytl , zwztl )
760	CALL wrk_dealloc( jpi, jpj, ztnwtl, ztnetl, ztswtl, ztsetl )
761	#if defined key_vvl
762	CALL wrk_dealloc( jpi, jpj, jpk, ze3f_tl )
763	#endif
764	IF ( kt == nitend .AND. .NOT. lk_vvl ) DEALLOCATE( ze3f_tl )
765	!
766	IF( nn_timing == 1 ) CALL timing_stop('vor_een_tan')
767	!
768	END SUBROUTINE vor_een_tan
769
770
771	SUBROUTINE dyn_vor_adj( kt )
772	!!----------------------------------------------------------------------
773	!!
774	!! ** Purpose of the direct routine:
775	!! compute the lateral ocean tracer physics.
776	!!
777	!! ** Action : - Update (ua,va) with the now vorticity term trend
778	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
779	!! and planetary vorticity trends) ('key_trddyn')
780	!!----------------------------------------------------------------------
781	!!
782	INTEGER, INTENT( in ) :: kt ! ocean time-step index
783	!!----------------------------------------------------------------------
784	!
785	IF( nn_timing == 1 ) CALL timing_start('dyn_vor_adj')
786	!
787	! ! vorticity term
788	SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend
789	!
790	CASE ( -1 ) ! esopa: test all possibility with control print
791	CALL vor_een_adj( kt, ntot, ua_ad, va_ad )
792	! CALL vor_mix_adj( kt )
793	CALL vor_ens_adj( kt, ntot, ua_ad, va_ad )
794	! CALL vor_ene_adj( kt, ntot, ua_ad, va_ad )
795	!
796	CASE ( 0 ) ! energy conserving scheme
797	CALL ctl_stop ('vor_ene_adj not available yet')
798	! CALL vor_ene_adj( kt, ntot, ua_ad, va_ad ) ! total vorticity
799	!
800	CASE ( 1 ) ! enstrophy conserving scheme
801	CALL vor_ens_adj( kt, ntot, ua_ad, va_ad ) ! total vorticity
802	!
803	CASE ( 2 ) ! mixed ene-ens scheme
804	CALL ctl_stop ('vor_mix_adj not available yet')
805	! CALL vor_mix_adj( kt ) ! total vorticity (mix=ens-ene)
806	!
807	CASE ( 3 ) ! energy and enstrophy conserving scheme
808	CALL vor_een_adj( kt, ntot, ua_ad, va_ad ) ! total vorticity
809	!
810	END SELECT
811	!
812	IF( nn_timing == 1 ) CALL timing_stop('dyn_vor_adj')
813	!
814	END SUBROUTINE dyn_vor_adj
815	SUBROUTINE vor_ens_adj( kt, kvor, pua_ad, pva_ad )
816	!!----------------------------------------------------------------------
817	!! * ROUTINE vor_ens_adj *
818	!!
819	!! ** Purpose of the direct routine:
820	!! Compute the now total vorticity trend and add it to
821	!! the general trend of the momentum equation.
822	!!
823	!! ** Method of the direct routine:
824	!! Trend evaluated using now fields (centered in time)
825	!! and the Sadourny (1975) flux FORM formulation : conserves the
826	!! potential enstrophy of a horizontally non-divergent flow. the
827	!! trend of the vorticity term is given by:
828	!! * s-coordinate (ln_sco=T), the e3. are inside the derivative:
829	!! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
830	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
831	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
832	!! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ]
833	!! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ]
834	!! Add this trend to the general momentum trend (ua,va):
835	!! (ua,va) = (ua,va) + ( voru , vorv )
836	!!
837	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
838	!! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
839	!! and planetary vorticity trends) ('key_trddyn')
840	!!
841	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
842	!!----------------------------------------------------------------------
843	INTEGER , INTENT(in ) :: kt ! ocean time-step index
844	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
845	! ! =nrvm (relative vorticity or metric)
846	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_ad ! total u-trend
847	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_ad ! total v-trend
848	!!
849	INTEGER :: ji, jj, jk ! dummy loop indices
850	REAL(wp) :: zfact1 ! temporary scalars
851	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! temporary 3D workspace
852	REAL(wp) :: zuav, zvau ! temporary scalars
853	REAL(wp) :: zuavad, zvauad ! temporary scalars
854	REAL(wp), POINTER, DIMENSION(:,:) :: zwxad, zwyad, zwzad ! temporary 3D workspace
855	!!----------------------------------------------------------------------
856	!
857	IF( nn_timing == 1 ) CALL timing_start('vor_ens_adj')
858	!
859	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
860	CALL wrk_alloc( jpi, jpj, zwxad, zwyad, zwzad )
861	!
862	IF( kt == nitend ) THEN
863	IF(lwp) WRITE(numout,*)
864	IF(lwp) WRITE(numout,*) 'dyn_vor_ens_adj : vorticity term: enstrophy conserving scheme'
865	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
866	ENDIF
867
868	! Local constant initialization
869	zfact1 = 0.5 * 0.25
870
871	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zwxad, zwyad, zwzad )
872	! ! ===============
873	DO jk = 1, jpkm1 ! Horizontal slab
874	! ! ===============
875	! Potential vorticity and horizontal fluxes
876	! -----------------------------------------
877	SELECT CASE( kvor ) ! vorticity considered
878	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
879	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
880	CASE ( 3 ) ! metric term
881	DO jj = 1, jpjm1
882	DO ji = 1, fs_jpim1 ! vector opt.
883	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
884	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) )&
885	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
886	END DO
887	END DO
888	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
889	CASE ( 5 ) ! total (coriolis + metric)
890	DO jj = 1, jpjm1
891	DO ji = 1, fs_jpim1 ! vector opt.
892	zwz(ji,jj) = ( ff (ji,jj) &
893	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
894	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
895	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
896	& )
897	END DO
898	END DO
899	END SELECT
900
901	IF( ln_sco ) THEN
902	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
903	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
904	zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
905	zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
906	zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
907	END DO
908	END DO
909	ELSE
910	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
911	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
912	zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
913	zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
914	END DO
915	END DO
916	ENDIF
917	! ===================
918	! Adjoint counterpart
919	! ===================
920	zuavad = 0.0_wp
921	zvauad = 0.0_wp
922	zwxad(:,:) = 0.0_wp
923	zwyad(:,:) = 0.0_wp
924	zwzad(:,:) = 0.0_wp
925	! Compute and add the vorticity term trend
926	! ----------------------------------------
927	DO jj = jpjm1, 2, -1
928	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
929	! Compute and add the vorticity term trend
930	! ----------------------------------------
931	!- Direct counterpart
932	zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
933	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
934	zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
935	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
936	!-
937	zvauad = zvauad + pva_ad(ji,jj,jk) * ( zwz(ji-1,jj) + zwz(ji,jj) )
938	zwzad(ji-1,jj) = zwzad(ji-1,jj) + pva_ad(ji,jj,jk) * zvau
939	zwzad(ji ,jj) = zwzad(ji ,jj) + pva_ad(ji,jj,jk) * zvau
940
941	zuavad = zuavad + pua_ad(ji,jj,jk) * ( zwz(ji,jj-1) + zwz(ji,jj) )
942	zwzad(ji,jj-1) = zwzad(ji,jj-1) + pua_ad(ji,jj,jk) * zuav
943	zwzad(ji,jj ) = zwzad(ji,jj ) + pua_ad(ji,jj,jk) * zuav
944
945	zwxad(ji-1,jj ) = zwxad(ji-1,jj ) - zvauad * zfact1 / e2v(ji,jj)
946	zwxad(ji-1,jj+1) = zwxad(ji-1,jj+1) - zvauad * zfact1 / e2v(ji,jj)
947	zwxad(ji ,jj ) = zwxad(ji ,jj ) - zvauad * zfact1 / e2v(ji,jj)
948	zwxad(ji ,jj+1) = zwxad(ji ,jj+1) - zvauad * zfact1 / e2v(ji,jj)
949	zvauad = 0.0_wp
950
951	zwyad(ji ,jj-1) = zwyad(ji ,jj-1) + zuavad * zfact1 / e1u(ji,jj)
952	zwyad(ji+1,jj-1) = zwyad(ji+1,jj-1) + zuavad * zfact1 / e1u(ji,jj)
953	zwyad(ji ,jj ) = zwyad(ji ,jj ) + zuavad * zfact1 / e1u(ji,jj)
954	zwyad(ji+1,jj ) = zwyad(ji+1,jj ) + zuavad * zfact1 / e1u(ji,jj)
955	zuavad = 0.0_wp
956	END DO
957	END DO
958	IF( ln_sco ) THEN
959	DO jj = jpj, 1, -1 ! caution: don't use (:,:) for this loop
960	DO ji = jpi, 1, -1 ! it causes optimization problems on NEC in auto-tasking
961	zwzad(ji,jj) = zwzad(ji,jj) / fse3f(ji,jj,jk)
962	un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + zwxad(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
963	vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + zwyad(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,jk)
964	zwxad(ji,jj) = 0.0_wp
965	zwyad(ji,jj) = 0.0_wp
966	END DO
967	END DO
968	ELSE
969	DO jj = jpj, 1, -1 ! caution: don't use (:,:) for this loop
970	DO ji = jpi, 1, -1 ! it causes optimization problems on NEC in auto-tasking
971	un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + e2u(ji,jj) * zwxad(ji,jj)
972	vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + e1v(ji,jj) * zwyad(ji,jj)
973	zwxad(ji,jj) = 0.0_wp
974	zwyad(ji,jj) = 0.0_wp
975	END DO
976	END DO
977	ENDIF
978	! Potential vorticity and horizontal fluxes
979	! -----------------------------------------
980	SELECT CASE( kvor ) ! vorticity considered
981	CASE ( 1 ) ! planetary vorticity (Coriolis)
982	zwzad(:,:) = 0.0_wp
983	CASE ( 2 ,4) ! relative vorticity
984	rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:)
985	zwzad(:,:) = 0.0_wp
986	CASE ( 3 ,5 ) ! metric term
987	DO jj = jpjm1, 1, -1
988	DO ji = fs_jpim1, 1, -1 ! vector opt.
989	vn_ad(ji+1,jj,jk) = vn_ad(ji+1,jj,jk) &
990	& + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) ) &
991	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
992	vn_ad(ji ,jj,jk) = vn_ad(ji ,jj,jk) &
993	& + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) ) &
994	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
995	un_ad(ji,jj+1,jk) = un_ad(ji,jj+1,jk) &
996	& - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) ) &
997	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
998	un_ad(ji,jj ,jk) = un_ad(ji,jj ,jk) &
999	& - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) ) &
1000	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
1001	zwzad(ji,jj) = 0.0_wp
1002	END DO
1003	END DO
1004	END SELECT
1005	! ! ===============
1006	END DO ! End of slab
1007	! ! ===============
1008	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
1009	CALL wrk_dealloc( jpi, jpj, zwxad, zwyad, zwzad )
1010	!
1011	IF( nn_timing == 1 ) CALL timing_stop('vor_ens_adj')
1012	!
1013	END SUBROUTINE vor_ens_adj
1014
1015
1016	SUBROUTINE vor_een_adj( kt, kvor, pua_ad, pva_ad )
1017	!!----------------------------------------------------------------------
1018	!! * ROUTINE vor_een_adj *
1019	!!
1020	!! ** Purpose : Compute the now total vorticity trend and add it to
1021	!! the general trend of the momentum equation.
1022	!!
1023	!! ** Method : Trend evaluated using now fields (centered in time)
1024	!! and the Arakawa and Lamb (19XX) flux form formulation : conserves
1025	!! both the horizontal kinetic energy and the potential enstrophy
1026	!! when horizontal divergence is zero.
1027	!! The trend of the vorticity term is given by:
1028	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
1029	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
1030	!! Add this trend to the general momentum trend (ua,va):
1031	!! (ua,va) = (ua,va) + ( voru , vorv )
1032	!!
1033	!! ** Action : - Update (ua,va) with the now vorticity term trend
1034	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
1035	!! and planetary vorticity trends) ('key_trddyn')
1036	!!
1037	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
1038	!!----------------------------------------------------------------------
1039	INTEGER , INTENT(in ) :: kt ! ocean time-step index
1040	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
1041	! ! =nrvm (relative vorticity or metric)
1042	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_ad ! total u-trend
1043	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_ad ! total v-trend
1044	!!
1045	INTEGER :: ji, jj, jk ! dummy loop indices
1046	INTEGER :: ierr
1047	REAL(wp) :: zfac12 ! temporary scalars
1048	REAL(wp) :: zuaad, zvaad ! temporary scalars
1049	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! temporary 2D workspace
1050	REAL(wp), POINTER, DIMENSION(:,:) :: ztnw, ztne, ztsw, ztse ! temporary 3D workspace
1051	REAL(wp), POINTER, DIMENSION(:,:) :: zwxad, zwyad, zwzad ! temporary 2D workspace
1052	REAL(wp), POINTER, DIMENSION(:,:) :: ztnwad, ztnead, ztswad, ztsead ! temporary 3D workspace
1053	#if defined key_vvl
1054	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3f_ad
1055	#else
1056	REAL(wp), POINTER, DIMENSION(:,:,:), SAVE :: ze3f_ad
1057	#endif
1058	!!----------------------------------------------------------------------
1059	!
1060	IF( nn_timing == 1 ) CALL timing_start('vor_een_adj')
1061	!
1062	CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz )
1063	CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
1064	CALL wrk_alloc( jpi, jpj, zwxad , zwyad , zwzad )
1065	CALL wrk_alloc( jpi, jpj, ztnwad, ztnead, ztswad, ztsead )
1066
1067	! local adjoint initailization
1068	zuaad = 0.0_wp ; zvaad = 0.0_wp
1069	zwx (:,:) = 0.0_wp ; zwy (:,:) = 0.0_wp ; zwz (:,:) = 0.0_wp
1070	ztnw(:,:) = 0.0_wp ; ztne(:,:) = 0.0_wp ; ztsw(:,:) = 0.0_wp ; ztse(:,:) = 0.0_wp
1071	zwxad (:,:) = 0.0_wp ; zwyad (:,:) = 0.0_wp ; zwzad (:,:) = 0.0_wp
1072	ztnwad(:,:) = 0.0_wp ; ztnead(:,:) = 0.0_wp ; ztswad(:,:) = 0.0_wp ; ztsead(:,:) = 0.0_wp
1073	ze3f_ad = 0._wp
1074
1075	#if defined key_vvl
1076	CALL wrk_alloc( jpi, jpj, jpk, ze3f_ad )
1077	#endif
1078
1079	IF( kt == nitend ) THEN
1080	IF(lwp) WRITE(numout,*)
1081	IF(lwp) WRITE(numout,*) 'dyn:vor_een_adj : vorticity term: energy and enstrophy conserving scheme'
1082	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
1083	IF( .NOT.lk_vvl ) THEN
1084	ALLOCATE( ze3f_ad(jpi,jpj,jpk) , STAT=ierr )
1085	IF( lk_mpp ) CALL mpp_sum ( ierr )
1086	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
1087	ze3f_ad = 0._wp
1088	ENDIF
1089	ENDIF
1090
1091	IF( kt == nitend .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t)
1092	DO jk = 1, jpk
1093	DO jj = 1, jpjm1
1094	DO ji = 1, jpim1
1095	ze3f_ad(ji,jj,jk) = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
1096	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) ) * 0.25_wp
1097	IF( ze3f_ad(ji,jj,jk) /= 0.0_wp ) ze3f_ad(ji,jj,jk) = 1.0_wp / ze3f_ad(ji,jj,jk)
1098	END DO
1099	END DO
1100	END DO
1101	CALL lbc_lnk( ze3f_ad, 'F', 1._wp )
1102	ENDIF
1103
1104	! Local constant initialization
1105	zfac12 = 1.0_wp / 12.0_wp
1106
1107	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
1108	! ! ===============
1109	DO jk = 1, jpkm1 ! Horizontal slab
1110	! ! ===============
1111
1112	! Potential vorticity and horizontal fluxes (Direct local variables init)
1113	! -----------------------------------------
1114	SELECT CASE( kvor ) ! vorticity considered
1115	CASE ( 1 )
1116	zwz(:,:) = ff(:,:) * ze3f_ad(:,:,jk) ! planetary vorticity (Coriolis)
1117	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) * ze3f_ad(:,:,jk) ! relative vorticity
1118	CASE ( 3 ) ! metric term
1119	DO jj = 1, jpjm1
1120	DO ji = 1, fs_jpim1 ! vector opt.
1121	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
1122	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) )&
1123	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_ad(ji,jj,jk)
1124	END DO
1125	END DO
1126	CALL lbc_lnk( zwz, 'F', 1._wp )
1127	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f_ad(:,:,jk) ! total (relative + planetary vorticity)
1128	CASE ( 5 ) ! total (coriolis + metric)
1129	DO jj = 1, jpjm1
1130	DO ji = 1, fs_jpim1 ! vector opt.
1131	zwz(ji,jj) = ( ff (ji,jj) &
1132	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
1133	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
1134	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
1135	& ) * ze3f_ad(ji,jj,jk)
1136	END DO
1137	END DO
1138	CALL lbc_lnk( zwz, 'F', 1._wp )
1139	END SELECT
1140
1141	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
1142	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
1143
1144	! Compute and add the vorticity term trend
1145	! ----------------------------------------
1146	jj=2
1147	ztne(1,:) = 0.0_wp ; ztnw(1,:) = 0.0_wp ; ztse(1,:) = 0.0_wp ; ztsw(1,:) = 0.0_wp
1148	DO ji = 2, jpi
1149	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1150	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1151	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1152	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1153	END DO
1154	DO jj = 3, jpj
1155	DO ji = fs_2, jpi ! vector opt.
1156	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1157	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1158	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1159	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1160	END DO
1161	END DO
1162
1163	! ===================
1164	! Adjoint counterpart
1165	! ===================
1166
1167	DO jj = jpjm1, 2, -1
1168	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
1169	zuaad = zuaad + pua_ad(ji,jj,jk)
1170	zvaad = zvaad + pva_ad(ji,jj,jk)
1171
1172	zvaad = - zvaad * zfac12 / e2v(ji,jj)
1173	ztswad(ji ,jj+1) = ztswad(ji ,jj+1) + zvaad * zwx (ji-1,jj+1)
1174	zwxad (ji-1,jj+1) = zwxad (ji-1,jj+1) + zvaad * ztsw(ji ,jj+1)
1175	ztsead(ji ,jj+1) = ztsead(ji ,jj+1) + zvaad * zwx (ji ,jj+1)
1176	zwxad (ji ,jj+1) = zwxad (ji ,jj+1) + zvaad * ztse(ji ,jj+1)
1177	ztnwad(ji ,jj ) = ztnwad(ji ,jj ) + zvaad * zwx (ji-1,jj )
1178	zwxad (ji-1,jj ) = zwxad (ji-1,jj ) + zvaad * ztnw(ji ,jj )
1179	ztnead(ji ,jj ) = ztnead(ji ,jj ) + zvaad * zwx (ji ,jj )
1180	zwxad (ji ,jj ) = zwxad (ji ,jj ) + zvaad * ztne(ji ,jj )
1181	zvaad = 0.0_wp
1182
1183	zuaad = zuaad * zfac12 / e1u(ji,jj)
1184	ztnead(ji ,jj ) = ztnead(ji ,jj ) + zuaad * zwy (ji ,jj )
1185	zwyad (ji ,jj ) = zwyad (ji ,jj ) + zuaad * ztne(ji ,jj )
1186	ztnwad(ji+1,jj ) = ztnwad(ji+1,jj ) + zuaad * zwy (ji+1,jj )
1187	zwyad (ji+1,jj ) = zwyad (ji+1,jj ) + zuaad * ztnw(ji+1,jj )
1188	ztsead(ji ,jj ) = ztsead(ji ,jj ) + zuaad * zwy (ji ,jj-1)
1189	zwyad (ji ,jj-1) = zwyad (ji ,jj-1) + zuaad * ztse(ji ,jj )
1190	ztswad(ji+1,jj ) = ztswad(ji+1,jj ) + zuaad * zwy (ji+1,jj-1)
1191	zwyad (ji+1,jj-1) = zwyad (ji+1,jj-1) + zuaad * ztsw(ji+1,jj )
1192	zuaad = 0.0_wp
1193	END DO
1194	END DO
1195	DO jj = jpj, 3, -1
1196	DO ji = jpi, fs_2, -1 ! vector opt.
1197	zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztswad(ji,jj)
1198	zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztswad(ji,jj)
1199	zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztswad(ji,jj)
1200	ztswad(ji ,jj ) = 0.0_wp
1201	zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztsead(ji,jj)
1202	zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztsead(ji,jj)
1203	zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztsead(ji,jj)
1204	ztsead(ji,jj) = 0.0_wp
1205	zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztnwad(ji,jj)
1206	zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnwad(ji,jj)
1207	zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnwad(ji,jj)
1208	ztnwad(ji ,jj ) = 0.0_wp
1209	zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnead(ji,jj)
1210	zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnead(ji,jj)
1211	zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztnead(ji,jj)
1212	ztnead(ji,jj) = 0.0_wp
1213	END DO
1214	END DO
1215	jj=2
1216	DO ji = jpi, 2, -1
1217	zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztswad(ji,jj)
1218	zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztswad(ji,jj)
1219	zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztswad(ji,jj)
1220	ztswad(ji,jj) = 0.0_wp
1221	zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztsead(ji,jj)
1222	zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztsead(ji,jj)
1223	zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztsead(ji,jj)
1224	ztsead(ji ,jj ) = 0.0_wp
1225	zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztnwad(ji,jj)
1226	zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnwad(ji,jj)
1227	zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnwad(ji,jj)
1228	ztnwad(ji ,jj ) = 0.0_wp
1229	zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnead(ji,jj)
1230	zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnead(ji,jj)
1231	zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztnead(ji,jj)
1232	ztnead(ji ,jj ) = 0.0_wp
1233	END DO
1234	ztnead(1,:) = 0.0_wp ; ztnwad(1,:) = 0.0_wp
1235	ztsead(1,:) = 0.0_wp ; ztswad(1,:) = 0.0_wp
1236
1237	vn_ad(:,:,jk) = vn_ad(:,:,jk) + zwyad(:,:) * e1v(:,:) * fse3v(:,:,jk)
1238	un_ad(:,:,jk) = un_ad(:,:,jk) + zwxad(:,:) * e2u(:,:) * fse3u(:,:,jk)
1239	zwyad(:,:) = 0.0_wp
1240	zwxad(:,:) = 0.0_wp
1241
1242	! Potential vorticity and horizontal fluxes
1243	! -----------------------------------------
1244	SELECT CASE( kvor ) ! vorticity considered
1245	CASE ( 1 )
1246	zwzad(:,:) = 0.0_wp
1247	CASE ( 2 )
1248	rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:) * ze3f_ad(:,:,jk)
1249	zwzad(:,:) = 0.0_wp
1250	CASE ( 3 )
1251	CALL lbc_lnk_adj( zwzad, 'F', 1._wp ) ! metric term
1252	DO jj = jpjm1, 1, -1
1253	DO ji = fs_jpim1, 1, -1 ! vector opt.
1254	zwzad(ji ,jj ) = zwzad(ji,jj) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_ad(ji,jj,jk)
1255	vn_ad(ji+1,jj ,jk) = vn_ad(ji+1,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) )
1256	vn_ad(ji ,jj ,jk) = vn_ad(ji ,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) )
1257	un_ad(ji ,jj+1,jk) = un_ad(ji ,jj+1,jk) - zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) )
1258	un_ad(ji ,jj ,jk) = un_ad(ji ,jj ,jk) - zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) )
1259	zwzad(ji ,jj ) = 0.0_wp
1260	END DO
1261	END DO
1262	CASE ( 4 )
1263	rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:) * ze3f_ad(:,:,jk)
1264	zwzad(:,:) = 0.0_wp
1265	CASE ( 5 )
1266	CALL lbc_lnk_adj( zwzad, 'F', 1._wp ) ! total (coriolis + metric)
1267	DO jj = jpjm1, 1, -1
1268	DO ji = fs_jpim1, 1, -1 ! vector opt.
1269	zwzad(ji ,jj ) = zwzad(ji,jj) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_ad(ji,jj,jk)
1270	vn_ad(ji+1,jj ,jk) = vn_ad(ji+1,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) )
1271	vn_tl(ji ,jj ,jk) = vn_tl(ji ,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) )
1272	un_ad(ji ,jj+1,jk) = un_ad(ji ,jj+1,jk) - zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) )
1273	un_ad(ji ,jj ,jk) = un_ad(ji ,jj ,jk) - zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) )
1274	zwzad(ji ,jj ) = 0.0_wp
1275	END DO
1276	END DO
1277	END SELECT
1278	! ! ===============
1279	END DO ! End of slab
1280	! ! ===============
1281	CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz )
1282	CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
1283	CALL wrk_dealloc( jpi, jpj, zwxad , zwyad , zwzad )
1284	CALL wrk_dealloc( jpi, jpj, ztnwad, ztnead, ztswad, ztsead )
1285
1286	#if defined key_vvl
1287	CALL wrk_dealloc( jpi, jpj, jpk, ze3f_ad )
1288	#endif
1289	IF ( kt == nit000 .AND. .NOT. lk_vvl ) DEALLOCATE( ze3f_ad )
1290
1291	IF( nn_timing == 1 ) CALL timing_stop('vor_een_adj')
1292	!
1293	END SUBROUTINE vor_een_adj
1294
1295	SUBROUTINE dyn_vor_init_tam
1296	!!---------------------------------------------------------------------
1297	!! * ROUTINE vor_ctl_tam *
1298	!!
1299	!! ** Purpose : Control the consistency between cpp options for
1300	!! tracer advection schemes
1301	!!----------------------------------------------------------------------
1302	INTEGER :: ioptio ! temporary integer
1303	INTEGER :: ji, jj, jk ! dummy loop indices
1304	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
1305	!!----------------------------------------------------------------------
1306
1307	REWIND ( numnam ) ! Read Namelist namdyn_vor : Vorticity scheme options
1308	READ ( numnam, namdyn_vor )
1309
1310	IF(lwp) THEN ! Namelist print
1311	WRITE(numout,*)
1312	WRITE(numout,*) 'dyn:vor_init_tam : vorticity term : read namelist and control the consistency'
1313	WRITE(numout,*) '~~~~~~~~~~~~~~~'
1314	WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme'
1315	WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
1316	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
1317	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
1318	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
1319	ENDIF
1320	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
1321	! at angles with three ocean points and one land point
1322	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
1323	DO jk = 1, jpk
1324	DO jj = 2, jpjm1
1325	DO ji = 2, jpim1
1326	IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
1327	fmask(ji,jj,jk) = 1._wp
1328	END DO
1329	END DO
1330	END DO
1331	!
1332	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
1333	!
1334	ENDIF
1335	!
1336	ioptio = 0 ! Control of vorticity scheme options
1337	IF( ln_dynvor_ene ) ioptio = ioptio + 1
1338	IF( ln_dynvor_ens ) ioptio = ioptio + 1
1339	IF( ln_dynvor_mix ) ioptio = ioptio + 1
1340	IF( ln_dynvor_een ) ioptio = ioptio + 1
1341	IF( lk_esopa ) ioptio = 1
1342
1343	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
1344
1345	! ! Set nvor (type of scheme for vorticity)
1346	IF( ln_dynvor_ene ) nvor = 0
1347	IF( ln_dynvor_ens ) nvor = 1
1348	IF( ln_dynvor_mix ) nvor = 2
1349	IF( ln_dynvor_een ) nvor = 3
1350	IF( lk_esopa ) nvor = -1
1351
1352	! ! Set ncor, nrvm, ntot (type of vorticity)
1353	IF(lwp) WRITE(numout,*)
1354	ncor = 1
1355	IF( ln_dynadv_vec ) THEN
1356	IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'
1357	nrvm = 2
1358	ntot = 4
1359	ELSE
1360	IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'
1361	nrvm = 3
1362	ntot = 5
1363	ENDIF
1364
1365	IF(lwp) THEN ! Print the choice
1366	WRITE(numout,*)
1367	IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme'
1368	IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme'
1369	IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme'
1370	IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme'
1371	IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options'
1372	ENDIF
1373	!
1374	END SUBROUTINE dyn_vor_init_tam
1375
1376	SUBROUTINE dyn_vor_adj_tst( kumadt )
1377	!!-----------------------------------------------------------------------
1378	!!
1379	!! * ROUTINE dyn_adv_adj_tst *
1380	!!
1381	!! ** Purpose : Test the adjoint routine.
1382	!!
1383	!! ** Method : Verify the scalar product
1384	!!
1385	!! ( L dx )^T W dy = dx^T L^T W dy
1386	!!
1387	!! where L = tangent routine
1388	!! L^T = adjoint routine
1389	!! W = diagonal matrix of scale factors
1390	!! dx = input perturbation (random field)
1391	!! dy = L dx
1392	!!
1393	!! ** Action : Separate tests are applied for the following dx and dy:
1394	!!
1395	!! 1) dx = ( SSH ) and dy = ( SSH )
1396	!!
1397	!! History :
1398	!! ! 08-08 (A. Vidard)
1399	!!-----------------------------------------------------------------------
1400	!! * Modules used
1401
1402	!! * Arguments
1403	INTEGER, INTENT(IN) :: &
1404	& kumadt ! Output unit
1405
1406	INTEGER :: &
1407	& ji, & ! dummy loop indices
1408	& jj, &
1409	& jk, &
1410	& jt
1411	INTEGER, DIMENSION(jpi,jpj) :: &
1412	& iseed_2d ! 2D seed for the random number generator
1413
1414	!! * Local declarations
1415	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1416	& zun_tlin, & ! Tangent input: now u-velocity
1417	& zvn_tlin, & ! Tangent input: now v-velocity
1418	& zrotn_tlin, & ! Tangent input: now rot
1419	& zun_adout, & ! Adjoint output: now u-velocity
1420	& zvn_adout, & ! Adjoint output: now v-velocity
1421	& zrotn_adout, & ! Adjoint output: now rot
1422	& zua_adout, & ! Tangent output: after u-velocity
1423	& zva_adout, & ! Tangent output: after v-velocity
1424	& zua_tlin, & ! Tangent output: after u-velocity
1425	& zva_tlin, & ! Tangent output: after v-velocity
1426	& zua_tlout, & ! Tangent output: after u-velocity
1427	& zva_tlout, & ! Tangent output: after v-velocity
1428	& zua_adin, & ! Tangent output: after u-velocity
1429	& zva_adin, & ! Tangent output: after v-velocity
1430	& zau, & ! 3D random field for rotn
1431	& zav, & ! 3D random field for rotn
1432	& znu, & ! 3D random field for u
1433	& znv ! 3D random field for v
1434	REAL(KIND=wp) :: &
1435	& zsp1, & ! scalar product involving the tangent routine
1436	& zsp1_1, & ! scalar product components
1437	& zsp1_2, &
1438	& zsp2, & ! scalar product involving the adjoint routine
1439	& zsp2_1, & ! scalar product components
1440	& zsp2_2, &
1441	& zsp2_3, &
1442	& zsp2_4, &
1443	& zsp2_5
1444	CHARACTER(LEN=14) :: cl_name
1445
1446	! Allocate memory
1447
1448	ALLOCATE( &
1449	& zun_tlin(jpi,jpj,jpk), &
1450	& zvn_tlin(jpi,jpj,jpk), &
1451	& zrotn_tlin(jpi,jpj,jpk), &
1452	& zun_adout(jpi,jpj,jpk), &
1453	& zvn_adout(jpi,jpj,jpk), &
1454	& zrotn_adout(jpi,jpj,jpk), &
1455	& zua_adout(jpi,jpj,jpk), &
1456	& zva_adout(jpi,jpj,jpk), &
1457	& zua_tlin(jpi,jpj,jpk), &
1458	& zva_tlin(jpi,jpj,jpk), &
1459	& zua_tlout(jpi,jpj,jpk), &
1460	& zva_tlout(jpi,jpj,jpk), &
1461	& zua_adin(jpi,jpj,jpk), &
1462	& zva_adin(jpi,jpj,jpk), &
1463	& zau(jpi,jpj,jpk), &
1464	& zav(jpi,jpj,jpk), &
1465	& znu(jpi,jpj,jpk), &
1466	& znv(jpi,jpj,jpk) &
1467	& )
1468
1469	! init ntot parameter
1470	CALL dyn_vor_init_tam ! initialisation & control of options
1471
1472	DO jt = 1, 2
1473	IF (jt == 1) nvor=1 ! enstrophy conserving scheme
1474	IF (jt == 2) nvor=3 ! energy and enstrophy conserving scheme
1475
1476	! Initialize rotn
1477	!AV: it calls cla_init the tries to reallocate already allocated arrays...nasty
1478	!AV CALL div_cur ( nit000 )
1479
1480	!==================================================================
1481	! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
1482	! dy = ( hdivb_tl, hdivn_tl )
1483	!==================================================================
1484
1485	!--------------------------------------------------------------------
1486	! Reset the tangent and adjoint variables
1487	!--------------------------------------------------------------------
1488
1489	zun_tlin(:,:,:) = 0.0_wp
1490	zvn_tlin(:,:,:) = 0.0_wp
1491	zrotn_tlin(:,:,:) = 0.0_wp
1492	zun_adout(:,:,:) = 0.0_wp
1493	zvn_adout(:,:,:) = 0.0_wp
1494	zrotn_adout(:,:,:) = 0.0_wp
1495	zua_tlout(:,:,:) = 0.0_wp
1496	zva_tlout(:,:,:) = 0.0_wp
1497	zua_adin(:,:,:) = 0.0_wp
1498	zva_adin(:,:,:) = 0.0_wp
1499	zua_adout(:,:,:) = 0.0_wp
1500	zva_adout(:,:,:) = 0.0_wp
1501	zua_tlin(:,:,:) = 0.0_wp
1502	zva_tlin(:,:,:) = 0.0_wp
1503	znu(:,:,:) = 0.0_wp
1504	znv(:,:,:) = 0.0_wp
1505	zau(:,:,:) = 0.0_wp
1506	zav(:,:,:) = 0.0_wp
1507
1508
1509	un_tl(:,:,:) = 0.0_wp
1510	vn_tl(:,:,:) = 0.0_wp
1511	ua_tl(:,:,:) = 0.0_wp
1512	va_tl(:,:,:) = 0.0_wp
1513	un_ad(:,:,:) = 0.0_wp
1514	vn_ad(:,:,:) = 0.0_wp
1515	ua_ad(:,:,:) = 0.0_wp
1516	va_ad(:,:,:) = 0.0_wp
1517	rotn_tl(:,:,:) = 0.0_wp
1518	rotn_ad(:,:,:) = 0.0_wp
1519
1520	!--------------------------------------------------------------------
1521	! Initialize the tangent input with random noise: dx
1522	!--------------------------------------------------------------------
1523
1524	CALL grid_random( znu, 'U', 0.0_wp, stdu )
1525	CALL grid_random( znv, 'V', 0.0_wp, stdv )
1526	CALL grid_random( zau, 'U', 0.0_wp, stdu )
1527	CALL grid_random( zav, 'V', 0.0_wp, stdv )
1528
1529	DO jk = 1, jpk
1530	DO jj = nldj, nlej
1531	DO ji = nldi, nlei
1532	zun_tlin(ji,jj,jk) = znu(ji,jj,jk)
1533	zvn_tlin(ji,jj,jk) = znv(ji,jj,jk)
1534	zua_tlin(ji,jj,jk) = zau(ji,jj,jk)
1535	zva_tlin(ji,jj,jk) = zav(ji,jj,jk)
1536	END DO
1537	END DO
1538	END DO
1539	un_tl(:,:,:) = zun_tlin(:,:,:)
1540	vn_tl(:,:,:) = zvn_tlin(:,:,:)
1541	ua_tl(:,:,:) = zua_tlin(:,:,:)
1542	va_tl(:,:,:) = zva_tlin(:,:,:)
1543
1544	! initialize rotn_tl with noise
1545	CALL div_cur_tan ( nit000 )
1546
1547	DO jk = 1, jpk
1548	DO jj = nldj, nlej
1549	DO ji = nldi, nlei
1550	zrotn_tlin(ji,jj,jk) = rotn_tl(ji,jj,jk)
1551	END DO
1552	END DO
1553	END DO
1554	rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
1555
1556
1557	IF (nvor == 1 ) CALL vor_ens_tan( nit000, ntot, ua_tl, va_tl )
1558	IF (nvor == 3 ) CALL vor_een_tan( nit000, ntot, ua_tl, va_tl )
1559	zua_tlout(:,:,:) = ua_tl(:,:,:)
1560	zva_tlout(:,:,:) = va_tl(:,:,:)
1561
1562	!--------------------------------------------------------------------
1563	! Initialize the adjoint variables: dy^* = W dy
1564	!--------------------------------------------------------------------
1565
1566	DO jk = 1, jpk
1567	DO jj = nldj, nlej
1568	DO ji = nldi, nlei
1569	zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
1570	& * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
1571	& * umask(ji,jj,jk)
1572	zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
1573	& * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
1574	& * vmask(ji,jj,jk)
1575	END DO
1576	END DO
1577	END DO
1578	!--------------------------------------------------------------------
1579	! Compute the scalar product: ( L dx )^T W dy
1580	!--------------------------------------------------------------------
1581
1582	zsp1_1 = DOT_PRODUCT( zua_tlout, zua_adin )
1583	zsp1_2 = DOT_PRODUCT( zva_tlout, zva_adin )
1584	zsp1 = zsp1_1 + zsp1_2
1585
1586	!--------------------------------------------------------------------
1587	! Call the adjoint routine: dx^* = L^T dy^*
1588	!--------------------------------------------------------------------
1589
1590	ua_ad(:,:,:) = zua_adin(:,:,:)
1591	va_ad(:,:,:) = zva_adin(:,:,:)
1592
1593
1594	IF (nvor == 1 ) CALL vor_ens_adj( nitend, ntot, ua_ad, va_ad )
1595	IF (nvor == 3 ) CALL vor_een_adj( nitend, ntot, ua_ad, va_ad )
1596	zun_adout(:,:,:) = un_ad(:,:,:)
1597	zvn_adout(:,:,:) = vn_ad(:,:,:)
1598	zrotn_adout(:,:,:) = rotn_ad(:,:,:)
1599	zua_adout(:,:,:) = ua_ad(:,:,:)
1600	zva_adout(:,:,:) = va_ad(:,:,:)
1601
1602	zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout )
1603	zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout )
1604	zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
1605	zsp2_4 = DOT_PRODUCT( zua_tlin, zua_adout )
1606	zsp2_5 = DOT_PRODUCT( zva_tlin, zva_adout )
1607	zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5
1608
1609	! Compare the scalar products
1610
1611	! 14 char:'12345678901234'
1612	IF (nvor == 1 ) cl_name = 'dynvor_adj ens'
1613	IF (nvor == 3 ) cl_name = 'dynvor_adj een'
1614
1615	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1616	END DO
1617
1618	DEALLOCATE( &
1619	& zun_tlin, &
1620	& zvn_tlin, &
1621	& zrotn_tlin, &
1622	& zun_adout, &
1623	& zvn_adout, &
1624	& zrotn_adout, &
1625	& zua_adout, &
1626	& zva_adout, &
1627	& zua_tlin, &
1628	& zva_tlin, &
1629	& zua_tlout, &
1630	& zva_tlout, &
1631	& zua_adin, &
1632	& zva_adin, &
1633	& zau, &
1634	& zav, &
1635	& znu, &
1636	& znv &
1637	& )
1638	END SUBROUTINE dyn_vor_adj_tst
1639	!!=============================================================================
1640	#endif
1641	END MODULE dynvor_tam

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

Download in other formats: