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.

step_tam.F90 in tags/TAM_v3_0/NEMOTAM/OPATAM_SRC – NEMO

Context Navigation

source: tags/TAM_v3_0/NEMOTAM/OPATAM_SRC/step_tam.F90 @ 8068

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

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

Download in other formats: