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_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

Context Navigation

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN/dynvor_tam.F90 @ 3611

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

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

Download in other formats: