Context Navigation

dynspg_flt_tam.F90 @ 2587

Last change on this file since 2587 was 2587, checked in by vidard, 13 years ago
refer to ticket #798
File size: 77.1 KB

Line
1	MODULE dynspg_flt_tam
2	!!----------------------------------------------------------------------
3	!! This software is governed by the CeCILL licence (Version 2)
4	!!----------------------------------------------------------------------
5	#if defined key_tam
6	!!======================================================================
7	!! * MODULE dynspg_flt_tam : TANGENT/ADJOINT OF MODULE dynspg_flt *
8	!!
9	!! Ocean dynamics: surface pressure gradient trend
10	!!
11	!!======================================================================
12	!! History of the direct module:
13	!! 8.0 ! 98-05 (G. Roullet) free surface
14	!! ! 98-10 (G. Madec, M. Imbard) release 8.2
15	!! 8.5 ! 02-08 (G. Madec) F90: Free form and module
16	!! " " ! 02-11 (C. Talandier, A-M Treguier) Open boundaries
17	!! 9.0 ! 04-08 (C. Talandier) New trends organization
18	!! " " ! 05-11 (V. Garnier) Surface pressure gradient organization
19	!! " " ! 06-07 (S. Masson) distributed restart using iom
20	!! History of the TAM module:
21	!! 9.0 ! 08-08 (A. Vidard) skeleton
22	!! - ! 08-08 (A. Weaver) original version
23	!!----------------------------------------------------------------------
24	# if defined key_dynspg_flt \|\| defined key_esopa
25	!!----------------------------------------------------------------------
26	!! 'key_dynspg_flt' filtered free surface
27	!!----------------------------------------------------------------------
28	!!----------------------------------------------------------------------
29	!! dyn_spg_flt_tan : update the momentum trend with the surface pressure
30	!! gradient in the filtered free surface case with
31	!! vector optimization (tangent routine)
32	!! dyn_spg_flt_adj : update the momentum trend with the surface pressure
33	!! gradient in the filtered free surface case with
34	!! vector optimization (adjoint routine)
35	!! dyn_spg_flt_adj_tst : Test of the adjoint routine
36	!!----------------------------------------------------------------------
37	USE par_kind , ONLY: & ! Precision variables
38	& wp
39	USE prtctl ! Print control
40	USE par_oce , ONLY: & ! Ocean space and time domain variables
41	& jpi, &
42	& jpj, &
43	& jpk, &
44	& jpim1, &
45	& jpjm1, &
46	& jpkm1, &
47	& jpr2di, &
48	& jpr2dj, &
49	& jpij, &
50	& jpiglo, &
51	& lk_vopt_loop
52	USE in_out_manager, ONLY: & ! I/O manager
53	& lwp, &
54	& numout, &
55	& nit000, &
56	& nitend, &
57	& ctl_stop, &
58	& ln_ctl, &
59	& ctmp1
60	USE phycst , ONLY: & ! Physical constants
61	& grav, &
62	& rauw
63	USE lib_mpp , ONLY: & ! distributed memory computing
64	& lk_mpp, &
65	& mpp_sum
66	USE mppsum , ONLY: & ! Summation of arrays across processors
67	& mpp_sum_indep
68	USE dom_oce , ONLY: & ! Ocean space and time domain
69	& e1u, &
70	& e2u, &
71	& e1v, &
72	& e2v, &
73	& e1t, &
74	& e2t, &
75	#if defined key_zco
76	& e3t_0, &
77	#else
78	& e3u, &
79	& e3v, &
80	#endif
81	& tmask, &
82	& umask, &
83	& vmask, &
84	& bmask, &
85	& rdt, &
86	& neuler, &
87	& n_cla, &
88	& mig, &
89	& mjg, &
90	& nldi, &
91	& nldj, &
92	& nlei, &
93	& nlej, &
94	& atfp, &
95	& atfp1, &
96	& lk_vvl
97	USE solver , ONLY: & ! Solvers
98	& sol_mat
99	USE sol_oce , ONLY: & ! Solver variables in memory
100	& nsolv, &
101	& ncut, &
102	& niter, &
103	& eps, &
104	& epsr, &
105	& rnorme, &
106	& res, &
107	& rnu, &
108	& c_solver_pt, &
109	& gcdprc
110	USE oce_tam , ONLY: & ! Tangent-linear and adjoint model variables
111	& ub_tl, &
112	& vb_tl, &
113	& ua_tl, &
114	& va_tl, &
115	& wn_tl, &
116	& ub_ad, &
117	& vb_ad, &
118	& ua_ad, &
119	& va_ad, &
120	& wn_ad, &
121	& spgu_tl, &
122	& spgv_tl, &
123	& spgu_ad, &
124	& spgv_ad, &
125	& sshb_tl, &
126	& sshn_tl, &
127	& sshb_ad, &
128	& sshn_ad
129	USE sbc_oce_tam , ONLY: & ! Surface BCs for tangent and adjoint model
130	& emp_tl, &
131	& emp_ad
132	#if defined key_obc
133	# if defined key_pomme_r025
134	USE obc_oce
135	USE obcdyn_tam
136	# else
137	Error, OBC not ready.
138	# endif
139	#endif
140	USE sol_oce , ONLY: & ! Solver arrays
141	& gcdmat
142	USE sol_oce_tam , ONLY: & ! Solver arrays for tangent and adjoint model
143	& gcx_tl, &
144	& gcxb_tl, &
145	& gcb_tl, &
146	& gcx_ad, &
147	& gcxb_ad, &
148	& gcb_ad
149	USE solsor_tam , ONLY: & ! SOR solver for tangent and adjoint model
150	& sol_sor_tan, &
151	& sol_sor_adj
152	USE solpcg_tam , ONLY: & ! PCG solver for tangent and adjoint model
153	& sol_pcg_tan, &
154	& sol_pcg_adj
155	USE solfet_tam , ONLY: & ! FETI solver for tangent and adjoint model
156	& sol_fet_tan, &
157	& sol_fet_adj
158	USE lbclnk , ONLY: & ! Lateral boundary conditions
159	& lbc_lnk, &
160	& lbc_lnk_e
161	USE lbclnk_tam , ONLY: & ! TAM lateral boundary conditions
162	& lbc_lnk_adj, &
163	& lbc_lnk_e_adj
164	USE gridrandom , ONLY: & ! Random Gaussian noise on grids
165	& grid_random
166	USE dotprodfld , ONLY: & ! Computes dot product for 3D and 2D fields
167	& dot_product
168	USE paresp , ONLY: & ! Normalized energy weights
169	& wesp_emp, &
170	& wesp_ssh
171	USE cla_dynspg_tam, ONLY: & ! Cross land advection on surface pressure
172	& dyn_spg_cla_tan, &
173	& dyn_spg_cla_adj
174	USE tstool_tam , ONLY: &
175	& stdu, &
176	& stdv, &
177	& stdw, &
178	& stdemp, &
179	& stdssh, &
180	& stdgc, &
181	& prntst_adj, &
182	& prntst_tlm
183
184	IMPLICIT NONE
185	PRIVATE
186
187	!! * Accessibility
188	PUBLIC dyn_spg_flt_tan, & ! routine called by step_tan.F90
189	& dyn_spg_flt_adj, & ! routine called by step_adj.F90
190	& dyn_spg_flt_adj_tst ! routine called by the tst.F90
191	#if defined key_tst_tlm
192	PUBLIC dyn_spg_flt_tlm_tst
193	#endif
194	!! * Substitutions
195	# include "domzgr_substitute.h90"
196	# include "vectopt_loop_substitute.h90"
197	!!----------------------------------------------------------------------
198	CONTAINS
199
200	SUBROUTINE dyn_spg_flt_tan( kt, kindic )
201	!!---------------------------------------------------------------------
202	!! * routine dyn_spg_flt_tan *
203	!!
204	!! ** Purpose of the direct routine:
205	!! Compute the now trend due to the surface pressure
206	!! gradient in case of filtered free surface formulation and add
207	!! it to the general trend of momentum equation.
208	!!
209	!! ** Method of the direct routine:
210	!! Filtered free surface formulation. The surface
211	!! pressure gradient is given by:
212	!! spgu = 1/rau0 d/dx(ps) = 1/e1u di( sshn + rnu btda )
213	!! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + rnu btda )
214	!! where sshn is the free surface elevation and btda is the after
215	!! of the free surface elevation
216	!! -1- compute the after sea surface elevation from the kinematic
217	!! surface boundary condition:
218	!! zssha = sshb + 2 rdt ( wn - emp )
219	!! Time filter applied on now sea surface elevation to avoid
220	!! the divergence of two consecutive time-steps and swap of free
221	!! surface arrays to start the next time step:
222	!! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
223	!! sshn = zssha
224	!! -2- evaluate the surface presure trend (including the addi-
225	!! tional force) in three steps:
226	!! a- compute the right hand side of the elliptic equation:
227	!! gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
228	!! where (spgu,spgv) are given by:
229	!! spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
230	!! - grav 2 rdt hu /e1u di[sshn + emp]
231	!! spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
232	!! - grav 2 rdt hv /e2v dj[sshn + emp]
233	!! and define the first guess from previous computation :
234	!! zbtd = btda
235	!! btda = 2 zbtd - btdb
236	!! btdb = zbtd
237	!! b- compute the relative accuracy to be reached by the
238	!! iterative solver
239	!! c- apply the solver by a call to sol... routine
240	!! -3- compute and add the free surface pressure gradient inclu-
241	!! ding the additional force used to stabilize the equation.
242	!!
243	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
244	!!
245	!! References : Roullet and Madec 1999, JGR.
246	!!---------------------------------------------------------------------
247	!! * Arguments
248	INTEGER, INTENT( IN ) :: &
249	& kt ! ocean time-step index
250	INTEGER, INTENT( OUT ) :: &
251	& kindic ! solver convergence flag (<0 if not converge)
252
253	!! * Local declarations
254	INTEGER :: &
255	& ji, & ! dummy loop indices
256	& jj, &
257	& jk
258	REAL(wp) :: &
259	& z2dt, & ! temporary scalars
260	& z2dtg, &
261	& zraur, &
262	& znugdt, &
263	& znurau, &
264	& zgcb, &
265	& zbtd, &
266	& ztdgu, &
267	& ztdgv
268	!! Variable volume
269	REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace
270	& zsshub, &
271	& zsshua, &
272	& zsshvb, &
273	& zsshva, &
274	& zssha
275	REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace
276	& zfse3ub, &
277	& zfse3ua, &
278	& zfse3vb, &
279	& zfse3va
280	!!----------------------------------------------------------------------
281	!
282	IF( kt == nit000 ) THEN
283
284	IF(lwp) WRITE(numout,*)
285	IF(lwp) WRITE(numout,*) 'dyn_spg_flt_tan : surface pressure gradient trend'
286	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ (free surface constant volume case)'
287
288	! set to zero free surface specific arrays
289	spgu_tl(:,:) = 0.0_wp ! surface pressure gradient (i-direction)
290	spgv_tl(:,:) = 0.0_wp ! surface pressure gradient (j-direction)
291
292	ENDIF
293
294	! Local constant initialization
295	IF ( neuler == 0 .AND. kt == nit000 ) THEN
296
297	z2dt = rdt ! time step: Euler if restart from rest
298	CALL sol_mat(kt) ! initialize matrix
299
300	ELSEIF ( neuler == 0 .AND. kt == nit000 + 1 ) THEN
301
302	z2dt = 2.0_wp * rdt ! time step: leap-frog
303	CALL sol_mat(kt) ! initialize matrix
304
305	ELSEIF ( neuler /= 0 .AND. kt == nit000 ) THEN
306
307	z2dt = 2.0_wp * rdt ! time step: leap-frog
308	CALL sol_mat(kt) ! initialize matrix
309
310	ELSE
311
312	z2dt = 2.0_wp * rdt ! time step: leap-frog
313
314	ENDIF
315
316	z2dtg = grav * z2dt
317	zraur = 1.0_wp / rauw
318	znugdt = rnu * grav * z2dt
319	znurau = znugdt * zraur
320
321	!! Explicit physics with thickness weighted updates
322	IF( lk_vvl ) THEN ! variable volume
323
324	CALL ctl_stop( 'dyn_spg_flt_tan: lk_vvl is not available' )
325
326	ELSE ! fixed volume
327
328	! Surface pressure gradient (now)
329	DO jj = 2, jpjm1
330	DO ji = fs_2, fs_jpim1 ! vector opt.
331	spgu_tl(ji,jj) = - grav * ( sshn_tl(ji+1,jj) - sshn_tl(ji,jj) ) / e1u(ji,jj)
332	spgv_tl(ji,jj) = - grav * ( sshn_tl(ji,jj+1) - sshn_tl(ji,jj) ) / e2v(ji,jj)
333	END DO
334	END DO
335
336	! Add the surface pressure trend to the general trend
337	DO jk = 1, jpkm1
338	DO jj = 2, jpjm1
339	DO ji = fs_2, fs_jpim1 ! vector opt.
340	ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + spgu_tl(ji,jj)
341	va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + spgv_tl(ji,jj)
342	END DO
343	END DO
344	END DO
345
346	! Evaluate the masked next velocity (effect of the additional force not included)
347	DO jk = 1, jpkm1
348	DO jj = 2, jpjm1
349	DO ji = fs_2, fs_jpim1 ! vector opt.
350	ua_tl(ji,jj,jk) = ( ub_tl(ji,jj,jk) + z2dt * ua_tl(ji,jj,jk) ) * umask(ji,jj,jk)
351	va_tl(ji,jj,jk) = ( vb_tl(ji,jj,jk) + z2dt * va_tl(ji,jj,jk) ) * vmask(ji,jj,jk)
352	END DO
353	END DO
354	END DO
355
356	ENDIF
357
358	#if defined key_obc
359	CALL obc_dyn_tan( kt ) ! Update velocities on each open boundary with the radiation algorithm
360	#endif
361	#if defined key_bdy
362	CALL ctl_stop( 'dyn_spg_flt_tan: key_bdy is not available' )
363	#endif
364	#if defined key_agrif
365	CALL ctl_stop( 'dyn_spg_flt_tan: key_agrif is not available' )
366	#endif
367	#if defined key_orca_r2
368	IF( n_cla == 1 ) CALL dyn_spg_cla_tan( kt ) ! Cross Land Advection (update (ua,va))
369	#endif
370
371	! compute the next vertically averaged velocity (effect of the additional force not included)
372	! ---------------------------------------------
373	DO jj = 2, jpjm1
374	DO ji = fs_2, fs_jpim1 ! vector opt.
375	spgu_tl(ji,jj) = 0.0_wp
376	spgv_tl(ji,jj) = 0.0_wp
377	END DO
378	END DO
379
380	! vertical sum
381	!CDIR NOLOOPCHG
382	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
383	DO jk = 1, jpkm1
384	DO ji = 1, jpij
385	spgu_tl(ji,1) = spgu_tl(ji,1) + fse3u(ji,1,jk) * ua_tl(ji,1,jk)
386	spgv_tl(ji,1) = spgv_tl(ji,1) + fse3v(ji,1,jk) * va_tl(ji,1,jk)
387	END DO
388	END DO
389	ELSE ! No vector opt.
390	DO jk = 1, jpkm1
391	DO jj = 2, jpjm1
392	DO ji = 2, jpim1
393	spgu_tl(ji,jj) = spgu_tl(ji,jj) + fse3u(ji,jj,jk) * ua_tl(ji,jj,jk)
394	spgv_tl(ji,jj) = spgv_tl(ji,jj) + fse3v(ji,jj,jk) * va_tl(ji,jj,jk)
395	END DO
396	END DO
397	END DO
398	ENDIF
399
400	! transport: multiplied by the horizontal scale factor
401	DO jj = 2, jpjm1
402	DO ji = fs_2, fs_jpim1 ! vector opt.
403	spgu_tl(ji,jj) = spgu_tl(ji,jj) * e2u(ji,jj)
404	spgv_tl(ji,jj) = spgv_tl(ji,jj) * e1v(ji,jj)
405	END DO
406	END DO
407
408	! Boundary conditions on (spgu_tl,spgv_tl)
409	CALL lbc_lnk( spgu_tl, 'U', -1.0_wp )
410	CALL lbc_lnk( spgv_tl, 'V', -1.0_wp )
411
412	IF( lk_vvl ) CALL ctl_stop( 'dyn_spg_flt_tan: lk_vvl is not available' )
413
414	! Right hand side of the elliptic equation and first guess
415	! -----------------------------------------------------------
416	DO jj = 2, jpjm1
417	DO ji = fs_2, fs_jpim1 ! vector opt.
418	! Divergence of the after vertically averaged velocity
419	zgcb = spgu_tl(ji,jj) - spgu_tl(ji-1,jj) &
420	& + spgv_tl(ji,jj) - spgv_tl(ji,jj-1)
421	gcb_tl(ji,jj) = gcdprc(ji,jj) * zgcb
422	! First guess of the after barotropic transport divergence
423	zbtd = gcx_tl(ji,jj)
424	gcx_tl (ji,jj) = 2.0_wp * zbtd - gcxb_tl(ji,jj)
425	gcxb_tl(ji,jj) = zbtd
426	END DO
427	END DO
428	! apply the lateral boundary conditions
429	IF( nsolv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) CALL lbc_lnk_e( gcb_tl, c_solver_pt, 1.0_wp )
430	#if defined key_agrif
431	CALL ctl_stop( 'dyn_spg_flt_tan: key_agrif is not available' )
432	#endif
433
434	! Relative precision (computation on one processor)
435	! ------------------
436	rnorme =0.0_wp
437	rnorme = SUM( gcb_tl(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb_tl(1:jpi,1:jpj) * bmask(:,:) )
438	IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain
439
440	epsr = eps * eps * rnorme
441	ncut = 0
442	! if rnorme is 0, the solution is 0, the solver is not called
443	IF( rnorme == 0.0_wp ) THEN
444	gcx_tl(:,:) = 0.0_wp
445	res = 0.0_wp
446	niter = 0
447	ncut = 999
448	ENDIF
449
450	! Evaluate the next transport divergence
451	! --------------------------------------
452	! Iterarive solver for the elliptic equation (except IF sol.=0)
453	! (output in gcx with boundary conditions applied)
454	kindic = 0
455	IF( ncut == 0 ) THEN
456	IF( nsolv == 1 ) THEN ! diagonal preconditioned conjuguate gradient
457	CALL sol_pcg_tan( kt, kindic )
458	ELSEIF( nsolv == 2 ) THEN ! successive-over-relaxation
459	CALL sol_sor_tan( kt, kindic )
460	ELSEIF( nsolv == 3 ) THEN ! FETI solver
461	CALL sol_fet_tan( kt, kindic )
462	ELSE ! e r r o r in nsolv namelist parameter
463	WRITE(ctmp1,*) ' ~~~~~~~~~~~ not = ', nsolv
464	CALL ctl_stop( ' dyn_spg_flt_tan : e r r o r, nsolv = 1, 2 or 3', ctmp1 )
465	ENDIF
466	ENDIF
467
468	! Transport divergence gradient multiplied by z2dt
469	! --------------------------------------------====
470	DO jj = 2, jpjm1
471	DO ji = fs_2, fs_jpim1 ! vector opt.
472	! trend of Transport divergence gradient
473	ztdgu = znugdt * ( gcx_tl(ji+1,jj ) - gcx_tl(ji,jj) ) / e1u(ji,jj)
474	ztdgv = znugdt * ( gcx_tl(ji ,jj+1) - gcx_tl(ji,jj) ) / e2v(ji,jj)
475	! multiplied by z2dt
476	#if defined key_obc
477	! caution : grad D = 0 along open boundaries
478	spgu_tl(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
479	spgv_tl(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
480	#elif defined key_bdy
481	CALL ctl_stop( 'dyn_spg_flt_tan: key_bdy is not available' )
482	#else
483	spgu_tl(ji,jj) = z2dt * ztdgu
484	spgv_tl(ji,jj) = z2dt * ztdgv
485	#endif
486	END DO
487	END DO
488
489	#if defined key_agrif
490	CALL ctl_stop( 'dyn_spg_flt_tan: key_agrif is not available' )
491	#endif
492
493	! 7. Add the trends multiplied by z2dt to the after velocity
494	! -----------------------------------------------------------
495	! ( c a u t i o n : (ua_tl,va_tl) here are the after velocity not the
496	! trend, the leap-frog time stepping will not
497	! be done in dynnxt_tan.F90 routine)
498	DO jk = 1, jpkm1
499	DO jj = 2, jpjm1
500	DO ji = fs_2, fs_jpim1 ! vector opt.
501	ua_tl(ji,jj,jk) = ( ua_tl(ji,jj,jk) + spgu_tl(ji,jj) ) * umask(ji,jj,jk)
502	va_tl(ji,jj,jk) = ( va_tl(ji,jj,jk) + spgv_tl(ji,jj) ) * vmask(ji,jj,jk)
503	END DO
504	END DO
505	END DO
506
507	! Sea surface elevation time stepping
508	! -----------------------------------
509	IF( lk_vvl ) THEN ! after free surface elevation
510	CALL ctl_stop( 'dyn_spg_flt_tan: lk_vvl is not available' )
511	ELSE
512	zssha(:,:) = sshb_tl(:,:) + z2dt * ( wn_tl(:,:,1) - emp_tl(:,:) * zraur ) * tmask(:,:,1)
513	ENDIF
514	#if defined key_obc
515	zssha(:,:) = zssha(:,:) * obctmsk(:,:)
516	CALL lbc_lnk(zssha,'T',1.) ! absolutly compulsory !! (jmm)
517	#endif
518
519	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler (forward) time stepping, no time filter
520	! swap of arrays
521	sshb_tl(:,:) = sshn_tl (:,:)
522	sshn_tl(:,:) = zssha(:,:)
523	ELSE ! Leap-frog time stepping and time filter
524	! time filter and array swap
525	sshb_tl(:,:) = atfp * ( sshb_tl(:,:) + zssha(:,:) ) + atfp1 * sshn_tl(:,:)
526	sshn_tl(:,:) = zssha(:,:)
527	ENDIF
528
529	! print sum trends (used for debugging)
530	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn_tl, clinfo1=' spg - ssh: ', mask1=tmask )
531	!
532
533	END SUBROUTINE dyn_spg_flt_tan
534
535	SUBROUTINE dyn_spg_flt_adj( kt, kindic )
536	!!----------------------------------------------------------------------
537	!! * routine dyn_spg_flt_adj *
538	!!
539	!! ** Purpose of the direct routine:
540	!! Compute the now trend due to the surface pressure
541	!! gradient in case of filtered free surface formulation and add
542	!! it to the general trend of momentum equation.
543	!!
544	!! ** Method of the direct routine:
545	!! Filtered free surface formulation. The surface
546	!! pressure gradient is given by:
547	!! spgu = 1/rau0 d/dx(ps) = 1/e1u di( sshn + rnu btda )
548	!! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + rnu btda )
549	!! where sshn is the free surface elevation and btda is the after
550	!! of the free surface elevation
551	!! -1- compute the after sea surface elevation from the kinematic
552	!! surface boundary condition:
553	!! zssha = sshb + 2 rdt ( wn - emp )
554	!! Time filter applied on now sea surface elevation to avoid
555	!! the divergence of two consecutive time-steps and swap of free
556	!! surface arrays to start the next time step:
557	!! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
558	!! sshn = zssha
559	!! -2- evaluate the surface presure trend (including the addi-
560	!! tional force) in three steps:
561	!! a- compute the right hand side of the elliptic equation:
562	!! gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
563	!! where (spgu,spgv) are given by:
564	!! spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
565	!! - grav 2 rdt hu /e1u di[sshn + emp]
566	!! spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
567	!! - grav 2 rdt hv /e2v dj[sshn + emp]
568	!! and define the first guess from previous computation :
569	!! zbtd = btda
570	!! btda = 2 zbtd - btdb
571	!! btdb = zbtd
572	!! b- compute the relative accuracy to be reached by the
573	!! iterative solver
574	!! c- apply the solver by a call to sol... routine
575	!! -3- compute and add the free surface pressure gradient inclu-
576	!! ding the additional force used to stabilize the equation.
577	!!
578	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
579	!!
580	!! References : Roullet and Madec 1999, JGR.
581	!!---------------------------------------------------------------------
582	!! * Arguments
583	INTEGER, INTENT( IN ) :: &
584	& kt ! ocean time-step index
585	INTEGER, INTENT( OUT ) :: &
586	& kindic ! solver convergence flag (<0 if not converge)
587
588	!! * Local declarations
589	INTEGER :: &
590	& ji, & ! dummy loop indices
591	& jj, &
592	& jk
593	REAL(wp) :: &
594	& z2dt, & ! temporary scalars
595	& z2dtg, &
596	& zraur, &
597	& znugdt, &
598	& znurau, &
599	& zgcb, &
600	& zbtd, &
601	& ztdgu, &
602	& ztdgv
603	!! Variable volume
604	REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace
605	& zsshub, &
606	& zsshua, &
607	& zsshvb, &
608	& zsshva, &
609	& zssha
610	REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace
611	& zfse3ub, &
612	& zfse3ua, &
613	& zfse3vb, &
614	& zfse3va
615	!!----------------------------------------------------------------------
616	!
617
618	! Local constant initialization
619	IF ( neuler == 0 .AND. kt == nit000 ) THEN
620
621	z2dt = rdt ! time step: Euler if restart from rest
622	CALL sol_mat(kt) ! initialize matrix
623
624	ELSEIF ( neuler == 0 .AND. kt == nit000 + 1 ) THEN
625
626	z2dt = 2.0_wp * rdt ! time step: leap-frog
627	CALL sol_mat(kt) ! reinitialize matrix
628
629	ELSEIF ( kt == nitend ) THEN
630
631	z2dt = 2.0_wp * rdt ! time step: leap-frog
632	CALL sol_mat(kt) ! reinitialize matrix
633
634	ELSEIF ( neuler /= 0 .AND. kt == nit000 ) THEN
635
636	z2dt = 2.0_wp * rdt ! time step: leap-frog
637	CALL sol_mat(kt) ! initialize matrix
638
639	ELSE
640
641	z2dt = 2.0_wp * rdt ! time step: leap-frog
642
643	ENDIF
644
645	z2dtg = grav * z2dt
646	zraur = 1.0_wp / rauw
647	znugdt = rnu * grav * z2dt
648	znurau = znugdt * zraur
649
650	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler (forward) time stepping, no time filter
651	! swap of arrays
652	zssha (:,:) = sshn_ad(:,:)
653	sshn_ad(:,:) = sshb_ad(:,:)
654	sshb_ad(:,:) = 0.0_wp
655	ELSE ! Leap-frog time stepping and time filter
656	! time filter and array swap
657	zssha (:,:) = sshn_ad(:,:) + atfp * sshb_ad(:,:)
658	sshn_ad(:,:) = atfp1 * sshb_ad(:,:)
659	sshb_ad(:,:) = atfp * sshb_ad(:,:)
660	ENDIF
661
662	! Sea surface elevation time stepping
663	! -----------------------------------
664	#if defined key_obc
665	CALL lbc_lnk_adj(zssha,'T',1.) ! absolutly compulsory !! (jmm)
666	zssha(:,:) = zssha(:,:) * obctmsk(:,:)
667	#endif
668	IF( lk_vvl ) THEN ! after free surface elevation
669	CALL ctl_stop( 'dyn_spg_flt_adj: lk_vvl is not available' )
670	ELSE
671	sshb_ad(:,:) = sshb_ad(:,:) + zssha(:,:)
672	zssha(:,:) = zssha(:,:) * z2dt * tmask(:,:,1)
673	wn_ad(:,:,1) = wn_ad(:,:,1) + zssha(:,:)
674	emp_ad(:,:) = emp_ad(:,:) - zraur * zssha(:,:)
675	zssha(:,:) = 0.0_wp
676	ENDIF
677
678	! 7. Add the trends multiplied by z2dt to the after velocity
679	! -----------------------------------------------------------
680	! ( c a u t i o n : (ua_ad,va_ad) here are the after velocity not the
681	! trend, the leap-frog time stepping will not
682	! be done in dynnxt_adj.F90 routine)
683	DO jk = jpkm1, 1, -1
684	DO jj = jpjm1, 2, -1
685	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
686	ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) * umask(ji,jj,jk)
687	va_ad(ji,jj,jk) = va_ad(ji,jj,jk) * vmask(ji,jj,jk)
688	spgu_ad(ji,jj) = spgu_ad(ji,jj) * umask(ji,jj,jk) + ua_ad(ji,jj,jk)
689	spgv_ad(ji,jj) = spgv_ad(ji,jj) * vmask(ji,jj,jk) + va_ad(ji,jj,jk)
690	END DO
691	END DO
692	END DO
693
694	#if defined key_agrif
695	CALL ctl_stop( 'dyn_spg_flt_adj: key_agrif is not available' )
696	#endif
697
698	! Transport divergence gradient multiplied by z2dt
699	! --------------------------------------------====
700	DO jj = jpjm1, 2, -1
701	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
702	! multiplied by z2dt
703	#if defined key_obc
704	! caution : grad D = 0 along open boundaries
705	ztdgu = z2dt * spgu_ad(ji,jj) * obcumask(ji,jj)
706	ztdgv = z2dt * spgv_ad(ji,jj) * obcvmask(ji,jj)
707	spgu_ad(ji,jj) = 0.0_wp
708	spgv_ad(ji,jj) = 0.0_wp
709	#elif defined key_bdy
710	CALL ctl_stop( 'dyn_spg_flt_adj: key_bdy is not available' )
711	#else
712	ztdgu = z2dt * spgu_ad(ji,jj)
713	ztdgv = z2dt * spgv_ad(ji,jj)
714	spgu_ad(ji,jj) = 0.0_wp
715	spgv_ad(ji,jj) = 0.0_wp
716	#endif
717	! trend of Transport divergence gradient
718	ztdgu = ztdgu * znugdt / e1u(ji,jj)
719	ztdgv = ztdgv * znugdt / e2v(ji,jj)
720	gcx_ad(ji ,jj ) = gcx_ad(ji ,jj ) - ztdgu - ztdgv
721	gcx_ad(ji ,jj+1) = gcx_ad(ji ,jj+1) + ztdgv
722	gcx_ad(ji+1,jj ) = gcx_ad(ji+1,jj ) + ztdgu
723	END DO
724	END DO
725
726	! Evaluate the next transport divergence
727	! --------------------------------------
728	! Iterative solver for the elliptic equation (except IF sol.=0)
729	! (output in gcx_ad with boundary conditions applied)
730	kindic = 0
731	ncut = 0 ! Force solver
732	IF( ncut == 0 ) THEN
733	IF( nsolv == 1 ) THEN ! diagonal preconditioned conjuguate gradient
734	CALL sol_pcg_adj( kt, kindic )
735	ELSEIF( nsolv == 2 ) THEN ! successive-over-relaxation
736	CALL sol_sor_adj( kt, kindic )
737	ELSEIF( nsolv == 3 ) THEN ! FETI solver
738	CALL sol_fet_adj( kt, kindic )
739	ELSE ! e r r o r in nsolv namelist parameter
740	WRITE(ctmp1,*) ' ~~~~~~~~~~~ not = ', nsolv
741	CALL ctl_stop( ' dyn_spg_flt_adj : e r r o r, nsolv = 2', ctmp1 )
742	ENDIF
743	ENDIF
744
745	#if defined key_agrif
746	CALL ctl_stop( 'dyn_spg_flt_adj: key_agrif is not available' )
747	#endif
748
749	! Right hand side of the elliptic equation and first guess
750	! -----------------------------------------------------------
751	! apply the lateral boundary conditions
752	IF( nsolv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) &
753	& CALL lbc_lnk_e_adj( gcb_ad, c_solver_pt, 1.0_wp )
754
755	DO jj = jpjm1, 2, -1
756	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
757	! First guess of the after barotropic transport divergence
758	zbtd = gcxb_ad(ji,jj) + 2.0_wp * gcx_ad(ji,jj)
759	gcxb_ad(ji,jj) = - gcx_ad(ji,jj)
760	gcx_ad (ji,jj) = zbtd
761	! Divergence of the after vertically averaged velocity
762	zgcb = gcb_ad(ji,jj) * gcdprc(ji,jj)
763	gcb_ad(ji,jj) = 0.0_wp
764	spgu_ad(ji-1,jj ) = spgu_ad(ji-1,jj ) - zgcb
765	spgu_ad(ji ,jj ) = spgu_ad(ji ,jj ) + zgcb
766	spgv_ad(ji ,jj-1) = spgv_ad(ji ,jj-1) - zgcb
767	spgv_ad(ji ,jj ) = spgv_ad(ji ,jj ) + zgcb
768	END DO
769	END DO
770
771	IF( lk_vvl ) CALL ctl_stop( 'dyn_spg_flt_adj: lk_vvl is not available' )
772
773	! Boundary conditions on (spgu_ad,spgv_ad)
774	CALL lbc_lnk_adj( spgu_ad, 'U', -1.0_wp )
775	CALL lbc_lnk_adj( spgv_ad, 'V', -1.0_wp )
776
777	! transport: multiplied by the horizontal scale factor
778	DO jj = jpjm1, 2, -1
779	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
780	spgu_ad(ji,jj) = spgu_ad(ji,jj) * e2u(ji,jj)
781	spgv_ad(ji,jj) = spgv_ad(ji,jj) * e1v(ji,jj)
782	END DO
783	END DO
784
785	! compute the next vertically averaged velocity (effect of the additional force not included)
786	! ---------------------------------------------
787
788	! vertical sum
789	!CDIR NOLOOPCHG
790	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
791	DO jk = jpkm1, 1, -1
792	DO ji = jpij, 1, -1
793	ua_ad(ji,1,jk) = ua_ad(ji,1,jk) + fse3u(ji,1,jk) * spgu_ad(ji,1)
794	va_ad(ji,1,jk) = va_ad(ji,1,jk) + fse3v(ji,1,jk) * spgv_ad(ji,1)
795	END DO
796	END DO
797	ELSE ! No vector opt.
798	DO jk = jpkm1, 1, -1
799	DO jj = jpjm1, 2, -1
800	DO ji = 2, jpim1
801	ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) + fse3u(ji,jj,jk) * spgu_ad(ji,jj)
802	va_ad(ji,jj,jk) = va_ad(ji,jj,jk) + fse3v(ji,jj,jk) * spgv_ad(ji,jj)
803	END DO
804	END DO
805	END DO
806	ENDIF
807
808	DO jj = jpjm1, 2, -1
809	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
810	spgu_ad(ji,jj) = 0.0_wp
811	spgv_ad(ji,jj) = 0.0_wp
812	END DO
813	END DO
814
815	#if defined key_orca_r2
816	IF( n_cla == 1 ) CALL dyn_spg_cla_adj( kt ) ! Cross Land Advection (update (ua,va))
817	#endif
818	#if defined key_agrif
819	CALL ctl_stop( 'dyn_spg_flt_adj: key_agrif is not available' )
820	#endif
821	#if defined key_bdy
822	CALL ctl_stop( 'dyn_spg_flt_adj: key_bdy is not available' )
823	#endif
824	#if defined key_obc
825	CALL obc_dyn_adj( kt ) ! Update velocities on each open boundary with the radiation algorithm
826	#endif
827
828	!! Explicit physics with thickness weighted updates
829	IF( lk_vvl ) THEN ! variable volume
830
831	CALL ctl_stop( 'dyn_spg_flt_adj: lk_vvl is not available' )
832
833	ELSE ! fixed volume
834
835	! Evaluate the masked next velocity (effect of the additional force not included)
836	DO jk = jpkm1, 1, -1
837	DO jj = jpjm1, 2, -1
838	DO ji = fs_2, fs_jpim1 ! vector opt.
839	ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) * umask(ji,jj,jk)
840	va_ad(ji,jj,jk) = va_ad(ji,jj,jk) * vmask(ji,jj,jk)
841	ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) + ua_ad(ji,jj,jk)
842	vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) + va_ad(ji,jj,jk)
843	ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) * z2dt
844	va_ad(ji,jj,jk) = va_ad(ji,jj,jk) * z2dt
845	END DO
846	END DO
847	END DO
848
849	! Add the surface pressure trend to the general trend
850	DO jk = jpkm1, 1, -1
851	DO jj = jpjm1, 2, -1
852	DO ji = fs_2, fs_jpim1 ! vector opt.
853	spgu_ad(ji,jj) = spgu_ad(ji,jj) + ua_ad(ji,jj,jk)
854	spgv_ad(ji,jj) = spgv_ad(ji,jj) + va_ad(ji,jj,jk)
855	END DO
856	END DO
857	END DO
858
859	! Surface pressure gradient (now)
860	DO jj = jpjm1, 2, -1
861	DO ji = fs_jpim1, fs_2, -1 ! vector opt.
862	spgu_ad(ji ,jj ) = spgu_ad(ji ,jj ) * grav / e1u(ji,jj)
863	spgv_ad(ji ,jj ) = spgv_ad(ji ,jj ) * grav / e2v(ji,jj)
864	sshn_ad(ji ,jj ) = sshn_ad(ji ,jj ) + spgu_ad(ji,jj)
865	sshn_ad(ji ,jj ) = sshn_ad(ji ,jj ) + spgv_ad(ji,jj)
866	sshn_ad(ji ,jj+1) = sshn_ad(ji ,jj+1) - spgv_ad(ji,jj)
867	sshn_ad(ji+1,jj ) = sshn_ad(ji+1,jj ) - spgu_ad(ji,jj)
868	spgu_ad(ji ,jj ) = 0.0_wp
869	spgv_ad(ji ,jj ) = 0.0_wp
870	END DO
871	END DO
872
873	ENDIF
874
875	IF( kt == nit000 ) THEN
876
877	IF(lwp) WRITE(numout,*)
878	IF(lwp) WRITE(numout,*) 'dyn_spg_flt_adj : surface pressure gradient trend'
879	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ (free surface constant volume case)'
880
881	! set to zero free surface specific arrays
882	spgu_ad(:,:) = 0.0_wp ! surface pressure gradient (i-direction)
883	spgv_ad(:,:) = 0.0_wp ! surface pressure gradient (j-direction)
884
885	ENDIF
886
887	END SUBROUTINE dyn_spg_flt_adj
888
889	SUBROUTINE dyn_spg_flt_adj_tst( kumadt )
890	!!-----------------------------------------------------------------------
891	!!
892	!! * ROUTINE dyn_spg_flt_adj_tst *
893	!!
894	!! ** Purpose : Test the adjoint routine.
895	!!
896	!! ** Method : Verify the scalar product
897	!!
898	!! ( L dx )^T W dy = dx^T L^T W dy
899	!!
900	!! where L = tangent routine
901	!! L^T = adjoint routine
902	!! W = diagonal matrix of scale factors
903	!! dx = input perturbation (random field)
904	!! dy = L dx
905	!!
906	!! ** Action :
907	!!
908	!! History :
909	!! ! 09-01 (A. Weaver)
910	!!-----------------------------------------------------------------------
911	!! * Modules used
912
913	!! * Arguments
914	INTEGER, INTENT(IN) :: &
915	& kumadt ! Output unit
916
917	!! * Local declarations
918	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
919	& zua_tlin, & ! Tangent input: ua_tl
920	& zva_tlin, & ! Tangent input: va_tl
921	& zub_tlin, & ! Tangent input: ub_tl
922	& zvb_tlin, & ! Tangent input: vb_tl
923	& zua_tlout, & ! Tangent output: ua_tl
924	& zva_tlout, & ! Tangent output: va_tl
925	#if defined key_obc
926	& zub_tlout, & ! Tangent output: ua_tl
927	& zvb_tlout, & ! Tangent output: va_tl
928	& zub_adin, & ! Tangent output: ua_tl
929	& zvb_adin, & ! Tangent output: va_tl
930	#endif
931	& zua_adin, & ! Adjoint input: ua_ad
932	& zva_adin, & ! Adjoint input: va_ad
933	& zua_adout, & ! Adjoint output: ua_ad
934	& zva_adout, & ! Adjoint output: va_ad
935	& zub_adout, & ! Adjoint oputput: ub_ad
936	& zvb_adout, & ! Adjoint output: vb_ad
937	& znu ! 3D random field for u
938	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
939	& zsshb_tlin, & ! Tangent input: sshb_tl
940	& zsshn_tlin, & ! Tangent input: sshn_tl
941	& zemp_tlin, & ! Tangent input: emp_tl
942	& zwn_tlin, & ! Tangent input: wn_tl
943	#if defined key_obc
944	& zemp_tlout, & ! Tangent input: emp_tl
945	& zwn_tlout, & ! Tangent input: wn_tl
946	& zemp_adin, & ! Adjoint output: emp_ad
947	#endif
948	& zsshb_tlout,& ! Tangent output: sshb_tl
949	& zsshn_tlout,& ! Tangent output: sshn_tl
950	& zsshb_adin, & ! Adjoint input: sshb_ad
951	& zsshn_adin, & ! Adjoint input: sshn_ad
952	& zsshb_adout,& ! Adjoint output: sshb_ad
953	& zsshn_adout,& ! Adjoint output: sshn_ad
954	& zemp_adout, & ! Adjoint output: emp_ad
955	& zwn_adout, & ! Adjoint output: wn_ad
956	#if defined key_obc
957	& zgcx_tlin, zgcxb_tlin, zgcb_tlin, zgcx_tlout, zgcxb_tlout, zgcb_tlout, &
958	& zgcx_adin, zgcxb_adin, zgcb_adin, zgcx_adout, zgcxb_adout, zgcb_adout, &
959	& zwn_adin, &
960	#endif
961	& znssh ! 2D random field for SSH
962	REAL(KIND=wp) :: &
963	& zsp1, & ! scalar product involving the tangent routine
964	& zsp2 ! scalar product involving the adjoint routine
965	INTEGER, DIMENSION(jpi,jpj) :: &
966	& iseed_2d ! 2D seed for the random number generator
967	INTEGER :: &
968	& indic, &
969	& istp
970	INTEGER :: &
971	& ji, &
972	& jj, &
973	& jk, &
974	& jstp
975	CHARACTER (LEN=14) :: &
976	& cl_name
977
978	! Allocate memory
979
980	ALLOCATE( &
981	& zua_tlin(jpi,jpj,jpk), &
982	& zva_tlin(jpi,jpj,jpk), &
983	& zub_tlin(jpi,jpj,jpk), &
984	& zvb_tlin(jpi,jpj,jpk), &
985	& zua_tlout(jpi,jpj,jpk), &
986	& zva_tlout(jpi,jpj,jpk), &
987	& zua_adin(jpi,jpj,jpk), &
988	& zva_adin(jpi,jpj,jpk), &
989	& zua_adout(jpi,jpj,jpk), &
990	& zva_adout(jpi,jpj,jpk), &
991	& zub_adout(jpi,jpj,jpk), &
992	& zvb_adout(jpi,jpj,jpk), &
993	& znu(jpi,jpj,jpk) &
994	& )
995	ALLOCATE( &
996	& zsshb_tlin(jpi,jpj), &
997	& zsshn_tlin(jpi,jpj), &
998	& zemp_tlin(jpi,jpj), &
999	& zwn_tlin(jpi,jpj), &
1000	& zsshb_tlout(jpi,jpj),&
1001	& zsshn_tlout(jpi,jpj),&
1002	& zsshb_adin(jpi,jpj), &
1003	& zsshn_adin(jpi,jpj), &
1004	& zsshb_adout(jpi,jpj),&
1005	& zsshn_adout(jpi,jpj),&
1006	& zemp_adout(jpi,jpj), &
1007	& zwn_adout(jpi,jpj), &
1008	& znssh(jpi,jpj) &
1009	& )
1010
1011	#if defined key_obc
1012	ALLOCATE( zgcx_tlin (jpi,jpj), zgcx_tlout (jpi,jpj), zgcx_adin (jpi,jpj), zgcx_adout (jpi,jpj), &
1013	zgcxb_tlin(jpi,jpj), zgcxb_tlout(jpi,jpj), zgcxb_adin(jpi,jpj), zgcxb_adout(jpi,jpj), &
1014	zgcb_tlin (jpi,jpj), zgcb_tlout (jpi,jpj), zgcb_adin (jpi,jpj), zgcb_adout (jpi,jpj), &
1015	zemp_tlout(jpi,jpj), zemp_adin (jpi,jpj), &
1016	zwn_tlout (jpi,jpj), zwn_adin (jpi,jpj) )
1017
1018	ALLOCATE ( zub_tlout(jpi,jpj,jpk), zvb_tlout(jpi,jpj,jpk), &
1019	zub_adin (jpi,jpj,jpk), zvb_adin (jpi,jpj,jpk) )
1020	#endif
1021
1022	!=========================================================================
1023	! dx = ( sshb_tl, sshn_tl, ub_tl, ua_tl, vb_tl, va_tl, wn_tl, emp_tl )
1024	! and dy = ( sshb_tl, sshn_tl, ua_tl, va_tl )
1025	!=========================================================================
1026
1027	! Test for time steps nit000 and nit000 + 1 (the matrix changes)
1028
1029	DO jstp = nit000, nit000 + 1
1030	istp = jstp
1031
1032	!--------------------------------------------------------------------
1033	! Reset the tangent and adjoint variables
1034	!--------------------------------------------------------------------
1035
1036	zub_tlin (:,:,:) = 0.0_wp
1037	zvb_tlin (:,:,:) = 0.0_wp
1038	zua_tlin (:,:,:) = 0.0_wp
1039	zva_tlin (:,:,:) = 0.0_wp
1040	zua_tlout(:,:,:) = 0.0_wp
1041	zva_tlout(:,:,:) = 0.0_wp
1042	zua_adin (:,:,:) = 0.0_wp
1043	zva_adin (:,:,:) = 0.0_wp
1044	zub_adout(:,:,:) = 0.0_wp
1045	zvb_adout(:,:,:) = 0.0_wp
1046	zua_adout(:,:,:) = 0.0_wp
1047	zva_adout(:,:,:) = 0.0_wp
1048
1049	zsshb_tlin (:,:) = 0.0_wp
1050	zsshn_tlin (:,:) = 0.0_wp
1051	zwn_tlin (:,:) = 0.0_wp
1052	zemp_tlin (:,:) = 0.0_wp
1053	zsshb_tlout(:,:) = 0.0_wp
1054	zsshn_tlout(:,:) = 0.0_wp
1055	zsshb_adin (:,:) = 0.0_wp
1056	zsshn_adin (:,:) = 0.0_wp
1057	zsshb_adout(:,:) = 0.0_wp
1058	zsshn_adout(:,:) = 0.0_wp
1059	zwn_adout (:,:) = 0.0_wp
1060	zemp_adout (:,:) = 0.0_wp
1061
1062	!--------------------------------------------------------------------
1063	! Initialize the tangent input with random noise: dx
1064	!--------------------------------------------------------------------
1065
1066	DO jj = 1, jpj
1067	DO ji = 1, jpi
1068	iseed_2d(ji,jj) = - ( 456953 + &
1069	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1070	END DO
1071	END DO
1072	CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu )
1073
1074	DO jk = 1, jpk
1075	DO jj = nldj, nlej
1076	DO ji = nldi, nlei
1077	zua_tlin(ji,jj,jk) = znu(ji,jj,jk)
1078	END DO
1079	END DO
1080	END DO
1081
1082	DO jj = 1, jpj
1083	DO ji = 1, jpi
1084	iseed_2d(ji,jj) = - ( 3434334 + &
1085	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1086	END DO
1087	END DO
1088	CALL grid_random( iseed_2d, znu, 'V', 0.0_wp, stdv )
1089
1090	DO jk = 1, jpk
1091	DO jj = nldj, nlej
1092	DO ji = nldi, nlei
1093	zva_tlin(ji,jj,jk) = znu(ji,jj,jk)
1094	END DO
1095	END DO
1096	END DO
1097
1098	DO jj = 1, jpj
1099	DO ji = 1, jpi
1100	iseed_2d(ji,jj) = - ( 935678 + &
1101	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1102	END DO
1103	END DO
1104	CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu )
1105
1106	DO jk = 1, jpk
1107	DO jj = nldj, nlej
1108	DO ji = nldi, nlei
1109	zub_tlin(ji,jj,jk) = znu(ji,jj,jk)
1110	END DO
1111	END DO
1112	END DO
1113
1114	DO jj = 1, jpj
1115	DO ji = 1, jpi
1116	iseed_2d(ji,jj) = - ( 1462578 + &
1117	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1118	END DO
1119	END DO
1120	CALL grid_random( iseed_2d, znu, 'V', 0.0_wp, stdv )
1121
1122	DO jk = 1, jpk
1123	DO jj = nldj, nlej
1124	DO ji = nldi, nlei
1125	zvb_tlin(ji,jj,jk) = znu(ji,jj,jk)
1126	END DO
1127	END DO
1128	END DO
1129
1130	DO jj = 1, jpj
1131	DO ji = 1, jpi
1132	iseed_2d(ji,jj) = - ( 285607 + &
1133	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1134	END DO
1135	END DO
1136	CALL grid_random( iseed_2d, znu, 'T', 0.0_wp, stdw )
1137
1138	DO jj = nldj, nlej
1139	DO ji = nldi, nlei
1140	zwn_tlin(ji,jj) = znu(ji,jj,1)
1141	END DO
1142	END DO
1143
1144	DO jj = 1, jpj
1145	DO ji = 1, jpi
1146	iseed_2d(ji,jj) = - ( 25674910 + &
1147	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1148	END DO
1149	END DO
1150	CALL grid_random( iseed_2d, znu, 'T', 0.0_wp, stdemp )
1151
1152	DO jj = nldj, nlej
1153	DO ji = nldi, nlei
1154	zemp_tlin(ji,jj) = znu(ji,jj,1)
1155	END DO
1156	END DO
1157
1158	DO jj = 1, jpj
1159	DO ji = 1, jpi
1160	iseed_2d(ji,jj) = - ( 93756381 + &
1161	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1162	END DO
1163	END DO
1164	CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh )
1165
1166	DO jj = nldj, nlej
1167	DO ji = nldi, nlei
1168	zsshb_tlin(ji,jj) = znssh(ji,jj)
1169	END DO
1170	END DO
1171
1172	DO jj = 1, jpj
1173	DO ji = 1, jpi
1174	iseed_2d(ji,jj) = - ( 12672456 + &
1175	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1176	END DO
1177	END DO
1178	CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh )
1179
1180	DO jj = nldj, nlej
1181	DO ji = nldi, nlei
1182	zsshn_tlin(ji,jj) = znssh(ji,jj)
1183	END DO
1184	END DO
1185
1186	#if defined key_obc
1187	zgcx_tlin (:,:) = ( zua_tlin(:,:,1) + zub_tlin(:,:,1) ) / 10.
1188	zgcxb_tlin (:,:) = ( zua_tlin(:,:,2) + zub_tlin(:,:,2) ) / 10.
1189	zgcb_tlin (:,:) = ( zua_tlin(:,:,3) + zub_tlin(:,:,3) ) / 10.
1190	#endif
1191	!--------------------------------------------------------------------
1192	! Call the tangent routine: dy = L dx
1193	!--------------------------------------------------------------------
1194
1195	ua_tl(:,:,:) = zua_tlin(:,:,:)
1196	va_tl(:,:,:) = zva_tlin(:,:,:)
1197	ub_tl(:,:,:) = zub_tlin(:,:,:)
1198	vb_tl(:,:,:) = zvb_tlin(:,:,:)
1199	wn_tl(:,:,1) = zwn_tlin (:,:)
1200	emp_tl (:,:) = zemp_tlin (:,:)
1201	sshb_tl(:,:) = zsshb_tlin(:,:)
1202	sshn_tl(:,:) = zsshn_tlin(:,:)
1203
1204	#if defined key_obc
1205	gcx_tl (:,:) = zgcx_tlin (:,:) ; gcxb_tl(:,:) = zgcxb_tlin(:,:)
1206	gcb_tl (:,:) = zgcb_tlin (:,:)
1207	#else
1208	gcx_tl (:,:) = 0.e0 ; gcxb_tl(:,:) = 0.e0
1209	gcb_tl (:,:) = 0.e0
1210	#endif
1211
1212	CALL dyn_spg_flt_tan( istp, indic )
1213
1214	zua_tlout(:,:,:) = ua_tl(:,:,:) ; zva_tlout(:,:,:) = va_tl(:,:,:)
1215	zsshb_tlout(:,:) = sshb_tl(:,:) ; zsshn_tlout(:,:) = sshn_tl(:,:)
1216	#if defined key_obc
1217	zub_tlout(:,:,:) = ub_tl(:,:,:) ; zvb_tlout(:,:,:) = vb_tl(:,:,:)
1218	zgcx_tlout (:,:) = gcx_tl (:,:) ; zgcxb_tlout(:,:) = gcxb_tl(:,:)
1219	zgcb_tlout (:,:) = gcb_tl (:,:)
1220	zwn_tlout (:,:) = wn_tl (:,:,1)
1221	zemp_tlout (:,:) = emp_tl (:,:)
1222	#endif
1223
1224	!--------------------------------------------------------------------
1225	! Initialize the adjoint variables: dy^* = W dy
1226	!--------------------------------------------------------------------
1227
1228	DO jk = 1, jpk
1229	DO jj = nldj, nlej
1230	DO ji = nldi, nlei
1231	zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
1232	& * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
1233	& * umask(ji,jj,jk)
1234	zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
1235	& * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
1236	& * vmask(ji,jj,jk)
1237	#if defined key_obc
1238	zub_adin(ji,jj,jk) = zub_tlout(ji,jj,jk) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
1239	zvb_adin(ji,jj,jk) = zvb_tlout(ji,jj,jk) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)
1240	#endif
1241	END DO
1242	END DO
1243	END DO
1244	DO jj = nldj, nlej
1245	DO ji = nldi, nlei
1246	zsshb_adin(ji,jj) = zsshb_tlout(ji,jj) &
1247	& * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
1248	& * tmask(ji,jj,1)
1249	zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) &
1250	& * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
1251	& * tmask(ji,jj,1)
1252	END DO
1253	END DO
1254
1255	#if defined key_obc
1256	DO jj = nldj, nlej
1257	DO ji = nldi, nlei
1258	zgcx_adin (ji,jj) = zgcx_tlout (ji,jj) * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,jk) * tmask(ji,jj,1)
1259	zgcxb_adin(ji,jj) = zgcxb_tlout(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,jk) * tmask(ji,jj,1)
1260	zgcb_adin (ji,jj) = zgcb_tlout (ji,jj) * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,jk) * tmask(ji,jj,1)
1261	zwn_adin (ji,jj) = zwn_tlout (ji,jj) * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,jk) * tmask(ji,jj,1)
1262	zemp_adin (ji,jj) = zemp_tlout (ji,jj) * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,jk) * tmask(ji,jj,1)
1263	END DO
1264	END DO
1265	#endif
1266
1267	!--------------------------------------------------------------------
1268	! Compute the scalar product: ( L dx )^T W dy
1269	!--------------------------------------------------------------------
1270
1271	zsp1 = DOT_PRODUCT( zua_tlout , zua_adin ) &
1272	& + DOT_PRODUCT( zva_tlout , zva_adin ) &
1273	#if defined key_obc
1274	& + DOT_PRODUCT( zgcx_tlout , zgcx_adin ) &
1275	& + DOT_PRODUCT( zgcxb_tlout, zgcxb_adin ) &
1276	& + DOT_PRODUCT( zgcb_tlout , zgcb_adin ) &
1277	& + DOT_PRODUCT( zub_tlout , zub_adin ) &
1278	& + DOT_PRODUCT( zvb_tlout , zvb_adin ) &
1279	& + DOT_PRODUCT( zwn_tlout , zwn_adin ) &
1280	& + DOT_PRODUCT( zemp_tlout , zemp_adin ) &
1281	#endif
1282	& + DOT_PRODUCT( zsshb_tlout, zsshb_adin ) &
1283	& + DOT_PRODUCT( zsshn_tlout, zsshn_adin )
1284
1285	!--------------------------------------------------------------------
1286	! Call the adjoint routine: dx^* = L^T dy^*
1287	!--------------------------------------------------------------------
1288
1289	ua_ad(:,:,:) = zua_adin(:,:,:)
1290	va_ad(:,:,:) = zva_adin(:,:,:)
1291	sshb_ad(:,:) = zsshb_adin(:,:)
1292	sshn_ad(:,:) = zsshn_adin(:,:)
1293
1294	#if defined key_obc
1295	gcx_ad (:,:) = zgcx_adin (:,:) ; gcxb_ad(:,:) = zgcxb_adin(:,:)
1296	gcb_ad (:,:) = zgcb_adin (:,:)
1297	ub_ad (:,:,:) = zub_adin (:,:,:) ; vb_ad (:,:,:) = zvb_adin (:,:,:)
1298	wn_ad (:,:,1) = zwn_adin (:,:) ; emp_ad (:,:) = zemp_adin (:,:)
1299	#else
1300	ub_ad(:,:,:) = 0.0_wp
1301	vb_ad(:,:,:) = 0.0_wp
1302	wn_ad(:,:,:) = 0.0_wp
1303	emp_ad (:,:) = 0.0_wp
1304	gcb_ad (:,:) = 0.0_wp
1305	gcx_ad (:,:) = 0.0_wp
1306	gcxb_ad(:,:) = 0.0_wp
1307	#endif
1308
1309	CALL dyn_spg_flt_adj( istp, indic )
1310
1311	zua_adout(:,:,:) = ua_ad(:,:,:)
1312	zva_adout(:,:,:) = va_ad(:,:,:)
1313	zub_adout(:,:,:) = ub_ad(:,:,:)
1314	zvb_adout(:,:,:) = vb_ad(:,:,:)
1315	zwn_adout (:,:) = wn_ad(:,:,1)
1316	zemp_adout (:,:) = emp_ad (:,:)
1317	zsshb_adout(:,:) = sshb_ad(:,:)
1318	zsshn_adout(:,:) = sshn_ad(:,:)
1319	#if defined key_obc
1320	zgcx_adout (:,:) = gcx_ad (:,:)
1321	zgcxb_adout(:,:) = gcxb_ad(:,:)
1322	zgcb_adout (:,:) = gcb_ad (:,:)
1323	#endif
1324
1325	!--------------------------------------------------------------------
1326	! Compute the scalar product: dx^T L^T W dy
1327	!--------------------------------------------------------------------
1328
1329	zsp2 = DOT_PRODUCT( zua_tlin , zua_adout ) &
1330	& + DOT_PRODUCT( zva_tlin , zva_adout ) &
1331	& + DOT_PRODUCT( zub_tlin , zub_adout ) &
1332	& + DOT_PRODUCT( zvb_tlin , zvb_adout ) &
1333	#if defined key_obc
1334	& + DOT_PRODUCT( zgcx_tlin , zgcx_adout ) &
1335	& + DOT_PRODUCT( zgcxb_tlin, zgcxb_adout ) &
1336	& + DOT_PRODUCT( zgcb_tlin , zgcb_adout ) &
1337	#endif
1338	& + DOT_PRODUCT( zsshb_tlin, zsshb_adout ) &
1339	& + DOT_PRODUCT( zsshn_tlin, zsshn_adout ) &
1340	& + DOT_PRODUCT( zwn_tlin , zwn_adout ) &
1341	& + DOT_PRODUCT( zemp_tlin , zemp_adout )
1342
1343	! Compare the scalar products
1344
1345	! 14 char:'12345678901234'
1346	IF ( istp == nit000 ) THEN
1347	cl_name = 'dyn_spg_flt T1'
1348	ELSEIF ( istp == nit000 + 1 ) THEN
1349	cl_name = 'dyn_spg_flt T2'
1350	END IF
1351	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1352
1353	END DO
1354
1355	! Deallocate memory
1356
1357	DEALLOCATE( &
1358	& zua_tlin, &
1359	& zva_tlin, &
1360	& zub_tlin, &
1361	& zvb_tlin, &
1362	& zua_tlout, &
1363	& zva_tlout, &
1364	& zua_adin, &
1365	& zva_adin, &
1366	& zua_adout, &
1367	& zva_adout, &
1368	& zub_adout, &
1369	& zvb_adout, &
1370	& znu &
1371	& )
1372	DEALLOCATE( &
1373	& zsshb_tlin, &
1374	& zsshn_tlin, &
1375	& zemp_tlin, &
1376	& zwn_tlin, &
1377	& zsshb_tlout,&
1378	& zsshn_tlout,&
1379	& zsshb_adin, &
1380	& zsshn_adin, &
1381	& zsshb_adout,&
1382	& zsshn_adout,&
1383	& zemp_adout, &
1384	& zwn_adout, &
1385	& znssh &
1386	& )
1387
1388	#if defined key_obc
1389	DEALLOCATE( zgcx_tlin , zgcx_tlout , zgcx_adin , zgcx_adout , &
1390	zgcxb_tlin, zgcxb_tlout, zgcxb_adin, zgcxb_adout, &
1391	zgcb_tlin , zgcb_tlout , zgcb_adin , zgcb_adout , &
1392	zemp_tlout, zemp_adin , &
1393	zwn_tlout , zwn_adin )
1394
1395	DEALLOCATE ( zub_tlout, zvb_tlout, zub_adin , zvb_adin )
1396	#endif
1397	END SUBROUTINE dyn_spg_flt_adj_tst
1398
1399	#if defined key_tst_tlm
1400	SUBROUTINE dyn_spg_flt_tlm_tst( kumadt )
1401	!!-----------------------------------------------------------------------
1402	!!
1403	!! * ROUTINE dyn_spg_flt_tlm_tst *
1404	!!
1405	!! ** Purpose : Test the tangent routine.
1406	!!
1407	!! ** Method : Verify the tangent with Taylor expansion
1408	!!
1409	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
1410	!!
1411	!! where L = tangent routine
1412	!! M = direct routine
1413	!! dx = input perturbation (random field)
1414	!! h = ration on perturbation
1415	!!
1416	!! In the tangent test we verify that:
1417	!! M(x+h*dx) - M(x)
1418	!! g(h) = ------------------ ---> 1 as h ---> 0
1419	!! L(h*dx)
1420	!! and
1421	!! g(h) - 1
1422	!! f(h) = ---------- ---> k (costant) as h ---> 0
1423	!! p
1424	!!
1425	!! History :
1426	!! ! 09-08 (A. Vigilant)
1427	!!-----------------------------------------------------------------------
1428	!! * Modules used
1429	USE opatam_tst_ini, ONLY: &
1430	& tlm_namrd
1431	USE dynspg_flt
1432	USE tamtrj ! writing out state trajectory
1433	USE par_tlm, ONLY: &
1434	& tlm_bch, &
1435	& cur_loop, &
1436	& h_ratio
1437	USE istate_mod
1438	USE gridrandom, ONLY: &
1439	& grid_rd_sd
1440	USE trj_tam
1441	USE oce , ONLY: & ! ocean dynamics and tracers variables
1442	& ua, va, ub, vb, &
1443	& un, vn, &
1444	& sshb, sshn, wn
1445	USE sbc_oce , ONLY: &
1446	& emp
1447	USE sol_oce , ONLY: & ! ocean dynamics and tracers variables
1448	& gcb, gcx
1449	USE tamctl, ONLY: & ! Control parameters
1450	& numtan, numtan_sc
1451	!! * Arguments
1452	INTEGER, INTENT(IN) :: &
1453	& kumadt ! Output unit
1454
1455	!! * Local declarations
1456	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1457	& zua_tlin, & ! Tangent input: ua_tl
1458	& zva_tlin, & ! Tangent input: va_tl
1459	& zub_tlin, & ! Tangent input: ub_tl
1460	& zvb_tlin, & !
1461	& zwn_tlin, & ! Tangent input: wn_tl
1462	& zua_out, & !
1463	& zva_out, & !
1464	& zua_wop, & !
1465	& zva_wop, &
1466	& z3r ! 3D field
1467
1468	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
1469	& zsshb_tlin, & ! Tangent input: sshb_tl
1470	& zsshn_tlin, & ! Tangent input: sshn_tl
1471	& zemp_tlin, & ! Tangent input: emp_tl
1472	& zsshb_out, & !
1473	& zsshn_out, & !
1474	& zsshb_wop, & !
1475	& zsshn_wop, &
1476	& z2r ! 2D field
1477	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
1478	& zgcb_tlin , & ! Tangent input
1479	& zgcx_tlin , & ! Tangent input
1480	& zgcb_out , & ! Direct output
1481	& zgcx_out , & ! Direct output
1482	& zgcb_wop , & ! Direct output without perturbation
1483	& zgcx_wop , & ! Direct output without perturbation
1484	& zr ! 3D random field
1485	REAL(KIND=wp) :: &
1486	& zsp1, zsp1_1, zsp1_2, zsp1_3, zsp1_4, & !
1487	& zsp2, zsp2_1, zsp2_2, zsp2_3, zsp2_4, & !
1488	& zsp3, zsp3_1, zsp3_2, zsp3_3, zsp3_4, & !
1489	& zzsp, zzsp_1, zzsp_2, zzsp_3, zzsp_4, &
1490	& gamma, &
1491	& zsp_5,zsp1_5, zsp2_5, zsp3_5, zsp4_5, &
1492	& zzsp_5, zsp_6, &
1493	& zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, &
1494	& zgsp6, zgsp7
1495	INTEGER :: &
1496	& indic, &
1497	& istp
1498	INTEGER :: &
1499	& ji, &
1500	& jj, &
1501	& jk
1502	CHARACTER(LEN=14) :: cl_name
1503	CHARACTER (LEN=128) :: file_out, file_wop, file_xdx
1504	CHARACTER (LEN=90) :: FMT
1505	REAL(KIND=wp), DIMENSION(100):: &
1506	& zscua, zscva, &
1507	& zscerrua, zscerrva, &
1508	& zscsshb, zscsshn, &
1509	& zscerrsshb, zscerrsshn
1510	REAL(KIND=wp), DIMENSION(jpi,jpj) :: &
1511	& zerrgcb, zerrgcx
1512	REAL(KIND=wp), DIMENSION(100):: &
1513	& zscgcb,zscgcx, &
1514	& zscerrgcb, zscerrgcx
1515	INTEGER, DIMENSION(100):: &
1516	& iiposgcb, ijposgcb, &
1517	& iiposgcx, ijposgcx
1518	INTEGER, DIMENSION(100):: &
1519	& iipossshb, iipossshn, iiposua, iiposva, &
1520	& ijpossshb, ijpossshn, ijposua, ijposva, &
1521	& ikposua, ikposva
1522	INTEGER:: &
1523	& ii, &
1524	& isamp=40, &
1525	& jsamp=40, &
1526	& ksamp=40, &
1527	& numsctlm
1528	REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: &
1529	& zerrua, zerrva
1530	REAL(KIND=wp), DIMENSION(jpi,jpj) :: &
1531	& zerrsshb, zerrsshn
1532	! Allocate memory
1533
1534	ALLOCATE( &
1535	& zua_tlin(jpi,jpj,jpk), &
1536	& zva_tlin(jpi,jpj,jpk), &
1537	& zub_tlin(jpi,jpj,jpk), &
1538	& zvb_tlin(jpi,jpj,jpk), &
1539	& zwn_tlin(jpi,jpj,jpk), &
1540	& zua_out( jpi,jpj,jpk), &
1541	& zva_out( jpi,jpj,jpk), &
1542	& zua_wop( jpi,jpj,jpk), &
1543	& zva_wop( jpi,jpj,jpk), &
1544	& z3r(jpi,jpj,jpk) &
1545	& )
1546	ALLOCATE( &
1547	& zsshb_tlin(jpi,jpj), &
1548	& zsshn_tlin(jpi,jpj), &
1549	& zemp_tlin(jpi,jpj), &
1550	& zsshb_out(jpi,jpj), &
1551	& zsshn_out(jpi,jpj), &
1552	& zsshb_wop(jpi,jpj), &
1553	& zsshn_wop(jpi,jpj), &
1554	& z2r(jpi,jpj) &
1555	& )
1556	ALLOCATE( &
1557	& zgcb_tlin( jpi,jpj), &
1558	& zgcx_tlin( jpi,jpj), &
1559	& zgcb_out ( jpi,jpj), &
1560	& zgcx_out ( jpi,jpj), &
1561	& zgcb_wop ( jpi,jpj), &
1562	& zgcx_wop ( jpi,jpj), &
1563	& zr( jpi,jpj) &
1564	& )
1565	!--------------------------------------------------------------------
1566	! Reset variables
1567	!--------------------------------------------------------------------
1568	zua_tlin( :,:,:) = 0.0_wp
1569	zva_tlin( :,:,:) = 0.0_wp
1570	zub_tlin( :,:,:) = 0.0_wp
1571	zvb_tlin( :,:,:) = 0.0_wp
1572	zwn_tlin( :,:,:) = 0.0_wp
1573	zsshb_tlin( :,:) = 0.0_wp
1574	zsshn_tlin( :,:) = 0.0_wp
1575	zemp_tlin( :,:) = 0.0_wp
1576	zua_out ( :,:,:) = 0.0_wp
1577	zva_out ( :,:,:) = 0.0_wp
1578	zsshb_out ( :,:) = 0.0_wp
1579	zsshn_out ( :,:) = 0.0_wp
1580	zua_wop ( :,:,:) = 0.0_wp
1581	zva_wop ( :,:,:) = 0.0_wp
1582	zsshb_wop( :,:) = 0.0_wp
1583	zsshn_wop( :,:) = 0.0_wp
1584
1585	zscua(:) = 0.0_wp
1586	zscva(:) = 0.0_wp
1587	zscerrua(:) = 0.0_wp
1588	zscerrva(:) = 0.0_wp
1589	zerrua(:,:,:) = 0.0_wp
1590	zerrva(:,:,:) = 0.0_wp
1591	zscsshb(:) = 0.0_wp
1592	zscsshn(:) = 0.0_wp
1593	zscerrsshb(:) = 0.0_wp
1594	zscerrsshn(:) = 0.0_wp
1595	zerrsshb(:,:) = 0.0_wp
1596	zerrsshn(:,:) = 0.0_wp
1597
1598	zgcb_tlin( :,:) = 0.0_wp
1599	zgcx_tlin( :,:) = 0.0_wp
1600	zgcb_out ( :,:) = 0.0_wp
1601	zgcx_out ( :,:) = 0.0_wp
1602	zgcb_wop ( :,:) = 0.0_wp
1603	zgcx_wop ( :,:) = 0.0_wp
1604	zr( :,:) = 0.0_wp
1605	!--------------------------------------------------------------------
1606	! Output filename Xn=F(X0)
1607	!--------------------------------------------------------------------
1608	CALL tlm_namrd
1609	gamma = h_ratio
1610	file_wop='trj_wop_dynspg'
1611	file_xdx='trj_xdx_dynspg'
1612	!--------------------------------------------------------------------
1613	! Initialize the tangent input with random noise: dx
1614	!--------------------------------------------------------------------
1615	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
1616	CALL grid_rd_sd( 456953, z3r, 'U', 0.0_wp, stdu)
1617	DO jk = 1, jpk
1618	DO jj = nldj, nlej
1619	DO ji = nldi, nlei
1620	zua_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1621	END DO
1622	END DO
1623	END DO
1624	CALL grid_rd_sd( 3434334, z3r, 'V', 0.0_wp, stdv)
1625	DO jk = 1, jpk
1626	DO jj = nldj, nlej
1627	DO ji = nldi, nlei
1628	zva_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1629	END DO
1630	END DO
1631	END DO
1632	CALL grid_rd_sd( 935678, z3r, 'U', 0.0_wp, stdu)
1633	DO jk = 1, jpk
1634	DO jj = nldj, nlej
1635	DO ji = nldi, nlei
1636	zub_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1637	END DO
1638	END DO
1639	END DO
1640	CALL grid_rd_sd( 1462578, z3r, 'V', 0.0_wp, stdv)
1641	DO jk = 1, jpk
1642	DO jj = nldj, nlej
1643	DO ji = nldi, nlei
1644	zvb_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1645	END DO
1646	END DO
1647	END DO
1648	CALL grid_rd_sd( 285607, z3r, 'T', 0.0_wp, stdw)
1649	DO jk = 1, jpk
1650	DO jj = nldj, nlej
1651	DO ji = nldi, nlei
1652	zwn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1653	END DO
1654	END DO
1655	END DO
1656	CALL grid_rd_sd( 25674910, z2r, 'T', 0.0_wp, stdemp)
1657	DO jj = nldj, nlej
1658	DO ji = nldi, nlei
1659	zemp_tlin(ji,jj) = z2r(ji,jj)
1660	END DO
1661	END DO
1662	CALL grid_rd_sd( 93756381, z2r,'T', 0.0_wp, stdssh)
1663	DO jj = nldj, nlej
1664	DO ji = nldi, nlei
1665	zsshb_tlin(ji,jj) = z2r(ji,jj)
1666	END DO
1667	END DO
1668	CALL grid_rd_sd( 12672456, z2r,'T', 0.0_wp, stdssh)
1669	DO jj = nldj, nlej
1670	DO ji = nldi, nlei
1671	zsshn_tlin(ji,jj) = z2r(ji,jj)
1672	END DO
1673	END DO
1674	CALL grid_rd_sd( 596035, zr, c_solver_pt, 0.0_wp, stdgc)
1675	DO jj = nldj, nlej
1676	DO ji = nldi, nlei
1677	zgcb_tlin(ji,jj) = zr(ji,jj)
1678	END DO
1679	END DO
1680	CALL grid_rd_sd( 264792, zr, c_solver_pt, 0.0_wp, stdgc)
1681	DO jj = nldj, nlej
1682	DO ji = nldi, nlei
1683	zgcx_tlin(ji,jj) = zr(ji,jj)
1684	END DO
1685	END DO
1686	ENDIF
1687
1688	!--------------------------------------------------------------------
1689	! Complete Init for Direct
1690	!-------------------------------------------------------------------
1691	CALL istate_p
1692	! *** initialize the reference trajectory
1693	! ------------
1694	CALL trj_rea( nit000-1, 1 )
1695	CALL trj_rea( nit000, 1 )
1696	ua(:,:,:)=un(:,:,:)
1697	va(:,:,:)=vn(:,:,:)
1698	ub(:,:,:)=un(:,:,:)
1699	vb(:,:,:)=vn(:,:,:)
1700	gcx (:,:) = ua(:,:,1) / 10.0_wp
1701	gcb (:,:) = ua(:,:,3) / 10.0_wp
1702
1703	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
1704	zua_tlin(:,:,:) = gamma * zua_tlin(:,:,:)
1705	ua(:,:,:) = ua(:,:,:) + zua_tlin(:,:,:)
1706
1707	zva_tlin(:,:,:) = gamma * zva_tlin(:,:,:)
1708	va(:,:,:) = va(:,:,:) + zva_tlin(:,:,:)
1709
1710	zub_tlin(:,:,:) = gamma * zub_tlin(:,:,:)
1711	ub(:,:,:) = ub(:,:,:) + zub_tlin(:,:,:)
1712
1713	zvb_tlin(:,:,:) = gamma * zvb_tlin(:,:,:)
1714	vb(:,:,:) = vb(:,:,:) + zvb_tlin(:,:,:)
1715
1716	zwn_tlin(:,:,:) = gamma * zwn_tlin(:,:,:)
1717	wn(:,:,:) = wn(:,:,:) + zwn_tlin(:,:,:)
1718
1719	zemp_tlin(:,:) = gamma * zemp_tlin(:,:)
1720	emp(:,:) = emp(:,:) + zemp_tlin(:,:)
1721
1722	zsshb_tlin(:,:) = gamma * zsshb_tlin(:,:)
1723	sshb(:,:) = sshb(:,:) + zsshb_tlin(:,:)
1724
1725	zsshn_tlin(:,:) = gamma * zsshn_tlin(:,:)
1726	sshn(:,:) = sshn(:,:) + zsshn_tlin(:,:)
1727
1728	zgcb_tlin(:,:) = gamma * zgcb_tlin(:,:)
1729	gcb(:,:) = gcb(:,:) + zgcb_tlin(:,:)
1730
1731	zgcx_tlin(:,:) = gamma * zgcx_tlin(:,:)
1732	gcx(:,:) = gcx(:,:) + zgcx_tlin(:,:)
1733	ENDIF
1734
1735	!--------------------------------------------------------------------
1736	! Compute the direct model F(X0,t=n) = Xn
1737	!--------------------------------------------------------------------
1738	IF ( tlm_bch /= 2 ) CALL dyn_spg_flt(nit000, indic)
1739	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
1740	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
1741	!--------------------------------------------------------------------
1742	! Compute the Tangent
1743	!--------------------------------------------------------------------
1744	IF ( tlm_bch == 2 ) THEN
1745	gcx_tl (:,:) = 0.0_wp
1746	gcxb_tl(:,:) = 0.0_wp
1747	gcb_tl (:,:) = 0.0_wp
1748	!--------------------------------------------------------------------
1749	! Initialize the tangent variables
1750	!--------------------------------------------------------------------
1751	CALL trj_rea( nit000-1, 1 )
1752	CALL trj_rea( nit000, 1 )
1753	ua_tl (:,:,:) = zua_tlin (:,:,:)
1754	va_tl (:,:,:) = zva_tlin (:,:,:)
1755	ub_tl (:,:,:) = zub_tlin (:,:,:)
1756	vb_tl (:,:,:) = zvb_tlin (:,:,:)
1757	wn_tl (:,:,:) = zwn_tlin (:,:,:)
1758	emp_tl (:,: ) = zemp_tlin (:,: )
1759	sshb_tl(:,: ) = zsshb_tlin(:,: )
1760	sshn_tl(:,: ) = zsshn_tlin(:,: )
1761	gcb_tl (:,:) = zgcb_tlin (:,:)
1762	gcx_tl (:,:) = zgcx_tlin (:,:)
1763
1764	CALL dyn_spg_flt_tan(nit000, indic)
1765	!--------------------------------------------------------------------
1766	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
1767	!--------------------------------------------------------------------
1768
1769	zsp2_1 = DOT_PRODUCT( ua_tl, ua_tl )
1770	zsp2_2 = DOT_PRODUCT( va_tl, va_tl )
1771	zsp2_3 = DOT_PRODUCT( sshb_tl, sshb_tl )
1772	zsp2_4 = DOT_PRODUCT( sshn_tl, sshn_tl )
1773	zsp_5 = DOT_PRODUCT( gcx_tl, gcx_tl )
1774	zsp_6 = DOT_PRODUCT( gcxb_tl, gcxb_tl )
1775
1776	zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp_5 + zsp_6
1777	!--------------------------------------------------------------------
1778	! Storing data
1779	!--------------------------------------------------------------------
1780	CALL trj_rd_spl(file_wop)
1781	zua_wop (:,:,:) = ua (:,:,:)
1782	zva_wop (:,:,:) = va (:,:,:)
1783	zsshb_wop(:,:) = sshb(:,:)
1784	zsshn_wop(:,:) = sshn(:,:)
1785	zgcx_wop (:,:) = gcx (:,:)
1786	CALL trj_rd_spl(file_xdx)
1787	zua_out (:,:,:) = ua (:,:,:)
1788	zva_out (:,:,:) = va (:,:,:)
1789	zsshb_out(:,:) = sshb(:,:)
1790	zsshn_out(:,:) = sshn(:,:)
1791	zgcx_out (:,:) = gcx (:,:)
1792	!--------------------------------------------------------------------
1793	! Compute the Linearization Error
1794	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
1795	! and
1796	! Compute the Linearization Error
1797	! En = Nn -TL(gamma.dX0, t0,tn)
1798	!--------------------------------------------------------------------
1799	! Warning: Here we re-use local variables z()_out and z()_wop
1800	ii=0
1801	DO jk = 1, jpk
1802	DO jj = 1, jpj
1803	DO ji = 1, jpi
1804	zua_out (ji,jj,jk) = zua_out (ji,jj,jk) - zua_wop (ji,jj,jk)
1805	zua_wop (ji,jj,jk) = zua_out (ji,jj,jk) - ua_tl (ji,jj,jk)
1806	IF ( ua_tl(ji,jj,jk) .NE. 0.0_wp ) &
1807	& zerrua(ji,jj,jk) = zua_out(ji,jj,jk)/ua_tl(ji,jj,jk)
1808	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1809	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1810	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1811	ii = ii+1
1812	iiposua(ii) = ji
1813	ijposua(ii) = jj
1814	ikposua(ii) = jk
1815	IF ( INT(umask(ji,jj,jk)) .NE. 0) THEN
1816	zscua (ii) = zua_wop(ji,jj,jk)
1817	zscerrua (ii) = ( zerrua(ji,jj,jk) - 1.0_wp ) /gamma
1818	ENDIF
1819	ENDIF
1820	END DO
1821	END DO
1822	END DO
1823	ii=0
1824	DO jk = 1, jpk
1825	DO jj = 1, jpj
1826	DO ji = 1, jpi
1827	zva_out (ji,jj,jk) = zva_out (ji,jj,jk) - zva_wop (ji,jj,jk)
1828	zva_wop (ji,jj,jk) = zva_out (ji,jj,jk) - va_tl (ji,jj,jk)
1829	IF ( va_tl(ji,jj,jk) .NE. 0.0_wp ) &
1830	& zerrva(ji,jj,jk) = zva_out(ji,jj,jk)/va_tl(ji,jj,jk)
1831	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1832	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1833	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1834	ii = ii+1
1835	iiposva(ii) = ji
1836	ijposva(ii) = jj
1837	ikposva(ii) = jk
1838	IF ( INT(vmask(ji,jj,jk)) .NE. 0) THEN
1839	zscva (ii) = zua_wop(ji,jj,jk)
1840	zscerrva (ii) = ( zerrva(ji,jj,jk) - 1.0_wp ) /gamma
1841	ENDIF
1842	ENDIF
1843	END DO
1844	END DO
1845	END DO
1846	ii=0
1847	DO jj = 1, jpj
1848	DO ji = 1, jpi
1849	zsshb_out (ji,jj) = zsshb_out (ji,jj) - zsshb_wop (ji,jj)
1850	zsshb_wop (ji,jj) = zsshb_out (ji,jj) - sshb_tl (ji,jj)
1851	IF ( sshb_tl(ji,jj) .NE. 0.0_wp ) &
1852	& zerrsshb(ji,jj) = zsshb_out(ji,jj)/sshb_tl(ji,jj)
1853	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1854	& (MOD(jj, jsamp) .EQ. 0) ) THEN
1855	ii = ii+1
1856	iipossshb(ii) = ji
1857	ijpossshb(ii) = jj
1858	IF ( INT(tmask(ji,jj,1)) .NE. 0) THEN
1859	zscsshb (ii) = zsshb_wop(ji,jj)
1860	zscerrsshb (ii) = ( zerrsshb(ji,jj) - 1.0_wp ) / gamma
1861	ENDIF
1862	ENDIF
1863	END DO
1864	END DO
1865	ii=0
1866	DO jj = 1, jpj
1867	DO ji = 1, jpi
1868	zsshn_out (ji,jj) = zsshn_out (ji,jj) - zsshn_wop (ji,jj)
1869	zsshn_wop (ji,jj) = zsshn_out (ji,jj) - sshn_tl (ji,jj)
1870	IF ( sshn_tl(ji,jj) .NE. 0.0_wp ) &
1871	& zerrsshn(ji,jj) = zsshn_out(ji,jj)/sshn_tl(ji,jj)
1872	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1873	& (MOD(jj, jsamp) .EQ. 0) ) THEN
1874	ii = ii+1
1875	iipossshn(ii) = ji
1876	ijpossshn(ii) = jj
1877	IF ( INT(tmask(ji,jj,1)) .NE. 0) THEN
1878	zscsshn (ii) = zsshn_wop(ji,jj)
1879	zscerrsshn (ii) = ( zerrsshn(ji,jj) - 1.0_wp ) /gamma
1880	ENDIF
1881	ENDIF
1882	END DO
1883	END DO
1884	ii=0
1885	DO jj = 1, jpj
1886	DO ji = 1, jpi
1887	zgcx_out (ji,jj) = zgcx_out (ji,jj) - zgcx_wop (ji,jj)
1888	zgcx_wop (ji,jj) = zgcx_out (ji,jj) - gcx_tl (ji,jj)
1889	IF ( gcx_tl(ji,jj) .NE. 0.0_wp ) zerrgcx(ji,jj) = zgcx_out(ji,jj)/gcx_tl(ji,jj)
1890	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1891	& (MOD(jj, jsamp) .EQ. 0) ) THEN
1892	ii = ii+1
1893	iiposgcx(ii) = ji
1894	ijposgcx(ii) = jj
1895	IF ( INT(tmask(ji,jj,1)) .NE. 0) THEN
1896	zscgcx (ii) = zgcx_wop(ji,jj)
1897	zscerrgcx (ii) = ( zerrgcx(ji,jj) - 1.0_wp ) / gamma
1898	ENDIF
1899	ENDIF
1900	END DO
1901	END DO
1902	zsp1_1 = DOT_PRODUCT( zua_out, zua_out )
1903	zsp1_2 = DOT_PRODUCT( zva_out, zva_out )
1904	zsp1_3 = DOT_PRODUCT( zsshb_out, zsshb_out )
1905	zsp1_4 = DOT_PRODUCT( zsshn_out, zsshn_out )
1906	zsp1_5 = DOT_PRODUCT( zgcx_out, zgcx_out )
1907	zsp1 = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4 + zsp1_5
1908	zsp3_1 = DOT_PRODUCT( zua_wop, zua_wop )
1909	zsp3_2 = DOT_PRODUCT( zva_wop, zva_wop )
1910	zsp3_3 = DOT_PRODUCT( zsshb_wop, zsshb_wop )
1911	zsp3_4 = DOT_PRODUCT( zsshn_wop, zsshn_wop )
1912	zsp3_5 = DOT_PRODUCT( zgcx_wop, zgcx_wop )
1913	zsp3 = zsp3_1 + zsp3_2 + zsp3_3 + zsp3_4 + zsp3_5
1914
1915	!--------------------------------------------------------------------
1916	! Print the linearization error En - norme 2
1917	!--------------------------------------------------------------------
1918	! 14 char:'12345678901234'
1919	cl_name = 'dynspg_tam:En '
1920	zzsp = SQRT(zsp3)
1921	zzsp_1 = SQRT(zsp3_1)
1922	zzsp_2 = SQRT(zsp3_2)
1923	zzsp_3 = SQRT(zsp3_3)
1924	zzsp_4 = SQRT(zsp3_4)
1925	zzsp_5 = SQRT(zsp3_5)
1926	zgsp5 = zzsp
1927	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1928	!--------------------------------------------------------------------
1929	! Compute TLM norm2
1930	!--------------------------------------------------------------------
1931	zzsp = SQRT(zsp2)
1932	zzsp_1 = SQRT(zsp2_1)
1933	zzsp_2 = SQRT(zsp2_2)
1934	zzsp_3 = SQRT(zsp2_3)
1935	zzsp_4 = SQRT(zsp2_4)
1936	zzsp_5 = SQRT(zsp2_5)
1937	zgsp4 = zzsp
1938	cl_name = 'dynspg_tam:Ln2'
1939	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1940	!--------------------------------------------------------------------
1941	! Print the linearization error Nn - norme 2
1942	!--------------------------------------------------------------------
1943	zzsp = SQRT(zsp1)
1944	zzsp_1 = SQRT(zsp1_1)
1945	zzsp_2 = SQRT(zsp1_2)
1946	zzsp_3 = SQRT(zsp1_3)
1947	zzsp_4 = SQRT(zsp1_4)
1948	zzsp_5 = SQRT(zsp1_5)
1949	cl_name = 'dynspg:Mhdx-Mx'
1950	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1951
1952
1953	zgsp3 = SQRT( zsp3/zsp2 )
1954	zgsp7 = zgsp3/gamma
1955	zgsp1 = zzsp
1956	zgsp2 = zgsp1 / zgsp4
1957	zgsp6 = (zgsp2 - 1.0_wp)/gamma
1958
1959	FMT = "(A8,2X,I4.4,2X,E6.1,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13)"
1960	WRITE(numtan,FMT) 'dspgflt ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7
1961	!--------------------------------------------------------------------
1962	! Unitary calculus
1963	!--------------------------------------------------------------------
1964	FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)"
1965	cl_name = 'dspgflt '
1966	IF(lwp) THEN
1967	DO ii=1, 100, 1
1968	IF ( zscua(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscua ', &
1969	& cur_loop, h_ratio, ii, iiposua(ii), ijposua(ii), zscua(ii)
1970	ENDDO
1971	DO ii=1, 100, 1
1972	IF ( zscva(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscva ', &
1973	& cur_loop, h_ratio, ii, iiposva(ii), ijposva(ii), zscva(ii)
1974	ENDDO
1975	DO ii=1, 100, 1
1976	IF ( zscsshb(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscsshb ', &
1977	& cur_loop, h_ratio, ii, iipossshb(ii), ijpossshb(ii), zscsshb(ii)
1978	ENDDO
1979	DO ii=1, 100, 1
1980	IF ( zscsshn(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscsshn ', &
1981	& cur_loop, h_ratio, ii, iipossshn(ii), ijpossshn(ii), zscsshn(ii)
1982	ENDDO
1983	DO ii=1, 100, 1
1984	IF ( zscerrua(ii) .NE. 0.0_wp ) WRITE(numsctlm,FMT) cl_name, 'zscerrua ', &
1985	& cur_loop, h_ratio, ii, iiposua(ii), ijposua(ii), zscerrua(ii)
1986	ENDDO
1987	DO ii=1, 100, 1
1988	IF ( zscerrva(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrva ', &
1989	& cur_loop, h_ratio, ii, iiposva(ii), ijposva(ii), zscerrva(ii)
1990	ENDDO
1991	DO ii=1, 100, 1
1992	IF ( zscerrsshb(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrsshb ', &
1993	& cur_loop, h_ratio, ii, iipossshb(ii), ijpossshb(ii), zscerrsshb(ii)
1994	ENDDO
1995	DO ii=1, 100, 1
1996	IF ( zscerrsshn(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrsshn ', &
1997	& cur_loop, h_ratio, ii, iipossshn(ii), ijpossshn(ii), zscerrsshn(ii)
1998	ENDDO
1999	! write separator
2000	WRITE(numtan_sc,"(A4)") '===='
2001	ENDIF
2002	ENDIF
2003
2004	DEALLOCATE( &
2005	& zemp_tlin, zwn_tlin, &
2006	& zsshb_tlin, zsshn_tlin, &
2007	& zua_tlin, zva_tlin, &
2008	& zua_out, zva_out, &
2009	& zua_wop, zva_wop, &
2010	& zsshb_out, zsshn_out, &
2011	& zsshb_wop, zsshn_wop, &
2012	& z3r, z2r &
2013	& )
2014	END SUBROUTINE dyn_spg_flt_tlm_tst
2015	#endif
2016	#endif
2017	#endif
2018	END MODULE dynspg_flt_tam

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

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

Context Navigation

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/DYN/dynspg_flt_tam.F90 @ 2587

Download in other formats: