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.

dynspg_flt_tam.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/DYN – NEMO

Context Navigation

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

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

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

Download in other formats: