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 NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC – NEMO

Context Navigation

source: NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC/step_tam.F90 @ 13432

Last change on this file since 13432 was 4589, checked in by pabouttier, 10 years ago
Remove useless options and cosmetic changes in step_tam.F90
Property svn:executable set to ``*
File size: 37.9 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	!! ! 07-12 (P.-A. Bouttier) Phasing with 3.4 version
31	!!----------------------------------------------------------------------
32	!! stp_tam : OPA system time-stepping (tangent linear)
33	!! stp_adj : OPA system time-stepping (adjoint)
34	!!----------------------------------------------------------------------
35
36	USE step_oce_tam
37	#if defined key_agrif
38	#error 'agrif not yet implemented in nemotam'
39	#endif
40
41	IMPLICIT NONE
42	PRIVATE
43
44	PUBLIC stp_tan, &
45	& stp_adj, & ! called by simvar.F90
46	& stp_adj_tst
47
48	!! * Substitutions
49	# include "domzgr_substitute.h90"
50	# include "zdfddm_substitute.h90"
51
52	CONTAINS
53	SUBROUTINE stp_tan( kstp )
54	!!----------------------------------------------------------------------
55	!! * ROUTINE stp_tan *
56	!!
57	!! ** Purpose of the direct routine:
58	!! - Time stepping of OPA (momentum and active tracer eqs.)
59	!! - Time stepping of LIM (dynamic and thermodynamic eqs.)
60	!! - Tme stepping of TRC (passive tracer eqs.)
61	!!
62	!! ** Method of the direct routine:
63	!! -1- Update forcings and data
64	!! -2- Update ocean physics
65	!! -3- Compute the t and s trends
66	!! -4- Update t and s
67	!! -5- Compute the momentum trends
68	!! -6- Update the horizontal velocity
69	!! -7- Compute the diagnostics variables (rd,N2, div,cur,w)
70	!! -8- Outputs and diagnostics
71	!!----------------------------------------------------------------------
72	!! * Arguments
73	INTEGER, INTENT( in ) :: kstp ! ocean time-step index
74
75	!! * local declarations
76	INTEGER :: indic ! error indicator if < 0
77	!! ---------------------------------------------------------------------
78
79	indic = 0 ! reset to no error condition
80	IF ( kstp == nit000 ) CALL tl_trj_wri(nit000-1)
81	IF ( kstp /= nit000 ) CALL day_tam( kstp, 0 ) ! Calendar (day was already called at nit000 in day_init)
82
83	CALL iom_setkt( kstp ) ! say to iom that we are at time step kstp
84	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
85	! Update data, open boundaries, surface boundary condition (including sea-ice)
86	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
87
88	! update 3D temperature data ... not needed in tangent
89
90	! update 3D salinity data ... not needed in tangent (to be investigated, see sbc_ssr)
91	CALL sbc_tan ( kstp ) ! Sea boundary condition (including sea-ice)
92
93	! Output the initial state and forcings ... not needed in tangent
94
95	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
96	! Ocean dynamics : ssh, wn, hdiv, rot !
97	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
98
99	CALL ssh_wzv_tan( kstp ) ! after ssh & vertical velocity
100
101	! saving direct variables ua,va, ta, sa before entering in tracer
102
103	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
104	! Ocean physics update (ua, va, ta, sa used as workspace)
105	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
106	CALL bn2_tan( tsb, tsb_tl, rn2b_tl ) ! now Brunt-Vaisala frequency
107	CALL bn2_tan( tsn, tsn_tl, rn2_tl ) ! now Brunt-Vaisala frequency
108
109	!-----------------------------------------------------------------------
110	! VERTICAL PHYSICS
111	!-----------------------------------------------------------------------
112	CALL zdf_bfr_tan ( kstp ) ! bottom friction...
113
114	! ! Vertical eddy viscosity and diffusivity coefficients
115	! Richardson number dependent Kz ... not available
116	! TKE closure scheme for Kz ... not available
117	! KPP closure scheme for Kz ... not available
118	! Constant Kz (reset avt, avm[uv] to the background value)...
119	! increase diffusivity at rivers mouths... not needed in tangent
120	! enhanced vertical eddy diffusivity ... not needed in tangent with lk_zdfcst_tan
121	! double diffusive mixing ... not needed in tangent with lk_zdfcst_tan
122	! mixed layer depth... not needed in tangent with lk_zdfcst_tan
123	!-----------------------------------------------------------------------
124	! LATERAL PHYSICS
125	!-----------------------------------------------------------------------
126	! N.B. ua, va, ta, sa arrays are used as workspace in this section
127	!-----------------------------------------------------------------------
128	! before slope of the lateral mixing... not needed in tangent with lk_zdfcst_tan
129	#if defined key_traldf_c2d
130	! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_tan
131	#endif
132	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
133	! diagnostics and outputs (ua, va, ta, sa used as workspace)
134	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
135	! not needed in tangent
136	#if defined key_top
137	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
138	! Passive Tracer Model
139	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
140	! N.B. ua, va, ta, sa arrays are used as workspace in this section
141	!-----------------------------------------------------------------------
142
143	! time-stepping... not needed in tangent for the time being
144
145	#endif
146
147	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
148	! Active tracers (ua, va used as workspace)
149	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
150	tsa_tl(:,:,:,:) = 0.0_wp ! set tracer trends to zero
151
152	! apply tracer assimilation increment ... not needed in tangent
153	CALL tra_sbc_tan( kstp ) ! surface boundary condition
154	!IF( ln_traqsr ) CALL tra_qsr_tan( kstp ) ! penetrative solar radiation qsr
155	IF( ln_trabbc ) CALL tra_bbc_tan( kstp ) ! bottom heat flux
156	IF( lk_trabbl ) CALL tra_bbl_tan( kstp ) ! diffusive bottom boundary layer scheme
157	!! advective (and/or diffusive) bottom boundary layer scheme ... currently not available
158	IF( ln_tradmp ) CALL tra_dmp_tan( kstp ) ! internal damping trends
159	CALL tra_adv_tan( kstp ) ! horizontal & vertical advection
160	CALL tra_ldf_tan( kstp ) ! lateral mixing
161	CALL tra_zdf_tan( kstp ) ! vertical mixing
162
163	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos)
164	! update the new (t,s) fields by non
165	! penetrative convective adjustment ... not available
166	CALL tra_nxt_tan( kstp ) ! tracer fields at next time step
167	CALL eos_tan( tsa, tsa_tl, & ! Time-filtered in situ density used in dynhpg module
168	& rhd_tl, rhop_tl )
169	IF( ln_zps ) CALL zps_hde_tan( kstp, jpts, tsa, tsa_tl, rhd_tl, &
170	& gtsu_tl, gru_tl, gtsv_tl, & ! Partial steps: time filtered hor. gradient
171	& grv_tl ) ! of t, s, rd at the bottom ocean level
172	ELSE ! centered hpg (default case)
173	CALL eos_tan( tsn, tsn_tl, & ! now (swap=before) in situ density for dynhpg module
174	& rhd_tl, rhop_tl )
175
176	IF( ln_zps ) CALL zps_hde_tan( kstp, jpts, tsn, tsn_tl, rhd_tl, &
177	& gtsu_tl, gru_tl, gtsv_tl, & ! Partial steps: time filtered hor. gradient
178	& grv_tl )
179	CALL tra_nxt_tan( kstp ) ! tracer fields at next time step
180	ENDIF
181
182	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
183	! Dynamics (ta, sa used as workspace)
184	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
185	ua_tl(:,:,:) = 0.0_wp ! set dynamics trends to zero
186	va_tl(:,:,:) = 0.0_wp
187
188	! apply dynamics assimilation increment ... not needed in tangent
189	CALL dyn_adv_tan( kstp ) ! advection (vector or flux form)
190	CALL dyn_vor_tan( kstp ) ! vorticity term including Coriolis
191	CALL dyn_ldf_tan( kstp ) ! lateral mixing
192	CALL dyn_hpg_tan( kstp ) ! horizontal gradient of Hydrostatic pressure
193	CALL dyn_bfr_tan( kstp ) ! bottom friction
194	CALL dyn_zdf_tan( kstp ) ! vertical diffusion
195	CALL dyn_spg_tan( kstp, indic ) ! surface pressure gradient
196	CALL dyn_nxt_tan( kstp ) ! lateral velocity at next time step
197	CALL ssh_nxt_tan( kstp ) ! sea surface height at next time step
198
199	! vertical mesh at next time step ... not available
200	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
201	! Computation of diagnostic variables
202	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
203	! N.B. ua, va, ta, sa arrays are used as workspace in this section
204	!-----------------------------------------------------------------------
205	CALL tl_trj_wri( kstp )
206	CALL trj_rea( kstp, 1) ! ... Read basic state trajectory at end of current step
207
208	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
209	! Control, and restarts
210	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
211	! N.B. ua, va, ta, sa arrays are used as workspace in this section
212	!-----------------------------------------------------------------------
213	! ! Time loop: control and print
214	!*B This fails with cgmod. To be revised CALL stp_ctl_tan( kstp, indic, 0 )
215	IF( indic < 0 ) CALL ctl_stop( 'step_tan: indic < 0' )
216	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
217	! diagnostics and outputs
218	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
219	! N.B. ua, va, ta, sa arrays are used as workspace in this section
220	!-----------------------------------------------------------------------
221	IF ( nstop == 0 ) THEN ! Diagnostics
222	! not needed in tangent
223	ENDIF
224	END SUBROUTINE stp_tan
225
226	SUBROUTINE stp_adj( kstp )
227	!!----------------------------------------------------------------------
228	!! * ROUTINE stp_adj *
229	!!
230	!! ** Purpose of the direct routine:
231	!! - Time stepping of OPA (momentum and active tracer eqs.)
232	!! - Time stepping of LIM (dynamic and thermodynamic eqs.)
233	!! - Tme stepping of TRC (passive tracer eqs.)
234	!!
235	!! ** Method of the direct routine:
236	!! -1- Update forcings and data
237	!! -2- Update ocean physics
238	!! -3- Compute the t and s trends
239	!! -4- Update t and s
240	!! -5- Compute the momentum trends
241	!! -6- Update the horizontal velocity
242	!! -7- Compute the diagnostics variables (rd,N2, div,cur,w)
243	!! -8- Outputs and diagnostics
244	!!----------------------------------------------------------------------
245	!! * Arguments
246	INTEGER, INTENT( in ) :: kstp ! ocean time-step index
247	!! * local declarations
248	INTEGER :: indic ! error indicator if < 0
249	!! ---------------------------------------------------------------------
250
251	indic = 1 ! reset to no error condition
252	CALL day_tam( kstp, 1 ) ! Calendar
253
254	! ... Read basic state trajectory at end of previous step
255	CALL trj_rea( ( kstp - 1 ), -1 )
256
257	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
258	! Dynamics (ta, sa used as workspace)
259	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
260	! apply dynamics assimilation increment ... not needed in adjoint
261	indic=0
262	! vertical mesh at next time step ... not available
263	CALL ssh_nxt_adj( kstp ) ! sea surface height at next time step
264	CALL dyn_nxt_adj( kstp ) ! lateral velocity at next time step
265	CALL dyn_spg_adj( kstp, indic ) ! surface pressure gradient
266	CALL dyn_zdf_adj( kstp ) ! vertical diffusion
267	CALL dyn_bfr_adj( kstp ) ! bottom friction
268	CALL dyn_hpg_adj( kstp ) ! horizontal gradient of Hydrostatic pressure
269	CALL dyn_ldf_adj( kstp ) ! lateral mixing
270	CALL dyn_vor_adj( kstp ) ! vorticity term including Coriolis
271	CALL dyn_adv_adj( kstp ) ! advection (vector or flux form)
272
273	ua_ad(:,:,:) = 0.0_wp ! set tracer trends to zero
274	va_ad(:,:,:) = 0.0_wp
275
276	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
277	! Active tracers (ua, va used as workspace)
278	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
279
280	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg
281	IF( ln_zps ) CALL zps_hde_adj( kstp, jpts, tsa, tsa_ad, & ! Partial steps: time filtered hor. gradient
282	& rhd_ad, gtsu_ad, gru_ad, gtsv_ad, & ! of t, s, rd at the bottom ocean level
283	& grv_ad )
284	CALL eos_adj( tsa, tsa_ad, rhd_ad, rhop_ad ) ! Time-filtered in situ density used in dynhpg module
285	CALL tra_nxt_adj( kstp ) ! tracer fields at next time step
286	ELSE ! centered hpg (default case)
287	CALL tra_nxt_adj( kstp ) ! tracer fields at next time step
288	IF( ln_zps ) CALL zps_hde_adj( kstp, jpts, tsn, tsn_ad, & ! Partial steps: time filtered hor. gradient
289	& rhd_ad, gtsu_ad, gru_ad, gtsv_ad, & ! of t, s, rd at the bottom ocean level
290	& grv_ad )
291	CALL eos_adj( tsn, tsn_ad, rhd_ad, rhop_ad ) ! now (swap=before) in situ density for dynhpg module
292	ENDIF
293
294	! update the new (t,s) fields by non
295	! penetrative convective adjustment ... not available
296
297	CALL tra_zdf_adj( kstp ) ! vertical mixing
298	CALL tra_ldf_adj( kstp ) ! lateral mixing
299	CALL tra_adv_adj( kstp ) ! horizontal & vertical advection
300	IF( ln_tradmp ) CALL tra_dmp_adj( kstp ) ! internal damping trends
301	!! advective (and/or diffusive) bottom boundary layer scheme ... currently not available
302	IF( lk_trabbl ) CALL tra_bbl_adj( kstp ) ! diffusive bottom boundary layer scheme
303	IF( ln_trabbc ) CALL tra_bbc_adj( kstp ) ! bottom heat flux
304
305	!IF( ln_traqsr ) CALL tra_qsr_adj( kstp ) ! penetrative solar radiation qsr
306
307	CALL tra_sbc_adj( kstp ) ! surface boundary condition
308
309	tsa_ad(:,:,:,:) = 0.0_wp ! set tracer trends to zero
310
311	! apply tracer assimilation increment ... not needed in adjoint
312	#if defined key_top
313	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
314	! Passive Tracer Model
315	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
316	! N.B. ua, va, ta, sa arrays are used as workspace in this section
317	!-----------------------------------------------------------------------
318
319	! time-stepping... not needed in adjoint for the time being
320
321	#endif
322	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
323	! Ocean physics update (ua, va, ta, sa used as workspace)
324	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
325	!-----------------------------------------------------------------------
326	! LATERAL PHYSICS
327	!-----------------------------------------------------------------------
328	#if defined key_traldf_c2d
329	! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_adj
330	#endif
331	! before slope of the lateral mixing... not needed in adjoint with lk_zdfcst_adj
332	!-----------------------------------------------------------------------
333	! VERTICAL PHYSICS
334	!-----------------------------------------------------------------------
335	CALL zdf_bfr_adj ( kstp )
336	!! N.B. ua, va, ta, sa arrays are used as workspace in this section
337	!!-----------------------------------------------------------------------
338	!! mixed layer depth... not needed in adjoint with lk_zdfcst_adj
339	!! bottom friction... not needed in adjoint with lk_zdfcst_adj
340	!! double diffusive mixing ... not needed in adjoint with lk_zdfcst_adj
341	!! enhanced vertical eddy diffusivity ... not needed in adjoint with lk_zdfcst_adj
342	!! ! increase diffusivity at rivers mouths... not needed in tangent
343	!! lk_zdfcst_adj: Constant Kz read from the reference trajectory
344	!! Constant Kz (reset avt, avm[uv] to the background value)... not available
345	!! KPP closure scheme for Kz ... not available
346	!! TKE closure scheme for Kz ... not available
347	!! Richardson number dependent Kz ... not available
348	!! ! Vertical eddy viscosity and diffusivity coefficients
349	CALL bn2_adj( tsn, tsn_ad, rn2_ad ) ! now Brunt-Vaisala frequency
350	CALL bn2_adj( tsb, tsb_ad, rn2b_ad ) ! now Brunt-Vaisala frequency
351	!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
352	!! Ocean dynamics : ssh, wn, hdiv, rot !
353	!!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
354	CALL ssh_wzv_adj( kstp ) ! after ssh & vertical velocity
355	!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
356	!! Update data, open boundaries and Forcings
357	!!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
358	!! Output the initial state and forcings ... not needed in adjoint
359	CALL sbc_adj ( kstp ) ! Sea boundary condition (including sea-ice)
360	! update 3D salinity data ... not needed in tangent
361	! update 3D temperature data ... not needed in adjoint
362
363	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
364	! Control, and restarts
365	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
366	! N.B. ua, va, ta, sa arrays are used as workspace in this section
367	!-----------------------------------------------------------------------
368	! ! Time loop: control and print
369	CALL stp_ctl_adj( kstp, indic, 0 )
370	IF( indic < 0 ) CALL ctl_stop( 'step_adj: indic < 0' )
371
372	! close input ocean restart file ... not needed in adjoint
373	! write output ocean restart file... not needed in adjoint
374	! write open boundary restart file... not needed in adjoint
375	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
376	! diagnostics and outputs
377	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
378	! N.B. ua, va, ta, sa arrays are used as workspace in this section
379	!-----------------------------------------------------------------------
380
381	IF ( nstop == 0 ) THEN ! Diagnostics
382	! not needed in adjoint
383	ENDIF
384
385	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
386	! Coupled mode
387	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
388
389	! coupled mode : field exchanges ... not available for the next 20 years
390	!
391	!
392	IF( nn_timing == 1 .AND. kstp == nit000 ) CALL timing_reset
393	!
394	END SUBROUTINE stp_adj
395
396	SUBROUTINE stp_adj_tst( kumadt )
397	!!-----------------------------------------------------------------------
398	!!
399	!! * ROUTINE example_adj_tst *
400	!!
401	!! ** Purpose : Test the adjoint routine.
402	!!
403	!! ** Method : Verify the scalar product
404	!!
405	!! ( L dx )^T W dy = dx^T L^T W dy
406	!!
407	!! where L = tangent routine
408	!! L^T = adjoint routine
409	!! W = diagonal matrix of scale factors
410	!! dx = input perturbation (random field)
411	!! dy = L dx
412	!!
413	!!
414	!! History :
415	!! ! 09-03 (F. Vigilant)
416	!!-----------------------------------------------------------------------
417	!! * Modules used
418	USE sbcssr_tam, ONLY: &
419	& qrp_tl, qrp_ad, &
420	& erp_tl, erp_ad
421	USE sbcfwb_tam, ONLY: &
422	& a_fwb_tl, & ! for before year.
423	& a_fwb_ad ! for before year.
424
425	!! * Arguments
426	INTEGER, INTENT(IN) :: &
427	& kumadt ! Output unit
428
429	!! * Local declarations
430	INTEGER :: &
431	& ji, & ! dummy loop indices
432	& jj, &
433	& jk, &
434	& istp
435	INTEGER, DIMENSION(jpi,jpj) :: &
436	& iseed_2d ! 2D seed for the random number generator
437	REAL(KIND=wp) :: &
438	& zsp1 , & ! scalar product involving the tangent routine
439	& zsp1_U , &
440	& zsp1_V , &
441	& zsp1_T , &
442	& zsp1_S , &
443	& zsp1_SSH, &
444	& zsp2 , & ! scalar product involving the adjoint routine
445	& zsp2_U , &
446	& zsp2_V , &
447	& zsp2_T , &
448	& zsp2_S , &
449	& zsp2_SSH
450	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
451	& zun_tlin , & ! Tangent input
452	& zvn_tlin , & ! Tangent input
453	& ztn_tlin , & ! Tangent input
454	& zsn_tlin , & ! Tangent input
455	& zun_tlout, & ! Tangent output
456	& zvn_tlout, & ! Tangent output
457	& ztn_tlout, & ! Tangent output
458	& zsn_tlout, & ! Tangent outpu
459	& zun_adin , & ! Adjoint input
460	& zvn_adin , & ! Adjoint input
461	& ztn_adin , & ! Adjoint input
462	& zsn_adin , & ! Adjoint input
463	& zun_adout, & ! Adjoint output
464	& zvn_adout, & ! Adjoint output
465	& ztn_adout, & ! Adjoint output
466	& zsn_adout, & ! Adjoint output
467	& z3r ! 3D random field
468	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
469	& zsshn_tlin , & ! Tangent input
470	& zsshn_tlout , & ! Tangent input
471	& zsshn_adin , & ! Adjoint output
472	& zsshn_adout , & ! Adjoint output
473	& z2r ! 2D random field
474
475	CHARACTER(LEN=14) :: &
476	& cl_name
477
478	INTEGER :: &
479	& jpert
480	INTEGER, PARAMETER :: jpertmax = 6
481
482	! Allocate memory
483
484	ALLOCATE( &
485	& zun_tlin( jpi,jpj,jpk), zvn_tlin( jpi,jpj,jpk), ztn_tlin( jpi,jpj,jpk), &
486	& zsn_tlin( jpi,jpj,jpk), zsshn_tlin( jpi,jpj ), zun_tlout( jpi,jpj,jpk), &
487	& zvn_tlout( jpi,jpj,jpk), ztn_tlout( jpi,jpj,jpk), zsn_tlout( jpi,jpj,jpk), &
488	& zsshn_tlout(jpi,jpj ), zun_adin( jpi,jpj,jpk), zvn_adin( jpi,jpj,jpk), &
489	& ztn_adin( jpi,jpj,jpk), zsn_adin( jpi,jpj,jpk), zsshn_adin(jpi,jpj ), &
490	& zun_adout( jpi,jpj,jpk), zvn_adout( jpi,jpj,jpk), ztn_adout( jpi,jpj,jpk), &
491	& zsn_adout( jpi,jpj,jpk), zsshn_adout(jpi,jpj ), z2r( jpi,jpj ), &
492	& z3r( jpi,jpj,jpk) &
493	& )
494
495	!==================================================================
496	! 1) dx = ( un_tl, vn_tl, tn_tl, sn_tl, sshn_tl ) and
497	! dy = ( ..... )
498	!==================================================================
499
500
501	DO jpert = jpertmax, 1, -1
502
503	!--------------------------------------------------------------------
504	! Reset the tangent and adjoint variables
505	!--------------------------------------------------------------------
506	zun_tlin (:,:,:) = 0.0_wp
507	zvn_tlin (:,:,:) = 0.0_wp
508	ztn_tlin (:,:,:) = 0.0_wp
509	zsn_tlin (:,:,:) = 0.0_wp
510	zsshn_tlin ( :,:) = 0.0_wp
511	zun_tlout (:,:,:) = 0.0_wp
512	zvn_tlout (:,:,:) = 0.0_wp
513	ztn_tlout (:,:,:) = 0.0_wp
514	zsn_tlout (:,:,:) = 0.0_wp
515	zsshn_tlout( :,:) = 0.0_wp
516	zun_adin (:,:,:) = 0.0_wp
517	zvn_adin (:,:,:) = 0.0_wp
518	ztn_adin (:,:,:) = 0.0_wp
519	zsn_adin (:,:,:) = 0.0_wp
520	zsshn_adin ( :,:) = 0.0_wp
521	zun_adout (:,:,:) = 0.0_wp
522	zvn_adout (:,:,:) = 0.0_wp
523	ztn_adout (:,:,:) = 0.0_wp
524	zsn_adout (:,:,:) = 0.0_wp
525	zsshn_adout( :,:) = 0.0_wp
526	z3r (:,:,:) = 0.0_wp
527	z2r ( :,:) = 0.0_wp
528
529	!--------------------------------------------------------------------
530	! Initialize the tangent input with random noise: dx
531	!--------------------------------------------------------------------
532
533	IF ( (jpert == 1) .OR. (jpert == jpertmax) ) THEN
534	CALL grid_random( z3r, 'U', 0.0_wp, stdu )
535	DO jk = 1, jpk
536	DO jj = nldj, nlej
537	DO ji = nldi, nlei
538	zun_tlin(ji,jj,jk) = z3r(ji,jj,jk)
539	END DO
540	END DO
541	END DO
542	ENDIF
543	IF ( (jpert == 2) .OR. (jpert == jpertmax) ) THEN
544	CALL grid_random( z3r, 'V', 0.0_wp, stdv )
545	DO jk = 1, jpk
546	DO jj = nldj, nlej
547	DO ji = nldi, nlei
548	zvn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
549	END DO
550	END DO
551	END DO
552	ENDIF
553	IF ( (jpert == 3) .OR. (jpert == jpertmax) ) THEN
554	CALL grid_random( z3r, 'T', 0.0_wp, stdt )
555	DO jk = 1, jpk
556	DO jj = nldj, nlej
557	DO ji = nldi, nlei
558	ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
559	END DO
560	END DO
561	END DO
562	ENDIF
563	IF ( (jpert == 4) .OR. (jpert == jpertmax) ) THEN
564	CALL grid_random( z3r, 'T', 0.0_wp, stds )
565	DO jk = 1, jpk
566	DO jj = nldj, nlej
567	DO ji = nldi, nlei
568	zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
569	END DO
570	END DO
571	END DO
572	ENDIF
573	IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN
574	CALL grid_random( z2r, 'T', 0.0_wp, stdssh )
575	DO jj = nldj, nlej
576	DO ji = nldi, nlei
577	zsshn_tlin(ji,jj) = z2r(ji,jj)
578	END DO
579	END DO
580	ENDIF
581
582	CALL oce_tam_init( 1 ) ! allocate/initialize tl variables
583	CALL sbc_oce_tam_init( 1 )
584	CALL sol_oce_tam_init( 1 )
585	#if defined key_tradmp
586	CALL trc_oce_tam_init( 1 )
587	strdmp_tl = 0.0_wp
588	ttrdmp_tl = 0.0_wp
589	#endif
590
591	!CALL oce_tam_init( 1 ) ! allocate/initialize tl variables
592	!CALL sbc_oce_tam_init( 1 )
593	!CALL sol_oce_tam_init( 1 )
594	!CALL trc_oce_tam_init( 1 )
595
596	qrp_tl = 0.0_wp
597	erp_tl = 0.0_wp
598
599	emp_tl(:,:) = 0.0_wp
600	a_fwb_tl = 0.0_wp
601
602	istp = nit000 - 1
603	CALL trj_rea( istp, 1 )
604	!--------------------------------------------------------------------
605	! Initialize the tangent variables: dy^* = W dy
606	!--------------------------------------------------------------------
607
608	un_tl (:,:,:) = zun_tlin (:,:,:)
609	vn_tl (:,:,:) = zvn_tlin (:,:,:)
610	tsn_tl (:,:,:,jp_tem) = ztn_tlin (:,:,:)
611	tsn_tl (:,:,:,jp_sal) = zsn_tlin (:,:,:)
612	sshn_tl( :,:) = zsshn_tlin ( :,:)
613
614	!CALL oce_tam_deallocate( 2 ) ! deallocate adj variables
615	!CALL sbc_oce_tam_deallocate( 2 )
616	!CALL sol_oce_tam_deallocate( 2 )
617
618	!-----------------------------------------------------------------------
619	! Initialization of the dynamics and tracer fields for the tangent
620	!-----------------------------------------------------------------------
621
622	CALL istate_init_tan
623
624	DO istp = nit000, nitend, 1
625	CALL stp_tan( istp )
626	END DO
627
628	zun_tlout ( :,:,:) = un_tl (:,:,:)
629	zvn_tlout ( :,:,:) = vn_tl (:,:,:)
630	ztn_tlout ( :,:,:) = tsn_tl (:,:,:,jp_tem)
631	zsn_tlout ( :,:,:) = tsn_tl (:,:,:,jp_sal)
632	zsshn_tlout( :,:) = sshn_tl ( :,:)
633
634	!--------------------------------------------------------------------
635	! Initialize the adjoint variables: dy^* = W dy
636	!--------------------------------------------------------------------
637
638	DO jk = 1, jpk
639	DO jj = nldj, nlej
640	DO ji = nldi, nlei
641	zun_adin(ji,jj,jk) = zun_tlout(ji,jj,jk) &
642	& * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
643	& * umask(ji,jj,jk) * wesp_u
644	zvn_adin(ji,jj,jk) = zvn_tlout(ji,jj,jk) &
645	& * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
646	& * vmask(ji,jj,jk) * wesp_u
647	ztn_adin(ji,jj,jk) = ztn_tlout(ji,jj,jk) &
648	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
649	& * tmask(ji,jj,jk) * wesp_t(jk)
650	zsn_adin(ji,jj,jk) = zsn_tlout(ji,jj,jk) &
651	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
652	& * tmask(ji,jj,jk) * wesp_s(jk)
653	END DO
654	END DO
655	END DO
656
657	DO jj = nldj, nlej
658	DO ji = nldi, nlei
659	zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) &
660	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
661	& * tmask(ji,jj,1) * wesp_ssh
662	END DO
663	END DO
664
665	!--------------------------------------------------------------------
666	! Compute the scalar product: ( L dx )^T W dy
667	!--------------------------------------------------------------------
668
669	zsp1_U = DOT_PRODUCT( zun_tlout , zun_adin )
670	zsp1_V = DOT_PRODUCT( zvn_tlout , zvn_adin )
671	zsp1_T = DOT_PRODUCT( ztn_tlout , ztn_adin )
672	zsp1_S = DOT_PRODUCT( zsn_tlout , zsn_adin )
673	zsp1_SSH = DOT_PRODUCT( zsshn_tlout, zsshn_adin )
674
675	zsp1 = zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH
676
677	!--------------------------------------------------------------------
678	! Call the adjoint routine: dx^* = L^T dy^*
679	!--------------------------------------------------------------------
680	CALL oce_tam_init( 2 ) ! allocate/initialize adj variables
681	CALL sbc_oce_tam_init( 2 )
682	CALL sol_oce_tam_init( 2 )
683	#if defined key_tradmp
684	CALL trc_oce_tam_init( 2 )
685	strdmp_ad = 0.0_wp
686	ttrdmp_ad = 0.0_wp
687	#endif
688	qrp_ad = 0.0_wp
689	erp_ad = 0.0_wp
690	emp_ad(:,:) = 0.0_wp
691
692	a_fwb_ad = 0.0_wp
693
694	un_ad (:,:,:) = zun_adin (:,:,:)
695	vn_ad (:,:,:) = zvn_adin (:,:,:)
696	tsn_ad (:,:,:,jp_tem) = ztn_adin (:,:,:)
697	tsn_ad (:,:,:,jp_sal) = zsn_adin (:,:,:)
698	sshn_ad( :,:) = zsshn_adin ( :,:)
699
700	!CALL oce_tam_deallocate( 1 ) !deallocate tl variables
701	!CALL sbc_oce_tam_deallocate( 1 )
702	!CALL sol_oce_tam_deallocate( 1 )
703	!CALL trc_oce_tam_deallocate( 1 )
704	istp = nitend
705	CALL trj_rea( istp, -1 )
706
707	DO istp = nitend, nit000, -1
708	CALL stp_adj(istp)
709	END DO
710
711	CALL istate_init_adj
712
713	zun_adout ( :,:,:) = un_ad (:,:,:)
714	zvn_adout ( :,:,:) = vn_ad (:,:,:)
715	ztn_adout ( :,:,:) = tsn_ad (:,:,:,jp_tem)
716	zsn_adout ( :,:,:) = tsn_ad (:,:,:,jp_sal)
717	zsshn_adout( :,:) = sshn_ad ( :,:)
718
719	zsp2_U = DOT_PRODUCT( zun_tlin , zun_adout )
720	zsp2_V = DOT_PRODUCT( zvn_tlin , zvn_adout )
721	zsp2_T = DOT_PRODUCT( ztn_tlin , ztn_adout )
722	zsp2_S = DOT_PRODUCT( zsn_tlin , zsn_adout )
723	zsp2_SSH = DOT_PRODUCT( zsshn_tlin, zsshn_adout )
724
725	zsp2 = zsp2_U + zsp2_V + zsp2_T + zsp2_S + zsp2_SSH
726
727	!CALL oce_tam_deallocate( 2 ) !deallocate adj variables
728	!CALL sbc_oce_tam_deallocate( 2 )
729	!CALL sol_oce_tam_deallocate( 2 )
730	!CALL trc_oce_tam_deallocate( 2 )
731
732	! 14 char:'12345678901234'
733	SELECT CASE (jpert)
734	CASE(1)
735	cl_name = 'step_adj U'
736	CASE(2)
737	cl_name = 'step_adj V'
738	CASE(3)
739	cl_name = 'step_adj T'
740	CASE(4)
741	cl_name = 'step_adj S'
742	CASE(5)
743	cl_name = 'step_adj SSH'
744	CASE(jpertmax)
745	cl_name = 'step_adj '
746	END SELECT
747	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
748	END DO
749
750	nitsor(:) = jp_it0adj ! restore nitsor to avoid non reproducible results with or without the tests
751
752	DEALLOCATE( &
753	& zun_tlin , zvn_tlin , ztn_tlin , &
754	& zsn_tlin , zsshn_tlin , zun_tlout , &
755	& zvn_tlout , ztn_tlout , zsn_tlout , &
756	& zsshn_tlout, zun_adin , zvn_adin , &
757	& ztn_adin , zsn_adin , zsshn_adin , &
758	& zun_adout , zvn_adout , ztn_adout , &
759	& zsn_adout , zsshn_adout, z2r , &
760	& z3r &
761	& )
762
763	END SUBROUTINE stp_adj_tst
764
765	!!----------------------------------------------------------------------
766
767	!!======================================================================
768	#endif
769	END MODULE step_tam

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

Download in other formats: