Context Navigation

step_tam.F90 @ 2587

Last change on this file since 2587 was 2587, checked in by vidard, 13 years ago
refer to ticket #798
Property svn:executable set to ``*
File size: 81.2 KB

Line
1	MODULE step_tam
2	#ifdef key_tam
3	!!======================================================================
4	!! * MODULE step_tam *
5	!! Time-stepping : manager of the adjoint ocean time stepping
6	!! Tangent and Adjoint module
7	!!======================================================================
8	!! History of
9	!! the direct: ! 91-03 () Original code
10	!! ! 91-11 (G. Madec)
11	!! ! 92-06 (M. Imbard) add a first output record
12	!! ! 96-04 (G. Madec) introduction of dynspg
13	!! ! 96-04 (M.A. Foujols) introduction of passive tracer
14	!! 8.0 ! 97-06 (G. Madec) new architecture of call
15	!! 8.2 ! 97-06 (G. Madec, M. Imbard, G. Roullet) free surface
16	!! 8.2 ! 99-02 (G. Madec, N. Grima) hpg implicit
17	!! 8.2 ! 00-07 (J-M Molines, M. Imbard) Open Bondary Conditions
18	!! 9.0 ! 02-06 (G. Madec) free form, suppress macro-tasking
19	!! " " ! 04-08 (C. Talandier) New trends organization
20	!! " " ! 05-01 (C. Ethe) Add the KPP closure scheme
21	!! " " ! 05-11 (V. Garnier) Surface pressure gradient organization
22	!! " " ! 05-11 (G. Madec) Reorganisation of tra and dyn calls
23	!! " " ! 06-01 (L. Debreu, C. Mazauric) Agrif implementation
24	!! " " ! 06-07 (S. Masson) restart using iom
25	!! " " " ! 07-04 (K. Mogensen, A. Weaver, M. Martin) Assimilation interface
26	!! History of the TAM
27	!! ! 08-06 (A. Vidard and A. Weaver) Tangent and Adjoint version of 9.0
28	!! ! 08-11 (A. Vidard) Nemo v3 update
29	!! ! 06-09 (F. Vigilant) Modified to split NEMOVAR / NEMOTAM
30	!!----------------------------------------------------------------------
31	!! stp_tam : OPA system time-stepping (tangent linear)
32	!! stp_adj : OPA system time-stepping (adjoint)
33	!!----------------------------------------------------------------------
34
35	USE par_kind , ONLY: & ! Precision variables
36	& wp
37	USE par_oce , ONLY: & ! Ocean space and time domain variables
38	& jpi, &
39	& jpj, &
40	& jpk, &
41	& jpiglo
42	USE oce , ONLY: & ! ocean dynamics and tracers variables
43	& tb, sb, tn, sn, ta, &
44	& un, vn, sshn, sshb, &
45	& sa, ub, vb, rn2, &
46	& ln_dynhpg_imp
47	USE zdfkpp , ONLY: &
48	& lk_zdfkpp ! KPP vertical mixing flag
49	USE dynspg_oce , ONLY: &
50	& lk_dynspg_rl ! Rigid-lid flag
51	USE dom_oce , ONLY: & ! ocean space and time domain variables
52	& n_cla, e1u, e2u, &
53	& e1v, e2v, e1t, e2t, &
54	# if defined key_vvl
55	& e3t_1, &
56	# else
57	# if defined key_zco
58	& e3t_0, e3w_0, &
59	# else
60	& e3t, e3u, e3v, e3w, &
61	# endif
62	# endif
63	& nldi, nldj, nlei, &
64	& nlej, lk_vvl, ln_zps,&
65	& narea, mig, mjg, &
66	& umask, vmask, tmask
67	USE in_out_manager, ONLY: & ! I/O manager
68	& numout, lwp, &
69	& nit000, nitend, &
70	& nstop, ctl_stop
71	USE daymod , ONLY: &
72	& adatrj
73	USE trabbc , ONLY: &! bottom boundary condition
74	& lk_trabbc
75	USE traqsr , ONLY: &
76	& ln_traqsr
77	! solar radiation penetration flag
78	USE asminc
79	USE asmbkg
80	USE oce_tam , ONLY: & ! Tangent linear and adjoint variables
81	& oce_tam_init, &
82	& un_tl, vn_tl, tn_tl, &
83	& sn_tl, ua_tl, va_tl, &
84	& ta_tl, sa_tl, ub_tl, &
85	& vb_tl, tb_tl, sb_tl, &
86	& un_ad, vn_ad, tn_ad, &
87	& sn_ad, ua_ad, va_ad, &
88	& ta_ad, sa_ad, ub_ad, &
89	& vb_ad, tb_ad, sb_ad, &
90	& sshn_tl, sshn_ad, &
91	#if defined key_obc
92	& sshb_ad, &
93	#endif
94	& rn2_tl, rn2_ad, &
95	& rhd_tl, rhop_tl, &
96	& rhd_ad, rhop_ad, &
97	& gtu_tl, gsu_tl, &
98	& gru_tl, gtv_tl, &
99	& gsv_tl, grv_tl, &
100	& gtu_ad, gsu_ad, &
101	& gru_ad, gtv_ad, &
102	& gsv_ad, grv_ad, &
103	& ssha_tl, sshb_tl, &
104	& rotb_tl, hdivb_tl, &
105	& rotn_tl, hdivn_tl, &
106	& wn_tl
107	USE lbclnk_tam
108	USE daymod_tam ! calendar (adjoint of day routine)
109	USE sbc_oce_tam
110	USE sbcmod_tam
111	USE traqsr_tam ! solar radiation penetration (adjoint of tra_qsr routine)
112	USE trasbc_tam ! surface boundary condition (adjoint of tra_sbc routine)
113	USE trabbc_tam ! bottom boundary condition (adjoint of tra_bbc routine)
114	USE tradmp_tam ! internal damping (adjoint of tra_dmp routine)
115	USE traadv_tam ! advection scheme control (adjoint of tra_adv_ctl routine)
116	USE traldf_tam ! lateral mixing (adjoint of tra_ldf routine)
117	USE cla_tam ! cross land advection (adjoint of tra_cla routine)
118	! zdfkpp ! KPP non-local tracer fluxes (adjoint of tra_kpp routine)
119	USE trazdf_tam ! vertical mixing (adjoint of tra_zdf routine)
120	USE tranxt_tam ! time-stepping (adjoint of tra_nxt routine)
121	USE eosbn2_tam ! equation of state (adjoint of eos_bn2 routine)
122
123	USE dynadv_tam ! advection (adjoint of dyn_adv routine)
124	USE dynvor_tam ! vorticity term (adjoint of dyn_vor routine)
125	USE dynhpg_tam ! hydrostatic pressure grad. (adjoint of dyn_hpg routine)
126	USE dynldf_tam ! lateral momentum diffusion (adjoint of dyn_ldf routine)
127	USE dynzdf_tam ! vertical diffusion (adjoint of dyn_zdf routine)
128	USE dynspg_tam ! surface pressure gradient (adjoint of dyn_spg routine)
129	USE dynnxt_tam ! time-stepping (adjoint of dyn_nxt routine)
130
131	! USE bdy_par_tam
132	! USE bdydta_tam
133
134	USE divcur_tam ! hor. divergence and curl (adjoint of div & cur routines)
135	USE cla_div_tam ! cross land: hor. divergence (adjoint of div_cla routine)
136	USE wzvmod_tam ! vertical velocity (adjoint of wzv routine)
137
138	!! USE zdfkpp_tam ! KPP vertical mixing
139	USE zdf_oce, ONLY : lk_zdfcst, avt, avt0, avmu, avmv, avm0, ln_zdfevd ! KPP vertical mixing
140	USE zdfddm, ONLY : & ! double diffusion mixing (zdf_ddm routine)
141	& lk_zdfddm, &
142	& zdf_ddm
143	USE zdfevd, only:zdf_evd ! double diffusion mixing (zdf_ddm routine)
144	USE zdfbfr, only:zdf_bfr ! double diffusion mixing (zdf_ddm routine)
145	USE zdfmxl, only:zdf_mxl ! double diffusion mixing (zdf_ddm routine)
146	USE eosbn2, ONLY: bn2
147	USE zpshde_tam ! partial step: hor. derivative (adjoint of zps_hde routine)
148
149	!! USE diaobs_tam ! obs-minus-model (adjoint of assimilation) (adjoint of dia_obs routine)
150
151	USE trj_tam
152
153	USE stpctl_tam ! time stepping control (adjoint of stp_ctl routine)
154
155	USE gridrandom, ONLY: &
156	! Random Gaussian noise on grids
157	& grid_random, &
158	& grid_rd_sd
159	USE dotprodfld, ONLY: & ! Computes dot product for 3D and 2D fields
160	& dot_product
161	USE tstool_tam, ONLY: &
162	& prntst_adj, &
163	& prntst_tlm, &
164	& stdemp, &
165	& stdu,stdv, &
166	& stdt, &
167	& stds, stdssh, &
168	& stdr
169
170	USE paresp, ONLY: &
171	! Weights for an energy-type scalar product
172	& wesp_t, &
173	& wesp_s, &
174	& wesp_u, &
175	& wesp_ssh
176
177	USE istate_tam !: Initial state setting (istate_init routine)
178	USE sol_oce_tam
179	USE trc_oce_tam
180
181	USE step
182	USE iom ! for dumping an abort state
183	#if defined key_obc
184	USE obc_par, ONLY : lk_obc
185	USE obcdta ! open boundary condition data (obc_dta routine)
186	#endif
187
188	#if defined key_agrif
189	#error 'agrif not yet implemented in nemotam'
190	#endif
191
192	IMPLICIT NONE
193	PRIVATE
194
195	PUBLIC stp_tan, &
196	& stp_adj, & ! called by simvar.F90
197	#if defined key_tst_tlm
198	& stp_tlm_tst, &
199	#endif
200	& stp_adj_tst
201
202	!! * Substitutions
203	# include "domzgr_substitute.h90"
204	# include "zdfddm_substitute.h90"
205
206	CONTAINS
207	#if defined key_agrif
208	#error 'agrif not yet implemented in nemotam'
209	! SUBROUTINE stp_tan( )
210	#else
211	SUBROUTINE stp_tan( kstp )
212	#endif
213	!!----------------------------------------------------------------------
214	!! * ROUTINE stp_tan *
215	!!
216	!! ** Purpose of the direct routine:
217	!! - Time stepping of OPA (momentum and active tracer eqs.)
218	!! - Time stepping of LIM (dynamic and thermodynamic eqs.)
219	!! - Tme stepping of TRC (passive tracer eqs.)
220	!!
221	!! ** Method of the direct routine:
222	!! -1- Update forcings and data
223	!! -2- Update ocean physics
224	!! -3- Compute the t and s trends
225	!! -4- Update t and s
226	!! -5- Compute the momentum trends
227	!! -6- Update the horizontal velocity
228	!! -7- Compute the diagnostics variables (rd,N2, div,cur,w)
229	!! -8- Outputs and diagnostics
230	!!----------------------------------------------------------------------
231	!! * Arguments
232	#if defined key_agrif
233	INTEGER :: kstp ! ocean time-step index
234	#else
235	INTEGER, INTENT( in ) :: kstp ! ocean time-step index
236	#endif
237
238	!! * local declarations
239	INTEGER :: indic ! error indicator if < 0
240	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zta_tmp, zsa_tmp
241	!! ---------------------------------------------------------------------
242
243	#if defined key_agrif
244	kstp = nit000 + Agrif_Nb_Step()
245	! IF ( Agrif_Root() .and. lwp) Write(,) '---'
246	! IF (lwp) Write(,) 'Grid Number',Agrif_Fixed(),' time step ',kstp
247	#endif
248	indic = 1 ! reset to no error condition
249
250	CALL day_tam( kstp, 0 ) ! Calendar
251
252	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
253	! Update data, open boundaries, surface boundary condition (including sea-ice)
254	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
255
256	! update 3D temperature data ... not needed in tangent
257
258	! update 3D salinity data ... not needed in tangent (to be investigated, see sbc_ssr)
259
260	CALL sbc_tan ( kstp ) ! Sea boundary condition (including sea-ice)
261
262	! IF( lk_obc ) CALL obc_dta( kstp ) ! ??? ! update dynamic and tracer data at open boundaries
263	! update dynamic and tracer data at open boundaries ... not needed for the time being, to investigate whenever we do obc.
264
265	! compute phase velocities at open boundaries ... not needed for the time being, to investigate whenever we do obc.
266	! IF( lk_bdy ) CALL bdy_dta_tan( kstp ) ! update dynamic and tracer data at unstructured open boundary
267
268	! Output the initial state and forcings ... not needed in tangent
269
270	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
271	! Ocean physics update
272	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
273	!-----------------------------------------------------------------------
274	! VERTICAL PHYSICS
275	!-----------------------------------------------------------------------
276	! N.B. ua, va, ta, sa arrays are used as workspace in this section
277	!-----------------------------------------------------------------------
278
279	CALL bn2_tan( tb, sb, tb_tl, sb_tl, rn2_tl ) ! before Brunt-Vaisala frequency
280
281	! ! Vertical eddy viscosity and diffusivity coefficients
282	! Richardson number dependent Kz ... not available
283	! TKE closure scheme for Kz ... not available
284	! KPP closure scheme for Kz ... not available
285
286	! Constant Kz (reset avt, avm[uv] to the background value)... not available
287
288	! lk_zdfcst_tan: Constant Kz read from the reference trajectory
289
290
291	! ! increase diffusivity at rivers mouths... not needed in tangent
292
293	! enhanced vertical eddy diffusivity ... not needed in tangent with lk_zdfcst_tan
294
295	! double diffusive mixing ... not needed in tangent with lk_zdfcst_tan
296
297	! bottom friction... not needed in tangent with lk_zdfcst_tan
298
299	! mixed layer depth... not needed in tangent with lk_zdfcst_tan
300
301
302	!-----------------------------------------------------------------------
303	! LATERAL PHYSICS
304	!-----------------------------------------------------------------------
305	! N.B. ua, va, ta, sa arrays are used as workspace in this section
306	!-----------------------------------------------------------------------
307
308	! before slope of the lateral mixing... not needed in tangent with lk_zdfcst_tan
309
310	#if defined key_traldf_c2d
311	! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_tan
312	#endif
313
314
315	#if defined key_top
316	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
317	! Passive Tracer Model
318	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
319	! N.B. ua, va, ta, sa arrays are used as workspace in this section
320	!-----------------------------------------------------------------------
321
322	! time-stepping... not needed in tangent for the time being
323
324	#endif
325	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
326	! Active tracers
327	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
328	! N.B. ua, va arrays are used as workspace in this section
329	!-----------------------------------------------------------------------
330
331	ta_tl(:,:,:) = 0.e0 ! set tracer trends to zero
332	sa_tl(:,:,:) = 0.e0
333
334	! apply tracer assimilation increment ... not needed in tangent
335	! apply bias directly to tn, sn, tb, sb ... not needed in tangent
336	CALL tra_sbc_tan( kstp ) ! surface boundary condition
337
338	IF( ln_traqsr ) CALL tra_qsr_tan( kstp ) ! penetrative solar radiation qsr
339
340	IF( lk_trabbc ) CALL tra_bbc_tan( kstp ) ! bottom heat flux
341
342	! diffusive bottom boundary layer scheme ... currently not available
343	! advective (and/or diffusive) bottom boundary layer scheme ... currently not available
344
345	IF( lk_tradmp ) CALL tra_dmp_tan( kstp ) ! internal damping trends
346
347	CALL tra_adv_tan( kstp ) ! horizontal & vertical advection
348
349	IF( n_cla == 1 ) CALL tra_cla_tan( kstp ) ! Cross Land Advection (Update Hor. advection)
350
351	! IF( lk_zdfkpp ) CALL tra_kpp_tan( kstp ) ! KPP non-local tracer fluxes
352
353	CALL tra_ldf_tan( kstp ) ! lateral mixing
354
355	#if defined key_agrif
356	! tracers sponge ... not available
357	#endif
358	CALL tra_zdf_tan( kstp ) ! vertical mixing
359
360	! update the new (t,s) fields by non
361	! penetrative convective adjustment ... not available
362
363	CALL tra_nxt_tan( kstp ) ! tracer fields at next time step
364
365	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg
366	CALL eos_tan( ta, sa, ta_tl, sa_tl, rhd_tl, rhop_tl ) ! Time-filtered in situ density used in dynhpg module
367	IF( ln_zps ) CALL zps_hde_tan( kstp, ta, sa, ta_tl, sa_tl, rhd_tl,& ! Partial steps: time filtered hor. gradient
368	& gtu_tl, gsu_tl, gru_tl, & ! of t, s, rd at the bottom ocean level
369	& gtv_tl, gsv_tl, grv_tl )
370	ELSE ! centered hpg (default case)
371	CALL eos_tan( tb, sb, tb_tl, sb_tl, rhd_tl, rhop_tl ) ! now (swap=before) in situ density for dynhpg module
372	IF( ln_zps ) CALL zps_hde_tan( kstp, tb, sb, tb_tl, sb_tl, rhd_tl,& ! Partial steps: now horizontal gradient
373	& gtu_tl, gsu_tl, gru_tl, & ! of t, s, rd at the bottom ocean level
374	& gtv_tl, gsv_tl, grv_tl )
375	ENDIF
376
377	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
378	! Dynamics
379	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
380	! N.B. ta, sa arrays are used as workspace in this section
381	!-----------------------------------------------------------------------
382
383	ua_tl(:,:,:) = 0.e0 ! set dynamics trends to zero
384	va_tl(:,:,:) = 0.e0
385
386	! apply dynamics assimilation increment ... not needed in tangent
387
388	CALL dyn_adv_tan( kstp ) ! advection (vector or flux form)
389	CALL dyn_vor_tan( kstp ) ! vorticity term including Coriolis
390	CALL dyn_ldf_tan( kstp ) ! lateral mixing
391	#if defined key_agrif
392	! momemtum sponge ... not available
393	#endif
394
395	CALL dyn_hpg_tan( kstp ) ! horizontal gradient of Hydrostatic pressure
396	CALL dyn_zdf_tan( kstp ) ! vertical diffusion
397
398	IF( lk_dynspg_rl ) THEN
399	! surface pressure gradient at open boundaries ... not available
400	ENDIF
401	indic=0
402	!i
403	CALL dyn_spg_tan( kstp, indic ) ! surface pressure gradient
404	CALL dyn_nxt_tan( kstp ) ! lateral velocity at next time step
405
406	! vertical mesh at next time step ... not available
407
408	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
409	! Computation of diagnostic variables
410	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
411	! N.B. ua, va, ta, sa arrays are used as workspace in this section
412	!-----------------------------------------------------------------------
413
414	CALL div_cur_tan( kstp ) ! Horizontal divergence & Relative vorticity
415
416	IF( n_cla == 1 ) CALL div_cla_tan( kstp ) ! Cross Land Advection (Update Hor. divergence)
417
418	CALL wzv_tan( kstp ) ! Vertical velocity
419
420	CALL trj_rea( kstp, 1) ! ... Read basic state trajectory at end of current step
421
422	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
423	! Control, and restarts
424	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
425	! N.B. ua, va, ta, sa arrays are used as workspace in this section
426	!-----------------------------------------------------------------------
427	! ! Time loop: control and print
428	CALL stp_ctl_tan( kstp, indic, 0 )
429	IF( indic < 0 ) CALL ctl_stop( 'step_tan: indic < 0' )
430
431	! close input ocean restart file ... not needed in tangent
432	! write output ocean restart file... not needed in tangent
433	! write open boundary restart file... not needed in tangent
434
435	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
436	! diagnostics and outputs
437	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
438	! N.B. ua, va, ta, sa arrays are used as workspace in this section
439	!-----------------------------------------------------------------------
440
441	IF ( nstop == 0 ) THEN ! Diagnostics
442	! drifting Floats... not needed in tangent
443	! trends: dynamics ... not needed in tangent
444	! trends: active tracers... not needed in tangent
445	! trends: Mixed-layer ... not needed in tangent
446	! trends: vorticity budget... not needed in tangent
447	! Surface pressure diagnostics... not needed in tangent
448	! Thermocline depth (20 degres isotherm depth)... not needed in tangent
449	! basin averaged diagnostics... not needed in tangent
450	! dynamical heigh diagnostics... not needed in tangent
451	! Fresh water budget diagnostics... not needed in tangent
452	!!! IF( lk_diaobs ) CALL dia_obs_tan( kstp ) ! obs-minus-model (assimilation) diagnostics NOTE: may be better off outside this module
453	! Poleward TRansports diagnostics... not needed in tangent
454
455
456	! ! Outputs
457	! ocean model: outputs... not needed in tangent
458
459	ENDIF
460
461	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
462	! Coupled mode
463	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
464
465	! coupled mode : field exchanges ... not available for the next 20 years
466	!
467	END SUBROUTINE stp_tan
468
469	#if defined key_agrif
470	#error 'agrif not yet implemented in nemotam'
471	! SUBROUTINE stp_adj( )
472	#else
473	SUBROUTINE stp_adj( kstp )
474	#endif
475	!!----------------------------------------------------------------------
476	!! * ROUTINE stp_adj *
477	!!
478	!! ** Purpose of the direct routine:
479	!! - Time stepping of OPA (momentum and active tracer eqs.)
480	!! - Time stepping of LIM (dynamic and thermodynamic eqs.)
481	!! - Tme stepping of TRC (passive tracer eqs.)
482	!!
483	!! ** Method of the direct routine:
484	!! -1- Update forcings and data
485	!! -2- Update ocean physics
486	!! -3- Compute the t and s trends
487	!! -4- Update t and s
488	!! -5- Compute the momentum trends
489	!! -6- Update the horizontal velocity
490	!! -7- Compute the diagnostics variables (rd,N2, div,cur,w)
491	!! -8- Outputs and diagnostics
492	!!----------------------------------------------------------------------
493	!! * Arguments
494	#if defined key_agrif
495	INTEGER :: kstp ! ocean time-step index
496	#else
497	INTEGER, INTENT( in ) :: kstp ! ocean time-step index
498	#endif
499
500	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zta_tmp, zsa_tmp
501
502	!! * local declarations
503	INTEGER :: indic ! error indicator if < 0
504	!! ---------------------------------------------------------------------
505
506	#if defined key_agrif
507	kstp = nit000 + Agrif_Nb_Step()
508	! IF ( Agrif_Root() .and. lwp) Write(,) '---'
509	! IF (lwp) Write(,) 'Grid Number',Agrif_Fixed(),' time step ',kstp
510	#endif
511	indic = 1 ! reset to no error condition
512
513	CALL day_tam( kstp, 0 ) ! Calendar
514
515	! ... Read basic state trajectory at end of previous step
516	CALL trj_rea( ( kstp - 1 ), -1 )
517
518	zta_tmp (:,:,:) = ta (:,:,:)
519	zsa_tmp (:,:,:) = sa (:,:,:)
520	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
521	! Computation of diagnostic variables
522	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
523	! N.B. ua, va, ta, sa arrays are used as workspace in this section
524	!-----------------------------------------------------------------------
525
526	! ocean surface freezing temperature ... not needed in adjoint
527
528	CALL wzv_adj( kstp ) ! Vertical velocity
529	IF( n_cla == 1 ) CALL div_cla_adj( kstp ) ! Cross Land Advection (Update Hor. divergence)
530
531
532	CALL div_cur_adj( kstp ) ! Horizontal divergence & Relative vorticity
533	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
534	! Dynamics
535	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
536	! N.B. ta, sa arrays are used as workspace in this section
537	!-----------------------------------------------------------------------
538	! apply dynamics assimilation increment ... not needed in adjoint
539	indic=0
540	!i bug lbc sur emp
541
542	! vertical mesh at next time step ... not available
543	CALL dyn_nxt_adj( kstp ) ! lateral velocity at next time step
544	CALL dyn_spg_adj( kstp, indic ) ! surface pressure gradient
545	!i
546	IF( lk_dynspg_rl ) THEN
547	! surface pressure gradient at open boundaries ... not available
548	ENDIF
549	CALL dyn_zdf_adj( kstp ) ! vertical diffusion
550	CALL dyn_hpg_adj( kstp ) ! horizontal gradient of Hydrostatic pressure
551	#if defined key_agrif
552	! momemtum sponge ... not available
553	#endif
554	CALL dyn_ldf_adj( kstp ) ! lateral mixing
555	CALL dyn_vor_adj( kstp ) ! vorticity term including Coriolis
556	CALL dyn_adv_adj( kstp ) ! advection (vector or flux form)
557
558	ta (:,:,:) = zta_tmp (:,:,:)
559	sa (:,:,:) = zsa_tmp (:,:,:)
560
561	ua_ad(:,:,:) = 0.e0 ! set tracer trends to zero
562	va_ad(:,:,:) = 0.e0
563
564	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg
565	IF( ln_zps ) CALL zps_hde_adj( kstp, ta, sa, ta_ad, sa_ad, rhd_ad,& ! Partial steps: time filtered hor. gradient
566	& gtu_ad, gsu_ad, gru_ad, & ! of t, s, rd at the bottom ocean level
567	& gtv_ad, gsv_ad, grv_ad )
568	CALL eos_adj( ta, sa, ta_ad, sa_ad, rhd_ad, rhop_ad ) ! Time-filtered in situ density used in dynhpg module
569	ELSE ! centered hpg (default case)
570	IF( ln_zps ) CALL zps_hde_adj( kstp, tb, sb, &
571	& tb_ad, sb_ad,rhd_ad,& ! Partial steps: now horizontal gradient
572	& gtu_ad, gsu_ad, gru_ad, & ! of t, s, rd at the bottom ocean level
573	& gtv_ad, gsv_ad, grv_ad )
574	CALL eos_adj( tb, sb, tb_ad, sb_ad, rhd_ad, rhop_ad ) ! now (swap=before) in situ density for dynhpg module
575	ENDIF
576
577	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
578	! Active tracers
579	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
580	! N.B. ua, va arrays are used as workspace in this section
581	!-----------------------------------------------------------------------
582
583	CALL tra_nxt_adj( kstp ) ! tracer fields at next time step
584	! update the new (t,s) fields by non
585	! penetrative convective adjustment ... not available
586	CALL tra_zdf_adj( kstp ) ! vertical mixing
587
588	#if defined key_agrif
589	! tracers sponge ... not available
590	#endif
591
592	CALL tra_ldf_adj( kstp ) ! lateral mixing
593
594	! IF( lk_zdfkpp ) CALL tra_kpp_adj( kstp ) ! KPP non-local tracer fluxes
595
596	IF( n_cla == 1 ) CALL tra_cla_adj( kstp ) ! Cross Land Advection (Update Hor. advection)
597
598	CALL tra_adv_adj( kstp ) ! horizontal & vertical advection
599
600	IF( lk_tradmp ) CALL tra_dmp_adj( kstp ) ! internal damping trends
601	! advective (and/or diffusive) bottom boundary layer scheme ... currently not available
602	! diffusive bottom boundary layer scheme ... currently not available
603
604	IF( lk_trabbc ) CALL tra_bbc_adj( kstp ) ! bottom heat flux
605
606	IF( ln_traqsr ) CALL tra_qsr_adj( kstp ) ! penetrative solar radiation qsr
607
608	CALL tra_sbc_adj( kstp ) ! surface boundary condition
609
610	ta_ad(:,:,:) = 0.e0 ! set tracer trends to zero
611	sa_ad(:,:,:) = 0.e0
612	! apply tracer assimilation increment ... not needed in adjoint
613	#if defined key_top
614	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
615	! Passive Tracer Model
616	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
617	! N.B. ua, va, ta, sa arrays are used as workspace in this section
618	!-----------------------------------------------------------------------
619
620	! time-stepping... not needed in adjoint for the time being
621
622	#endif
623	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
624	! Ocean physics update
625	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
626	!-----------------------------------------------------------------------
627	! LATERAL PHYSICS
628	!-----------------------------------------------------------------------
629	! N.B. ua, va, ta, sa arrays are used as workspace in this section
630	!-----------------------------------------------------------------------
631	#if defined key_traldf_c2d
632	! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_adj
633	#endif
634	! before slope of the lateral mixing... not needed in adjoint with lk_zdfcst_adj
635	!-----------------------------------------------------------------------
636	! VERTICAL PHYSICS
637	!-----------------------------------------------------------------------
638	! N.B. ua, va, ta, sa arrays are used as workspace in this section
639	!-----------------------------------------------------------------------
640	! mixed layer depth... not needed in adjoint with lk_zdfcst_adj
641	! bottom friction... not needed in adjoint with lk_zdfcst_adj
642	! double diffusive mixing ... not needed in adjoint with lk_zdfcst_adj
643	! enhanced vertical eddy diffusivity ... not needed in adjoint with lk_zdfcst_adj
644	! ! increase diffusivity at rivers mouths... not needed in tangent
645	! lk_zdfcst_adj: Constant Kz read from the reference trajectory
646	! Constant Kz (reset avt, avm[uv] to the background value)... not available
647	! KPP closure scheme for Kz ... not available
648	! TKE closure scheme for Kz ... not available
649	! Richardson number dependent Kz ... not available
650	! ! Vertical eddy viscosity and diffusivity coefficients
651	CALL bn2_adj( tb, sb, tb_ad, sb_ad, rn2_ad ) ! before Brunt-Vaisala frequency
652
653	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
654	! Update data, open boundaries and Forcings
655	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
656	! Output the initial state and forcings ... not needed in adjoint
657	! IF( lk_bdy ) CALL bdy_dta_adj( kstp ) ! update dynamic and tracer data at unstructured open boundary
658	! compute phase velocities at open boundaries ... not needed for the time being, to investigate whenever we do obc.
659	! update dynamic and tracer data at open boundaries ... not needed for the time being, to investigate whenever we do obc.
660
661	CALL sbc_adj ( kstp ) ! Sea boundary condition (including sea-ice)
662	! update 3D salinity data ... not needed in tangent
663	! update 3D temperature data ... not needed in adjoint
664
665	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
666	! Control, and restarts
667	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
668	! N.B. ua, va, ta, sa arrays are used as workspace in this section
669	!-----------------------------------------------------------------------
670	! ! Time loop: control and print
671	CALL stp_ctl_adj( kstp, indic, 0 )
672	IF( indic < 0 ) CALL ctl_stop( 'step_adj: indic < 0' )
673
674	! close input ocean restart file ... not needed in adjoint
675	! write output ocean restart file... not needed in adjoint
676	! write open boundary restart file... not needed in adjoint
677
678	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
679	! diagnostics and outputs
680	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
681	! N.B. ua, va, ta, sa arrays are used as workspace in this section
682	!-----------------------------------------------------------------------
683
684	IF ( nstop == 0 ) THEN ! Diagnostics
685	! drifting Floats... not needed in adjoint
686	! trends: dynamics ... not needed in adjoint
687	! trends: active tracers... not needed in adjoint
688	! trends: Mixed-layer ... not needed in adjoint
689	! trends: vorticity budget... not needed in adjoint
690	! Surface pressure diagnostics... not needed in adjoint
691	! Thermocline depth (20 degres isotherm depth)... not needed in adjoint
692	! basin averaged diagnostics... not needed in adjoint
693	! dynamical heigh diagnostics... not needed in adjoint
694	! Fresh water budget diagnostics... not needed in adjoint
695	!!! IF( lk_diaobs ) CALL dia_obs_adj( kstp ) ! obs-minus-model (assimilation) diagnostics NOTE: may be better off outside this module
696	! Poleward TRansports diagnostics... not needed in adjoint
697
698
699	! ! Outputs
700	! ocean model: outputs... not needed in adjoint
701
702	ENDIF
703
704	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
705	! Coupled mode
706	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
707
708	! coupled mode : field exchanges ... not available for the next 20 years
709	!
710	END SUBROUTINE stp_adj
711
712	SUBROUTINE stp_adj_tst( kumadt )
713	!!-----------------------------------------------------------------------
714	!!
715	!! * ROUTINE example_adj_tst *
716	!!
717	!! ** Purpose : Test the adjoint routine.
718	!!
719	!! ** Method : Verify the scalar product
720	!!
721	!! ( L dx )^T W dy = dx^T L^T W dy
722	!!
723	!! where L = tangent routine
724	!! L^T = adjoint routine
725	!! W = diagonal matrix of scale factors
726	!! dx = input perturbation (random field)
727	!! dy = L dx
728	!!
729	!!
730	!! History :
731	!! ! 09-03 (F. Vigilant)
732	!!-----------------------------------------------------------------------
733	!! * Modules used
734	USE sbcssr_tam, ONLY: &
735	& qrp_tl, qrp_ad, &
736	& erp_tl, erp_ad
737	USE sbcfwb_tam, ONLY: &
738	& a_fwb_tl, & ! for before year.
739	& a_fwb_ad ! for before year.
740
741	!! * Arguments
742	INTEGER, INTENT(IN) :: &
743	& kumadt ! Output unit
744
745	!! * Local declarations
746	INTEGER :: &
747	& ji, & ! dummy loop indices
748	& jj, &
749	& jk, &
750	& istp
751	INTEGER, DIMENSION(jpi,jpj) :: &
752	& iseed_2d ! 2D seed for the random number generator
753	REAL(KIND=wp) :: &
754	& zsp1 , & ! scalar product involving the tangent routine
755	& zsp1_U , &
756	& zsp1_V , &
757	& zsp1_T , &
758	& zsp1_S , &
759	& zsp1_SSH, &
760	& zsp2 , & ! scalar product involving the adjoint routine
761	& zsp2_U , &
762	& zsp2_V , &
763	& zsp2_T , &
764	& zsp2_S , &
765	& zsp2_SSH
766	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
767	& zun_tlin , & ! Tangent input
768	& zvn_tlin , & ! Tangent input
769	& ztn_tlin , & ! Tangent input
770	& zsn_tlin , & ! Tangent input
771	& zun_tlout, & ! Tangent output
772	& zvn_tlout, & ! Tangent output
773	& ztn_tlout, & ! Tangent output
774	& zsn_tlout, & ! Tangent outpu
775	& zun_adin , & ! Adjoint input
776	& zvn_adin , & ! Adjoint input
777	& ztn_adin , & ! Adjoint input
778	& zsn_adin , & ! Adjoint input
779	#if defined key_obc
780	& zub_adin , & ! Adjoint input
781	& zvb_adin , & ! Adjoint input
782	& ztb_adin , & ! Adjoint input
783	& zsb_adin , & ! Adjoint input
784	& zub_tlout , & ! Adjoint input
785	& zvb_tlout , & ! Adjoint input
786	& ztb_tlout , & ! Adjoint input
787	& zsb_tlout , & ! Adjoint input
788	#endif
789	& zun_adout, & ! Adjoint output
790	& zvn_adout, & ! Adjoint output
791	& ztn_adout, & ! Adjoint output
792	& zsn_adout, & ! Adjoint output
793	& z3r ! 3D random field
794	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
795	& zsshn_tlin , & ! Tangent input
796	& zsshn_tlout , & ! Tangent input
797	& zsshn_adin , & ! Adjoint output
798	& zsshn_adout , & ! Adjoint output
799	#if defined key_obc
800	& zsshb_tlout , & ! Tangent input
801	& zsshb_adin , & ! Tangent input
802	#endif
803	& z2r ! 2D random field
804
805	CHARACTER(LEN=14) :: &
806	& cl_name
807
808	INTEGER :: &
809	& jpert
810	INTEGER, PARAMETER :: jpertmax = 6
811
812	! Allocate memory
813
814	ALLOCATE( &
815	& zun_tlin( jpi,jpj,jpk), zvn_tlin( jpi,jpj,jpk), ztn_tlin( jpi,jpj,jpk), &
816	& zsn_tlin( jpi,jpj,jpk), zsshn_tlin( jpi,jpj ), zun_tlout( jpi,jpj,jpk), &
817	& zvn_tlout( jpi,jpj,jpk), ztn_tlout( jpi,jpj,jpk), zsn_tlout( jpi,jpj,jpk), &
818	& zsshn_tlout(jpi,jpj ), zun_adin( jpi,jpj,jpk), zvn_adin( jpi,jpj,jpk), &
819	& ztn_adin( jpi,jpj,jpk), zsn_adin( jpi,jpj,jpk), zsshn_adin(jpi,jpj ), &
820	& zun_adout( jpi,jpj,jpk), zvn_adout( jpi,jpj,jpk), ztn_adout( jpi,jpj,jpk), &
821	& zsn_adout( jpi,jpj,jpk), zsshn_adout(jpi,jpj ), z2r( jpi,jpj ), &
822	& z3r( jpi,jpj,jpk) &
823	& )
824
825	#if defined key_obc
826	ALLOCATE( zub_adin (jpi,jpj,jpk), zvb_adin (jpi,jpj,jpk) , &
827	& ztb_adin (jpi,jpj,jpk), zsb_adin (jpi,jpj,jpk) , &
828	& zub_tlout (jpi,jpj,jpk), zvb_tlout (jpi,jpj,jpk) , &
829	& ztb_tlout (jpi,jpj,jpk), zsb_tlout (jpi,jpj,jpk) , &
830	& zsshb_tlout(jpi,jpj) , zsshb_adin(jpi,jpj) )
831	#endif
832	!==================================================================
833	! 1) dx = ( un_tl, vn_tl, tn_tl, sn_tl, sshn_tl ) and
834	! dy = ( ..... )
835	!==================================================================
836
837
838	DO jpert = 1, jpertmax
839
840	!--------------------------------------------------------------------
841	! Reset the tangent and adjoint variables
842	!--------------------------------------------------------------------
843	zun_tlin (:,:,:) = 0.0_wp
844	zvn_tlin (:,:,:) = 0.0_wp
845	ztn_tlin (:,:,:) = 0.0_wp
846	zsn_tlin (:,:,:) = 0.0_wp
847	zsshn_tlin ( :,:) = 0.0_wp
848	zun_tlout (:,:,:) = 0.0_wp
849	zvn_tlout (:,:,:) = 0.0_wp
850	ztn_tlout (:,:,:) = 0.0_wp
851	zsn_tlout (:,:,:) = 0.0_wp
852	zsshn_tlout( :,:) = 0.0_wp
853	zun_adin (:,:,:) = 0.0_wp
854	zvn_adin (:,:,:) = 0.0_wp
855	ztn_adin (:,:,:) = 0.0_wp
856	zsn_adin (:,:,:) = 0.0_wp
857	zsshn_adin ( :,:) = 0.0_wp
858	zun_adout (:,:,:) = 0.0_wp
859	zvn_adout (:,:,:) = 0.0_wp
860	ztn_adout (:,:,:) = 0.0_wp
861	zsn_adout (:,:,:) = 0.0_wp
862	zsshn_adout( :,:) = 0.0_wp
863	z3r (:,:,:) = 0.0_wp
864	z2r ( :,:) = 0.0_wp
865
866	!--------------------------------------------------------------------
867	! Initialize the tangent input with random noise: dx
868	!--------------------------------------------------------------------
869
870	IF ( (jpert == 1) .OR. (jpert == jpertmax) ) THEN
871	DO jj = 1, jpj
872	DO ji = 1, jpi
873	iseed_2d(ji,jj) = - ( 596035 + &
874	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
875	END DO
876	END DO
877	CALL grid_random( iseed_2d, z3r, 'U', 0.0_wp, stdu )
878	DO jk = 1, jpk
879	DO jj = nldj, nlej
880	DO ji = nldi, nlei
881	zun_tlin(ji,jj,jk) = z3r(ji,jj,jk)
882	END DO
883	END DO
884	END DO
885	ENDIF
886	IF ( (jpert == 2) .OR. (jpert == jpertmax) ) THEN
887	DO jj = 1, jpj
888	DO ji = 1, jpi
889	iseed_2d(ji,jj) = - ( 371836 + &
890	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
891	END DO
892	END DO
893	CALL grid_random( iseed_2d, z3r, 'V', 0.0_wp, stdv )
894	DO jk = 1, jpk
895	DO jj = nldj, nlej
896	DO ji = nldi, nlei
897	zvn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
898	END DO
899	END DO
900	END DO
901	ENDIF
902	IF ( (jpert == 3) .OR. (jpert == jpertmax) ) THEN
903	DO jj = 1, jpj
904	DO ji = 1, jpi
905	iseed_2d(ji,jj) = - ( 284035 + &
906	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
907	END DO
908	END DO
909	CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stdt )
910	DO jk = 1, jpk
911	DO jj = nldj, nlej
912	DO ji = nldi, nlei
913	ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
914	END DO
915	END DO
916	END DO
917	ENDIF
918	IF ( (jpert == 4) .OR. (jpert == jpertmax) ) THEN
919	DO jj = 1, jpj
920	DO ji = 1, jpi
921	iseed_2d(ji,jj) = - ( 471426 + &
922	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
923	END DO
924	END DO
925	CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stds )
926	DO jk = 1, jpk
927	DO jj = nldj, nlej
928	DO ji = nldi, nlei
929	zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
930	END DO
931	END DO
932	END DO
933	ENDIF
934	IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN
935	DO jj = 1, jpj
936	DO ji = 1, jpi
937	iseed_2d(ji,jj) = - ( 712035 + &
938	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
939	END DO
940	END DO
941	CALL grid_random( iseed_2d, z2r, 'T', 0.0_wp, stdssh )
942	DO jj = nldj, nlej
943	DO ji = nldi, nlei
944	zsshn_tlin(ji,jj) = z2r(ji,jj)
945	END DO
946	END DO
947	ENDIF
948	IF (jpert>0) THEN
949	CALL oce_tam_init( 0 ) ! kindic = 0 allocate/initialize tl and adj variables
950	CALL sbc_oce_tam_init( 0 ) ! kindic = 0 allocate/initialize tl and adj variables
951	CALL sol_oce_tam_init( 1 )
952	CALL sol_oce_tam_init( 2 )
953	CALL trc_oce_tam_init( 0 )
954
955	qrp_tl = 0.0_wp
956	qrp_ad = 0.0_wp
957	erp_tl = 0.0_wp
958	erp_ad = 0.0_wp
959
960	emp_tl = 0.e0
961	emp_ad = 0.e0
962
963	a_fwb_tl = 0.0_wp
964	a_fwb_ad = 0.0_wp
965
966	istp = nit000 - 1
967	CALL trj_rea( istp, 1 )
968	!--------------------------------------------------------------------
969	! Initialize the tangent variables: dy^* = W dy
970	!--------------------------------------------------------------------
971
972	un_tl (:,:,:) = zun_tlin (:,:,:)
973	vn_tl (:,:,:) = zvn_tlin (:,:,:)
974	tn_tl (:,:,:) = ztn_tlin (:,:,:)
975	sn_tl (:,:,:) = zsn_tlin (:,:,:)
976	sshn_tl( :,:) = zsshn_tlin ( :,:)
977
978	#if defined key_pomme_r025
979	IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN
980	!DO ji = 1, jpi
981	! DO jj = 1, jpj
982	! sshn_tl(ji,jj) = cos( (2.rpi)(FLOAT(ji)-0.5)/FLOAT(jpi) + rpi/2. ) &
983	! * sin( (2.rpi)(FLOAT(jj)-0.5)/FLOAT(jpj) ) / 100.
984	! END DO
985	!END DO
986
987	DO ji = 1, jpi
988	DO jj = 1, jpj
989	zsshn_tlin(ji,jj) = exp( -((float(ji)-float(jpi)/2.)/(float(jpi)/5.))**2 &
990	-((float(jj)-float(jpj)/2.)/(float(jpj)/5.))**2 ) / 100. &
991	* tmask(ji,jj,1)
992	END DO
993	END DO
994	sshn_tl(:,:) = zsshn_tlin(:,:)
995	ENDIF
996	#endif
997
998	!-----------------------------------------------------------------------
999	! Initialization of the dynamics and tracer fields for the tangent
1000	!-----------------------------------------------------------------------
1001
1002	CALL istate_init_tan
1003
1004	DO istp = nit000, nitend, 1
1005	CALL stp_tan( istp )
1006	END DO
1007
1008
1009	zun_tlout ( :,:,:) = un_tl (:,:,:)
1010	zvn_tlout ( :,:,:) = vn_tl (:,:,:)
1011	ztn_tlout ( :,:,:) = tn_tl (:,:,:)
1012	zsn_tlout ( :,:,:) = sn_tl (:,:,:)
1013	zsshn_tlout( :,:) = sshn_tl ( :,:)
1014
1015	#if defined key_obc
1016	zub_tlout (:,:,:) = ub_tl (:,:,:)
1017	zvb_tlout (:,:,:) = vb_tl (:,:,:)
1018	ztb_tlout (:,:,:) = tb_tl (:,:,:)
1019	zsb_tlout (:,:,:) = sb_tl (:,:,:)
1020	zsshb_tlout(:,:) = sshb_tl(:,:)
1021	#endif
1022
1023	!--------------------------------------------------------------------
1024	! Initialize the adjoint variables: dy^* = W dy
1025	!--------------------------------------------------------------------
1026
1027	DO jk = 1, jpk
1028	DO jj = nldj, nlej
1029	DO ji = nldi, nlei
1030	zun_adin(ji,jj,jk) = zun_tlout(ji,jj,jk) &
1031	& * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
1032	& * umask(ji,jj,jk) * wesp_u
1033	#if defined key_obc
1034	zub_adin(ji,jj,jk) = zub_tlout(ji,jj,jk) &
1035	& * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
1036	& * umask(ji,jj,jk) * wesp_u
1037	#endif
1038	END DO
1039	END DO
1040	END DO
1041
1042	DO jk = 1, jpk
1043	DO jj = nldj, nlej
1044	DO ji = nldi, nlei
1045	zvn_adin(ji,jj,jk) = zvn_tlout(ji,jj,jk) &
1046	& * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
1047	& * vmask(ji,jj,jk) * wesp_u
1048	#if defined key_obc
1049	zvb_adin(ji,jj,jk) = zvb_tlout(ji,jj,jk) &
1050	& * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
1051	& * vmask(ji,jj,jk) * wesp_u
1052	#endif
1053	END DO
1054	END DO
1055	END DO
1056
1057	DO jk = 1, jpk
1058	DO jj = nldj, nlej
1059	DO ji = nldi, nlei
1060	ztn_adin(ji,jj,jk) = ztn_tlout(ji,jj,jk) &
1061	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1062	& * tmask(ji,jj,jk) * wesp_t(jk)
1063	#if defined key_obc
1064	ztb_adin(ji,jj,jk) = ztb_tlout(ji,jj,jk) &
1065	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1066	& * tmask(ji,jj,jk) * wesp_t(jk)
1067	#endif
1068	END DO
1069	END DO
1070	END DO
1071
1072	DO jk = 1, jpk
1073	DO jj = nldj, nlej
1074	DO ji = nldi, nlei
1075	zsn_adin(ji,jj,jk) = zsn_tlout(ji,jj,jk) &
1076	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1077	& * tmask(ji,jj,jk) * wesp_s(jk)
1078	#if defined key_obc
1079	zsb_adin(ji,jj,jk) = zsb_tlout(ji,jj,jk) &
1080	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1081	& * tmask(ji,jj,jk) * wesp_s(jk)
1082	#endif
1083	END DO
1084	END DO
1085	END DO
1086
1087	DO jj = nldj, nlej
1088	DO ji = nldi, nlei
1089	zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) &
1090	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
1091	& * tmask(ji,jj,1) * wesp_ssh
1092	#if defined key_obc
1093	zsshb_adin(ji,jj) = zsshb_tlout(ji,jj) &
1094	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
1095	& * tmask(ji,jj,1) * wesp_ssh
1096	#endif
1097	END DO
1098	END DO
1099
1100	!--------------------------------------------------------------------
1101	! Compute the scalar product: ( L dx )^T W dy
1102	!--------------------------------------------------------------------
1103
1104
1105	zsp1_U = DOT_PRODUCT( zun_tlout , zun_adin )
1106	zsp1_V = DOT_PRODUCT( zvn_tlout , zvn_adin )
1107	zsp1_T = DOT_PRODUCT( ztn_tlout , ztn_adin )
1108	zsp1_S = DOT_PRODUCT( zsn_tlout , zsn_adin )
1109	zsp1_SSH = DOT_PRODUCT( zsshn_tlout, zsshn_adin )
1110
1111	zsp1 = zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH
1112
1113	#if defined key_obc
1114	zsp1_U = DOT_PRODUCT( zub_tlout , zub_adin )
1115	zsp1_V = DOT_PRODUCT( zvb_tlout , zvb_adin )
1116	zsp1_T = DOT_PRODUCT( ztb_tlout , ztb_adin )
1117	zsp1_S = DOT_PRODUCT( zsb_tlout , zsb_adin )
1118	zsp1_SSH = DOT_PRODUCT( zsshb_tlout, zsshb_adin )
1119
1120	zsp1 = zsp1 + ( zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH )
1121	#endif
1122	!--------------------------------------------------------------------
1123	! Call the adjoint routine: dx^* = L^T dy^*
1124	!--------------------------------------------------------------------
1125
1126	un_ad (:,:,:) = zun_adin (:,:,:)
1127	vn_ad (:,:,:) = zvn_adin (:,:,:)
1128	tn_ad (:,:,:) = ztn_adin (:,:,:)
1129	sn_ad (:,:,:) = zsn_adin (:,:,:)
1130	sshn_ad( :,:) = zsshn_adin ( :,:)
1131
1132	#if defined key_obc
1133	ub_ad (:,:,:) = zub_adin (:,:,:)
1134	vb_ad (:,:,:) = zvb_adin (:,:,:)
1135	tb_ad (:,:,:) = ztb_adin (:,:,:)
1136	sb_ad (:,:,:) = zsb_adin (:,:,:)
1137	sshb_ad(:,:) = zsshb_adin (:,:)
1138	#endif
1139
1140	DO istp = nitend, nit000, -1
1141	CALL stp_adj(istp)
1142	END DO
1143
1144	CALL istate_init_adj
1145
1146	zun_adout ( :,:,:) = un_ad (:,:,:)
1147	zvn_adout ( :,:,:) = vn_ad (:,:,:)
1148	ztn_adout ( :,:,:) = tn_ad (:,:,:)
1149	zsn_adout ( :,:,:) = sn_ad (:,:,:)
1150	zsshn_adout( :,:) = sshn_ad ( :,:)
1151
1152	zsp2_U = DOT_PRODUCT( zun_tlin , zun_adout )
1153	zsp2_V = DOT_PRODUCT( zvn_tlin , zvn_adout )
1154	zsp2_T = DOT_PRODUCT( ztn_tlin , ztn_adout )
1155	zsp2_S = DOT_PRODUCT( zsn_tlin , zsn_adout )
1156	zsp2_SSH = DOT_PRODUCT( zsshn_tlin, zsshn_adout )
1157
1158	zsp2 = zsp2_U + zsp2_V + zsp2_T + zsp2_S + zsp2_SSH
1159
1160	! 14 char:'12345678901234'
1161	SELECT CASE (jpert)
1162	CASE(1)
1163	cl_name = 'step_adj U'
1164	CASE(2)
1165	cl_name = 'step_adj V'
1166	CASE(3)
1167	cl_name = 'step_adj T'
1168	CASE(4)
1169	cl_name = 'step_adj S'
1170	CASE(5)
1171	cl_name = 'step_adj SSH'
1172	CASE(jpertmax)
1173	cl_name = 'step_adj '
1174	END SELECT
1175	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1176	endif
1177	END DO
1178
1179	DEALLOCATE( &
1180	& zun_tlin , zvn_tlin , ztn_tlin , &
1181	& zsn_tlin , zsshn_tlin , zun_tlout , &
1182	& zvn_tlout , ztn_tlout , zsn_tlout , &
1183	& zsshn_tlout, zun_adin , zvn_adin , &
1184	& ztn_adin , zsn_adin , zsshn_adin , &
1185	& zun_adout , zvn_adout , ztn_adout , &
1186	& zsn_adout , zsshn_adout, z2r , &
1187	& z3r &
1188	& )
1189
1190	END SUBROUTINE stp_adj_tst
1191	#if defined key_tst_tlm
1192	SUBROUTINE stp_tlm_tst( kumadt )
1193	!!-----------------------------------------------------------------------
1194	!!
1195	!! * ROUTINE example_tl_tst *
1196	!!
1197	!! ** Purpose : Test the tangent linear routine.
1198	!!
1199	!! ** Method : Verify the tangent with Taylor expansion
1200	!!
1201	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
1202	!!
1203	!! where L = tangent routine
1204	!! M = direct routine
1205	!! dx = input perturbation (random field)
1206	!! h = ration on perturbation
1207	!!
1208	!! In the tangent test we verify that:
1209	!! M(x+h*dx) - M(x)
1210	!! g(h) = ------------------ ---> 1 as h ---> 0
1211	!! L(h*dx)
1212	!! and
1213	!! g(h) - 1
1214	!! f(h) = ---------- ---> k (costant) as h ---> 0
1215	!! p
1216	!!
1217	!! History :
1218	!! ! 09-06 (F. Vigilant)
1219	!!-----------------------------------------------------------------------
1220	!! * Modules used
1221	USE opatam_tst_ini, ONLY: &
1222	& tlm_namrd
1223	USE sbcssr_tam, ONLY: &
1224	& qrp_tl,erp_tl
1225	USE sbcfwb_tam, ONLY: &
1226	& a_fwb_tl ! for before year.
1227	USE step_c1d ! Time stepping loop for the 1D configuration
1228	USE tamtrj ! writing out state trajectory
1229	USE lib_mpp, ONLY: & ! distributed memory computing
1230	& lk_mpp, &
1231	& mpp_max
1232	USE c1d, ONLY: & ! 1D configuration
1233	& lk_c1d
1234	USE par_tlm, ONLY: &
1235	& tlm_bch, &
1236	& cur_loop, &
1237	& h_ratio
1238	USE istate_mod
1239	USE wzvmod ! vertical velocity
1240	USE sbc_oce , ONLY: & ! ocean dynamics and tracers variables
1241	& emp, &
1242	& emps
1243	USE tamctl, ONLY: & ! Control parameters
1244	& numtan, numtan_sc
1245	USE step
1246	!! * Arguments
1247	INTEGER, INTENT(IN) :: &
1248	& kumadt ! Output unit
1249
1250	!! * Local declarations
1251	INTEGER :: &
1252	& ji, & ! dummy loop indices
1253	& jj, &
1254	& jk, &
1255	& istp, &
1256	& pgamma, & ! control the number of pertubation loops
1257	& igamma, &
1258	& nit_loop
1259	INTEGER, DIMENSION(jpi,jpj) :: &
1260	& iseed_2d ! 2D seed for the random number generator
1261	REAL(KIND=wp) :: &
1262	& zsp1, zsp2, zsp3, zzsp, & ! scalar product
1263	& zsp1_U, zsp1_V, zsp1_T, zsp1_S, zsp1_SSH, &
1264	& zsp2_U, zsp2_V, zsp2_T, zsp2_S, zsp2_SSH, &
1265	& zsp3_U, zsp3_V, zsp3_T, zsp3_S, zsp3_SSH, &
1266	& zzsp_U, zzsp_V, zzsp_T, zzsp_S, zzsp_SSH, &
1267	& gamma, &
1268	& zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, &
1269	& zgsp6, zgsp7
1270	REAL(KIND=wp) :: &
1271	& zgsp1_U, zgsp1_V, zgsp1_T, zgsp1_S, zgsp1_SSH, &
1272	& zgsp4_U, zgsp4_V, zgsp4_T, zgsp4_S, zgsp4_SSH, &
1273	& zgsp5_U, zgsp5_V, zgsp5_T, zgsp5_S, zgsp5_SSH
1274	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1275	& ta_tmp , &
1276	& sa_tmp , &
1277	& zun_tlin , & ! Pertubation input
1278	& zvn_tlin , & ! Pertubation input
1279	& ztn_tlin , & ! Pertubation input
1280	& zsn_tlin , & ! Pertubation input
1281	& zun_out , & ! Direct output
1282	& zvn_out , & ! Direct output
1283	& ztn_out , & ! Direct output
1284	& zsn_out , & ! Direct output
1285	& zun_wop , & ! Direct without pertubation output
1286	& zvn_wop , & ! Direct without pertubation output
1287	& ztn_wop , & ! Direct without pertubation output
1288	& zsn_wop , & ! Direct without pertubation output
1289	& z3r ! 3D random field
1290	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
1291	& emp_tmp, &
1292	& emps_tmp, &
1293	& sshn_tmp , &
1294	& zsshn_tlin , & ! Pertubation input
1295	& zsshn_out , & ! Direct output
1296	& zsshn_wop , & ! Direct without pertubation output
1297	& z2r ! 2D random field
1298
1299	CHARACTER(LEN=14) :: cl_name
1300	CHARACTER (LEN=128) :: file_out, file_wop, file_wop2, file_xdx
1301	CHARACTER (LEN=90) :: FMT
1302
1303	REAL(KIND=wp), DIMENSION(100):: &
1304	& zscu , zscv , zsct , zscs , zscssh, &
1305	& zscerru, zscerrv, zscerrt, zscerrs, zscerrssh
1306	INTEGER, DIMENSION(100):: &
1307	& iiposu, iiposv, iipost, iiposs, iiposssh, &
1308	& ijposu, ijposv, ijpost, ijposs, ijposssh, &
1309	& ikposu, ikposv, ikpost, ikposs
1310	INTEGER:: &
1311	& ii, &
1312	& isamp=40, &
1313	& jsamp=40, &
1314	& ksamp=10, &
1315	& numsctlm
1316	REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: &
1317	& zerru, zerrv, zerrt, zerrs ! Pertubation input
1318	REAL(KIND=wp), DIMENSION(jpi,jpj) :: &
1319	& zerrssh ! Pertubation input
1320	INTEGER :: &
1321	& jdir, jdirmax=6 ! Perturbation direction (U, V, T, S, SSH, all)
1322
1323	! Allocate memory
1324
1325	ALLOCATE( &
1326	& zun_tlin( jpi,jpj,jpk), &
1327	& zvn_tlin( jpi,jpj,jpk), &
1328	& ztn_tlin( jpi,jpj,jpk), &
1329	& zsn_tlin( jpi,jpj,jpk), &
1330	& zsshn_tlin( jpi,jpj ), &
1331	& zun_out( jpi,jpj,jpk), &
1332	& zvn_out( jpi,jpj,jpk), &
1333	& ztn_out( jpi,jpj,jpk), &
1334	& zsn_out( jpi,jpj,jpk), &
1335	& zsshn_out( jpi,jpj ), &
1336	& zun_wop( jpi,jpj,jpk), &
1337	& zvn_wop( jpi,jpj,jpk), &
1338	& ztn_wop( jpi,jpj,jpk), &
1339	& zsn_wop( jpi,jpj,jpk), &
1340	& zsshn_wop( jpi,jpj ), &
1341	& z2r( jpi,jpj ), &
1342	& z3r( jpi,jpj,jpk), &
1343	& emp_tmp( jpi,jpj ), &
1344	& emps_tmp( jpi,jpj ), &
1345	& ta_tmp( jpi,jpj,jpk), &
1346	& sa_tmp( jpi,jpj,jpk), &
1347	& sshn_tmp( jpi,jpj ) &
1348	& )
1349
1350	!--------------------------------------------------------------------
1351	! Reset the scalar
1352	!--------------------------------------------------------------------
1353	zscu (:) = 0.0_wp
1354	zscv (:) = 0.0_wp
1355	zsct (:) = 0.0_wp
1356	zscs (:) = 0.0_wp
1357	zscssh (:) = 0.0_wp
1358
1359	iiposu (:) = 0
1360	iiposv (:) = 0
1361	iipost (:) = 0
1362	iiposs (:) = 0
1363	iiposssh (:) = 0
1364	ijposu (:) = 0
1365	ijposv (:) = 0
1366	ijpost (:) = 0
1367	ijposs (:) = 0
1368	ijposssh (:) = 0
1369	ikposu (:) = 0
1370	ikposv (:) = 0
1371	ikpost (:) = 0
1372	ikposs (:) = 0
1373
1374	zscerru (:) = 0.0_wp
1375	zscerrv (:) = 0.0_wp
1376	zscerrt (:) = 0.0_wp
1377	zscerrs (:) = 0.0_wp
1378	zscerrssh (:) = 0.0_wp
1379
1380	zerru (:,:,:) = 0.0_wp
1381	zerrv (:,:,:) = 0.0_wp
1382	zerrt (:,:,:) = 0.0_wp
1383	zerrs (:,:,:) = 0.0_wp
1384	zerrssh (:,:) = 0.0_wp
1385
1386	!--------------------------------------------------------------------
1387	! Reset the tangent variables
1388	!--------------------------------------------------------------------
1389	zun_tlin ( :,:,:) = 0.0_wp
1390	zvn_tlin ( :,:,:) = 0.0_wp
1391	ztn_tlin ( :,:,:) = 0.0_wp
1392	zsn_tlin ( :,:,:) = 0.0_wp
1393	zsshn_tlin ( :,:) = 0.0_wp
1394	z3r ( :,:,:) = 0.0_wp
1395	z2r ( :,:) = 0.0_wp
1396
1397	!--------------------------------------------------------------------
1398	! Output filename Xn=F(X0)
1399	!--------------------------------------------------------------------
1400	CALL tlm_namrd
1401	gamma = h_ratio
1402	file_wop='trj_wop_step'
1403	file_xdx='trj_xdx_step'
1404	!--------------------------------------------------------------------
1405	! Initialize the tangent input with random noise: dx
1406	!--------------------------------------------------------------------
1407	jdir = jdirmax
1408	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
1409	IF (( jdir .EQ. 1 ) .OR. ( jdir .EQ. jdirmax ) ) THEN
1410	CALL grid_rd_sd( 596035, z3r, 'U', 0.0_wp, stdu)
1411	DO jk = 1, jpk
1412	DO jj = nldj, nlej
1413	DO ji = nldi, nlei
1414	zun_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1415	END DO
1416	END DO
1417	END DO
1418
1419	ENDIF
1420	IF (( jdir .EQ. 2 ) .OR. ( jdir .EQ. jdirmax ) ) THEN
1421	CALL grid_rd_sd( 371836, z3r, 'V', 0.0_wp, stdv)
1422	DO jk = 1, jpk
1423	DO jj = nldj, nlej
1424	DO ji = nldi, nlei
1425	zvn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1426	END DO
1427	END DO
1428	END DO
1429
1430	ENDIF
1431	IF (( jdir .EQ. 3 ) .OR. ( jdir .EQ. jdirmax ) ) THEN
1432	CALL grid_rd_sd( 284035, z3r, 'T', 0.0_wp, stdt)
1433	DO jk = 1, jpk
1434	DO jj = nldj, nlej
1435	DO ji = nldi, nlei
1436	ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1437	END DO
1438	END DO
1439	END DO
1440
1441	ENDIF
1442	IF (( jdir .EQ. 4 ) .OR. ( jdir .EQ. jdirmax ) ) THEN
1443	CALL grid_rd_sd( 471426, z3r, 'T', 0.0_wp, stds)
1444	DO jk = 1, jpk
1445	DO jj = nldj, nlej
1446	DO ji = nldi, nlei
1447	zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
1448	END DO
1449	END DO
1450	END DO
1451
1452	ENDIF
1453	IF (( jdir .EQ. 5 ) .OR. ( jdir .EQ. jdirmax ) ) THEN
1454	CALL grid_rd_sd( 712035, z2r,'T', 0.0_wp, stdssh)
1455	DO jj = nldj, nlej
1456	DO ji = nldi, nlei
1457	zsshn_tlin(ji,jj) = z2r(ji,jj)
1458	END DO
1459	END DO
1460
1461	ENDIF
1462	ENDIF
1463	IF ( tlm_bch /= 2 ) CALL istate_p
1464
1465	!--------------------------------------------------------------------
1466	! Open File to save scalar result
1467	!--------------------------------------------------------------------
1468
1469	! check that all process are still there... If some process have an error,
1470	! they will never enter in step and other processes will wait until the end of the cpu time!
1471	IF( lk_mpp ) CALL mpp_max( nstop )
1472
1473	istp = nit000
1474	IF ( tlm_bch /= 2 ) THEN
1475	IF( lk_c1d ) THEN ! 1D configuration (no AGRIF zoom)
1476	DO WHILE ( istp <= nitend .AND. nstop == 0 )
1477	CALL stp_c1d( istp )
1478	istp = istp + 1
1479	END DO
1480	ENDIF
1481	ENDIF
1482
1483	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
1484	zun_tlin(:,:,:) = gamma * zun_tlin(:,:,:)
1485	un(:,:,:) = un(:,:,:) + zun_tlin(:,:,:)
1486
1487	zvn_tlin(:,:,:) = gamma * zvn_tlin(:,:,:)
1488	vn(:,:,:) = vn(:,:,:) + zvn_tlin(:,:,:)
1489
1490	ztn_tlin(:,:,:) = gamma * ztn_tlin(:,:,:)
1491	tn(:,:,:) = tn(:,:,:) + ztn_tlin(:,:,:)
1492
1493	zsn_tlin(:,:,:) = gamma * zsn_tlin(:,:,:)
1494	sn(:,:,:) = sn(:,:,:) + zsn_tlin(:,:,:)
1495
1496	zsshn_tlin(:,:) = gamma * zsshn_tlin(:,:)
1497	sshn (:,:) = sshn (:,:) + zsshn_tlin(:,:)
1498	ENDIF
1499
1500	IF( .NOT. lk_vvl ) CALL wzv(nit000)
1501
1502	!--------------------------------------------------------------------
1503	! Compute the direct model F(X0,t=n) = Xn
1504	!--------------------------------------------------------------------
1505	IF ( tlm_bch /= 2 ) THEN
1506	DO istp = nit000, nitend, 1
1507	CALL stp( istp )
1508	END DO
1509	ENDIF
1510	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
1511	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
1512	!--------------------------------------------------------------------
1513	! Compute the Tangent
1514	!--------------------------------------------------------------------
1515	IF ( tlm_bch == 2 ) THEN
1516
1517	!--------------------------------------------------------------------
1518	! Initialize the tangent variables: dy^* = W dy
1519	!--------------------------------------------------------------------
1520	qrp_tl = 0.0_wp
1521	#if defined key_tradmp
1522	strdmp_tl = 0.0_wp
1523	ttrdmp_tl = 0.0_wp
1524	#endif
1525	a_fwb_tl = 0.0_wp
1526
1527	CALL trj_rea( nit000-1, 1 )
1528
1529	un_tl (:,:,:) = zun_tlin (:,:,:)
1530	vn_tl (:,:,:) = zvn_tlin (:,:,:)
1531	tn_tl (:,:,:) = ztn_tlin (:,:,:)
1532	sn_tl (:,:,:) = zsn_tlin (:,:,:)
1533	sshn_tl( :,:) = zsshn_tlin ( :,:)
1534
1535
1536	!-----------------------------------------------------------------------
1537	! Initialization of the dynamics and tracer fields for the tangent
1538	!-----------------------------------------------------------------------
1539
1540	CALL istate_init_tan
1541
1542	DO istp = nit000, nitend, 1
1543	CALL stp_tan( istp )
1544	!CALL stp_tan_cpd( istp )
1545	END DO
1546
1547	!--------------------------------------------------------------------
1548	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
1549	!--------------------------------------------------------------------
1550
1551	zsp2_U = DOT_PRODUCT( un_tl , un_tl )
1552	zsp2_V = DOT_PRODUCT( vn_tl , vn_tl )
1553	zsp2_T = DOT_PRODUCT( tn_tl , tn_tl )
1554	zsp2_S = DOT_PRODUCT( sn_tl , sn_tl )
1555	zsp2_SSH = DOT_PRODUCT( sshn_tl, sshn_tl )
1556
1557	zsp2 = zsp2_U + zsp2_V + zsp2_T + zsp2_S + zsp2_SSH
1558
1559	!--------------------------------------------------------------------
1560	! Storing data
1561	!--------------------------------------------------------------------
1562	CALL trj_rd_spl(file_wop)
1563	zun_wop (:,:,:) = un (:,:,:)
1564	zvn_wop (:,:,:) = vn (:,:,:)
1565	ztn_wop (:,:,:) = tn (:,:,:)
1566	zsn_wop (:,:,:) = sn (:,:,:)
1567	zsshn_wop (:,: ) = sshn (:,: )
1568	CALL trj_rd_spl(file_xdx)
1569	zun_out (:,:,:) = un (:,:,:)
1570	zvn_out (:,:,:) = vn (:,:,:)
1571	ztn_out (:,:,:) = tn (:,:,:)
1572	zsn_out (:,:,:) = sn (:,:,:)
1573	zsshn_out (:,: ) = sshn (:,: )
1574	!--------------------------------------------------------------------
1575	! Compute the Linearization Error
1576	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
1577	! and
1578	! Compute the Linearization Error
1579	! En = Nn -TL(gamma.dX0, t0,tn)
1580	!--------------------------------------------------------------------
1581	! Warning: Here we re-use local variables z()_out and z()_wop
1582	ii=0
1583	DO jk = 1, jpk
1584	DO jj = 1, jpj
1585	DO ji = 1, jpi
1586	zun_out (ji,jj,jk) = zun_out (ji,jj,jk) - zun_wop (ji,jj,jk)
1587	zun_wop (ji,jj,jk) = zun_out (ji,jj,jk) - un_tl (ji,jj,jk)
1588	IF ( un_tl(ji,jj,jk) .NE. 0.0_wp ) &
1589	& zerru(ji,jj,jk) = zun_out(ji,jj,jk)/un_tl(ji,jj,jk)
1590	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1591	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1592	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1593	ii = ii+1
1594	iiposu(ii) = ji
1595	ijposu(ii) = jj
1596	ikposu(ii) = jk
1597	IF ( INT(umask(ji,jj,jk)) .NE. 0) THEN
1598	zscu (ii) = zun_out(ji,jj,jk)
1599	zscerru (ii) = ( zerru(ji,jj,jk) - 1.0_wp ) / gamma
1600	ENDIF
1601	ENDIF
1602	END DO
1603	END DO
1604	END DO
1605	ii=0
1606	DO jk = 1, jpk
1607	DO jj = 1, jpj
1608	DO ji = 1, jpi
1609	zvn_out (ji,jj,jk) = zvn_out (ji,jj,jk) - zvn_wop (ji,jj,jk)
1610	zvn_wop (ji,jj,jk) = zvn_out (ji,jj,jk) - vn_tl (ji,jj,jk)
1611	IF ( vn_tl(ji,jj,jk) .NE. 0.0_wp ) &
1612	& zerrv(ji,jj,jk) = zvn_out(ji,jj,jk)/vn_tl(ji,jj,jk)
1613	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1614	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1615	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1616	ii = ii+1
1617	iiposv(ii) = ji
1618	ijposv(ii) = jj
1619	ikposv(ii) = jk
1620	IF ( INT(vmask(ji,jj,jk)) .NE. 0) THEN
1621	zscv (ii) = zvn_out (ji,jj,jk)
1622	zscerrv (ii) = ( zerrv(ji,jj,jk) - 1.0_wp ) / gamma
1623	ENDIF
1624	ENDIF
1625	END DO
1626	END DO
1627	END DO
1628	ii=0
1629	DO jk = 1, jpk
1630	DO jj = 1, jpj
1631	DO ji = 1, jpi
1632	ztn_out (ji,jj,jk) = ztn_out (ji,jj,jk) - ztn_wop (ji,jj,jk)
1633	ztn_wop (ji,jj,jk) = ztn_out (ji,jj,jk) - tn_tl (ji,jj,jk)
1634	IF ( tn_tl(ji,jj,jk) .NE. 0.0_wp ) &
1635	& zerrt(ji,jj,jk) = ztn_out(ji,jj,jk)/tn_tl(ji,jj,jk)
1636	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1637	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1638	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1639	ii = ii+1
1640	iipost(ii) = ji
1641	ijpost(ii) = jj
1642	ikpost(ii) = jk
1643	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
1644	zsct (ii) = ztn_out (ji,jj,jk)
1645	zscerrt (ii) = ( zerrt(ji,jj,jk) - 1.0_wp ) / gamma
1646	ENDIF
1647	ENDIF
1648	END DO
1649	END DO
1650	END DO
1651	ii=0
1652	DO jk = 1, jpk
1653	DO jj = 1, jpj
1654	DO ji = 1, jpi
1655	zsn_out (ji,jj,jk) = zsn_out (ji,jj,jk) - zsn_wop (ji,jj,jk)
1656	zsn_wop (ji,jj,jk) = zsn_out (ji,jj,jk) - sn_tl (ji,jj,jk)
1657	IF ( sn_tl(ji,jj,jk) .NE. 0.0_wp ) &
1658	& zerrs(ji,jj,jk) = zsn_out(ji,jj,jk)/sn_tl(ji,jj,jk)
1659	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1660	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1661	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1662	ii = ii+1
1663	iiposs(ii) = ji
1664	ijposs(ii) = jj
1665	ikposs(ii) = jk
1666	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
1667	zscs (ii) = zsn_out (ji,jj,jk)
1668	zscerrs (ii) = ( zerrs(ji,jj,jk) - 1.0_wp ) / gamma
1669	ENDIF
1670	ENDIF
1671	END DO
1672	END DO
1673	END DO
1674	ii=0
1675	jk=1
1676	DO jj = 1, jpj
1677	DO ji = 1, jpi
1678	zsshn_out (ji,jj ) = zsshn_out (ji,jj ) - zsshn_wop (ji,jj )
1679	zsshn_wop (ji,jj ) = zsshn_out (ji,jj ) - sshn_tl (ji,jj )
1680	IF ( sshn_tl(ji,jj) .NE. 0.0_wp ) &
1681	& zerrssh(ji,jj) = zsshn_out(ji,jj)/sshn_tl(ji,jj)
1682	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1683	& (MOD(jj, jsamp) .EQ. 0) ) THEN
1684	ii = ii+1
1685	iiposssh(ii) = ji
1686	ijposssh(ii) = jj
1687	IF ( INT(tmask(ji,jj,jk)) .NE. 0 ) THEN
1688	zscssh (ii) = zsshn_out (ji,jj)
1689	zscerrssh (ii) = ( zerrssh(ji,jj) - 1.0_wp ) / gamma
1690	ENDIF
1691	ENDIF
1692	END DO
1693	END DO
1694
1695	zsp1_U = DOT_PRODUCT( zun_out , zun_out )
1696	zsp1_V = DOT_PRODUCT( zvn_out , zvn_out )
1697	zsp1_T = DOT_PRODUCT( ztn_out , ztn_out )
1698	zsp1_S = DOT_PRODUCT( zsn_out , zsn_out )
1699	zsp1_SSH = DOT_PRODUCT( zsshn_out, zsshn_out )
1700
1701	zsp1 = zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH
1702
1703	zsp3_U = DOT_PRODUCT( zun_wop , zun_wop )
1704	zsp3_V = DOT_PRODUCT( zvn_wop , zvn_wop )
1705	zsp3_T = DOT_PRODUCT( ztn_wop , ztn_wop )
1706	zsp3_S = DOT_PRODUCT( zsn_wop , zsn_wop )
1707	zsp3_SSH = DOT_PRODUCT( zsshn_wop, zsshn_wop )
1708
1709	zsp3 = zsp3_U + zsp3_V + zsp3_T + zsp3_S + zsp3_SSH
1710
1711	!--------------------------------------------------------------------
1712	! Print the linearization error En - norme 2
1713	!--------------------------------------------------------------------
1714	! 14 char:'12345678901234'
1715	cl_name = 'step_tam:En '
1716	zzsp = SQRT(zsp3)
1717	zzsp_U = SQRT(zsp3_U)
1718	zzsp_V = SQRT(zsp3_V)
1719	zzsp_T = SQRT(zsp3_T)
1720	zzsp_S = SQRT(zsp3_S)
1721	zzsp_SSH = SQRT(zsp3_SSH)
1722	zgsp5 = zzsp
1723	zgsp5_U=zzsp_U
1724	zgsp5_V=zzsp_V
1725	zgsp5_T=zzsp_T
1726	zgsp5_S=zzsp_S
1727	zgsp5_SSH=zzsp_SSH
1728	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1729
1730	CALL prntst_tlm( 'U ', kumadt, zzsp_U , h_ratio )
1731	CALL prntst_tlm( 'V ', kumadt, zzsp_V , h_ratio )
1732	CALL prntst_tlm( 'T ', kumadt, zzsp_T , h_ratio )
1733	CALL prntst_tlm( 'S ', kumadt, zzsp_S , h_ratio )
1734	CALL prntst_tlm( 'SSH ', kumadt, zzsp_SSH, h_ratio )
1735
1736	!--------------------------------------------------------------------
1737	! Compute TLM norm2
1738	!--------------------------------------------------------------------
1739	zzsp = SQRT(zsp2)
1740	zzsp_U = SQRT(zsp2_U)
1741	zzsp_V = SQRT(zsp2_V)
1742	zzsp_T = SQRT(zsp2_T)
1743	zzsp_S = SQRT(zsp2_S)
1744	zzsp_SSH = SQRT(zsp2_SSH)
1745	zgsp4 = zzsp
1746	zgsp4_U=zzsp_U
1747	zgsp4_V=zzsp_V
1748	zgsp4_T=zzsp_T
1749	zgsp4_S=zzsp_S
1750	zgsp4_SSH=zzsp_SSH
1751	cl_name = 'step_tam:Ln2 '
1752	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1753
1754	CALL prntst_tlm( 'U ', kumadt, zzsp_U , h_ratio )
1755	CALL prntst_tlm( 'V ', kumadt, zzsp_V , h_ratio )
1756	CALL prntst_tlm( 'T ', kumadt, zzsp_T , h_ratio )
1757	CALL prntst_tlm( 'S ', kumadt, zzsp_S , h_ratio )
1758	CALL prntst_tlm( 'SSH ', kumadt, zzsp_SSH, h_ratio )
1759
1760	!--------------------------------------------------------------------
1761	! Print the linearization error Nn - norme 2
1762	!--------------------------------------------------------------------
1763	! 14 char:'12345678901234'
1764	cl_name = 'step_tam:Nn '
1765	zzsp = SQRT(zsp1)
1766	zzsp_U = SQRT(zsp1_U)
1767	zzsp_V = SQRT(zsp1_V)
1768	zzsp_T = SQRT(zsp1_T)
1769	zzsp_S = SQRT(zsp1_S)
1770	zzsp_SSH = SQRT(zsp1_SSH)
1771	! zzsp = zzsp_U + zzsp_V + zzsp_T + zzsp_S + zzsp_SSH
1772	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1773
1774	CALL prntst_tlm( 'U ', kumadt, zzsp_U , h_ratio )
1775	CALL prntst_tlm( 'V ', kumadt, zzsp_V , h_ratio )
1776	CALL prntst_tlm( 'T ', kumadt, zzsp_T , h_ratio )
1777	CALL prntst_tlm( 'S ', kumadt, zzsp_S , h_ratio )
1778	CALL prntst_tlm( 'SSH ', kumadt, zzsp_SSH, h_ratio )
1779
1780	zgsp3 = SQRT( zsp3/zsp2 )
1781	zgsp7 = zgsp3/gamma
1782	zgsp1 = zzsp
1783	zgsp1_U=zzsp_U
1784	zgsp1_V=zzsp_V
1785	zgsp1_T=zzsp_T
1786	zgsp1_S=zzsp_S
1787	zgsp1_SSH=zzsp_SSH
1788	zgsp2 = zgsp1 / zgsp4
1789	zgsp6 = (zgsp2 - 1.0_wp)/gamma
1790
1791	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)"
1792
1793	WRITE(numtan,FMT) 'step ', cur_loop, h_ratio, zgsp1, zgsp1_T ,zgsp4_T, zgsp5_T,zgsp1_S,zgsp4_S,zgsp5_S!,
1794	!--------------------------------------------------------------------
1795	! Unitary calculus
1796	!--------------------------------------------------------------------
1797	FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)"
1798	cl_name = 'step '
1799	IF (lwp) THEN
1800	DO ii=1, 100, 1
1801	IF ( zscu(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscu ', &
1802	& cur_loop, h_ratio, iiposu(ii), ijposu(ii), ikposu(ii), zscu(ii)
1803	ENDDO
1804	DO ii=1, 100, 1
1805	IF ( zscv(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscv ', &
1806	& cur_loop, h_ratio, iiposv(ii), ijposv(ii), ikposv(ii), zscv(ii)
1807	ENDDO
1808	DO ii=1, 100, 1
1809	IF ( zsct(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zsct ', &
1810	& cur_loop, h_ratio, iipost(ii), ijpost(ii), ikpost(ii), zsct(ii)
1811	ENDDO
1812	DO ii=1, 100, 1
1813	IF ( zscs(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscs ', &
1814	& cur_loop, h_ratio, iiposs(ii), ijposs(ii), ikposs(ii), zscs(ii)
1815	ENDDO
1816	DO ii=1, 100, 1
1817	IF ( zscssh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscssh ', &
1818	& cur_loop, h_ratio, iiposssh(ii), ijposssh(ii), ijposssh(ii), zscssh(ii)
1819	ENDDO
1820
1821	DO ii=1, 100, 1
1822	IF ( zscerru(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerru ', &
1823	& cur_loop, h_ratio, iiposu(ii), ijposu(ii), ikposu(ii), zscerru(ii)
1824	ENDDO
1825	DO ii=1, 100, 1
1826	IF ( zscerrv(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrv ', &
1827	& cur_loop, h_ratio, iiposv(ii), ijposv(ii), ikposv(ii), zscerrv(ii)
1828	ENDDO
1829	DO ii=1, 100, 1
1830	IF ( zscerrt(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrt ', &
1831	& cur_loop, h_ratio, iipost(ii), ijpost(ii), ikpost(ii), zscerrt(ii)
1832	ENDDO
1833	DO ii=1, 100, 1
1834	IF ( zscerrs(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrs ', &
1835	& cur_loop, h_ratio, iiposs(ii), ijposs(ii), ikposs(ii), zscerrs(ii)
1836	ENDDO
1837	DO ii=1, 100, 1
1838	IF ( zscerrssh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerssh', &
1839	& cur_loop, h_ratio, iiposssh(ii), ijposssh(ii), ijposssh(ii), zscerrssh(ii)
1840	ENDDO
1841	! write separator
1842	WRITE(numtan_sc,"(A4)") '===='
1843	ENDIF
1844
1845	ENDIF
1846
1847	DEALLOCATE( &
1848	& zun_tlin , &
1849	& zvn_tlin , &
1850	& ztn_tlin , &
1851	& zsn_tlin , &
1852	& zsshn_tlin , &
1853	& zun_out , &
1854	& zvn_out , &
1855	& ztn_out , &
1856	& zsn_out , &
1857	& zsshn_out , &
1858	& zun_wop , &
1859	& zvn_wop , &
1860	& ztn_wop , &
1861	& zsn_wop , &
1862	& zsshn_wop , &
1863	& z2r , &
1864	& z3r &
1865	& )
1866
1867	END SUBROUTINE stp_tlm_tst
1868
1869	!!----------------------------------------------------------------------
1870
1871	SUBROUTINE asm_trj_wop_tasa_wri
1872	!!-----------------------------------------------------------------------
1873	!!
1874	!! * ROUTINE asm_trj_wop_tasa_wri *
1875	!!
1876	!! ** Purpose : Write to file the model state trajectory for use with
1877	!! 4D-Var.
1878	!!
1879	!! ** Method :
1880	!!
1881	!! ** Action :
1882	!!
1883	!! History :
1884	!! ! 09-07 (F. Vigilant)
1885	!!-----------------------------------------------------------------------
1886	!! *Module udes
1887	USE iom
1888	USE oce , ONLY: & ! ocean dynamics and tracers variables
1889	& tn, sn
1890	!! * Arguments
1891	!! * Local declarations
1892	INTEGER :: &
1893	& inum ! File unit number
1894	CHARACTER (LEN=50) :: &
1895	& filename
1896	! Define the output file
1897	filename='trj_wop_tasa'
1898	WRITE(filename, FMT='(A,A)' ) TRIM( filename ), '.nc'
1899	filename = TRIM( filename )
1900	CALL iom_open( filename, inum, ldwrt = .TRUE., kiolib = jprstlib)
1901
1902	! Output trajectory fields
1903	CALL iom_rstput( 1, 1, inum, 'tn' , tn )
1904	CALL iom_rstput( 1, 1, inum, 'sn' , sn )
1905	CALL iom_close( inum )
1906	END SUBROUTINE asm_trj_wop_tasa_wri
1907
1908	SUBROUTINE asm_trj_wop_tasa_rd
1909	!!-----------------------------------------------------------------------
1910	!!
1911	!! * ROUTINE asm_trj_wop_tasa_rd *
1912	!!
1913	!! ** Purpose : Read file the model state trajectory for use with
1914	!! 4D-Var.
1915	!!
1916	!! ** Method :
1917	!!
1918	!! ** Action :
1919	!!
1920	!! History :
1921	!! ! 09-07 (F. Vigilant)
1922	!!-----------------------------------------------------------------------
1923	!! *Module udes
1924	USE iom ! I/O module
1925	USE oce , ONLY: & ! ocean dynamics and tracers variables
1926	& tn, sn
1927	!! * Arguments
1928	!! * Local declarations
1929	INTEGER :: &
1930	& inum ! File unit number
1931	CHARACTER (LEN=50) :: &
1932	& filename
1933	! Define the output file
1934	filename='trj_wop_tasa'
1935	WRITE(filename, FMT='(A,A)' ) TRIM( filename ), '.nc'
1936	filename = TRIM( filename )
1937	CALL iom_open( filename, inum)
1938
1939	! Output trajectory fields
1940	CALL iom_get( inum, jpdom_data, 'tn' , tn, 1 )
1941	CALL iom_get( inum, jpdom_data, 'sn' , sn, 1 )
1942	CALL iom_close( inum )
1943	END SUBROUTINE asm_trj_wop_tasa_rd
1944
1945
1946	!!======================================================================
1947	#endif
1948	#endif
1949	END MODULE step_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/step_tam.F90 @ 2587

Download in other formats: